OOPS
Minimizer.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_MINIMIZER_H_
12 #define OOPS_ASSIMILATION_MINIMIZER_H_
13 
14 #include <map>
15 #include <string>
16 #include <boost/noncopyable.hpp>
17 
18 #include "eckit/config/Configuration.h"
19 #include "oops/assimilation/ControlIncrement.h"
20 #include "oops/assimilation/ControlVariable.h"
21 #include "oops/assimilation/CostFunction.h"
22 #include "oops/assimilation/DualVector.h"
23 #include "oops/assimilation/HMatrix.h"
24 #include "oops/assimilation/HtMatrix.h"
25 #include "oops/interface/Increment.h"
26 #include "oops/interface/State.h"
27 #include "oops/util/dot_product.h"
28 #include "oops/util/formats.h"
29 #include "oops/util/Logger.h"
30 #include "oops/util/PrintAdjTest.h"
31 
32 namespace oops {
33 
34 // -----------------------------------------------------------------------------
35 
36 /// A Minimizer knows how to minimize a cost function
37 template<typename MODEL, typename OBS> class Minimizer : private boost::noncopyable {
38  typedef CostFunction<MODEL, OBS> CostFct_;
39  typedef ControlIncrement<MODEL, OBS> CtrlInc_;
40  typedef DualVector<MODEL, OBS> Dual_;
41  typedef HMatrix<MODEL, OBS> H_;
42  typedef HtMatrix<MODEL, OBS> Ht_;
43  typedef State<MODEL> State_;
44 
45  public:
46  explicit Minimizer(const CostFct_ & J): J_(J), outerIteration_(0) {}
47  virtual ~Minimizer() {}
48  ControlIncrement<MODEL, OBS> * minimize(const eckit::Configuration &);
49  virtual const std::string classname() const = 0;
50 
51  private:
52  virtual ControlIncrement<MODEL, OBS> * doMinimize(const eckit::Configuration &) = 0;
53 
54  void adjTests(const eckit::Configuration &);
55  void adjModelTest(const Ht_ &, const H_ &);
56  void adjObsTest(const Ht_ &, const H_ &);
57 
58  void tlmTests(const eckit::Configuration &);
59  void tlmApproxTest(const H_ &);
60  void tlmTaylorTest(const H_ &);
61 
62  void tlmPropagTest(const eckit::Configuration & config, const CtrlInc_ &);
63 
64  const CostFct_ & J_;
65  int outerIteration_;
66 };
67 
68 // -----------------------------------------------------------------------------
69 
70 template<typename MODEL, typename OBS>
71 ControlIncrement<MODEL, OBS> *
72 Minimizer<MODEL, OBS>::minimize(const eckit::Configuration & config) {
73  // TLM tests
74  this->tlmTests(config);
75 
76  // ADJ tests
77  this->adjTests(config);
78 
79  // Minimize
80  ControlIncrement<MODEL, OBS> * dx = this->doMinimize(config);
81 
82  // TLM propagation test
83  this->tlmPropagTest(config, *dx);
84 
85  // Update outer loop counter
86  outerIteration_++;
87 
88  return dx;
89 }
90 
91 // -----------------------------------------------------------------------------
92 
93 template<typename MODEL, typename OBS>
94 void Minimizer<MODEL, OBS>::tlmTests(const eckit::Configuration & config) {
95 // Tangent Linear tests
96 
97  if (config.has("online diagnostics")) {
98  const eckit::LocalConfiguration onlineDiag(config, "online diagnostics");
99  bool runTlmTaylorTest = onlineDiag.getBool("tlm taylor test", false);
100  bool runTlmApproxTest = onlineDiag.getBool("tlm approx test", false);
101 
102  const H_ H(J_);
103  const Ht_ Ht(J_);
104 
105  // Online tests
106  // TLM linear approximation test
107  if (outerIteration_ == 0 && runTlmApproxTest) this->tlmApproxTest(H);
108 
109  // TLM Taylor test
110  if (outerIteration_ == 0 && runTlmTaylorTest) this->tlmTaylorTest(H);
111  }
112 }
113 
114 // -----------------------------------------------------------------------------
115 
116 template<typename MODEL, typename OBS>
117 void Minimizer<MODEL, OBS>::adjTests(const eckit::Configuration & config) {
118 // Adjoint tests
119 
120  if (config.has("online diagnostics")) {
121  const eckit::LocalConfiguration onlineDiag(config, "online diagnostics");
122  bool runAdjTlmTest = onlineDiag.getBool("adj tlm test", false);
123  bool runAdjObsTest = onlineDiag.getBool("adj obs test", false);
124 
125  const H_ H(J_);
126  const Ht_ Ht(J_);
127 
128  // Online tests
129  // Model adjoint test
130  if (runAdjTlmTest) this->adjModelTest(Ht, H);
131 
132  // Obs adjoint test
133  if (runAdjObsTest) this->adjObsTest(Ht, H);
134  }
135 }
136 
137 // -----------------------------------------------------------------------------
138 
139 template<typename MODEL, typename OBS>
140 void Minimizer<MODEL, OBS>::tlmApproxTest(const H_ & H) {
141 /* TL approx test:
142  calculate and store:
143  M(x + dx), M(x), M'(dx)
144  where dx is initialized using randomization
145 
146  TO DO:
147  write to file depending on requriements */
148 
149  Log::info() << std::endl
150  << "TLM Linear Approximation Test - starting: " << outerIteration_
151  << std::endl << std::endl;
152 
153 // initialize incremets
154  CtrlInc_ dx(J_.jb());
155 
156 // randomize increments
157  J_.jb().randomize(dx);
158 
159 // run NL for unperturbed initial condition
160  PostProcessor<State_> pp;
161  ControlVariable<MODEL, OBS> mxx(J_.jb().getBackground());
162  J_.runNL(mxx, pp);
163 
164 // run NL for perturbed initial condition
165  ControlVariable<MODEL, OBS> mpertxx(J_.jb().getBackground());
166  J_.addIncrement(mpertxx, dx);
167  J_.runNL(mpertxx, pp);
168 
169 // run TL for the perturbation
170  Dual_ dummy;
171  CtrlInc_ mdx(dx);
172  H.multiply(mdx, dummy);
173 
174 // print log
175  Log::info() << std::endl << "TLM Linear Approximation Test: done." << std::endl;
176 }
177 
178 // -----------------------------------------------------------------------------
179 
180 template<typename MODEL, typename OBS>
181 void Minimizer<MODEL, OBS>::tlmPropagTest(const eckit::Configuration & config,
182  const CtrlInc_ & dx) {
183 /* TL propagation test:
184  calculate and store:
185  M'(dx)
186  where dx comes from minimization stopped at requested iteration
187 
188  TO DO:
189  write to file depending on requriements */
190 
191  if (config.has("online diagnostics")) {
192  const eckit::LocalConfiguration onlineDiag(config, "online diagnostics");
193  bool runTlmPropagTest = onlineDiag.getBool("tlm propagation test", false);
194 
195  if (runTlmPropagTest) {
196  // print log
197  Log::info() << "TLM Propagation Test - starting: " << outerIteration_
198  << std::endl << std::endl;
199 
200  // construct propagation matrix
201  const H_ H(J_);
202 
203  // run TL for the input perturbation: M'(dx)
204  Dual_ dummy;
205  CtrlInc_ mdx(dx);
206  H.multiply(mdx, dummy);
207 
208  // print log
209  Log::info() << std::endl << "TLM Propagation Test: done." << std::endl;
210  }
211  }
212 }
213 
214 // -----------------------------------------------------------------------------
215 
216 template<typename MODEL, typename OBS>
217 void Minimizer<MODEL, OBS>::tlmTaylorTest(const H_ & H) {
218 /* TL Taylor test:
219  ||M(x + p dx) - M(x)|| / ||M'(p dx)|| -> 1
220  where p is a scalar damping factor and dx is initialized
221  using randomization */
222 
223  Log::info() << std::endl
224  << "TLM Taylor Test: " << outerIteration_
225  << std::endl;
226 
227 // initialize incremets
228  CtrlInc_ dx(J_.jb());
229 
230 // randomize increments
231  J_.jb().randomize(dx);
232 
233 // run NL for unperturbed initial condition: M(x)
234  PostProcessor<State_> pp;
235  ControlVariable<MODEL, OBS> mxx(J_.jb().getBackground());
236  J_.runNL(mxx, pp);
237 
238 // run TL for the un-damped perturbation: M'(dx)
239  Dual_ dummy;
240  CtrlInc_ mdx(dx);
241  H.multiply(mdx, dummy);
242 
243 // loop over decreasing increments
244  for (unsigned int jj = 0; jj < 14; ++jj) {
245  // ||p M'(dx)||
246  CtrlInc_ pmdx(mdx);
247  pmdx *= 1./pow(10.0, jj);
248  double denom = sqrt(dot_product(pmdx, pmdx));
249 
250  // run perturbed NL: M(x+pdx)
251  ControlVariable<MODEL, OBS> mpertxx(J_.jb().getBackground());
252  CtrlInc_ pdx(dx);
253  pdx *= 1./pow(10.0, jj);
254  J_.addIncrement(mpertxx, pdx);
255  J_.runNL(mpertxx, pp);
256 
257  // ||M(x+pdx) - M(x)||
258  CtrlInc_ diff_nl(mdx, false);
259  diff_nl.state().diff(mpertxx.state(), mxx.state());
260  double nom = sqrt(dot_product(diff_nl, diff_nl));
261 
262  // print results
263  Log::info() << std::endl
264  << "TLM Taylor test: p = " << std::setw(8) << 1./pow(10.0, jj)
265  << ", ||M(x) - M(x+p dx)|| / ||p M'(dx)|| = "
266  << util::full_precision(1.+std::abs(1.-nom/denom))
267  << std::endl;
268  }
269  Log::info() << std::endl;
270 }
271 
272 // -----------------------------------------------------------------------------
273 
274 template<typename MODEL, typename OBS>
275 void Minimizer<MODEL, OBS>::adjModelTest(const Ht_ & Ht,
276  const H_ & H) {
277 /* Perform the adjoint test of the linear model
278  <M dx1, dx2> - <dx1, Mt dx2>)/<M dx1, dx2>
279  where dx1 and dx2 are increments obtained through randomization */
280 
281 // Model adjoint test
282  CtrlInc_ dx1(J_.jb());
283  CtrlInc_ dx2(J_.jb());
284 
285 // randomize increments
286  J_.jb().randomize(dx1);
287  J_.jb().randomize(dx2);
288 
289 /* zero out modvar and obsvar
290  since only model is being tested
291 */
292  dx1.obsVar().zero();
293  dx1.modVar().zero();
294  dx2.obsVar().zero();
295  dx2.modVar().zero();
296 
297 // run TL
298  Dual_ dummy;
299  CtrlInc_ mdx1(dx1);
300  H.multiply(mdx1, dummy, false);
301 
302 // run ADJ
303  dummy.zero();
304  CtrlInc_ mtdx2(dx2);
305  mtdx2.state().updateTime(mdx1.state().validTime() - dx1.state().validTime());
306  Ht.multiply(dummy, mtdx2, false);
307 
308 // calculate FWD < M dx1, dx2 >
309  double adj_tst_fwd = dot_product(mdx1, dx2);
310 
311 // calculate BWD < dx1, Mt dx2 >
312  double adj_tst_bwd = dot_product(dx1, mtdx2);
313 
314 // print results
315  Log::info() << "Model Adjoint Test: " << outerIteration_ << std::endl
316  << util::PrintAdjTest(adj_tst_fwd, adj_tst_bwd, "M")
317  << std::endl << std::endl;
318 }
319 
320 // -----------------------------------------------------------------------------
321 
322 template<typename MODEL, typename OBS>
323 void Minimizer<MODEL, OBS>::adjObsTest(const Ht_ & Ht,
324  const H_ & H) {
325 /* Perform the adjoint test of the linear observation operator
326  (<H dx, dy> - <dx, Ht dy>)/<H dx, dy>
327  where dx is an increment obtained through randomization and
328  dy is a DualVector obtained by first creating an increment
329  through randomization and then applying a linear observation
330  operator */
331 
332 // Model adjoint test
333  CtrlInc_ dx1(J_.jb());
334  CtrlInc_ dx2(J_.jb());
335  CtrlInc_ hthdx2(J_.jb());
336  Dual_ hdx1;
337  Dual_ hdx2;
338 
339 // randomize increments
340  J_.jb().randomize(dx1);
341  J_.jb().randomize(dx2);
342 
343 // run TL
344  H.multiply(dx1, hdx1, true);
345  H.multiply(dx2, hdx2, true);
346 
347 // run ADJ
348  J_.zeroAD(hthdx2);
349  Ht.multiply(hdx2, hthdx2, true);
350 
351 // calculate FWD < H dx1, hdx2 >
352  double adj_tst_fwd = dot_product(hdx1, hdx2);
353 
354 // calculate BWD < dx1, Ht hdx2 >
355  double adj_tst_bwd = dot_product(dx1, hthdx2);
356 
357 // print results
358  Log::info() << "Obs Adjoint Test: " << outerIteration_ << std::endl
359  << util::PrintAdjTest(adj_tst_fwd, adj_tst_bwd, "H")
360  << std::endl << std::endl;
361 }
362 
363 // -----------------------------------------------------------------------------
364 
365 /// Minimizer Factory
366 template <typename MODEL, typename OBS>
367 class MinFactory {
368  typedef CostFunction<MODEL, OBS> CostFct_;
369  public:
370  static Minimizer<MODEL, OBS> * create(const eckit::Configuration &, const CostFct_ &);
371  virtual ~MinFactory() = default;
372  protected:
373  explicit MinFactory(const std::string &);
374  private:
375  virtual Minimizer<MODEL, OBS> * make(const eckit::Configuration &, const CostFct_ &) = 0;
376  static std::map < std::string, MinFactory<MODEL, OBS> * > & getMakers() {
377  static std::map < std::string, MinFactory<MODEL, OBS> * > makers_;
378  return makers_;
379  }
380 };
381 
382 // -----------------------------------------------------------------------------
383 
384 template<class MODEL, class OBS, class FCT>
385 class MinMaker : public MinFactory<MODEL, OBS> {
386  typedef CostFunction<MODEL, OBS> CostFct_;
387  virtual Minimizer<MODEL, OBS> * make(const eckit::Configuration & conf, const CostFct_ & J) {
388  return new FCT(conf, J);
389  }
390  public:
391  explicit MinMaker(const std::string & name) : MinFactory<MODEL, OBS>(name) {}
392 };
393 
394 // =============================================================================
395 
396 template <typename MODEL, typename OBS>
397 MinFactory<MODEL, OBS>::MinFactory(const std::string & name) {
398  if (getMakers().find(name) != getMakers().end()) {
399  throw std::runtime_error(name + " already registered in minimizer factory.");
400  }
401  getMakers()[name] = this;
402 }
403 
404 // -----------------------------------------------------------------------------
405 
406 template <typename MODEL, typename OBS>
407 Minimizer<MODEL, OBS>* MinFactory<MODEL, OBS>::create(const eckit::Configuration & config,
408  const CostFct_ & J) {
409  std::string id = config.getString("algorithm");
410  Log::info() << "Minimizer algorithm=" << id << std::endl;
411  typename std::map<std::string, MinFactory<MODEL, OBS>*>::iterator j = getMakers().find(id);
412  if (j == getMakers().end()) {
413  throw std::runtime_error(id + " does not exist in minimizer factory.");
414  }
415  return (*j).second->make(config, J);
416 }
417 
418 // -----------------------------------------------------------------------------
419 
420 } // namespace oops
421 
422 #endif // OOPS_ASSIMILATION_MINIMIZER_H_
oops::Minimizer::Dual_
DualVector< MODEL, OBS > Dual_
Definition: Minimizer.h:40
oops::HtMatrix::multiply
void multiply(const DualVector< MODEL, OBS > &dy, ControlIncrement< MODEL, OBS > &dx, const bool idModel=false) const
Definition: HtMatrix.h:40
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
HtMatrix.h
oops::ControlIncrement::obsVar
ObsAuxIncrs_ & obsVar()
Get augmented observation control variable.
Definition: ControlIncrement.h:94
oops::Minimizer::adjModelTest
void adjModelTest(const Ht_ &, const H_ &)
Definition: Minimizer.h:275
oops::MinFactory::MinFactory
MinFactory(const std::string &)
Definition: Minimizer.h:397
oops::Minimizer::Ht_
HtMatrix< MODEL, OBS > Ht_
Definition: Minimizer.h:42
oops::HtMatrix
The matrix.
Definition: HtMatrix.h:33
oops::Minimizer
A Minimizer knows how to minimize a cost function.
Definition: Minimizer.h:37
oops::Minimizer::adjObsTest
void adjObsTest(const Ht_ &, const H_ &)
Definition: Minimizer.h:323
oops::ControlVariable
Control variable.
Definition: ControlVariable.h:41
CostFunction.h
HMatrix.h
oops::HMatrix
The matrix.
Definition: HMatrix.h:33
oops::MinFactory::CostFct_
CostFunction< MODEL, OBS > CostFct_
Definition: Minimizer.h:368
oops::ControlIncrement::modVar
ModelAuxIncr_ & modVar()
Get augmented model control variable.
Definition: ControlIncrement.h:90
oops::DualVector
Container of dual space vectors for all terms of the cost function.
Definition: DualVector.h:34
oops::Minimizer::adjTests
void adjTests(const eckit::Configuration &)
Definition: Minimizer.h:117
oops::Minimizer::CtrlInc_
ControlIncrement< MODEL, OBS > CtrlInc_
Definition: Minimizer.h:39
oops::Minimizer::Minimizer
Minimizer(const CostFct_ &J)
Definition: Minimizer.h:46
oops::ControlIncrement
Definition: ControlIncrement.h:46
oops::Increment::diff
void diff(const State_ &, const State_ &)
Interactions with State.
Definition: oops/interface/Increment.h:196
oops::MinMaker::MinMaker
MinMaker(const std::string &name)
Definition: Minimizer.h:391
oops::Minimizer::tlmApproxTest
void tlmApproxTest(const H_ &)
Definition: Minimizer.h:140
oops::Minimizer::H_
HMatrix< MODEL, OBS > H_
Definition: Minimizer.h:41
oops::Increment::validTime
const util::DateTime validTime() const
Time.
Definition: oops/interface/Increment.h:72
oops::Minimizer::classname
virtual const std::string classname() const =0
oops::Minimizer::~Minimizer
virtual ~Minimizer()
Definition: Minimizer.h:47
oops::Minimizer::outerIteration_
int outerIteration_
Definition: Minimizer.h:65
oops::DualVector::zero
void zero()
Definition: DualVector.h:199
oops::MinFactory::create
static Minimizer< MODEL, OBS > * create(const eckit::Configuration &, const CostFct_ &)
Definition: Minimizer.h:407
oops::Minimizer::tlmTests
void tlmTests(const eckit::Configuration &)
Definition: Minimizer.h:94
oops::MinFactory::~MinFactory
virtual ~MinFactory()=default
oops::MinFactory::make
virtual Minimizer< MODEL, OBS > * make(const eckit::Configuration &, const CostFct_ &)=0
oops::MinMaker
Definition: Minimizer.h:385
oops::MinFactory
Minimizer Factory.
Definition: Minimizer.h:367
oops::Minimizer::J_
const CostFct_ & J_
Definition: Minimizer.h:64
oops::Minimizer::minimize
ControlIncrement< MODEL, OBS > * minimize(const eckit::Configuration &)
Definition: Minimizer.h:72
oops::ModelAuxIncrement::zero
void zero()
Definition: oops/interface/ModelAuxIncrement.h:143
oops::MinFactory::getMakers
static std::map< std::string, MinFactory< MODEL, OBS > * > & getMakers()
Definition: Minimizer.h:376
oops::Minimizer::tlmPropagTest
void tlmPropagTest(const eckit::Configuration &config, const CtrlInc_ &)
Definition: Minimizer.h:181
oops::ControlIncrement::state
Increment_ & state()
Get state control variable.
Definition: ControlIncrement.h:86
oops::MinMaker::make
virtual Minimizer< MODEL, OBS > * make(const eckit::Configuration &conf, const CostFct_ &J)
Definition: Minimizer.h:387
oops::Increment::updateTime
void updateTime(const util::Duration &dt)
Definition: oops/interface/Increment.h:73
oops::State
Encapsulates the model state.
Definition: CostJbState.h:28
oops::HMatrix::multiply
void multiply(CtrlInc_ &dx, DualVector< MODEL, OBS > &dy, const bool idModel=false) const
Definition: HMatrix.h:41
oops::PostProcessor
Control model post processing.
Definition: PostProcessor.h:30
oops::Minimizer::State_
State< MODEL > State_
Definition: Minimizer.h:43
State.h
ControlIncrement.h
oops::Minimizer::doMinimize
virtual ControlIncrement< MODEL, OBS > * doMinimize(const eckit::Configuration &)=0
oops::Minimizer::tlmTaylorTest
void tlmTaylorTest(const H_ &)
Definition: Minimizer.h:217
ControlVariable.h
oops::ObsAuxIncrements::zero
void zero()
Definition: ObsAuxIncrements.h:143
oops::CostFunction
Cost Function.
Definition: CostFunction.h:53
oops::Minimizer::CostFct_
CostFunction< MODEL, OBS > CostFct_
Definition: Minimizer.h:38
DualVector.h
oops::ControlVariable::state
State_ & state()
Get state control variable.
Definition: ControlVariable.h:69
oops::MinMaker::CostFct_
CostFunction< MODEL, OBS > CostFct_
Definition: Minimizer.h:386
Increment.h