OOPS
test/interface/LinearModel.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 TEST_INTERFACE_LINEARMODEL_H_
12 #define TEST_INTERFACE_LINEARMODEL_H_
13 
14 #include <algorithm>
15 #include <cmath>
16 #include <iomanip>
17 #include <iostream>
18 #include <limits>
19 #include <memory>
20 #include <string>
21 #include <vector>
22 
23 #define ECKIT_TESTING_SELF_REGISTER_CASES 0
24 
25 #include <boost/noncopyable.hpp>
26 
27 #include "eckit/config/LocalConfiguration.h"
28 #include "eckit/testing/Test.h"
29 #include "oops/base/instantiateCovarFactory.h"
30 #include "oops/base/ModelSpaceCovarianceBase.h"
31 #include "oops/base/PostProcessor.h"
32 #include "oops/base/PostProcessorTLAD.h"
33 #include "oops/base/TrajectorySaver.h"
34 #include "oops/base/Variables.h"
35 #include "oops/generic/instantiateTlmFactory.h"
36 #include "oops/interface/Geometry.h"
37 #include "oops/interface/Increment.h"
38 #include "oops/interface/LinearModel.h"
39 #include "oops/interface/Model.h"
40 #include "oops/interface/ModelAuxControl.h"
41 #include "oops/interface/ModelAuxIncrement.h"
42 #include "oops/interface/State.h"
43 #include "oops/mpi/mpi.h"
44 #include "oops/runs/Test.h"
45 #include "oops/util/DateTime.h"
46 #include "oops/util/dot_product.h"
47 #include "oops/util/Duration.h"
48 #include "test/TestEnvironment.h"
49 
50 namespace test {
51 
52 // =============================================================================
53 
54 template <typename MODEL> class LinearModelFixture : private boost::noncopyable {
55  typedef oops::Geometry<MODEL> Geometry_;
56  typedef oops::Increment<MODEL> Increment_;
57  typedef oops::LinearModel<MODEL> LinearModel_;
58  typedef oops::Model<MODEL> Model_;
59  typedef oops::ModelAuxControl<MODEL> ModelAux_;
60  typedef oops::ModelAuxIncrement<MODEL> ModelAuxIncr_;
61  typedef oops::State<MODEL> State_;
62  typedef oops::ModelSpaceCovarianceBase<MODEL> Covariance_;
63 
64  public:
65  static const eckit::Configuration & test() {return *getInstance().test_;}
66  static const Geometry_ & resol() {return *getInstance().resol_;}
67  static const oops::Variables & ctlvars() {return *getInstance().ctlvars_;}
68  static const util::DateTime & time() {return *getInstance().time_;}
69  static const Covariance_ & covariance() {return *getInstance().B_;}
70  static const Model_ & model() {return *getInstance().model_;}
71  static const State_ & xref() {return *getInstance().xref_;}
72  static const ModelAux_ & bias() {return *getInstance().bias_;}
73  static const ModelAuxIncr_ & dbias() {return *getInstance().dbias_;}
74  static const LinearModel_ & tlm() {return *getInstance().tlm_;}
75  static void reset() {
76  getInstance().tlm_.reset();
77  getInstance().B_.reset();
78  getInstance().time_.reset();
79  getInstance().xref_.reset();
80  getInstance().model_.reset();
81  getInstance().dbias_.reset();
82  getInstance().bias_.reset();
83  getInstance().ctlvars_.reset();
84  getInstance().resol_.reset();
85  getInstance().test_.reset();
86  }
87 
88  private:
89  static LinearModelFixture<MODEL>& getInstance() {
90  static LinearModelFixture<MODEL> theLinearModelFixture;
91  return theLinearModelFixture;
92  }
93 
94  LinearModelFixture<MODEL>() {
95  test_.reset(new eckit::LocalConfiguration(TestEnvironment::config(), "linear model test"));
96  const util::Duration len(test_->getString("forecast length"));
97 
98  const eckit::LocalConfiguration resolConfig(TestEnvironment::config(), "geometry");
99  resol_.reset(new Geometry_(resolConfig, oops::mpi::world()));
100 
101  ctlvars_.reset(new oops::Variables(TestEnvironment::config(), "analysis variables"));
102 
103  const eckit::LocalConfiguration biasConf(TestEnvironment::config(), "model aux control");
104  bias_.reset(new ModelAux_(*resol_, biasConf));
105 
106  const eckit::LocalConfiguration dbiasConf =
107  TestEnvironment::config().getSubConfiguration("model aux error");
108  dbias_.reset(new ModelAuxIncr_(*resol_, dbiasConf));
109 
110  const eckit::LocalConfiguration nlConf(TestEnvironment::config(), "model");
111  model_.reset(new Model_(*resol_, nlConf));
112 
113  const eckit::LocalConfiguration iniConf(TestEnvironment::config(), "initial condition");
114  xref_.reset(new State_(*resol_, iniConf));
115  time_.reset(new util::DateTime(xref_->validTime()));
116 
117 // Create a covariance matrix
118  oops::instantiateCovarFactory<MODEL>();
119  const eckit::LocalConfiguration covar(TestEnvironment::config(), "background error");
120  B_.reset(oops::CovarianceFactory<MODEL>::create(covar, *resol_, *ctlvars_, *xref_, *xref_));
121 
122 // Linear model configuration
123  tlConf_.reset(new eckit::LocalConfiguration(TestEnvironment::config(), "linear model"));
124 
125 // Setup trajectory for TL and AD
126  oops::instantiateTlmFactory<MODEL>();
127  oops::PostProcessor<State_> post;
128  oops::PostProcessorTLAD<MODEL> pptraj;
129  tlm_.reset(new LinearModel_(*resol_, *tlConf_));
130  post.enrollProcessor(new oops::TrajectorySaver<MODEL>(*tlConf_, *resol_, *bias_, tlm_, pptraj));
131  State_ xx(*xref_);
132  model_->forecast(xx, *bias_, len, post);
133  }
134 
135  ~LinearModelFixture<MODEL>() {}
136 
137  std::unique_ptr<const eckit::LocalConfiguration> test_;
138  std::unique_ptr<const eckit::LocalConfiguration> tlConf_;
139  std::unique_ptr<const Geometry_> resol_;
140  std::unique_ptr<const util::DateTime> time_;
141  std::unique_ptr<const oops::Variables> ctlvars_;
142  std::unique_ptr<const State_> xref_;
143  std::unique_ptr<const Model_> model_;
144  std::unique_ptr<const ModelAux_> bias_;
145  std::unique_ptr<const ModelAuxIncr_> dbias_;
146  std::unique_ptr<const Covariance_> B_;
147  std::shared_ptr<LinearModel_> tlm_;
148 };
149 
150 // =============================================================================
151 
152 template <typename MODEL> void testLinearModelConstructor() {
153  typedef LinearModelFixture<MODEL> Test_;
154 
155  const util::Duration zero(0);
156  EXPECT(Test_::tlm().timeResolution() > zero);
157 }
158 
159 // -----------------------------------------------------------------------------
160 
161 template <typename MODEL> void testLinearModelZeroLength() {
162  typedef LinearModelFixture<MODEL> Test_;
163  typedef oops::Increment<MODEL> Increment_;
164  typedef oops::ModelAuxIncrement<MODEL> ModelAuxIncr_;
165 
166  const util::DateTime vt(Test_::time());
167  const util::Duration zero(0);
168 
169  Increment_ dxref(Test_::resol(), Test_::ctlvars(), vt);
170  Test_::covariance().randomize(dxref);
171  ModelAuxIncr_ daux(Test_::dbias());
172  const double ininorm = dxref.norm();
173  EXPECT(ininorm > 0.0);
174 
175  Increment_ dx(Test_::resol(), Test_::ctlvars(), vt);
176  Increment_ dxm(Test_::resol(), Test_::tlm().variables(), vt);
177  dxm = dxref;
178  Test_::tlm().forecastTL(dxm, daux, zero);
179  dx = dxm;
180  EXPECT(dx.validTime() == vt);
181  EXPECT(dx.norm() == ininorm);
182 
183  dxm.zero();
184  dxm = dxref;
185  Test_::tlm().forecastAD(dxm, daux, zero);
186  dx = dxm;
187  EXPECT(dx.validTime() == vt);
188  EXPECT(dx.norm() == ininorm);
189 }
190 
191 // -----------------------------------------------------------------------------
192 
193 template <typename MODEL> void testLinearModelZeroPert() {
194  typedef LinearModelFixture<MODEL> Test_;
195  typedef oops::Increment<MODEL> Increment_;
196  typedef oops::ModelAuxIncrement<MODEL> ModelAuxIncr_;
197 
198  const util::Duration len(Test_::test().getString("forecast length"));
199  const util::DateTime t1(Test_::time());
200  const util::DateTime t2(t1 + len);
201  EXPECT(t2 > t1);
202 
203  Increment_ dx(Test_::resol(), Test_::tlm().variables(), t1);
204  ModelAuxIncr_ daux(Test_::dbias());
205 
206  dx.zero();
207  daux.zero();
208  EXPECT(dx.norm() == 0.0);
209  Test_::tlm().forecastTL(dx, daux, len);
210  EXPECT(dx.validTime() == t2);
211  EXPECT(dx.norm() == 0.0);
212 
213  dx.zero();
214  daux.zero();
215  EXPECT(dx.norm() == 0.0);
216  Test_::tlm().forecastAD(dx, daux, len);
217  EXPECT(dx.validTime() == t1);
218  EXPECT(dx.norm() == 0.0);
219 }
220 
221 // -----------------------------------------------------------------------------
222 
223 template <typename MODEL> void testLinearModelLinearity() {
224  typedef LinearModelFixture<MODEL> Test_;
225  typedef oops::Increment<MODEL> Increment_;
226  typedef oops::ModelAuxIncrement<MODEL> ModelAuxIncr_;
227 
228  const util::Duration len(Test_::test().getString("forecast length"));
229  const util::DateTime t1(Test_::time());
230  const util::DateTime t2(t1 + len);
231  EXPECT(t2 > t1);
232  const double zz = 3.1415;
233 
234  Increment_ dx1(Test_::resol(), Test_::tlm().variables(), t1);
235  Test_::covariance().randomize(dx1);
236  ModelAuxIncr_ daux1(Test_::dbias());
237  EXPECT(dx1.norm() > 0.0);
238 
239  Increment_ dx2(dx1);
240  ModelAuxIncr_ daux2(daux1);
241 
242  Test_::tlm().forecastTL(dx1, daux1, len);
243  EXPECT(dx1.validTime() == t2);
244  dx1 *= zz;
245  daux1 *= zz;
246 
247  dx2 *= zz;
248  daux2 *= zz;
249  Test_::tlm().forecastTL(dx2, daux2, len);
250  EXPECT(dx2.validTime() == t2);
251 
252  const double tol = Test_::test().getDouble("tolerance AD");
253  EXPECT(oops::is_close(dx1.norm(), dx2.norm(), tol));
254 }
255 
256 // -----------------------------------------------------------------------------
257 
258 template <typename MODEL> void testLinearApproximation() {
259  typedef LinearModelFixture<MODEL> Test_;
260  typedef oops::Increment<MODEL> Increment_;
261  typedef oops::State<MODEL> State_;
262 
263  const util::Duration len(Test_::test().getString("forecast length"));
264  const util::DateTime t1(Test_::time());
265  const util::DateTime t2(t1 + len);
266  EXPECT(t2 > t1);
267 
268  Increment_ dx0(Test_::resol(), Test_::tlm().variables(), t1);
269  Test_::covariance().randomize(dx0);
270  EXPECT(dx0.norm() > 0.0);
271 
272  Increment_ dx(dx0);
273  Test_::tlm().forecastTL(dx, Test_::dbias(), len);
274  const double dxnorm = dx.norm();
275 
276  oops::PostProcessor<State_> post;
277  State_ xx0(Test_::xref());
278  Test_::model().forecast(xx0, Test_::bias(), len, post);
279 
280  const unsigned int ntest = Test_::test().getInt("iterations TL");
281  double zz = 1.0;
282  if (Test_::test().has("first multiplier TL")) {
283  zz = Test_::test().getDouble("first multiplier TL");
284  }
285 
286  std::vector<double> errors;
287  for (unsigned int jtest = 0; jtest < ntest; ++jtest) {
288  State_ xx(Test_::xref());
289  Increment_ pert(dx0);
290  pert *= zz;
291  xx += pert;
292  Test_::model().forecast(xx, Test_::bias(), len, post);
293 
294  Increment_ diff(Test_::resol(), Test_::tlm().variables(), t2);
295  diff.diff(xx, xx0);
296  const double difnorm = diff.norm();
297  const double err = zz * dxnorm / difnorm;
298  Increment_ derr(dx);
299  derr *= zz;
300  derr -= diff;
301  const double errnorm = derr.norm();
302  errors.push_back(errnorm / difnorm);
303  oops::Log::test() << "TL error = " << std::setprecision(16) << err
304  << ", relative error = " << errnorm / difnorm << std::endl;
305  zz /= 10.0;
306  }
307 
308 // Analyze results
309  const double approx = *std::min_element(errors.begin(), errors.end());
310  oops::Log::test() << "Test TL min error = " << approx << std::endl;
311  const double tol = Test_::test().getDouble("tolerance TL");
312  EXPECT(approx < tol);
313 }
314 
315 // -----------------------------------------------------------------------------
316 
317 template <typename MODEL> void testLinearModelAdjoint() {
318  typedef LinearModelFixture<MODEL> Test_;
319  typedef oops::Increment<MODEL> Increment_;
320  typedef oops::ModelAuxIncrement<MODEL> ModelAuxIncr_;
321 
322  const util::Duration len(Test_::test().getString("forecast length"));
323  const util::DateTime t1(Test_::time());
324  const util::DateTime t2(t1 + len);
325  EXPECT(t2 > t1);
326 
327  Increment_ dx11(Test_::resol(), Test_::tlm().variables(), t1);
328  Test_::covariance().randomize(dx11);
329  ModelAuxIncr_ daux1(Test_::dbias());
330  EXPECT(dx11.norm() > 0.0);
331  Increment_ dx12(dx11);
332  Test_::tlm().forecastTL(dx12, daux1, len);
333  EXPECT(dx12.norm() > 0.0);
334  Increment_ dx22(Test_::resol(), Test_::tlm().variables(), t2);
335  Test_::covariance().randomize(dx22);
336  ModelAuxIncr_ daux2(Test_::dbias());
337  EXPECT(dx22.norm() > 0.0);
338  Increment_ dx21(dx22);
339  Test_::tlm().forecastAD(dx21, daux2, len);
340  EXPECT(dx21.norm() > 0.0);
341 
342  EXPECT(dx11.norm() != dx22.norm());
343  EXPECT(dx11.validTime() == t1);
344  EXPECT(dx21.validTime() == t1);
345  EXPECT(dx12.validTime() == t2);
346  EXPECT(dx22.validTime() == t2);
347 
348  const double dot1 = dot_product(dx11, dx21);
349  const double dot2 = dot_product(dx12, dx22);
350  const double tol = Test_::test().getDouble("tolerance AD");
351  EXPECT(oops::is_close(dot1, dot2, tol));
352 }
353 
354 // =============================================================================
355 
356 template <typename MODEL>
357 class LinearModel : public oops::Test {
358  public:
359  LinearModel() {}
360  virtual ~LinearModel() {LinearModelFixture<MODEL>::reset();}
361 
362  private:
363  std::string testid() const override {return "test::LinearModel<" + MODEL::name() + ">";}
364 
365  void register_tests() const override {
366  std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
367 
368  ts.emplace_back(CASE("interface/GeometryIterator/testLinearModelConstructor")
369  { testLinearModelConstructor<MODEL>(); });
370  ts.emplace_back(CASE("interface/GeometryIterator/testLinearModelZeroLength")
371  { testLinearModelZeroLength<MODEL>(); });
372  ts.emplace_back(CASE("interface/GeometryIterator/testLinearModelZeroPert")
373  { testLinearModelZeroPert<MODEL>(); });
374  ts.emplace_back(CASE("interface/GeometryIterator/testLinearModelLinearity")
375  { testLinearModelLinearity<MODEL>(); });
376  ts.emplace_back(CASE("interface/GeometryIterator/testLinearApproximation")
377  { testLinearApproximation<MODEL>(); });
378  ts.emplace_back(CASE("interface/GeometryIterator/testLinearModelAdjoint")
379  { testLinearModelAdjoint<MODEL>(); });
380  }
381 
382  void clear() const override {}
383 };
384 
385 // =============================================================================
386 
387 } // namespace test
388 
389 #endif // TEST_INTERFACE_LINEARMODEL_H_
test::LinearModelFixture::ModelAuxIncr_
oops::ModelAuxIncrement< MODEL > ModelAuxIncr_
Definition: test/interface/LinearModel.h:60
test::LinearModelFixture::ctlvars_
std::unique_ptr< const oops::Variables > ctlvars_
Definition: test/interface/LinearModel.h:141
test::LinearModelFixture::Increment_
oops::Increment< MODEL > Increment_
Definition: test/interface/LinearModel.h:56
test::testLinearModelZeroLength
void testLinearModelZeroLength()
Definition: test/interface/LinearModel.h:161
test::LinearModelFixture::model_
std::unique_ptr< const Model_ > model_
Definition: test/interface/LinearModel.h:143
instantiateTlmFactory.h
test::LinearModelFixture::getInstance
static LinearModelFixture< MODEL > & getInstance()
Definition: test/interface/LinearModel.h:89
test::LinearModelFixture::xref
static const State_ & xref()
Definition: test/interface/LinearModel.h:71
TrajectorySaver.h
test::LinearModelFixture::test_
std::unique_ptr< const eckit::LocalConfiguration > test_
Definition: test/interface/LinearModel.h:137
ModelSpaceCovarianceBase.h
oops::CovarianceFactory
Covariance Factory.
Definition: ModelSpaceCovarianceBase.h:93
test::LinearModelFixture::dbias
static const ModelAuxIncr_ & dbias()
Definition: test/interface/LinearModel.h:73
test::testLinearModelAdjoint
void testLinearModelAdjoint()
Definition: test/interface/LinearModel.h:317
mpi.h
test::LinearModelFixture::test
static const eckit::Configuration & test()
Definition: test/interface/LinearModel.h:65
oops::ModelAuxControl
Definition: oops/interface/ModelAuxControl.h:35
test::LinearModelFixture::reset
static void reset()
Definition: test/interface/LinearModel.h:75
test::CASE
CASE("test_linearmodelparameterswrapper_valid_name")
Definition: LinearModelFactory.cc:22
test::LinearModelFixture::time_
std::unique_ptr< const util::DateTime > time_
Definition: test/interface/LinearModel.h:140
test
Definition: LinearModelFactory.cc:20
test::LinearModelFixture::model
static const Model_ & model()
Definition: test/interface/LinearModel.h:70
test::LinearModelFixture::dbias_
std::unique_ptr< const ModelAuxIncr_ > dbias_
Definition: test/interface/LinearModel.h:145
test::testLinearModelZeroPert
void testLinearModelZeroPert()
Definition: test/interface/LinearModel.h:193
LinearModel.h
test::testLinearModelLinearity
void testLinearModelLinearity()
Definition: test/interface/LinearModel.h:223
test::testLinearApproximation
void testLinearApproximation()
Definition: test/interface/LinearModel.h:258
test::LinearModelFixture
Definition: test/interface/LinearModel.h:54
test::LinearModel::clear
void clear() const override
Definition: test/interface/LinearModel.h:382
oops::PostProcessor::enrollProcessor
void enrollProcessor(PostBase_ *pp)
Definition: PostProcessor.h:38
test::LinearModelFixture::xref_
std::unique_ptr< const State_ > xref_
Definition: test/interface/LinearModel.h:142
test::LinearModelFixture::bias_
std::unique_ptr< const ModelAux_ > bias_
Definition: test/interface/LinearModel.h:144
Test.h
test::LinearModel::testid
std::string testid() const override
Definition: test/interface/LinearModel.h:363
test::LinearModelFixture::ModelAux_
oops::ModelAuxControl< MODEL > ModelAux_
Definition: test/interface/LinearModel.h:59
test::LinearModelFixture::Covariance_
oops::ModelSpaceCovarianceBase< MODEL > Covariance_
Definition: test/interface/LinearModel.h:62
test::LinearModelFixture::resol
static const Geometry_ & resol()
Definition: test/interface/LinearModel.h:66
test::LinearModelFixture::resol_
std::unique_ptr< const Geometry_ > resol_
Definition: test/interface/LinearModel.h:139
test::TestEnvironment::config
static const eckit::Configuration & config()
Definition: TestEnvironment.h:40
test::LinearModelFixture::ctlvars
static const oops::Variables & ctlvars()
Definition: test/interface/LinearModel.h:67
PostProcessor.h
test::LinearModelFixture::B_
std::unique_ptr< const Covariance_ > B_
Definition: test/interface/LinearModel.h:146
ModelAuxIncrement.h
test::LinearModelFixture::Model_
oops::Model< MODEL > Model_
Definition: test/interface/LinearModel.h:58
test::LinearModel::LinearModel
LinearModel()
Definition: test/interface/LinearModel.h:359
test::LinearModelFixture::LinearModel_
oops::LinearModel< MODEL > LinearModel_
Definition: test/interface/LinearModel.h:57
test::LinearModelFixture::bias
static const ModelAux_ & bias()
Definition: test/interface/LinearModel.h:72
oops::PostProcessorTLAD
Control model post processing.
Definition: PostProcessorTLAD.h:33
test::LinearModelFixture::tlm_
std::shared_ptr< LinearModel_ > tlm_
Definition: test/interface/LinearModel.h:147
oops::Geometry
Geometry class used in oops; subclass of interface class above.
Definition: oops/interface/Geometry.h:189
Model.h
instantiateCovarFactory.h
test::LinearModelFixture::Geometry_
oops::Geometry< MODEL > Geometry_
Definition: test/interface/LinearModel.h:55
TestEnvironment.h
test::LinearModelFixture::tlm
static const LinearModel_ & tlm()
Definition: test/interface/LinearModel.h:74
test::LinearModelFixture::time
static const util::DateTime & time()
Definition: test/interface/LinearModel.h:68
test::LinearModelFixture::State_
oops::State< MODEL > State_
Definition: test/interface/LinearModel.h:61
test::LinearModel::register_tests
void register_tests() const override
Definition: test/interface/LinearModel.h:365
oops::ModelAuxIncrement
Definition: oops/interface/ModelAuxIncrement.h:39
oops_variables_mod::has
logical function has(this, var)
Definition: variables_mod.F90:140
test::testLinearModelConstructor
void testLinearModelConstructor()
Definition: test/interface/LinearModel.h:152
oops::State
Encapsulates the model state.
Definition: CostJbState.h:28
oops::PostProcessor
Control model post processing.
Definition: PostProcessor.h:30
oops::Test
Definition: Test.h:39
State.h
oops::mpi::world
const eckit::mpi::Comm & world()
Default communicator with all MPI tasks (ie MPI_COMM_WORLD)
Definition: oops/mpi/mpi.cc:22
oops::Variables
Definition: oops/base/Variables.h:23
test::LinearModelFixture::covariance
static const Covariance_ & covariance()
Definition: test/interface/LinearModel.h:69
oops::ModelSpaceCovarianceBase
Definition: ModelSpaceCovarianceBase.h:61
oops::Increment
Increment Class: Difference between two states.
Definition: CostJbState.h:27
oops::TrajectorySaver
Save trajectory during forecast run.
Definition: TrajectorySaver.h:32
test::LinearModelFixture::tlConf_
std::unique_ptr< const eckit::LocalConfiguration > tlConf_
Definition: test/interface/LinearModel.h:138
PostProcessorTLAD.h
oops::LinearModel
Encapsulates the linear forecast model.
Definition: oops/interface/LinearModel.h:65
Variables.h
oops::Model
Encapsulates the nonlinear forecast model Note: to see methods that need to be implemented in the for...
Definition: oops/interface/Model.h:54
ModelAuxControl.h
test::LinearModel
Definition: test/interface/LinearModel.h:357
Geometry.h
Increment.h
test::LinearModel::~LinearModel
virtual ~LinearModel()
Definition: test/interface/LinearModel.h:360