OOPS
oops/interface/ObsVector.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2009-2016 ECMWF.
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  * In applying this licence, ECMWF does not waive the privileges and immunities
7  * granted to it by virtue of its status as an intergovernmental organisation nor
8  * does it submit to any jurisdiction.
9  */
10 
11 #ifndef OOPS_INTERFACE_OBSVECTOR_H_
12 #define OOPS_INTERFACE_OBSVECTOR_H_
13 
14 #include <Eigen/Dense>
15 #include <math.h>
16 #include <memory>
17 #include <ostream>
18 #include <string>
19 
20 #include "oops/interface/ObsDataVector.h"
21 #include "oops/interface/ObsSpace.h"
22 #include "oops/util/gatherPrint.h"
23 #include "oops/util/Logger.h"
24 #include "oops/util/ObjectCounter.h"
25 #include "oops/util/Printable.h"
26 #include "oops/util/Timer.h"
27 
28 namespace eckit {
29  class Configuration;
30 }
31 
32 namespace util {
33  class DateTime;
34 }
35 
36 namespace oops {
37 
38 // -----------------------------------------------------------------------------
39 
40 template <typename OBS>
41 class ObsVector : public util::Printable,
42  private util::ObjectCounter<ObsVector<OBS> > {
43  typedef typename OBS::ObsVector ObsVector_;
44 
45  public:
46  static const std::string classname() {return "oops::ObsVector";}
47 
48  explicit ObsVector(const ObsSpace<OBS> &, const std::string name = "", const bool fail = true);
49  explicit ObsVector(const ObsVector &);
50  ObsVector(const ObsSpace<OBS> &, const ObsVector &);
51  ~ObsVector();
52 
53 /// Interfacing
54  ObsVector_ & obsvector() {return *data_;}
55  const ObsVector_ & obsvector() const {return *data_;}
56 
57 // Linear algebra
58  ObsVector & operator = (const ObsVector &);
59  ObsVector & operator*= (const double &);
60  ObsVector & operator+= (const ObsVector &);
61  ObsVector & operator-= (const ObsVector &);
62  ObsVector & operator*= (const ObsVector &);
63  ObsVector & operator/= (const ObsVector &);
64 
65 /// Pack into an Eigen vector (excluding vector elements that are masked out)
66  Eigen::VectorXd packEigen() const;
67 
68  void zero();
69  void axpy(const double &, const ObsVector &);
70  void invert();
71  void random();
72  double dot_product_with(const ObsVector &) const;
73  double rms() const;
74 /// Mask out elements of the vector where the passed in flags are > 0
75  void mask(const ObsDataVector<OBS, int> &);
76 
77 // I/O
78  void save(const std::string &) const;
79  void read(const std::string &);
80 
81  unsigned int nobs() const;
82 
83  private:
84  void print(std::ostream &) const;
85  std::unique_ptr<ObsVector_> data_;
86  const eckit::mpi::Comm & commTime_;
87 };
88 
89 // -----------------------------------------------------------------------------
90 template <typename OBS>
91 ObsVector<OBS>::ObsVector(const ObsSpace<OBS> & os, const std::string name,
92  const bool fail): data_(), commTime_(os.timeComm()) {
93  Log::trace() << "ObsVector<OBS>::ObsVector starting " << name << std::endl;
94  util::Timer timer(classname(), "ObsVector");
95 
96  data_.reset(new ObsVector_(os.obsspace(), name, fail));
97 
98  Log::trace() << "ObsVector<OBS>::ObsVector done" << std::endl;
99 }
100 // -----------------------------------------------------------------------------
101 template <typename OBS>
102 ObsVector<OBS>::ObsVector(const ObsVector & other): data_(), commTime_(other.commTime_) {
103  Log::trace() << "ObsVector<OBS>::ObsVector starting" << std::endl;
104  util::Timer timer(classname(), "ObsVector");
105 
106  data_.reset(new ObsVector_(*other.data_));
107 
108  Log::trace() << "ObsVector<OBS>::ObsVector done" << std::endl;
109 }
110 // -----------------------------------------------------------------------------
111 template <typename OBS>
112 ObsVector<OBS>::ObsVector(const ObsSpace<OBS> & os, const ObsVector & other)
113  : data_(), commTime_(os.timeComm())
114 {
115  Log::trace() << "ObsVector<OBS>::ObsVector starting" << std::endl;
116  util::Timer timer(classname(), "ObsVector");
117 
118  data_.reset(new ObsVector_(os.obsspace(), other.obsvector()));
119 
120  Log::trace() << "ObsVector<OBS>::ObsVector done" << std::endl;
121 }
122 // -----------------------------------------------------------------------------
123 template <typename OBS>
124 ObsVector<OBS>::~ObsVector() {
125  Log::trace() << "ObsVector<OBS>::~ObsVector starting" << std::endl;
126  util::Timer timer(classname(), "~ObsVector");
127 
128  data_.reset();
129 
130  Log::trace() << "ObsVector<OBS>::~ObsVector done" << std::endl;
131 }
132 // -----------------------------------------------------------------------------
133 template <typename OBS>
134 ObsVector<OBS> & ObsVector<OBS>::operator=(const ObsVector & rhs) {
135  Log::trace() << "ObsVector<OBS>::operator= starting" << std::endl;
136  util::Timer timer(classname(), "operator=");
137 
138  *data_ = *rhs.data_;
139 
140  Log::trace() << "ObsVector<OBS>::operator= done" << std::endl;
141  return *this;
142 }
143 // -----------------------------------------------------------------------------
144 template <typename OBS>
145 ObsVector<OBS> & ObsVector<OBS>::operator*=(const double & zz) {
146  Log::trace() << "ObsVector<OBS>::operator*= starting" << std::endl;
147  util::Timer timer(classname(), "operator*=");
148 
149  *data_ *= zz;
150 
151  Log::trace() << "ObsVector<OBS>::operator*= done" << std::endl;
152  return *this;
153 }
154 // -----------------------------------------------------------------------------
155 template <typename OBS>
156 ObsVector<OBS> & ObsVector<OBS>::operator+=(const ObsVector & rhs) {
157  Log::trace() << "ObsVector<OBS>::operator+= starting" << std::endl;
158  util::Timer timer(classname(), "operator+=");
159 
160  *data_ += *rhs.data_;
161 
162  Log::trace() << "ObsVector<OBS>::operator+= done" << std::endl;
163  return *this;
164 }
165 // -----------------------------------------------------------------------------
166 template <typename OBS>
167 ObsVector<OBS> & ObsVector<OBS>::operator-=(const ObsVector & rhs) {
168  Log::trace() << "ObsVector<OBS>::operator-= starting" << std::endl;
169  util::Timer timer(classname(), "operator-=");
170 
171  *data_ -= *rhs.data_;
172 
173  Log::trace() << "ObsVector<OBS>::operator-= done" << std::endl;
174  return *this;
175 }
176 // -----------------------------------------------------------------------------
177 template <typename OBS>
178 ObsVector<OBS> & ObsVector<OBS>::operator*=(const ObsVector & rhs) {
179  Log::trace() << "ObsVector<OBS>::operator*= starting" << std::endl;
180  util::Timer timer(classname(), "operator*=");
181 
182  *data_ *= *rhs.data_;
183 
184  Log::trace() << "ObsVector<OBS>::operator*= done" << std::endl;
185  return *this;
186 }
187 // -----------------------------------------------------------------------------
188 template <typename OBS>
189 ObsVector<OBS> & ObsVector<OBS>::operator/=(const ObsVector & rhs) {
190  Log::trace() << "ObsVector<OBS>::operator/= starting" << std::endl;
191  util::Timer timer(classname(), "operator/=");
192 
193  *data_ /= *rhs.data_;
194 
195  Log::trace() << "ObsVector<OBS>::operator/= done" << std::endl;
196  return *this;
197 }
198 // -----------------------------------------------------------------------------
199 template <typename OBS>
200 void ObsVector<OBS>::zero() {
201  Log::trace() << "ObsVector<OBS>::zero starting" << std::endl;
202  util::Timer timer(classname(), "zero");
203 
204  data_->zero();
205 
206  Log::trace() << "ObsVector<OBS>::zero done" << std::endl;
207 }
208 // -----------------------------------------------------------------------------
209 template <typename OBS>
210 void ObsVector<OBS>::axpy(const double & zz, const ObsVector & rhs) {
211  Log::trace() << "ObsVector<OBS>::axpy starting" << std::endl;
212  util::Timer timer(classname(), "axpy");
213 
214  data_->axpy(zz, *rhs.data_);
215 
216  Log::trace() << "ObsVector<OBS>::axpy done" << std::endl;
217 }
218 // -----------------------------------------------------------------------------
219 template <typename OBS>
220 void ObsVector<OBS>::invert() {
221  Log::trace() << "ObsVector<OBS>::invert starting" << std::endl;
222  util::Timer timer(classname(), "invert");
223 
224  data_->invert();
225 
226  Log::trace() << "ObsVector<OBS>::invert done" << std::endl;
227 }
228 // -----------------------------------------------------------------------------
229 template <typename OBS>
230 void ObsVector<OBS>::random() {
231  Log::trace() << "ObsVector<OBS>::random starting" << std::endl;
232  util::Timer timer(classname(), "random");
233 
234  data_->random();
235 
236  Log::trace() << "ObsVector<OBS>::random done" << std::endl;
237 }
238 // -----------------------------------------------------------------------------
239 template <typename OBS>
240 double ObsVector<OBS>::dot_product_with(const ObsVector & other) const {
241  Log::trace() << "ObsVector<OBS>::dot_product starting" << std::endl;
242  util::Timer timer(classname(), "dot_product");
243 
244  double zz = data_->dot_product_with(*other.data_);
245  commTime_.allReduceInPlace(zz, eckit::mpi::Operation::SUM);
246 
247  Log::trace() << "ObsVector<OBS>::dot_product done" << std::endl;
248  return zz;
249 }
250 // -----------------------------------------------------------------------------
251 template <typename OBS>
252 void ObsVector<OBS>::mask(const ObsDataVector<OBS, int> & qc) {
253  Log::trace() << "ObsVector<OBS>::mask starting" << std::endl;
254  util::Timer timer(classname(), "mask");
255  data_->mask(qc.obsdatavector());
256  Log::trace() << "ObsVector<OBS>::mask done" << std::endl;
257 }
258 // -----------------------------------------------------------------------------
259 template <typename OBS>
260 double ObsVector<OBS>::rms() const {
261  Log::trace() << "ObsVector<OBS>::rms starting" << std::endl;
262  util::Timer timer(classname(), "rms");
263 
264  double zz = 0.0;
265  size_t ntot = this->nobs();
266  if (ntot > 0) {
267  zz = data_->rms();
268  double zzz = zz * zz * static_cast<double>(data_->nobs());
269  commTime_.allReduceInPlace(zzz, eckit::mpi::Operation::SUM);
270  zzz /= static_cast<double>(ntot);
271  zz = std::sqrt(zzz);
272  }
273 
274  Log::trace() << "ObsVector<OBS>::rms done" << std::endl;
275  return zz;
276 }
277 // -----------------------------------------------------------------------------
278 template <typename OBS>
279 unsigned int ObsVector<OBS>::nobs() const {
280  int nobs = data_->nobs();
281  commTime_.allReduceInPlace(nobs, eckit::mpi::Operation::SUM);
282  return nobs;
283 }
284 // -----------------------------------------------------------------------------
285 template <typename OBS>
286 void ObsVector<OBS>::print(std::ostream & os) const {
287  Log::trace() << "ObsVector<OBS>::print starting" << std::endl;
288  util::Timer timer(classname(), "print");
289  if (commTime_.size() > 1) {
290  gatherPrint(os, *data_, commTime_);
291  } else {
292  os << *data_;
293  }
294  Log::trace() << "ObsVector<OBS>::print done" << std::endl;
295 }
296 // -----------------------------------------------------------------------------
297 template <typename OBS>
298 void ObsVector<OBS>::save(const std::string & name) const {
299  Log::trace() << "ObsVector<OBS>::save starting " << name << std::endl;
300  util::Timer timer(classname(), "save");
301 
302  data_->save(name);
303 
304  Log::trace() << "ObsVector<OBS>::save done" << std::endl;
305 }
306 // -----------------------------------------------------------------------------
307 template <typename OBS>
308 Eigen::VectorXd ObsVector<OBS>::packEigen() const {
309  Log::trace() << "ObsVector<OBS>::packEigen starting " << std::endl;
310  util::Timer timer(classname(), "packEigen");
311 
312  Eigen::VectorXd vec = data_->packEigen();
313  ASSERT(vec.size() == nobs());
314 
315  Log::trace() << "ObsVector<OBS>::packEigen done" << std::endl;
316  return vec;
317 }
318 // -----------------------------------------------------------------------------
319 template <typename OBS>
320 void ObsVector<OBS>::read(const std::string & name) {
321  Log::trace() << "ObsVector<OBS>::read starting " << name << std::endl;
322  util::Timer timer(classname(), "read");
323 
324  data_->read(name);
325 
326  Log::trace() << "ObsVector<OBS>::read done" << std::endl;
327 }
328 // -----------------------------------------------------------------------------
329 
330 } // namespace oops
331 
332 #endif // OOPS_INTERFACE_OBSVECTOR_H_
oops::ObsVector::dot_product_with
double dot_product_with(const ObsVector &) const
Definition: oops/interface/ObsVector.h:240
oops::ObsVector::obsvector
ObsVector_ & obsvector()
Interfacing.
Definition: oops/interface/ObsVector.h:54
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
oops::ObsVector::read
void read(const std::string &)
Definition: oops/interface/ObsVector.h:320
oops::ObsVector::classname
static const std::string classname()
Definition: oops/interface/ObsVector.h:46
oops::ObsVector::ObsVector_
OBS::ObsVector ObsVector_
Definition: oops/interface/ObsVector.h:43
oops::ObsVector::operator*=
ObsVector & operator*=(const double &)
Definition: oops/interface/ObsVector.h:145
oops::ObsVector::nobs
unsigned int nobs() const
Definition: oops/interface/ObsVector.h:279
oops::ObsVector::~ObsVector
~ObsVector()
Definition: oops/interface/ObsVector.h:124
oops::ObsSpace
Definition: oops/interface/ObsSpace.h:42
oops::ObsVector::rms
double rms() const
Definition: oops/interface/ObsVector.h:260
oops::ObsSpace::obsspace
ObsSpace_ & obsspace() const
Interfacing.
Definition: oops/interface/ObsSpace.h:61
oops::ObsVector::data_
std::unique_ptr< ObsVector_ > data_
Definition: oops/interface/ObsVector.h:85
oops::ObsVector::operator-=
ObsVector & operator-=(const ObsVector &)
Definition: oops/interface/ObsVector.h:167
oops::ObsVector::operator/=
ObsVector & operator/=(const ObsVector &)
Definition: oops/interface/ObsVector.h:189
oops::ObsVector
Definition: oops/interface/ObsSpace.h:36
oops::ObsDataVector::obsdatavector
ObsDataVec_ & obsdatavector()
Interfacing.
Definition: ObsDataVector.h:39
oops::ObsVector::random
void random()
Definition: oops/interface/ObsVector.h:230
eckit
Definition: FieldL95.h:22
oops::ObsVector::obsvector
const ObsVector_ & obsvector() const
Definition: oops/interface/ObsVector.h:55
oops::ObsVector::zero
void zero()
Definition: oops/interface/ObsVector.h:200
oops::ObsVector::operator+=
ObsVector & operator+=(const ObsVector &)
Definition: oops/interface/ObsVector.h:156
oops::ObsVector::operator=
ObsVector & operator=(const ObsVector &)
Definition: oops/interface/ObsVector.h:134
oops::ObsVector::ObsVector
ObsVector(const ObsSpace< OBS > &, const std::string name="", const bool fail=true)
Definition: oops/interface/ObsVector.h:91
oops::ObsVector::save
void save(const std::string &) const
Definition: oops/interface/ObsVector.h:298
ObsSpace.h
oops::ObsVector::axpy
void axpy(const double &, const ObsVector &)
Definition: oops/interface/ObsVector.h:210
oops::ObsDataVector
Definition: ObsDataVector.h:28
oops::ObsVector::print
void print(std::ostream &) const
Definition: oops/interface/ObsVector.h:286
oops::ObsVector::invert
void invert()
Definition: oops/interface/ObsVector.h:220
oops::ObsVector::packEigen
Eigen::VectorXd packEigen() const
Pack into an Eigen vector (excluding vector elements that are masked out)
Definition: oops/interface/ObsVector.h:308
oops::ObsVector::mask
void mask(const ObsDataVector< OBS, int > &)
Mask out elements of the vector where the passed in flags are > 0.
Definition: oops/interface/ObsVector.h:252
oops::ObsVector::commTime_
const eckit::mpi::Comm & commTime_
Definition: oops/interface/ObsVector.h:86
util
Definition: ObservationL95.h:32
ObsDataVector.h