UFO
ObsBiasIncrement.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2018-2019 UCAR
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  */
7 
8 #include "ufo/ObsBiasIncrement.h"
9 
10 #include <iomanip>
11 #include <set>
12 
13 #include "eckit/mpi/Comm.h"
14 
15 #include "ioda/ObsSpace.h"
16 #include "ioda/ObsVector.h"
17 
18 #include "oops/util/IntSetParser.h"
19 #include "oops/util/Logger.h"
20 #include "oops/util/missingValues.h"
21 
22 #include "ufo/ObsBias.h"
23 
24 namespace ufo {
25 
26 // -----------------------------------------------------------------------------
27 
28 ObsBiasIncrement::ObsBiasIncrement(const ioda::ObsSpace & odb, const eckit::Configuration & conf)
29  : predbases_(0), jobs_(0), odb_(odb), conf_(conf) {
30  oops::Log::trace() << "ObsBiasIncrement::create starting." << std::endl;
31 
32  // Get the jobs(channels)
33  if (conf_.has("obs bias.jobs")) {
34  const std::set<int> jobs = oops::parseIntSet(conf_.getString("obs bias.jobs"));
35  jobs_.assign(jobs.begin(), jobs.end());
36  }
37 
38  // Predictor factory
39  if (conf_.has("obs bias.predictors")) {
40  std::vector<eckit::LocalConfiguration> confs;
41  conf_.get("obs bias.predictors", confs);
42  typedef std::unique_ptr<PredictorBase> predictor;
43  for (std::size_t j = 0; j < confs.size(); ++j) {
44  predbases_.push_back(predictor(PredictorFactory::create(confs[j], jobs_)));
45  prednames_.push_back(predbases_[j]->name());
46  }
47  }
48 
49  // initialize bias coefficient perturbations
50  biascoeffsinc_.resize(prednames_.size()*jobs_.size());
51  std::fill(biascoeffsinc_.begin(), biascoeffsinc_.end(), 0.0);
52 
53  oops::Log::trace() << "ObsBiasIncrement::create done." << std::endl;
54 }
55 
56 // -----------------------------------------------------------------------------
57 
58 ObsBiasIncrement::ObsBiasIncrement(const ObsBiasIncrement & other, const bool copy)
59  : odb_(other.odb_), conf_(other.conf_), predbases_(other.predbases_),
60  prednames_(other.prednames_), jobs_(other.jobs_) {
61  oops::Log::trace() << "ObsBiasIncrement::copy ctor starting" << std::endl;
62 
63  // initialize bias coefficient perturbations
64  biascoeffsinc_.resize(prednames_.size()*jobs_.size());
65  std::fill(biascoeffsinc_.begin(), biascoeffsinc_.end(), 0.0);
66 
67  // Copy the bias model coeff data
68  if (biascoeffsinc_.size() > 0) *this = other;
69 
70  oops::Log::trace() << "ObsBiasIncrement::copy ctor done." << std::endl;
71 }
72 
73 // -----------------------------------------------------------------------------
74 
75 ObsBiasIncrement::ObsBiasIncrement(const ObsBiasIncrement & other,
76  const eckit::Configuration & conf)
77  : odb_(other.odb_), conf_(conf), predbases_(), prednames_(), jobs_() {
78  oops::Log::trace() << "ObsBiasIncrement::copy ctor starting." << std::endl;
79  // Get the jobs(channels)
80  if (conf_.has("obs bias.jobs")) {
81  const std::set<int> jobs = oops::parseIntSet(conf_.getString("obs bias.jobs"));
82  jobs_.assign(jobs.begin(), jobs.end());
83  }
84 
85  // Predictor factory
86  if (conf_.has("obs bias.predictors")) {
87  std::vector<eckit::LocalConfiguration> confs;
88  conf_.get("obs bias.predictors", confs);
89  typedef std::unique_ptr<PredictorBase> predictor;
90  for (std::size_t j = 0; j < confs.size(); ++j) {
91  predbases_.push_back(predictor(PredictorFactory::create(confs[j], jobs_)));
92  prednames_.push_back(predbases_[j]->name());
93  }
94  }
95 
96  // initialize bias coefficient perturbations
97  biascoeffsinc_.resize(prednames_.size()*jobs_.size());
98  std::fill(biascoeffsinc_.begin(), biascoeffsinc_.end(), 0.0);
99 
100  // Copy the data
101  if (biascoeffsinc_.size() > 0) *this = other;
102 
103  oops::Log::trace() << "ObsBiasIncrement::copy ctor done." << std::endl;
104 }
105 
106 // -----------------------------------------------------------------------------
107 
108 void ObsBiasIncrement::diff(const ObsBias & b1, const ObsBias & b2) {
109  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
110  biascoeffsinc_[jj]= b1[jj] - b2[jj];
111 }
112 
113 // -----------------------------------------------------------------------------
114 
115 void ObsBiasIncrement::zero() {
116  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
117  biascoeffsinc_[jj]= 0.0;
118 }
119 
120 // -----------------------------------------------------------------------------
121 
122 ObsBiasIncrement & ObsBiasIncrement::operator=(const ObsBiasIncrement & rhs) {
123  if (rhs) {
124  predbases_.clear();
125  jobs_.clear();
126 
127  predbases_ = rhs.predbases_;
128  jobs_ = rhs.jobs_;
129  biascoeffsinc_ = rhs.biascoeffsinc_;
130  }
131  return *this;
132 }
133 
134 // -----------------------------------------------------------------------------
135 
136 ObsBiasIncrement & ObsBiasIncrement::operator+=(const ObsBiasIncrement & rhs) {
137  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
138  biascoeffsinc_[jj] += rhs[jj];
139  return *this;
140 }
141 
142 // -----------------------------------------------------------------------------
143 
144 ObsBiasIncrement & ObsBiasIncrement::operator-=(const ObsBiasIncrement & rhs) {
145  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
146  biascoeffsinc_[jj] -= rhs[jj];
147  return *this;
148 }
149 
150 // -----------------------------------------------------------------------------
151 
152 ObsBiasIncrement & ObsBiasIncrement::operator*=(const double fact) {
153  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
154  biascoeffsinc_[jj] *= fact;
155  return *this;
156 }
157 
158 // -----------------------------------------------------------------------------
159 
160 void ObsBiasIncrement::axpy(const double fact, const ObsBiasIncrement & rhs) {
161  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj)
162  biascoeffsinc_[jj] += fact * rhs[jj];
163 }
164 
165 // -----------------------------------------------------------------------------
166 
167 double ObsBiasIncrement::dot_product_with(const ObsBiasIncrement & rhs) const {
168  double zz = 0.0;
169  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj) {
170  zz += biascoeffsinc_[jj] * rhs[jj];
171  }
172  return zz;
173 }
174 
175 // -----------------------------------------------------------------------------
176 
177 double ObsBiasIncrement::norm() const {
178  double zz = 0.0;
179  for (std::size_t jj = 0; jj < biascoeffsinc_.size(); ++jj) {
180  zz += biascoeffsinc_[jj] * biascoeffsinc_[jj];
181  }
182  if (biascoeffsinc_.size() > 0) zz = std::sqrt(zz/biascoeffsinc_.size());
183  return zz;
184 }
185 
186 // -----------------------------------------------------------------------------
187 
188 void ObsBiasIncrement::computeObsBiasTL(const GeoVaLs & geovals,
189  const std::vector<ioda::ObsVector> & predData,
190  ioda::ObsVector & ybiasinc) const {
191  oops::Log::trace() << "ObsBiasIncrement::computeObsBiasTL starts." << std::endl;
192 
193  if (serialSize() > 0) {
194  const std::size_t nlocs = ybiasinc.nlocs();
195  const std::size_t npreds = prednames_.size();
196  const std::size_t njobs = jobs_.size();
197 
198  ASSERT(biascoeffsinc_.size() == npreds*njobs);
199  ASSERT(predData.size() == npreds);
200  ASSERT(ybiasinc.nvars() == njobs);
201 
202  /* predData memory layout (npreds X nlocs X njobs)
203  * Loc 0 1 2 3
204  * --------------------------
205  * pred1 Chan1 | 0 3 6 9
206  * Chan2 | 1 4 7 10
207  * .... | 2 5 8 11
208  *
209  * pred2 Chan1 |12 15 18 21
210  * Chan2 |13 16 19 22
211  * .... |14 17 20 23
212  */
213 
214  /* ybiasinc memory layout (nlocs X njobs)
215  * ch1 ch2 ch3 ch4
216  * Loc --------------------------
217  * 0 | 0 1 2 3
218  * 1 | 4 5 6 7
219  * 2 | 8 9 10 11
220  * 3 |12 13 14 15
221  * 4 |16 17 18 19
222  * ...|
223  */
224 
225  ybiasinc.zero();
226 
227  // map bias coeffs to eigen matrix (read only)
228  Eigen::Map<const Eigen::MatrixXd> coeffs(biascoeffsinc_.data(), npreds, njobs);
229 
230  // For each channel: ( nlocs X 1 ) = ( nlocs X npreds ) * ( npreds X 1 )
231  for (std::size_t jch = 0; jch < njobs; ++jch) {
232  for (std::size_t jp = 0; jp < npreds; ++jp) {
233  const double beta = coeffs(jp, jch);
234  /// axpy
235  for (std::size_t jl = 0; jl < nlocs; ++jl) {
236  ybiasinc[jl*njobs+jch] += predData[jp][jl*njobs+jch] * beta;
237  }
238  }
239  }
240  }
241 
242  oops::Log::trace() << "ObsBiasIncrement::computeObsBiasTL done." << std::endl;
243 }
244 
245 // -----------------------------------------------------------------------------
246 
247 void ObsBiasIncrement::computeObsBiasAD(GeoVaLs & geovals,
248  const std::vector<ioda::ObsVector> & predData,
249  const ioda::ObsVector & ybiasinc) {
250  oops::Log::trace() << "ObsBiasIncrement::computeObsBiasAD starts." << std::endl;
251 
252  if (this->serialSize() > 0) {
253  const std::size_t nlocs = ybiasinc.nlocs();
254  const std::size_t npreds = prednames_.size();
255  const std::size_t njobs = jobs_.size();
256 
257  ASSERT(biascoeffsinc_.size() == npreds*njobs);
258  ASSERT(predData.size() == npreds);
259  ASSERT(ybiasinc.nvars() == njobs);
260 
261  // map bias coeffs to eigen matrix (writable)
262  Eigen::Map<Eigen::MatrixXd> coeffs(biascoeffsinc_.data(), npreds, njobs);
263 
264  std::size_t indx;
265  Eigen::VectorXd tmp = Eigen::VectorXd::Zero(nlocs, 1);
266  for (std::size_t jch = 0; jch < njobs; ++jch) {
267  for (std::size_t jl = 0; jl < nlocs; ++jl) {
268  indx = jl*njobs+jch;
269  if (ybiasinc[indx] != util::missingValue(ybiasinc[indx])) {
270  tmp(jl) = ybiasinc[indx];
271  }
272  }
273  // For each channel: ( npreds X 1 ) = ( npreds X nlocs ) * ( nlocs X 1 )
274  for (std::size_t jp = 0; jp < npreds; ++jp) {
275  for (std::size_t jl = 0; jl < nlocs; ++jl) {
276  coeffs(jp, jch) += predData[jp][jl*njobs+jch] * tmp(jl);
277  }
278  }
279 
280  // zero out for next job
281  tmp.setConstant(0.0);
282  }
283 
284  // Sum across the processros
285  odb_.comm().allReduceInPlace(biascoeffsinc_.begin(), biascoeffsinc_.end(), eckit::mpi::sum());
286  }
287 
288  oops::Log::trace() << "ObsBiasIncrement::computeAD done." << std::endl;
289 }
290 
291 // -----------------------------------------------------------------------------
292 
293 void ObsBiasIncrement::print(std::ostream & os) const {
294  if (this->serialSize() > 0) {
295  // map bias coeffs to eigen matrix (writable)
296  Eigen::Map<const Eigen::MatrixXd>
297  coeffs(biascoeffsinc_.data(), prednames_.size(), jobs_.size());
298 
299  os << "ObsBiasIncrement::print " << std::endl;
300  os << "---------------------------------------------------------------" << std::endl;
301  auto njobs = jobs_.size();
302  for (std::size_t p = 0; p < prednames_.size(); ++p) {
303  os << std::fixed << std::setw(20) << prednames_[p]
304  << ": Min= " << std::setw(15) << std::setprecision(8)
305  << coeffs.row(p).minCoeff()
306  << ", Max= " << std::setw(15) << std::setprecision(8)
307  << coeffs.row(p).maxCoeff()
308  << ", Norm= " << std::setw(15) << std::setprecision(8)
309  << coeffs.row(p).norm()
310  << std::endl;
311  }
312  os << "---------------------------------------------------------------" << std::endl;
313  os << "ObsBiasIncrement::print done" << std::endl;
314  }
315 }
316 
317 // -----------------------------------------------------------------------------
318 
319 } // namespace ufo
ufo::ObsBiasIncrement::biascoeffsinc_
std::vector< double > biascoeffsinc_
Definition: ObsBiasIncrement.h:88
ufo::ObsBiasIncrement::odb_
const ioda::ObsSpace & odb_
Definition: ObsBiasIncrement.h:85
ufo::ObsBiasIncrement::norm
double norm() const
Definition: ObsBiasIncrement.cc:177
ufo::ObsBiasIncrement::ObsBiasIncrement
ObsBiasIncrement(const ioda::ObsSpace &, const eckit::Configuration &)
Definition: ObsBiasIncrement.cc:28
ObsBiasIncrement.h
ObsBias.h
ufo::ObsBiasIncrement::operator+=
ObsBiasIncrement & operator+=(const ObsBiasIncrement &)
Definition: ObsBiasIncrement.cc:136
ufo::ObsBiasIncrement::operator*=
ObsBiasIncrement & operator*=(const double)
Definition: ObsBiasIncrement.cc:152
ufo::ObsBiasIncrement::print
void print(std::ostream &) const
Definition: ObsBiasIncrement.cc:293
ufo::ObsBiasIncrement::diff
void diff(const ObsBias &, const ObsBias &)
Definition: ObsBiasIncrement.cc:108
ufo::ObsBiasIncrement::axpy
void axpy(const double, const ObsBiasIncrement &)
Definition: ObsBiasIncrement.cc:160
ufo::ObsBiasIncrement::computeObsBiasTL
void computeObsBiasTL(const GeoVaLs &, const std::vector< ioda::ObsVector > &, ioda::ObsVector &) const
Definition: ObsBiasIncrement.cc:188
ufo
Definition: RunCRTM.h:27
ufo::ObsBiasIncrement::predbases_
std::vector< std::shared_ptr< PredictorBase > > predbases_
Definition: ObsBiasIncrement.h:89
ufo::ObsBias
Class to handle observation bias parameters.
Definition: ObsBias.h:44
ufo::ObsBiasIncrement::jobs_
std::vector< int > jobs_
Definition: ObsBiasIncrement.h:91
ufo::ObsBiasIncrement::operator=
ObsBiasIncrement & operator=(const ObsBiasIncrement &)
Definition: ObsBiasIncrement.cc:122
ufo::ObsBiasIncrement::conf_
const eckit::LocalConfiguration conf_
Definition: ObsBiasIncrement.h:86
ufo::ObsBiasIncrement::operator-=
ObsBiasIncrement & operator-=(const ObsBiasIncrement &)
Definition: ObsBiasIncrement.cc:144
ufo::GeoVaLs
GeoVaLs: geophysical values at locations.
Definition: src/ufo/GeoVaLs.h:39
ufo::ObsBiasIncrement
Definition: ObsBiasIncrement.h:39
ufo::ObsBiasIncrement::serialSize
std::size_t serialSize() const
Definition: ObsBiasIncrement.h:75
ufo::PredictorFactory::create
static PredictorBase * create(const eckit::Configuration &, const std::vector< int > &)
Definition: PredictorBase.cc:39
ufo::ObsBiasIncrement::computeObsBiasAD
void computeObsBiasAD(GeoVaLs &, const std::vector< ioda::ObsVector > &, const ioda::ObsVector &)
Definition: ObsBiasIncrement.cc:247
ufo::ObsBiasIncrement::prednames_
std::vector< std::string > prednames_
Definition: ObsBiasIncrement.h:90
ufo::ObsBiasIncrement::dot_product_with
double dot_product_with(const ObsBiasIncrement &) const
Definition: ObsBiasIncrement.cc:167
ufo::ObsBiasIncrement::zero
void zero()
Definition: ObsBiasIncrement.cc:115
conf
Definition: conf.py:1