OOPS
CostJbTotal.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_ASSIMILATION_COSTJBTOTAL_H_
12 #define OOPS_ASSIMILATION_COSTJBTOTAL_H_
13 
14 #include <limits>
15 #include <memory>
16 
17 #include "eckit/config/LocalConfiguration.h"
18 #include "oops/assimilation/ControlIncrement.h"
19 #include "oops/assimilation/ControlVariable.h"
20 #include "oops/assimilation/CostJbState.h"
21 #include "oops/base/ObsAuxCovariances.h"
22 #include "oops/base/ObsSpaces.h"
23 #include "oops/interface/Geometry.h"
24 #include "oops/interface/ModelAuxCovariance.h"
25 #include "oops/interface/State.h"
26 #include "oops/util/Logger.h"
27 
28 namespace oops {
29  template<typename MODEL> class JqTermTLAD;
30 
31 // -----------------------------------------------------------------------------
32 
33 /// Total Jb cost function for all components of the control variable.
34 
35 template<typename MODEL, typename OBS> class CostJbTotal {
36  typedef ControlIncrement<MODEL, OBS> CtrlInc_;
37  typedef ControlVariable<MODEL, OBS> CtrlVar_;
38  typedef State<MODEL> State_;
39  typedef CostJbState<MODEL> JbState_;
40  typedef JqTermTLAD<MODEL> JqTermTLAD_;
41  typedef Geometry<MODEL> Geometry_;
42  typedef ModelAuxCovariance<MODEL> ModelAuxCovar_;
43  typedef ObsAuxCovariances<OBS> ObsAuxCovars_;
44  typedef ObsSpaces<OBS> ObsSpaces_;
45 
46  public:
47 /// Construct \f$ J_b\f$.
48  CostJbTotal(const CtrlVar_ &, JbState_ *, const eckit::Configuration &,
49  const Geometry_ &, const ObsSpaces_ & odb);
50 
51 /// Destructor
52  ~CostJbTotal() {}
53 
54 /// Initialize before nonlinear model integration.
55  void initialize(const CtrlVar_ &) const;
56  JqTermTLAD_ * initializeTraj(const CtrlVar_ &, const Geometry_ &,
57  const eckit::Configuration &);
58 
59 /// Finalize computation after nonlinear model integration.
60  double finalize(const CtrlVar_ &) const;
61  void finalizeTraj(JqTermTLAD_ *);
62 
63 /// Initialize before starting the TL run.
64  JqTermTLAD_ * initializeTL() const;
65  void finalizeTL(JqTermTLAD_ *, const CtrlInc_ &, CtrlInc_ &) const;
66 
67 /// Initialize before starting the AD run.
68  JqTermTLAD_ * initializeAD(CtrlInc_ &, const CtrlInc_ &) const;
69  void finalizeAD(JqTermTLAD_ *) const;
70 
71 /// Multiply by covariance matrix and its inverse.
72  void multiplyB(const CtrlInc_ &, CtrlInc_ &) const;
73  void multiplyBinv(const CtrlInc_ &, CtrlInc_ &) const;
74 
75 /// Randomize
76  void randomize(CtrlInc_ &) const;
77 
78 /// Add Jb gradient at first guess.
79  void addGradientFG(CtrlInc_ &) const;
80  void addGradientFG(CtrlInc_ &, CtrlInc_ &) const;
81 
82 /// Return background.
83  const CtrlVar_ & getBackground() const {return xb_;}
84 
85 /// Return first guess \f$ x_0-x_b\f$.
86  const CtrlInc_ & getFirstGuess() const {return *dxFG_;}
87 
88 /// Jb terms for ControlIncrement constructor.
89  const Geometry_ & resolution() const {return *resol_;}
90  const JbState_ & jbState() const {return *jb_;}
91  const ModelAuxCovar_ & jbModBias() const {return jbModBias_;}
92  const ObsAuxCovars_ & jbObsBias() const {return jbObsBias_;}
93  const util::DateTime & windowBegin() const {return windowBegin_;}
94  const util::DateTime & windowEnd() const {return windowEnd_;}
95 
96  private:
97  double evaluate(const CtrlInc_ &) const;
98 
99  const CtrlVar_ & xb_;
100  std::unique_ptr<JbState_> jb_;
101  ModelAuxCovar_ jbModBias_;
102  ObsAuxCovars_ jbObsBias_;
103 
104  const CtrlVar_ mutable * fg_;
105 
106 /// First guess increment \f$x_0-x_b\f$ or more generally \f$ x_i-M(x_{i-1})\f$.
107  std::unique_ptr<CtrlInc_> dxFG_;
108 
109 /// Inner loop resolution
110  std::unique_ptr<Geometry_> resol_;
111  const util::DateTime windowBegin_;
112  const util::DateTime windowEnd_;
113 };
114 
115 // =============================================================================
116 
117 template<typename MODEL, typename OBS>
118 CostJbTotal<MODEL, OBS>::CostJbTotal(const CtrlVar_ & xb, JbState_ * jb,
119  const eckit::Configuration & conf,
120  const Geometry_ & resol, const ObsSpaces_ & odb)
121  : xb_(xb), jb_(jb),
122  jbModBias_(conf.getSubConfiguration("model aux error"), resol),
123  jbObsBias_(odb, conf), dxFG_(), resol_(),
124  windowBegin_(conf.getString("window begin")),
125  windowEnd_(windowBegin_ + util::Duration(conf.getString("window length")))
126 {
127  Log::trace() << "CostJbTotal contructed." << std::endl;
128 }
129 
130 // -----------------------------------------------------------------------------
131 
132 template<typename MODEL, typename OBS>
133 void CostJbTotal<MODEL, OBS>::initialize(const CtrlVar_ & fg) const {
134  Log::trace() << "CostJbTotal::initialize start" << std::endl;
135  fg_ = &fg;
136  Log::trace() << "CostJbTotal::initialize done" << std::endl;
137 }
138 
139 // -----------------------------------------------------------------------------
140 
141 template<typename MODEL, typename OBS>
142 double CostJbTotal<MODEL, OBS>::finalize(const CtrlVar_ & mx) const {
143  Log::trace() << "CostJbTotal::finalize start" << std::endl;
144  ASSERT(fg_);
145  CtrlInc_ dx(*this);
146 
147 // Compute x_0 - x_b for Jb (and Jq is present)
148  jb_->computeIncrement(xb_.state(), fg_->state(), mx.state(), dx.state());
149 
150 // Model and Obs biases
151  dx.modVar().diff(fg_->modVar(), xb_.modVar());
152  dx.obsVar().diff(fg_->obsVar(), xb_.obsVar());
153 
154 // Print increment
155  Log::info() << "CostJb: FG-BG" << dx << std::endl;
156 
157 // Compute Jb value
158  double zjb = this->evaluate(dx);
159  Log::trace() << "CostJbTotal::finalize done" << std::endl;
160  return zjb;
161 }
162 
163 // -----------------------------------------------------------------------------
164 
165 template<typename MODEL, typename OBS>
166 JqTermTLAD<MODEL> * CostJbTotal<MODEL, OBS>::initializeTraj(const CtrlVar_ & fg,
167  const Geometry_ & resol,
168  const eckit::Configuration & inner) {
169  Log::trace() << "CostJbTotal::initializeTraj start" << std::endl;
170  fg_ = &fg;
171  resol_.reset(new Geometry_(resol));
172 // Linearize terms
173  jb_->linearize(fg.state(), *resol_);
174  jbModBias_.linearize(fg.modVar(), *resol_);
175  jbObsBias_.linearize(fg.obsVar(), inner);
176 
177  JqTermTLAD_ * jqlin = jb_->initializeJqTLAD();
178 
179  Log::trace() << "CostJbTotal::initializeTraj done" << std::endl;
180  return jqlin;
181 }
182 
183 // -----------------------------------------------------------------------------
184 
185 template<typename MODEL, typename OBS>
186 void CostJbTotal<MODEL, OBS>::finalizeTraj(JqTermTLAD_ * jqlin) {
187  Log::trace() << "CostJbTotal::finalizeTraj start" << std::endl;
188  ASSERT(fg_);
189 // Compute and save first guess increment.
190  dxFG_.reset(new CtrlInc_(*this));
191 
192 // Compute x_0 - x_b for Jb (and Jq is present)
193  const State_ * mx = &fg_->state();
194  if (jqlin) mx = &jqlin->getMxi();
195  jb_->computeIncrement(xb_.state(), fg_->state(), *mx, dxFG_->state());
196 
197 // Model and Obs biases
198  dxFG_->modVar().diff(fg_->modVar(), xb_.modVar());
199  dxFG_->obsVar().diff(fg_->obsVar(), xb_.obsVar());
200 
201 // Print increment
202  Log::info() << "CostJb: FG-BG" << *dxFG_ << std::endl;
203  Log::trace() << "CostJbTotal::finalizeTraj done" << std::endl;
204 }
205 
206 // -----------------------------------------------------------------------------
207 
208 template<typename MODEL, typename OBS>
209 double CostJbTotal<MODEL, OBS>::evaluate(const CtrlInc_ & dx) const {
210  Log::trace() << "CostJbTotal::evaluate start" << std::endl;
211  CtrlInc_ gg(*this);
212  this->multiplyBinv(dx, gg);
213 
214  double zjb = 0.0;
215  double zz = 0.5 * dot_product(dx.state(), gg.state());
216  Log::info() << "CostJb : Nonlinear Jb State = " << zz << std::endl;
217  zjb += zz;
218  zz = 0.5 * dot_product(dx.modVar(), gg.modVar());
219  Log::info() << "CostJb : Nonlinear Jb Model Aux = " << zz << std::endl;
220  zjb += zz;
221  zz = 0.5 * dot_product(dx.obsVar(), gg.obsVar());
222  Log::info() << "CostJb : Nonlinear Jb Obs Aux = " << zz << std::endl;
223  zjb += zz;
224 
225  Log::info() << "CostJb : Nonlinear Jb = " << zjb << std::endl;
226 
227 // Get rid of very small values for test
228  double ztest = zjb;
229  if (zjb >= 0.0 && zjb <= std::numeric_limits<double>::epsilon()) ztest = 0.0;
230  Log::test() << "CostJb : Nonlinear Jb = " << ztest << std::endl;
231 
232  Log::trace() << "CostJbTotal::evaluate done" << std::endl;
233  return zjb;
234 }
235 
236 // -----------------------------------------------------------------------------
237 
238 template<typename MODEL, typename OBS>
239 void CostJbTotal<MODEL, OBS>::addGradientFG(CtrlInc_ & grad) const {
240  Log::trace() << "CostJbTotal::addGradientFG 1 start" << std::endl;
241  CtrlInc_ gg(grad, false);
242  this->multiplyBinv(*dxFG_, gg);
243  grad += gg;
244  Log::trace() << "CostJbTotal::addGradientFG 1 done" << std::endl;
245 }
246 
247 // -----------------------------------------------------------------------------
248 
249 template<typename MODEL, typename OBS>
250 void CostJbTotal<MODEL, OBS>::addGradientFG(CtrlInc_ & grad, CtrlInc_ & gradJb) const {
251  Log::trace() << "CostJbTotal::addGradientFG 2 start" << std::endl;
252  jb_->addGradient(dxFG_->state(), grad.state(), gradJb.state());
253  grad.modVar() += gradJb.modVar();
254  grad.obsVar() += gradJb.obsVar();
255  Log::trace() << "CostJbTotal::addGradientFG 2 done" << std::endl;
256 }
257 
258 // -----------------------------------------------------------------------------
259 
260 template<typename MODEL, typename OBS>
261 JqTermTLAD<MODEL> * CostJbTotal<MODEL, OBS>::initializeTL() const {
262  Log::trace() << "CostJbTotal::initializeTL start" << std::endl;
263  JqTermTLAD_ * jqtl = jb_->initializeJqTL();
264  Log::trace() << "CostJbTotal::initializeTL done" << std::endl;
265  return jqtl;
266 }
267 
268 // -----------------------------------------------------------------------------
269 
270 template<typename MODEL, typename OBS>
271 void CostJbTotal<MODEL, OBS>::finalizeTL(JqTermTLAD_ * jqtl, const CtrlInc_ & bgns,
272  CtrlInc_ & dx) const {
273  Log::trace() << "CostJbTotal::finalizeTL start" << std::endl;
274  dx = bgns;
275  if (jqtl) jqtl->computeModelErrorTL(dx.state());
276  Log::trace() << "CostJbTotal::finalizeTL done" << std::endl;
277 }
278 
279 // -----------------------------------------------------------------------------
280 
281 template<typename MODEL, typename OBS>
282 JqTermTLAD<MODEL> * CostJbTotal<MODEL, OBS>::initializeAD(CtrlInc_ & bgns,
283  const CtrlInc_ & dx) const {
284  Log::trace() << "CostJbTotal::initializeAD start" << std::endl;
285  JqTermTLAD_ * jqad = jb_->initializeJqAD(dx.state());
286  bgns += dx;
287  Log::trace() << "CostJbTotal::initializeAD done" << std::endl;
288  return jqad;
289 }
290 
291 // -----------------------------------------------------------------------------
292 
293 template<typename MODEL, typename OBS>
294 void CostJbTotal<MODEL, OBS>::finalizeAD(JqTermTLAD_ * jqad) const {
295  Log::trace() << "CostJbTotal::finalizeAD start" << std::endl;
296  if (jqad) jqad->clear();
297  Log::trace() << "CostJbTotal::finalizeAD done" << std::endl;
298 }
299 
300 // -----------------------------------------------------------------------------
301 
302 template<typename MODEL, typename OBS>
303 void CostJbTotal<MODEL, OBS>::multiplyB(const CtrlInc_ & dxin, CtrlInc_ & dxout) const {
304  Log::trace() << "CostJbTotal::multiplyB start" << std::endl;
305  jb_->Bmult(dxin.state(), dxout.state());
306  jbModBias_.multiply(dxin.modVar(), dxout.modVar());
307  jbObsBias_.multiply(dxin.obsVar(), dxout.obsVar());
308  Log::trace() << "CostJbTotal::multiplyB done" << std::endl;
309 }
310 
311 // -----------------------------------------------------------------------------
312 
313 template<typename MODEL, typename OBS>
314 void CostJbTotal<MODEL, OBS>::multiplyBinv(const CtrlInc_ & dxin, CtrlInc_ & dxout) const {
315  Log::trace() << "CostJbTotal::multiplyBinv start" << std::endl;
316  jb_->Bminv(dxin.state(), dxout.state());
317  jbModBias_.inverseMultiply(dxin.modVar(), dxout.modVar());
318  jbObsBias_.inverseMultiply(dxin.obsVar(), dxout.obsVar());
319  Log::trace() << "CostJbTotal::multiplyBinv done" << std::endl;
320 }
321 
322 // -----------------------------------------------------------------------------
323 
324 template<typename MODEL, typename OBS>
325 void CostJbTotal<MODEL, OBS>::randomize(CtrlInc_ & dx) const {
326  Log::trace() << "CostJbTotal::randomize start" << std::endl;
327  jb_->randomize(dx.state());
328  jbModBias_.randomize(dx.modVar());
329  jbObsBias_.randomize(dx.obsVar());
330  Log::trace() << "CostJbTotal::randomize done" << std::endl;
331 }
332 
333 // -----------------------------------------------------------------------------
334 
335 } // namespace oops
336 
337 #endif // OOPS_ASSIMILATION_COSTJBTOTAL_H_
oops::CostJbTotal::multiplyB
void multiplyB(const CtrlInc_ &, CtrlInc_ &) const
Multiply by covariance matrix and its inverse.
Definition: CostJbTotal.h:303
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
oops::CostJbTotal::multiplyBinv
void multiplyBinv(const CtrlInc_ &, CtrlInc_ &) const
Definition: CostJbTotal.h:314
oops::ControlVariable::obsVar
ObsAuxCtrls_ & obsVar()
Get augmented observation control variable.
Definition: ControlVariable.h:77
oops::CostJbTotal::JbState_
CostJbState< MODEL > JbState_
Definition: CostJbTotal.h:39
oops::ControlIncrement::obsVar
ObsAuxIncrs_ & obsVar()
Get augmented observation control variable.
Definition: ControlIncrement.h:94
oops::CostJbTotal::ObsSpaces_
ObsSpaces< OBS > ObsSpaces_
Definition: CostJbTotal.h:44
oops::CostJbTotal::jbState
const JbState_ & jbState() const
Definition: CostJbTotal.h:90
oops::CostJbTotal::finalizeTL
void finalizeTL(JqTermTLAD_ *, const CtrlInc_ &, CtrlInc_ &) const
Definition: CostJbTotal.h:271
oops::ControlVariable
Control variable.
Definition: ControlVariable.h:41
ObsSpaces.h
oops::CostJbTotal::CtrlInc_
ControlIncrement< MODEL, OBS > CtrlInc_
Definition: CostJbTotal.h:36
oops::ControlIncrement::modVar
ModelAuxIncr_ & modVar()
Get augmented model control variable.
Definition: ControlIncrement.h:90
oops::CostJbTotal::fg_
const CtrlVar_ * fg_
Definition: CostJbTotal.h:104
oops::CostJbTotal::initialize
void initialize(const CtrlVar_ &) const
Initialize before nonlinear model integration.
Definition: CostJbTotal.h:133
oops::CostJbTotal::resol_
std::unique_ptr< Geometry_ > resol_
Inner loop resolution.
Definition: CostJbTotal.h:110
oops::JqTermTLAD::computeModelErrorTL
void computeModelErrorTL(Increment_ &)
Definition: JqTermTLAD.h:141
oops::CostJbState
Jb Cost Function Base Class.
Definition: CostJbState.h:39
oops::ControlIncrement
Definition: ControlIncrement.h:46
oops::JqTermTLAD
Definition: CostJb3D.h:29
ModelAuxCovariance.h
oops::ControlVariable::modVar
ModelAux_ & modVar()
Get augmented model control variable.
Definition: ControlVariable.h:73
oops::CostJbTotal::finalizeTraj
void finalizeTraj(JqTermTLAD_ *)
Definition: CostJbTotal.h:186
oops::CostJbTotal::jb_
std::unique_ptr< JbState_ > jb_
Definition: CostJbTotal.h:100
oops::ModelAuxIncrement::diff
void diff(const ModelAuxControl_ &, const ModelAuxControl_ &)
Linear algebra operators.
Definition: oops/interface/ModelAuxIncrement.h:135
oops::CostJbTotal::windowBegin_
const util::DateTime windowBegin_
Definition: CostJbTotal.h:111
oops::JqTermTLAD::getMxi
State_ & getMxi() const
Definition: JqTermTLAD.h:122
oops::ObsAuxIncrements::diff
void diff(const ObsAuxControls_ &, const ObsAuxControls_ &)
Linear algebra operators.
Definition: ObsAuxIncrements.h:133
oops::State::state
State_ & state()
Interfacing.
Definition: oops/interface/State.h:56
oops::CostJbTotal::windowEnd_
const util::DateTime windowEnd_
Definition: CostJbTotal.h:112
oops::CostJbTotal::resolution
const Geometry_ & resolution() const
Jb terms for ControlIncrement constructor.
Definition: CostJbTotal.h:89
oops::CostJbTotal::jbModBias
const ModelAuxCovar_ & jbModBias() const
Definition: CostJbTotal.h:91
oops::JqTermTLAD::clear
void clear()
Definition: JqTermTLAD.h:37
oops::CostJbTotal::windowBegin
const util::DateTime & windowBegin() const
Definition: CostJbTotal.h:93
oops::CostJbTotal::randomize
void randomize(CtrlInc_ &) const
Randomize.
Definition: CostJbTotal.h:325
CostJbState.h
oops::CostJbTotal::xb_
const CtrlVar_ & xb_
Definition: CostJbTotal.h:99
oops::CostJbTotal::jbObsBias_
ObsAuxCovars_ jbObsBias_
Definition: CostJbTotal.h:102
oops::Geometry
Geometry class used in oops; subclass of interface class above.
Definition: oops/interface/Geometry.h:189
oops::CostJbTotal::ObsAuxCovars_
ObsAuxCovariances< OBS > ObsAuxCovars_
Definition: CostJbTotal.h:43
oops::CostJbTotal::State_
State< MODEL > State_
Definition: CostJbTotal.h:38
oops::ControlIncrement::state
Increment_ & state()
Get state control variable.
Definition: ControlIncrement.h:86
oops::CostJbTotal::evaluate
double evaluate(const CtrlInc_ &) const
Definition: CostJbTotal.h:209
oops::CostJbTotal::jbObsBias
const ObsAuxCovars_ & jbObsBias() const
Definition: CostJbTotal.h:92
oops::CostJbTotal::CtrlVar_
ControlVariable< MODEL, OBS > CtrlVar_
Definition: CostJbTotal.h:37
oops::CostJbTotal::initializeAD
JqTermTLAD_ * initializeAD(CtrlInc_ &, const CtrlInc_ &) const
Initialize before starting the AD run.
Definition: CostJbTotal.h:282
oops::CostJbTotal::addGradientFG
void addGradientFG(CtrlInc_ &) const
Add Jb gradient at first guess.
Definition: CostJbTotal.h:239
oops::ModelAuxCovariance
Definition: oops/interface/ModelAuxCovariance.h:36
oops::CostJbTotal::JqTermTLAD_
JqTermTLAD< MODEL > JqTermTLAD_
Definition: CostJbTotal.h:40
oops::CostJbTotal::jbModBias_
ModelAuxCovar_ jbModBias_
Definition: CostJbTotal.h:101
oops::CostJbTotal::Geometry_
Geometry< MODEL > Geometry_
Definition: CostJbTotal.h:41
ObsAuxCovariances.h
oops::State
Encapsulates the model state.
Definition: CostJbState.h:28
oops::CostJbTotal::initializeTraj
JqTermTLAD_ * initializeTraj(const CtrlVar_ &, const Geometry_ &, const eckit::Configuration &)
Definition: CostJbTotal.h:166
oops::CostJbTotal::getFirstGuess
const CtrlInc_ & getFirstGuess() const
Return first guess .
Definition: CostJbTotal.h:86
State.h
ControlIncrement.h
oops::ObsSpaces
Definition: ObsSpaces.h:41
oops::CostJbTotal::windowEnd
const util::DateTime & windowEnd() const
Definition: CostJbTotal.h:94
oops::CostJbTotal::finalize
double finalize(const CtrlVar_ &) const
Finalize computation after nonlinear model integration.
Definition: CostJbTotal.h:142
oops::CostJbTotal::finalizeAD
void finalizeAD(JqTermTLAD_ *) const
Definition: CostJbTotal.h:294
oops::CostJbTotal::dxFG_
std::unique_ptr< CtrlInc_ > dxFG_
First guess increment or more generally .
Definition: CostJbTotal.h:107
oops::CostJbTotal::CostJbTotal
CostJbTotal(const CtrlVar_ &, JbState_ *, const eckit::Configuration &, const Geometry_ &, const ObsSpaces_ &odb)
Construct .
Definition: CostJbTotal.h:118
util
Definition: ObservationL95.h:32
ControlVariable.h
oops::CostJbTotal::~CostJbTotal
~CostJbTotal()
Destructor.
Definition: CostJbTotal.h:52
oops::CostJbTotal::ModelAuxCovar_
ModelAuxCovariance< MODEL > ModelAuxCovar_
Definition: CostJbTotal.h:42
oops::CostJbTotal::getBackground
const CtrlVar_ & getBackground() const
Return background.
Definition: CostJbTotal.h:83
oops::ControlVariable::state
State_ & state()
Get state control variable.
Definition: ControlVariable.h:69
Geometry.h
oops::CostJbTotal::initializeTL
JqTermTLAD_ * initializeTL() const
Initialize before starting the TL run.
Definition: CostJbTotal.h:261
oops::ObsAuxCovariances
Definition: ObsAuxCovariances.h:32