OOPS
CostFunction.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_COSTFUNCTION_H_
12 #define OOPS_ASSIMILATION_COSTFUNCTION_H_
13 
14 #include <map>
15 #include <memory>
16 #include <string>
17 #include <vector>
18 #include <boost/noncopyable.hpp>
19 #include <boost/ptr_container/ptr_vector.hpp>
20 
21 #include "eckit/config/LocalConfiguration.h"
22 #include "eckit/mpi/Comm.h"
23 #include "oops/assimilation/ControlIncrement.h"
24 #include "oops/assimilation/ControlVariable.h"
25 #include "oops/assimilation/CostJbState.h"
26 #include "oops/assimilation/CostJbTotal.h"
27 #include "oops/assimilation/CostJo.h"
28 #include "oops/assimilation/CostTermBase.h"
29 #include "oops/assimilation/DualVector.h"
30 #include "oops/assimilation/JqTermTLAD.h"
31 #include "oops/base/PostProcessor.h"
32 #include "oops/base/PostProcessorTLAD.h"
33 #include "oops/base/Variables.h"
34 #include "oops/interface/Geometry.h"
35 #include "oops/interface/Increment.h"
36 #include "oops/interface/State.h"
37 #include "oops/mpi/mpi.h"
38 #include "oops/util/DateTime.h"
39 #include "oops/util/dot_product.h"
40 #include "oops/util/Duration.h"
41 #include "oops/util/Logger.h"
42 
43 namespace oops {
44 
45 // -----------------------------------------------------------------------------
46 
47 /// Cost Function
48 /*!
49  * The CostFunction defines and manages the computation of all the terms
50  * of the variational data assimilation cost function.
51  */
52 
53 template<typename MODEL, typename OBS> class CostFunction : private boost::noncopyable {
54  typedef ControlIncrement<MODEL, OBS> CtrlInc_;
55  typedef ControlVariable<MODEL, OBS> CtrlVar_;
56  typedef CostJbTotal<MODEL, OBS> JbTotal_;
57  typedef CostTermBase<MODEL, OBS> CostBase_;
58  typedef JqTermTLAD<MODEL> JqTermTLAD_;
59  typedef Geometry<MODEL> Geometry_;
60  typedef State<MODEL> State_;
61  typedef Increment<MODEL> Increment_;
62 
63  public:
64  explicit CostFunction(const eckit::Configuration &);
65  virtual ~CostFunction() {}
66 
67  double evaluate(const CtrlVar_ &,
68  const eckit::Configuration & config = eckit::LocalConfiguration(),
69  PostProcessor<State_> post = PostProcessor<State_>() );
70  double linearize(const CtrlVar_ &, const eckit::Configuration &,
71  PostProcessor<State_> post = PostProcessor<State_>() );
72 
73  virtual void runTLM(CtrlInc_ &, PostProcessorTLAD<MODEL> &,
74  PostProcessor<Increment_> post = PostProcessor<Increment_>(),
75  const bool idModel = false) const = 0;
76  virtual void runADJ(CtrlInc_ &, PostProcessorTLAD<MODEL> &,
77  PostProcessor<Increment_> post = PostProcessor<Increment_>(),
78  const bool idModel = false) const = 0;
79  virtual void zeroAD(CtrlInc_ &) const = 0;
80 
81  virtual void runNL(CtrlVar_ &, PostProcessor<State_>&) const = 0;
82 
83  void addIncrement(CtrlVar_ &, const CtrlInc_ &,
84  PostProcessor<Increment_> post = PostProcessor<Increment_>() ) const;
85  void resetLinearization();
86 
87 /// Compute cost function gradient at first guess (without Jb).
88  void computeGradientFG(CtrlInc_ &) const;
89 
90 /// Access \f$ J_b\f$
91  const JbTotal_ & jb() const {return *jb_;}
92 /// Access terms of the cost function other than \f$ J_b\f$
93  const CostBase_ & jterm(const unsigned ii) const {return jterms_[ii];}
94  unsigned nterms() const {return jterms_.size();}
95  double getCostJb() const {return costJb_;}
96  double getCostJoJc() const {return costJoJc_;}
97 
98  protected:
99  void setupTerms(const eckit::Configuration &);
100  void setupTerms(const eckit::Configuration &, const State_ &); // generic 1d-var
101  const CtrlVar_ & background() const {return *xb_;}
102 
103  private:
104  virtual void addIncr(CtrlVar_ &, const CtrlInc_ &, PostProcessor<Increment_>&) const = 0;
105 
106  virtual CostJbState<MODEL> * newJb(const eckit::Configuration &, const Geometry_ &,
107  const CtrlVar_ &) const = 0;
108  virtual CostJo<MODEL, OBS> * newJo(const eckit::Configuration &) const = 0;
109  virtual CostTermBase<MODEL, OBS> * newJc(const eckit::Configuration &,
110  const Geometry_ &) const = 0;
111  virtual void doLinearize(const Geometry_ &, const eckit::Configuration &,
112  const CtrlVar_ &, const CtrlVar_ &,
113  PostProcessor<State_> &, PostProcessorTLAD<MODEL> &) = 0;
114  virtual const Geometry_ & geometry() const = 0;
115 
116 // Data members
117  std::unique_ptr<const CtrlVar_> xb_;
118  std::unique_ptr<JbTotal_> jb_;
119  boost::ptr_vector<CostBase_> jterms_;
120 
121  mutable double costJb_;
122  mutable double costJoJc_;
123 };
124 
125 // -----------------------------------------------------------------------------
126 
127 /// Cost Function Factory
128 template <typename MODEL, typename OBS>
129 class CostFactory {
130  public:
131  static CostFunction<MODEL, OBS> * create(const eckit::Configuration &, const eckit::mpi::Comm &);
132  virtual ~CostFactory() = default;
133 
134  protected:
135  explicit CostFactory(const std::string &);
136  private:
137  virtual CostFunction<MODEL, OBS> * make(const eckit::Configuration &,
138  const eckit::mpi::Comm &) = 0;
139  static std::map < std::string, CostFactory<MODEL, OBS> * > & getMakers() {
140  static std::map < std::string, CostFactory<MODEL, OBS> * > makers_;
141  return makers_;
142  }
143 };
144 
145 template<class MODEL, class OBS, class FCT>
146 class CostMaker : public CostFactory<MODEL, OBS> {
147  private:
148  CostFunction<MODEL, OBS> * make(const eckit::Configuration & config,
149  const eckit::mpi::Comm & comm) override
150  {return new FCT(config, comm);}
151  public:
152  explicit CostMaker(const std::string & name) : CostFactory<MODEL, OBS>(name) {}
153 };
154 
155 // =============================================================================
156 
157 // Factory
158 // -----------------------------------------------------------------------------
159 
160 template <typename MODEL, typename OBS>
161 CostFactory<MODEL, OBS>::CostFactory(const std::string & name) {
162  if (getMakers().find(name) != getMakers().end()) {
163  throw std::runtime_error(name + " already registered in cost function factory.");
164  }
165  getMakers()[name] = this;
166 }
167 
168 // -----------------------------------------------------------------------------
169 
170 template <typename MODEL, typename OBS>
171 CostFunction<MODEL, OBS>* CostFactory<MODEL, OBS>::create(const eckit::Configuration & config,
172  const eckit::mpi::Comm & comm) {
173  std::string id = config.getString("cost type");
174  Log::trace() << "Variational Assimilation Type=" << id << std::endl;
175  typename std::map<std::string, CostFactory<MODEL, OBS>*>::iterator j = getMakers().find(id);
176  if (j == getMakers().end()) {
177  throw std::runtime_error(id + " does not exist in cost function factory.");
178  }
179  Log::trace() << "CostFactory::create found cost function type" << std::endl;
180  return (*j).second->make(config, comm);
181 }
182 
183 // -----------------------------------------------------------------------------
184 // Cost Function
185 // -----------------------------------------------------------------------------
186 
187 template<typename MODEL, typename OBS>
188 CostFunction<MODEL, OBS>::CostFunction(const eckit::Configuration & config)
189  : jb_(), jterms_()
190 {}
191 
192 // -----------------------------------------------------------------------------
193 
194 template<typename MODEL, typename OBS>
195 void CostFunction<MODEL, OBS>::setupTerms(const eckit::Configuration & config) {
196  Log::trace() << "CostFunction::setupTerms start" << std::endl;
197 
198 // Jo
199  CostJo<MODEL, OBS> * jo = this->newJo(config);
200  jterms_.push_back(jo);
201  Log::trace() << "CostFunction::setupTerms Jo added" << std::endl;
202 
203 // Jb
204  xb_.reset(new CtrlVar_(config, this->geometry(), jo->obspaces()));
205  jb_.reset(new JbTotal_(*xb_, this->newJb(config, this->geometry(), *xb_),
206  config, this->geometry(), jo->obspaces()));
207  Log::trace() << "CostFunction::setupTerms Jb added" << std::endl;
208 
209 // Other constraints
210  std::vector<eckit::LocalConfiguration> jcs;
211  config.get("constraints", jcs);
212  for (size_t jj = 0; jj < jcs.size(); ++jj) {
213  CostTermBase<MODEL, OBS> * jc = this->newJc(jcs[jj], this->geometry());
214  jterms_.push_back(jc);
215  }
216  Log::trace() << "CostFunction::setupTerms done" << std::endl;
217 }
218 
219 // -----------------------------------------------------------------------------
220 // This setup terms method has been written for the generic 1d-var which is
221 // under development in UFO
222 // -----------------------------------------------------------------------------
223 template<typename MODEL, typename OBS>
224 void CostFunction<MODEL, OBS>::setupTerms(const eckit::Configuration & config,
225  const State_ & statein) {
226  Log::trace() << "CostFunction::setupTerms start" << std::endl;
227 
228 // Jo
229  CostJo<MODEL, OBS> * jo = this->newJo(config);
230  jterms_.push_back(jo);
231  Log::trace() << "CostFunction::setupTerms Jo added" << std::endl;
232 
233 // Jb
234  xb_.reset(new CtrlVar_(config, statein, jo->obspaces()));
235  jb_.reset(new JbTotal_(*xb_, this->newJb(config, this->geometry(), *xb_),
236  config, this->geometry(), jo->obspaces()));
237  Log::trace() << "CostFunction::setupTerms Jb added" << std::endl;
238 
239 // Other constraints
240  std::vector<eckit::LocalConfiguration> jcs;
241  config.get("constraints", jcs);
242  for (size_t jj = 0; jj < jcs.size(); ++jj) {
243  CostTermBase<MODEL, OBS> * jc = this->newJc(jcs[jj], this->geometry());
244  jterms_.push_back(jc);
245  }
246  Log::trace() << "CostFunction::setupTerms done" << std::endl;
247 }
248 
249 // -----------------------------------------------------------------------------
250 
251 template<typename MODEL, typename OBS>
252 double CostFunction<MODEL, OBS>::evaluate(const CtrlVar_ & fguess,
253  const eckit::Configuration & config,
254  PostProcessor<State_> post) {
255  Log::trace() << "CostFunction::evaluate start" << std::endl;
256 // Setup terms of cost function
257  PostProcessor<State_> pp(post);
258  jb_->initialize(fguess);
259  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
260  pp.enrollProcessor(jterms_[jj].initialize(fguess, config));
261  }
262 
263 // Run NL model
264  CtrlVar_ mfguess(fguess);
265  this->runNL(mfguess, pp);
266 
267 // Cost function value
268  double zzz = 0.0;
269  costJb_ = jb_->finalize(mfguess);
270  zzz += costJb_;
271  costJoJc_ = 0.0;
272  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
273  costJoJc_ += jterms_[jj].finalize();
274  }
275  zzz += costJoJc_;
276  Log::test() << "CostFunction: Nonlinear J = " << zzz << std::endl;
277  Log::trace() << "CostFunction::evaluate done" << std::endl;
278  return zzz;
279 }
280 
281 // -----------------------------------------------------------------------------
282 
283 template<typename MODEL, typename OBS>
284 double CostFunction<MODEL, OBS>::linearize(const CtrlVar_ & fguess,
285  const eckit::Configuration & innerConf,
286  PostProcessor<State_> post) {
287  Log::trace() << "CostFunction::linearize start" << std::endl;
288 // Inner loop resolution
289  const eckit::LocalConfiguration resConf(innerConf, "geometry");
290  const Geometry_ lowres(resConf, this->geometry().getComm(), this->geometry().timeComm());
291 
292 // Setup trajectory for terms of cost function
293  PostProcessorTLAD<MODEL> pptraj;
294  JqTermTLAD_ * jq = jb_->initializeTraj(fguess, lowres, innerConf);
295  pptraj.enrollProcessor(jq);
296  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
297  pptraj.enrollProcessor(jterms_[jj].initializeTraj(fguess, lowres, innerConf));
298  }
299 
300 // Specific linearization if needed (including TLM)
301  this->doLinearize(lowres, innerConf, *xb_, fguess, post, pptraj);
302 
303 // Run NL model
304  double zzz = this->evaluate(fguess, innerConf, post);
305 
306 // Finalize trajectory setup
307  jb_->finalizeTraj(jq);
308  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
309  jterms_[jj].finalizeTraj();
310  }
311 
312  Log::trace() << "CostFunction::linearize done" << std::endl;
313  return zzz;
314 }
315 
316 // -----------------------------------------------------------------------------
317 
318 template<typename MODEL, typename OBS>
319 void CostFunction<MODEL, OBS>::computeGradientFG(CtrlInc_ & grad) const {
320  Log::trace() << "CostFunction::computeGradientFG start" << std::endl;
321  PostProcessor<Increment_> pp;
322  PostProcessorTLAD<MODEL> costad;
323  this->zeroAD(grad);
324 
325  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
326  std::shared_ptr<const GeneralizedDepartures> tmp(jterms_[jj].newGradientFG());
327  costad.enrollProcessor(jterms_[jj].setupAD(tmp, grad));
328  }
329 
330  this->runADJ(grad, costad, pp);
331  Log::info() << "CostFunction::computeGradientFG: gradient:" << grad << std::endl;
332  Log::trace() << "CostFunction::computeGradientFG done" << std::endl;
333 }
334 
335 // -----------------------------------------------------------------------------
336 
337 template<typename MODEL, typename OBS>
338 void CostFunction<MODEL, OBS>::addIncrement(CtrlVar_ & xx, const CtrlInc_ & dx,
339  PostProcessor<Increment_> post) const {
340  Log::trace() << "CostFunction::addIncrement start" << std::endl;
341  Log::info() << "CostFunction::addIncrement: First guess:" << xx << std::endl;
342  Log::info() << "CostFunction::addIncrement: Increment:" << dx << std::endl;
343 
344  xx.obsVar() += dx.obsVar();
345  xx.modVar() += dx.modVar();
346  this->addIncr(xx, dx, post);
347 
348  Log::info() << "CostFunction::addIncrement: Analysis: " << xx << std::endl;
349  Log::test() << "CostFunction::addIncrement: Analysis: " << xx << std::endl;
350  Log::trace() << "CostFunction::addIncrement done" << std::endl;
351 }
352 
353 // -----------------------------------------------------------------------------
354 
355 template<typename MODEL, typename OBS>
356 void CostFunction<MODEL, OBS>::resetLinearization() {
357  Log::trace() << "CostFunction::resetLinearization start" << std::endl;
358  for (unsigned jj = 0; jj < jterms_.size(); ++jj) {
359  jterms_[jj].resetLinearization();
360  }
361  Log::trace() << "CostFunction::resetLinearization done" << std::endl;
362 }
363 
364 // -----------------------------------------------------------------------------
365 
366 } // namespace oops
367 
368 #endif // OOPS_ASSIMILATION_COSTFUNCTION_H_
oops
The namespace for the main oops code.
Definition: ErrorCovarianceL95.cc:22
oops::CostFunction::zeroAD
virtual void zeroAD(CtrlInc_ &) const =0
oops::CostFunction::~CostFunction
virtual ~CostFunction()
Definition: CostFunction.h:65
oops::ControlVariable::obsVar
ObsAuxCtrls_ & obsVar()
Get augmented observation control variable.
Definition: ControlVariable.h:77
oops::CostMaker::make
CostFunction< MODEL, OBS > * make(const eckit::Configuration &config, const eckit::mpi::Comm &comm) override
Definition: CostFunction.h:148
oops::ControlIncrement::obsVar
ObsAuxIncrs_ & obsVar()
Get augmented observation control variable.
Definition: ControlIncrement.h:94
oops::CostTermBase
Base Class for Cost Function Terms.
Definition: CostTermBase.h:37
oops::CostFunction::jterms_
boost::ptr_vector< CostBase_ > jterms_
Definition: CostFunction.h:119
oops::CostFunction::linearize
double linearize(const CtrlVar_ &, const eckit::Configuration &, PostProcessor< State_ > post=PostProcessor< State_ >())
Definition: CostFunction.h:284
oops::ControlVariable
Control variable.
Definition: ControlVariable.h:41
oops::CostFunction::nterms
unsigned nterms() const
Definition: CostFunction.h:94
oops::CostFunction::background
const CtrlVar_ & background() const
Definition: CostFunction.h:101
oops::CostFunction::JqTermTLAD_
JqTermTLAD< MODEL > JqTermTLAD_
Definition: CostFunction.h:58
mpi.h
oops::CostFunction::costJoJc_
double costJoJc_
Definition: CostFunction.h:122
oops::ControlIncrement::modVar
ModelAuxIncr_ & modVar()
Get augmented model control variable.
Definition: ControlIncrement.h:90
oops::CostFunction::runTLM
virtual void runTLM(CtrlInc_ &, PostProcessorTLAD< MODEL > &, PostProcessor< Increment_ > post=PostProcessor< Increment_ >(), const bool idModel=false) const =0
oops::CostFunction::jterm
const CostBase_ & jterm(const unsigned ii) const
Access terms of the cost function other than .
Definition: CostFunction.h:93
oops::CostFunction::evaluate
double evaluate(const CtrlVar_ &, const eckit::Configuration &config=eckit::LocalConfiguration(), PostProcessor< State_ > post=PostProcessor< State_ >())
Definition: CostFunction.h:252
oops::CostJbState
Jb Cost Function Base Class.
Definition: CostJbState.h:39
oops::ControlIncrement
Definition: ControlIncrement.h:46
oops::CostFunction::addIncrement
void addIncrement(CtrlVar_ &, const CtrlInc_ &, PostProcessor< Increment_ > post=PostProcessor< Increment_ >()) const
Definition: CostFunction.h:338
oops::CostFunction::CtrlVar_
ControlVariable< MODEL, OBS > CtrlVar_
Definition: CostFunction.h:55
oops::CostFunction::CtrlInc_
ControlIncrement< MODEL, OBS > CtrlInc_
Definition: CostFunction.h:54
oops::JqTermTLAD
Definition: CostJb3D.h:29
oops::CostFunction::resetLinearization
void resetLinearization()
Definition: CostFunction.h:356
CostJbTotal.h
oops::CostFunction::jb_
std::unique_ptr< JbTotal_ > jb_
Definition: CostFunction.h:118
oops::ControlVariable::modVar
ModelAux_ & modVar()
Get augmented model control variable.
Definition: ControlVariable.h:73
oops::CostFunction::jb
const JbTotal_ & jb() const
Access .
Definition: CostFunction.h:91
oops::CostFunction::computeGradientFG
void computeGradientFG(CtrlInc_ &) const
Compute cost function gradient at first guess (without Jb).
Definition: CostFunction.h:319
CostTermBase.h
oops::CostFunction::getCostJoJc
double getCostJoJc() const
Definition: CostFunction.h:96
oops::CostFactory
Cost Function Factory.
Definition: CostFunction.h:129
oops::CostFunction::newJc
virtual CostTermBase< MODEL, OBS > * newJc(const eckit::Configuration &, const Geometry_ &) const =0
oops::CostFunction::addIncr
virtual void addIncr(CtrlVar_ &, const CtrlInc_ &, PostProcessor< Increment_ > &) const =0
oops::CostFunction::CostBase_
CostTermBase< MODEL, OBS > CostBase_
Definition: CostFunction.h:57
oops::CostFunction::getCostJb
double getCostJb() const
Definition: CostFunction.h:95
oops::PostProcessor::enrollProcessor
void enrollProcessor(PostBase_ *pp)
Definition: PostProcessor.h:38
oops::CostMaker::CostMaker
CostMaker(const std::string &name)
Definition: CostFunction.h:152
oops::CostFunction::Geometry_
Geometry< MODEL > Geometry_
Definition: CostFunction.h:59
oops::CostFunction::xb_
std::unique_ptr< const CtrlVar_ > xb_
Definition: CostFunction.h:117
oops::PostProcessorTLAD::enrollProcessor
void enrollProcessor(PostBaseTLAD_ *pp)
Definition: PostProcessorTLAD.h:43
oops::CostFunction::State_
State< MODEL > State_
Definition: CostFunction.h:60
CostJbState.h
oops::CostJo::obspaces
const ObsSpaces_ & obspaces() const
Definition: CostJo.h:116
oops::CostFunction::CostFunction
CostFunction(const eckit::Configuration &)
Definition: CostFunction.h:188
PostProcessor.h
JqTermTLAD.h
oops::CostFunction::runADJ
virtual void runADJ(CtrlInc_ &, PostProcessorTLAD< MODEL > &, PostProcessor< Increment_ > post=PostProcessor< Increment_ >(), const bool idModel=false) const =0
oops::CostFactory::getMakers
static std::map< std::string, CostFactory< MODEL, OBS > * > & getMakers()
Definition: CostFunction.h:139
oops::CostFunction::newJb
virtual CostJbState< MODEL > * newJb(const eckit::Configuration &, const Geometry_ &, const CtrlVar_ &) const =0
oops::PostProcessorTLAD
Control model post processing.
Definition: PostProcessorTLAD.h:33
oops::Geometry
Geometry class used in oops; subclass of interface class above.
Definition: oops/interface/Geometry.h:189
oops::CostFunction::JbTotal_
CostJbTotal< MODEL, OBS > JbTotal_
Definition: CostFunction.h:56
oops::CostFactory::~CostFactory
virtual ~CostFactory()=default
oops::CostJbTotal
Control variable increment.
Definition: ControlIncrement.h:44
oops::CostFunction::newJo
virtual CostJo< MODEL, OBS > * newJo(const eckit::Configuration &) const =0
oops::State
Encapsulates the model state.
Definition: CostJbState.h:28
oops::CostFactory::CostFactory
CostFactory(const std::string &)
Definition: CostFunction.h:161
oops::CostFactory::make
virtual CostFunction< MODEL, OBS > * make(const eckit::Configuration &, const eckit::mpi::Comm &)=0
oops::PostProcessor
Control model post processing.
Definition: PostProcessor.h:30
State.h
oops::CostFunction::runNL
virtual void runNL(CtrlVar_ &, PostProcessor< State_ > &) const =0
ControlIncrement.h
CostJo.h
oops::CostFunction::setupTerms
void setupTerms(const eckit::Configuration &)
Definition: CostFunction.h:195
oops::CostFunction::geometry
virtual const Geometry_ & geometry() const =0
oops::CostMaker
Definition: CostFunction.h:146
oops::Increment
Increment Class: Difference between two states.
Definition: CostJbState.h:27
ControlVariable.h
oops::CostFunction
Cost Function.
Definition: CostFunction.h:53
PostProcessorTLAD.h
DualVector.h
oops::CostFunction::costJb_
double costJb_
Definition: CostFunction.h:121
oops::CostFactory::create
static CostFunction< MODEL, OBS > * create(const eckit::Configuration &, const eckit::mpi::Comm &)
Definition: CostFunction.h:171
oops::CostFunction::Increment_
Increment< MODEL > Increment_
Definition: CostFunction.h:61
oops::CostFunction::doLinearize
virtual void doLinearize(const Geometry_ &, const eckit::Configuration &, const CtrlVar_ &, const CtrlVar_ &, PostProcessor< State_ > &, PostProcessorTLAD< MODEL > &)=0
Variables.h
oops::CostJo
Jo Cost Function.
Definition: CostJo.h:54
Geometry.h
Increment.h