IODA Bundle
test/interface/LinearObsOperator.h
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2017-2018 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 #ifndef TEST_INTERFACE_LINEAROBSOPERATOR_H_
9 #define TEST_INTERFACE_LINEAROBSOPERATOR_H_
10 
11 #include <memory>
12 #include <string>
13 #include <vector>
14 
15 #define ECKIT_TESTING_SELF_REGISTER_CASES 0
16 
17 
18 #include "eckit/testing/Test.h"
19 #include "oops/base/Variables.h"
20 #include "oops/interface/LinearObsOperator.h"
21 #include "oops/interface/ObsAuxControl.h"
22 #include "oops/interface/ObsAuxCovariance.h"
23 #include "oops/interface/ObsAuxIncrement.h"
24 #include "oops/interface/ObsDiagnostics.h"
25 #include "oops/interface/ObsOperator.h"
26 #include "oops/runs/Test.h"
27 #include "oops/util/dot_product.h"
28 #include "oops/util/Expect.h"
29 #include "oops/util/Logger.h"
30 #include "test/interface/ObsTestsFixture.h"
31 #include "test/TestEnvironment.h"
32 
33 namespace test {
34 
35 const char *expectConstructorToThrow = "expect constructor to throw exception with message";
36 const char *expectSetTrajectoryToThrow = "expect setTrajectory to throw exception with message";
37 const char *expectSimulateObsToThrow = "expect simulateObs to throw exception with message";
38 const char *expectSimulateObsTLToThrow = "expect simulateObsTL to throw exception with message";
39 const char *expectSimulateObsADToThrow = "expect simulateObsAD to throw exception with message";
40 
41 // -----------------------------------------------------------------------------
42 /// \brief tests constructor and print method
43 template <typename OBS> void testConstructor() {
44  typedef ObsTestsFixture<OBS> Test_;
45  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
46 
47  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
48  const eckit::LocalConfiguration & conf = Test_::config(jj);
49  // Use ObsOperator section of yaml unless LinearObsOperator is specified
50  std::string confname = "obs operator";
51  if (conf.has("linear obs operator")) confname = "linear obs operator";
52  eckit::LocalConfiguration linobsopconf(conf, confname);
53 
54  if (!Test_::config(jj).has(expectConstructorToThrow)) {
55  std::unique_ptr<LinearObsOperator_> linobsop(
56  new LinearObsOperator_(Test_::obspace()[jj], linobsopconf));
57  EXPECT(linobsop.get());
58  oops::Log::test() << "Testing LinearObsOperator: " << *linobsop << std::endl;
59  linobsop.reset();
60  EXPECT(!linobsop.get());
61  } else {
62  // The constructor is expected to throw an exception containing the specified string.
63  const std::string expectedMessage = Test_::config(jj).getString(expectConstructorToThrow);
64  EXPECT_THROWS_MSG(LinearObsOperator_(Test_::obspace()[jj], linobsopconf),
65  expectedMessage.c_str());
66  }
67  }
68 }
69 
70 // -----------------------------------------------------------------------------
71 
72 template <typename OBS> void testLinearity() {
73  typedef ObsTestsFixture<OBS> Test_;
74  typedef oops::GeoVaLs<OBS> GeoVaLs_;
75  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
76  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
77  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
78  typedef oops::ObsOperator<OBS> ObsOperator_;
79  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
80  typedef oops::ObsVector<OBS> ObsVector_;
81 
82  const double zero = 0.0;
83  const double coef = 3.14;
84  const double tol = 1.0e-11;
85 
86  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
87  const eckit::LocalConfiguration & conf = Test_::config(jj);
88  if (conf.has(expectConstructorToThrow) ||
89  conf.has(expectSetTrajectoryToThrow) ||
90  conf.has(expectSimulateObsToThrow) ||
91  conf.has(expectSimulateObsTLToThrow))
92  continue;
93 
94  // initialize observation operator (set variables requested from the model,
95  // variables simulated by the observation operator, other init)
96  eckit::LocalConfiguration obsopconf(conf, "obs operator");
97  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
98  // initialize TL/AD observation operator (set model variables for Jacobian),
99  // other init)
100  // Use ObsOperator section of yaml unless LinearObsOperator is specified
101  std::string confname = "obs operator";
102  if (conf.has("linear obs operator")) confname = "linear obs operator";
103  eckit::LocalConfiguration linobsopconf(conf, confname);
104  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
105 
106  // initialize obs bias
107  eckit::LocalConfiguration biasconf = conf.getSubConfiguration("obs bias");
108  typename ObsAuxCtrl_::Parameters_ biasparams;
109  biasparams.validateAndDeserialize(biasconf);
110  const ObsAuxCtrl_ ybias(Test_::obspace()[jj], biasparams);
111  ObsAuxIncr_ ybinc(Test_::obspace()[jj], biasparams);
112 
113  // read geovals from the file
114  const eckit::LocalConfiguration gconf(conf, "geovals");
115  oops::Variables hopvars = hop.requiredVars();
116  hopvars += ybias.requiredVars();
117  const GeoVaLs_ gval(gconf, Test_::obspace()[jj], hopvars);
118 
119  // initialize Obs. Bias Covariance
120  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], biasparams);
121 
122  // set trajectory for TL/AD to be the geovals from the file
123  hoptl.setTrajectory(gval, ybias);
124 
125  // create obsvector
126  ObsVector_ dy1(Test_::obspace()[jj]);
127 
128  // create geovals
129  GeoVaLs_ dx(gconf, Test_::obspace()[jj], hoptl.requiredVars());
130 
131  // test rms(H * (dx, ybinc)) = 0, when dx = 0
132  dx.zero();
133  ybinc.zero();
134  hoptl.simulateObsTL(dx, dy1, ybinc);
135  EXPECT(dy1.rms() == zero);
136 
137  // test rms(H * (dx, ybinc)) > 0, when dx is random
138  dx.random();
139  Bobsbias.randomize(ybinc);
140  hoptl.simulateObsTL(dx, dy1, ybinc);
141  EXPECT(dy1.rms() > zero);
142 
143  // test k * H * (dx, ybinc) ~ H * (k*dx, k*ybinc)
144  dy1 *= coef;
145  dx *= coef;
146  ybinc *= coef;
147  ObsVector_ dy2(Test_::obspace()[jj]);
148  hoptl.simulateObsTL(dx, dy2, ybinc);
149 
150  dy2 -= dy1;
151  EXPECT(dy2.rms() / dy1.rms() < tol);
152  }
153 }
154 
155 // -----------------------------------------------------------------------------
156 
157 template <typename OBS> void testAdjoint() {
158  typedef ObsTestsFixture<OBS> Test_;
159  typedef oops::GeoVaLs<OBS> GeoVaLs_;
160  typedef oops::ObsOperator<OBS> ObsOperator_;
161  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
162  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
163  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
164  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
165  typedef oops::ObsVector<OBS> ObsVector_;
166 
167  const double zero = 0.0;
168 
169  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
170  const eckit::LocalConfiguration & conf = Test_::config(jj);
171  if (conf.has(expectConstructorToThrow) ||
172  conf.has(expectSetTrajectoryToThrow) ||
173  conf.has(expectSimulateObsToThrow) ||
174  conf.has(expectSimulateObsTLToThrow) ||
175  conf.has(expectSimulateObsADToThrow))
176  continue;
177 
178  // initialize observation operator (set variables requested from the model,
179  // variables simulated by the observation operator, other init)
180  eckit::LocalConfiguration obsopconf(conf, "obs operator");
181  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
182  // initialize TL/AD observation operator (set model variables for Jacobian),
183  // other init)
184  // Use ObsOperator section of yaml unless LinearObsOperator is specified
185  std::string confname = "obs operator";
186  if (conf.has("linear obs operator")) confname = "linear obs operator";
187  eckit::LocalConfiguration linobsopconf(conf, confname);
188  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
189 
190  const double tol = conf.getDouble("linear obs operator test.tolerance AD");
191  // initialize bias correction
192  eckit::LocalConfiguration biasconf = conf.getSubConfiguration("obs bias");
193  typename ObsAuxCtrl_::Parameters_ biasparams;
194  biasparams.validateAndDeserialize(biasconf);
195  const ObsAuxCtrl_ ybias(Test_::obspace()[jj], biasparams);
196  ObsAuxIncr_ ybinc1(Test_::obspace()[jj], biasparams); // TL
197  ObsAuxIncr_ ybinc2(Test_::obspace()[jj], biasparams); // AD
198 
199  // initialize Obs. Bias Covariance
200  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], biasparams);
201 
202  // read geovals from the file
203  eckit::LocalConfiguration gconf(conf, "geovals");
204  oops::Variables hopvars = hop.requiredVars();
205  hopvars += ybias.requiredVars();
206  const GeoVaLs_ gval(gconf, Test_::obspace()[jj], hopvars);
207 
208  // set TL/AD trajectory to the geovals from the file
209  hoptl.setTrajectory(gval, ybias);
210 
211  ObsVector_ dy1(Test_::obspace()[jj]);
212  ObsVector_ dy2(Test_::obspace()[jj]);
213  GeoVaLs_ dx1(gconf, Test_::obspace()[jj], hoptl.requiredVars());
214  GeoVaLs_ dx2(gconf, Test_::obspace()[jj], hoptl.requiredVars());
215 
216  // calculate dy1 = H (dx1, ybinc1) (with random dx1, and random ybinc1)
217  dx1.random();
218  EXPECT(dot_product(dx1, dx1) > zero); // BOOST_REQUIRE
219  Bobsbias.randomize(ybinc1);
220  hoptl.simulateObsTL(dx1, dy1, ybinc1);
221  EXPECT(dot_product(dy1, dy1) > zero);
222 
223  // calculate (dx2, ybinc2) = HT dy2 (with random dy2)
224  dy2.random();
225  EXPECT(dot_product(dy2, dy2) > zero); // BOOST_REQUIRE
226  dx2.zero();
227  ybinc2.zero();
228  hoptl.simulateObsAD(dx2, dy2, ybinc2);
229  EXPECT(dot_product(dx2, dx2) > zero);
230 
231  const double zz1 = dot_product(dx1, dx2) + dot_product(ybinc1, ybinc2);
232  const double zz2 = dot_product(dy1, dy2);
233 
234  oops::Log::test() << "Adjoint test result: (<x,HTy>-<Hx,y>)/<Hx,y> = "
235  << (zz1-zz2)/zz2 << std::endl;
236 
237  EXPECT(zz1 != zero);
238  EXPECT(zz2 != zero);
239  EXPECT(oops::is_close(zz1, zz2, tol));
240  }
241 }
242 
243 // -----------------------------------------------------------------------------
244 
245 template <typename OBS> void testTangentLinear() {
246  // Test ||(hop(x+alpha*dx)-hop(x)) - hoptl(alpha*dx)|| < tol
247  typedef ObsTestsFixture<OBS> Test_;
248  typedef oops::GeoVaLs<OBS> GeoVaLs_;
249  typedef oops::ObsDiagnostics<OBS> ObsDiags_;
250  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
251  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
252  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
253  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
254  typedef oops::ObsOperator<OBS> ObsOperator_;
255  typedef oops::ObsVector<OBS> ObsVector_;
256 
257  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
258  const eckit::LocalConfiguration & conf = Test_::config(jj);
259  if (conf.has(expectConstructorToThrow) ||
260  conf.has(expectSetTrajectoryToThrow) ||
261  conf.has(expectSimulateObsToThrow) ||
262  conf.has(expectSimulateObsTLToThrow))
263  continue;
264 
265  // initialize observation operator (set variables requested from the model,
266  // variables simulated by the observation operator, other init)
267  eckit::LocalConfiguration obsopconf(conf, "obs operator");
268  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
269  // initialize TL/AD observation operator (set model variables for Jacobian),
270  // other init)
271  // Use ObsOperator section of yaml unless LinearObsOperator is specified
272  std::string confname = "obs operator";
273  if (conf.has("linear obs operator")) confname = "linear obs operator";
274  eckit::LocalConfiguration linobsopconf(conf, confname);
275  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
276 
277  const double tol = conf.getDouble("linear obs operator test.tolerance TL");
278  const double alpha = conf.getDouble("linear obs operator test.coef TL", 0.1);
279  const int iter = conf.getInt("linear obs operator test.iterations TL", 1);
280 
281  // initialize obs bias from file
282  eckit::LocalConfiguration biasconf = conf.getSubConfiguration("obs bias");
283  typename ObsAuxCtrl_::Parameters_ biasparams;
284  biasparams.validateAndDeserialize(biasconf);
285  const ObsAuxCtrl_ ybias0(Test_::obspace()[jj], biasparams);
286  ObsAuxCtrl_ ybias(Test_::obspace()[jj], biasparams);
287 
288  // initialize Obs. Bias Covariance
289  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], biasparams);
290 
291  // read geovals from the file
292  const eckit::LocalConfiguration gconf(conf, "geovals");
293  oops::Variables hopvars = hop.requiredVars();
294  hopvars += ybias0.requiredVars();
295  const GeoVaLs_ x0(gconf, Test_::obspace()[jj], hopvars);
296  GeoVaLs_ x(gconf, Test_::obspace()[jj], hopvars);
297 
298  // set TL trajectory to the geovals and the bias coeff. from the files
299  hoptl.setTrajectory(x0, ybias0);
300 
301  // create obsvectors
302  ObsVector_ y1(Test_::obspace()[jj]);
303  ObsVector_ y2(Test_::obspace()[jj]);
304  ObsVector_ y3(Test_::obspace()[jj]);
305 
306  // create obsdatavector to hold diags
307  oops::Variables diagvars;
308  diagvars += ybias0.requiredHdiagnostics();
309  ObsDiags_ ydiag(Test_::obspace()[jj], hop.locations(), diagvars);
310 
311  // y1 = hop(x0, ybias0)
312  hop.simulateObs(x0, y1, ybias0, ydiag);
313 
314  // randomize dx and ybinc
315  GeoVaLs_ dx(gconf, Test_::obspace()[jj], hoptl.requiredVars());
316  dx.random();
317  ObsAuxIncr_ ybinc(Test_::obspace()[jj], biasparams);
318  Bobsbias.randomize(ybinc);
319 
320  // scale dx by x0
321  dx *= x0;
322 
323  for (int jter = 0; jter < iter; ++jter) {
324  // x = x0 + alpha*dx
325  dx *= alpha;
326  x = x0;
327  x += dx;
328  // ybias = ybias0 + alpha*ybinc
329  ybinc *= alpha;
330  ybias = ybias0;
331  ybias += ybinc;
332 
333  // y2 = hop(x0+alpha*dx, ybias0+alpha*ybinc)
334  hop.simulateObs(x, y2, ybias, ydiag);
335  y2 -= y1;
336  // y3 = hoptl(alpha*dx, alpha*ybinc)
337  hoptl.simulateObsTL(dx, y3, ybinc);
338  y2 -= y3;
339 
340  double test_norm = y2.rms();
341  oops::Log::test() << "Iter:" << jter << " ||(h(x+alpha*dx)-h(x)-h'*(alpha*dx))||="
342  << test_norm << std::endl;
343  }
344  EXPECT(y2.rms() < tol);
345  }
346 }
347 
348 // -----------------------------------------------------------------------------
349 
350 template <typename OBS> void testException() {
351  typedef ObsTestsFixture<OBS> Test_;
352  typedef oops::GeoVaLs<OBS> GeoVaLs_;
353  typedef oops::ObsOperator<OBS> ObsOperator_;
354  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
355  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
356  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
357  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
358  typedef oops::ObsVector<OBS> ObsVector_;
359 
360  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
361  const eckit::LocalConfiguration & conf = Test_::config(jj);
362  if (conf.has(expectConstructorToThrow))
363  continue;
364 
365  // Set up objects prior to throwing exceptions.
366  eckit::LocalConfiguration obsopconf(conf, "obs operator");
367  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
368  std::string confname = "obs operator";
369  if (conf.has("linear obs operator")) confname = "linear obs operator";
370  eckit::LocalConfiguration linobsopconf(conf, confname);
371  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
372  eckit::LocalConfiguration biasconf = conf.getSubConfiguration("obs bias");
373  typename ObsAuxCtrl_::Parameters_ biasparams;
374  biasparams.validateAndDeserialize(biasconf);
375  const ObsAuxCtrl_ ybias(Test_::obspace()[jj], biasparams);
376  ObsAuxIncr_ ybinc(Test_::obspace()[jj], biasparams);
377  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], biasparams);
378  eckit::LocalConfiguration gconf(conf, "geovals");
379  oops::Variables hopvars = hop.requiredVars();
380  hopvars += ybias.requiredVars();
381  const GeoVaLs_ gval(gconf, Test_::obspace()[jj], hopvars);
382  oops::Variables diagvars;
383  diagvars += ybias.requiredHdiagnostics();
384 
385  if (Test_::config(jj).has(expectSetTrajectoryToThrow)) {
386  // The setTrajectory method is expected to throw an exception
387  // containing the specified string.
388  const std::string expectedMessage =
389  Test_::config(jj).getString(expectSetTrajectoryToThrow);
390  EXPECT_THROWS_MSG(hoptl.setTrajectory(gval, ybias),
391  expectedMessage.c_str());
392  // Do not continue further because setTrajectory must be run
393  // before simulateObsTL and simulateObsAD.
394  continue;
395  }
396 
397  if (Test_::config(jj).has(expectSimulateObsTLToThrow)) {
398  hoptl.setTrajectory(gval, ybias);
399  ObsVector_ dy1(Test_::obspace()[jj]);
400  GeoVaLs_ dx1(gconf, Test_::obspace()[jj], hoptl.requiredVars());
401  dx1.random();
402  Bobsbias.randomize(ybinc);
403  // The simulateObsTL method is expected to throw an exception
404  // containing the specified string.
405  const std::string expectedMessage =
406  Test_::config(jj).getString(expectSimulateObsTLToThrow);
407  EXPECT_THROWS_MSG(hoptl.simulateObsTL(dx1, dy1, ybinc),
408  expectedMessage.c_str());
409  }
410 
411  if (Test_::config(jj).has(expectSimulateObsADToThrow)) {
412  hoptl.setTrajectory(gval, ybias);
413  ObsVector_ dy2(Test_::obspace()[jj]);
414  GeoVaLs_ dx2(gconf, Test_::obspace()[jj], hoptl.requiredVars());
415  Bobsbias.randomize(ybinc);
416  dy2.random();
417  dx2.zero();
418  ybinc.zero();
419  // The simulateObsAD method is expected to throw an exception
420  // containing the specified string.
421  const std::string expectedMessage =
422  Test_::config(jj).getString(expectSimulateObsADToThrow);
423  EXPECT_THROWS_MSG(hoptl.simulateObsAD(dx2, dy2, ybinc),
424  expectedMessage.c_str());
425  }
426  }
427 }
428 
429 // -----------------------------------------------------------------------------
430 
431 template <typename OBS>
432 class LinearObsOperator : public oops::Test {
433  typedef ObsTestsFixture<OBS> Test_;
434 
435  public:
436  LinearObsOperator() {}
437  virtual ~LinearObsOperator() {}
438 
439  private:
440  std::string testid() const override {return "test::LinearObsOperator<" + OBS::name() + ">";}
441 
442  void register_tests() const override {
443  std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
444 
445  ts.emplace_back(CASE("interface/LinearObsOperator/testConstructor")
446  { testConstructor<OBS>(); });
447  ts.emplace_back(CASE("interface/LinearObsOperator/testLinearity")
448  { testLinearity<OBS>(); });
449  ts.emplace_back(CASE("interface/LinearObsOperator/testTangentLinear")
450  { testTangentLinear<OBS>(); });
451  ts.emplace_back(CASE("interface/LinearObsOperator/testAdjoint")
452  { testAdjoint<OBS>(); });
453  ts.emplace_back(CASE("interface/LinearObsOperator/testException")
454  { testException<OBS>(); });
455  }
456 
457  void clear() const override {
458  Test_::reset();
459  }
460 };
461 
462 // -----------------------------------------------------------------------------
463 
464 } // namespace test
465 
466 #endif // TEST_INTERFACE_LINEAROBSOPERATOR_H_
test
program test
Definition: OpenMPTryFlag.f90:2
test::LinearObsOperator::testid
std::string testid() const override
Definition: test/interface/LinearObsOperator.h:440
03-VariablesIntro.name
name
Definition: 03-VariablesIntro.py:74
oops_variables_mod::has
logical function has(this, var)
Definition: variables_mod.F90:164
test::expectConstructorToThrow
const char * expectConstructorToThrow
Definition: test/interface/LinearObsOperator.h:35
test::expectSimulateObsTLToThrow
const char * expectSimulateObsTLToThrow
Definition: test/interface/LinearObsOperator.h:38
test::expectSetTrajectoryToThrow
const char * expectSetTrajectoryToThrow
Definition: test/interface/LinearObsOperator.h:36
test::expectSimulateObsToThrow
const char * expectSimulateObsToThrow
Definition: test/interface/LinearObsOperator.h:37
test::expectSimulateObsADToThrow
const char * expectSimulateObsADToThrow
Definition: test/interface/LinearObsOperator.h:39
test::CASE
CASE("assimilation/FullGMRES/FullGMRES")
Definition: test/assimilation/FullGMRES.h:56
test::testConstructor
void testConstructor()
Tests creation and destruction of ObsErrorCovariances.
Definition: ObsErrorCovariance.h:32