OOPS
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/Logger.h"
29 #include "test/interface/ObsTestsFixture.h"
30 #include "test/TestEnvironment.h"
31 
32 namespace test {
33 
34 // -----------------------------------------------------------------------------
35 
36 template <typename OBS> void testConstructor() {
37  typedef ObsTestsFixture<OBS> Test_;
38  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
39 
40  std::vector<eckit::LocalConfiguration> conf;
41  TestEnvironment::config().get("observations", conf);
42 
43  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
44  // Use ObsOperator section of yaml unless LinearObsOperator is specified
45  std::string confname = "obs operator";
46  if (conf[jj].has("linear obs operator")) confname = "linear obs operator";
47  eckit::LocalConfiguration linobsopconf(conf[jj], confname);
48  std::unique_ptr<LinearObsOperator_> ov(
49  new LinearObsOperator_(Test_::obspace()[jj], linobsopconf));
50  EXPECT(ov.get());
51 
52  ov.reset();
53  EXPECT(!ov.get());
54  }
55 }
56 
57 // -----------------------------------------------------------------------------
58 
59 template <typename OBS> void testLinearity() {
60  typedef ObsTestsFixture<OBS> Test_;
61  typedef oops::GeoVaLs<OBS> GeoVaLs_;
62  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
63  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
64  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
65  typedef oops::ObsOperator<OBS> ObsOperator_;
66  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
67  typedef oops::ObsVector<OBS> ObsVector_;
68 
69  const double zero = 0.0;
70  const double coef = 3.14;
71  const double tol = 1.0e-11;
72  std::vector<eckit::LocalConfiguration> conf;
73  TestEnvironment::config().get("observations", conf);
74 
75  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
76  // initialize observation operator (set variables requested from the model,
77  // variables simulated by the observation operator, other init)
78  eckit::LocalConfiguration obsopconf(conf[jj], "obs operator");
79  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
80  // initialize TL/AD observation operator (set model variables for Jacobian),
81  // other init)
82  // Use ObsOperator section of yaml unless LinearObsOperator is specified
83  std::string confname = "obs operator";
84  if (conf[jj].has("linear obs operator")) confname = "linear obs operator";
85  eckit::LocalConfiguration linobsopconf(conf[jj], confname);
86  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
87 
88  // initialize obs bias
89  const ObsAuxCtrl_ ybias(Test_::obspace()[jj], conf[jj]);
90  ObsAuxIncr_ ybinc(Test_::obspace()[jj], conf[jj]);
91 
92  // read geovals from the file
93  const eckit::LocalConfiguration gconf(conf[jj], "geovals");
94  oops::Variables hopvars = hop.requiredVars();
95  hopvars += ybias.requiredVars();
96  const GeoVaLs_ gval(gconf, Test_::obspace()[jj], hopvars);
97 
98  // initialize Obs. Bias Covariance
99  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], conf[jj]);
100 
101  // set trajectory for TL/AD to be the geovals from the file
102  hoptl.setTrajectory(gval, ybias);
103 
104  // create obsvector
105  ObsVector_ dy1(Test_::obspace()[jj]);
106 
107  // create geovals
108  GeoVaLs_ dx(gconf, Test_::obspace()[jj], hoptl.requiredVars());
109 
110  // test rms(H * (dx, ybinc)) = 0, when dx = 0
111  dx.zero();
112  ybinc.zero();
113  hoptl.simulateObsTL(dx, dy1, ybinc);
114  EXPECT(dy1.rms() == zero);
115 
116  // test rms(H * (dx, ybinc)) > 0, when dx is random
117  dx.random();
118  Bobsbias.randomize(ybinc);
119  hoptl.simulateObsTL(dx, dy1, ybinc);
120  EXPECT(dy1.rms() > zero);
121 
122  // test k * H * (dx, ybinc) ~ H * (k*dx, k*ybinc)
123  dy1 *= coef;
124  dx *= coef;
125  ybinc *= coef;
126  ObsVector_ dy2(Test_::obspace()[jj]);
127  hoptl.simulateObsTL(dx, dy2, ybinc);
128 
129  dy2 -= dy1;
130  EXPECT(dy2.rms() / dy1.rms() < tol);
131  }
132 }
133 
134 // -----------------------------------------------------------------------------
135 
136 template <typename OBS> void testAdjoint() {
137  typedef ObsTestsFixture<OBS> Test_;
138  typedef oops::GeoVaLs<OBS> GeoVaLs_;
139  typedef oops::ObsOperator<OBS> ObsOperator_;
140  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
141  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
142  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
143  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
144  typedef oops::ObsVector<OBS> ObsVector_;
145 
146  const double zero = 0.0;
147  std::vector<eckit::LocalConfiguration> conf;
148  TestEnvironment::config().get("observations", conf);
149 
150  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
151  // initialize observation operator (set variables requested from the model,
152  // variables simulated by the observation operator, other init)
153  eckit::LocalConfiguration obsopconf(conf[jj], "obs operator");
154  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
155  // initialize TL/AD observation operator (set model variables for Jacobian),
156  // other init)
157  // Use ObsOperator section of yaml unless LinearObsOperator is specified
158  std::string confname = "obs operator";
159  if (conf[jj].has("linear obs operator")) confname = "linear obs operator";
160  eckit::LocalConfiguration linobsopconf(conf[jj], confname);
161  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
162 
163  const double tol = conf[jj].getDouble("linear obs operator test.tolerance AD");
164  // initialize bias correction
165  const ObsAuxCtrl_ ybias(Test_::obspace()[jj], conf[jj]);
166  ObsAuxIncr_ ybinc1(Test_::obspace()[jj], conf[jj]); // TL
167  ObsAuxIncr_ ybinc2(Test_::obspace()[jj], conf[jj]); // AD
168 
169  // initialize Obs. Bias Covariance
170  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], conf[jj]);
171 
172  // read geovals from the file
173  eckit::LocalConfiguration gconf(conf[jj], "geovals");
174  oops::Variables hopvars = hop.requiredVars();
175  hopvars += ybias.requiredVars();
176  const GeoVaLs_ gval(gconf, Test_::obspace()[jj], hopvars);
177 
178  // set TL/AD trajectory to the geovals from the file
179  hoptl.setTrajectory(gval, ybias);
180 
181  ObsVector_ dy1(Test_::obspace()[jj]);
182  ObsVector_ dy2(Test_::obspace()[jj]);
183  GeoVaLs_ dx1(gconf, Test_::obspace()[jj], hoptl.requiredVars());
184  GeoVaLs_ dx2(gconf, Test_::obspace()[jj], hoptl.requiredVars());
185 
186  // calculate dy1 = H (dx1, ybinc1) (with random dx1, and random ybinc1)
187  dx1.random();
188  EXPECT(dot_product(dx1, dx1) > zero); // BOOST_REQUIRE
189  Bobsbias.randomize(ybinc1);
190  hoptl.simulateObsTL(dx1, dy1, ybinc1);
191  EXPECT(dot_product(dy1, dy1) > zero);
192 
193  // calculate (dx2, ybinc2) = HT dy2 (with random dy2)
194  dy2.random();
195  EXPECT(dot_product(dy2, dy2) > zero); // BOOST_REQUIRE
196  dx2.zero();
197  ybinc2.zero();
198  hoptl.simulateObsAD(dx2, dy2, ybinc2);
199  EXPECT(dot_product(dx2, dx2) > zero);
200 
201  const double zz1 = dot_product(dx1, dx2) + dot_product(ybinc1, ybinc2);
202  const double zz2 = dot_product(dy1, dy2);
203 
204  oops::Log::info() << "Adjoint test result: (<x,HTy>-<Hx,y>)/<Hx,y> = "
205  << (zz1-zz2)/zz2 << std::endl;
206 
207  EXPECT(zz1 != zero);
208  EXPECT(zz2 != zero);
209  EXPECT(oops::is_close(zz1, zz2, tol));
210  }
211 }
212 
213 // -----------------------------------------------------------------------------
214 
215 template <typename OBS> void testTangentLinear() {
216  // Test ||(hop(x+alpha*dx)-hop(x)) - hoptl(alpha*dx)|| < tol
217  typedef ObsTestsFixture<OBS> Test_;
218  typedef oops::GeoVaLs<OBS> GeoVaLs_;
219  typedef oops::ObsDiagnostics<OBS> ObsDiags_;
220  typedef oops::ObsAuxControl<OBS> ObsAuxCtrl_;
221  typedef oops::ObsAuxIncrement<OBS> ObsAuxIncr_;
222  typedef oops::ObsAuxCovariance<OBS> ObsAuxCov_;
223  typedef oops::LinearObsOperator<OBS> LinearObsOperator_;
224  typedef oops::ObsOperator<OBS> ObsOperator_;
225  typedef oops::ObsVector<OBS> ObsVector_;
226 
227  std::vector<eckit::LocalConfiguration> conf;
228  TestEnvironment::config().get("observations", conf);
229 
230  for (std::size_t jj = 0; jj < Test_::obspace().size(); ++jj) {
231  // initialize observation operator (set variables requested from the model,
232  // variables simulated by the observation operator, other init)
233  eckit::LocalConfiguration obsopconf(conf[jj], "obs operator");
234  ObsOperator_ hop(Test_::obspace()[jj], obsopconf);
235  // initialize TL/AD observation operator (set model variables for Jacobian),
236  // other init)
237  // Use ObsOperator section of yaml unless LinearObsOperator is specified
238  std::string confname = "obs operator";
239  if (conf[jj].has("linear obs operator")) confname = "linear obs operator";
240  eckit::LocalConfiguration linobsopconf(conf[jj], confname);
241  LinearObsOperator_ hoptl(Test_::obspace()[jj], linobsopconf);
242 
243  const double tol = conf[jj].getDouble("linear obs operator test.tolerance TL");
244  const double alpha = conf[jj].getDouble("linear obs operator test.coef TL", 0.1);
245  const int iter = conf[jj].getInt("linear obs operator test.iterations TL", 1);
246 
247  // initialize obs bias from file
248  const ObsAuxCtrl_ ybias0(Test_::obspace()[jj], conf[jj]);
249  ObsAuxCtrl_ ybias(Test_::obspace()[jj], conf[jj]);
250 
251  // initialize Obs. Bias Covariance
252  const ObsAuxCov_ Bobsbias(Test_::obspace()[jj], conf[jj]);
253 
254  // read geovals from the file
255  const eckit::LocalConfiguration gconf(conf[jj], "geovals");
256  oops::Variables hopvars = hop.requiredVars();
257  hopvars += ybias0.requiredVars();
258  const GeoVaLs_ x0(gconf, Test_::obspace()[jj], hopvars);
259  GeoVaLs_ x(gconf, Test_::obspace()[jj], hopvars);
260 
261  // set TL trajectory to the geovals and the bias coeff. from the files
262  hoptl.setTrajectory(x0, ybias0);
263 
264  // create obsvectors
265  ObsVector_ y1(Test_::obspace()[jj]);
266  ObsVector_ y2(Test_::obspace()[jj]);
267  ObsVector_ y3(Test_::obspace()[jj]);
268 
269  // create obsdatavector to hold diags
270  oops::Variables diagvars;
271  diagvars += ybias0.requiredHdiagnostics();
272  ObsDiags_ ydiag(Test_::obspace()[jj],
273  hop.locations(Test_::obspace()[jj].windowStart(),
274  Test_::obspace()[jj].windowEnd()),
275  diagvars);
276 
277  // y1 = hop(x0, ybias0)
278  hop.simulateObs(x0, y1, ybias0, ydiag);
279 
280  // randomize dx and ybinc
281  GeoVaLs_ dx(gconf, Test_::obspace()[jj], hoptl.requiredVars());
282  dx.random();
283  ObsAuxIncr_ ybinc(Test_::obspace()[jj], conf[jj]);
284  Bobsbias.randomize(ybinc);
285 
286  // scale dx by x0
287  dx *= x0;
288 
289  for (int jter = 0; jter < iter; ++jter) {
290  // x = x0 + alpha*dx
291  dx *= alpha;
292  x = x0;
293  x += dx;
294  // ybias = ybias0 + alpha*ybinc
295  ybinc *= alpha;
296  ybias = ybias0;
297  ybias += ybinc;
298 
299  // y2 = hop(x0+alpha*dx, ybias0+alpha*ybinc)
300  hop.simulateObs(x, y2, ybias, ydiag);
301  y2 -= y1;
302  // y3 = hoptl(alpha*dx, alpha*ybinc)
303  hoptl.simulateObsTL(dx, y3, ybinc);
304  y2 -= y3;
305 
306  double test_norm = y2.rms();
307  oops::Log::info() << "Iter:" << jter << " ||(h(x+alpha*dx)-h(x)-h'*(alpha*dx))||="
308  << test_norm << std::endl;
309  }
310  EXPECT(y2.rms() < tol);
311  }
312 }
313 
314 // -----------------------------------------------------------------------------
315 
316 template <typename OBS>
317 class LinearObsOperator : public oops::Test {
318  typedef ObsTestsFixture<OBS> Test_;
319  public:
320  LinearObsOperator() {}
321  virtual ~LinearObsOperator() {}
322  private:
323  std::string testid() const override {return "test::LinearObsOperator<" + OBS::name() + ">";}
324 
325  void register_tests() const override {
326  std::vector<eckit::testing::Test>& ts = eckit::testing::specification();
327 
328  ts.emplace_back(CASE("interface/GeometryIterator/testConstructor")
329  { testConstructor<OBS>(); });
330  ts.emplace_back(CASE("interface/GeometryIterator/testLinearity")
331  { testLinearity<OBS>(); });
332  ts.emplace_back(CASE("interface/GeometryIterator/testTangentLinear")
333  { testTangentLinear<OBS>(); });
334  ts.emplace_back(CASE("interface/GeometryIterator/testAdjoint")
335  { testAdjoint<OBS>(); });
336  }
337 
338  void clear() const override {
339  Test_::reset();
340  }
341 };
342 
343 // -----------------------------------------------------------------------------
344 
345 } // namespace test
346 
347 #endif // TEST_INTERFACE_LINEAROBSOPERATOR_H_
test::ObsTestsFixture
Definition: ObsTestsFixture.h:29
test::LinearObsOperator::clear
void clear() const override
Definition: test/interface/LinearObsOperator.h:338
ObsAuxIncrement.h
test::testAdjoint
void testAdjoint()
Definition: test/interface/LinearObsOperator.h:136
ObsAuxCovariance.h
ObsAuxControl.h
test::ObsTestsFixture::reset
static void reset()
Definition: ObsTestsFixture.h:37
test::LinearObsOperator::Test_
ObsTestsFixture< OBS > Test_
Definition: test/interface/LinearObsOperator.h:318
oops::ObsAuxIncrement
Definition: oops/interface/ObsAuxIncrement.h:38
test::CASE
CASE("test_linearmodelparameterswrapper_valid_name")
Definition: LinearModelFactory.cc:22
ObsOperator.h
test
Definition: LinearModelFactory.cc:20
oops::ObsAuxControl
Definition: oops/interface/ObsAuxControl.h:35
oops::ObsVector
Definition: oops/interface/ObsSpace.h:36
oops::ObsAuxCovariance
Definition: oops/interface/ObsAuxCovariance.h:36
ObsDiagnostics.h
LinearObsOperator.h
oops::ObsOperator
Definition: oops/interface/ObsOperator.h:38
Test.h
test::TestEnvironment::config
static const eckit::Configuration & config()
Definition: TestEnvironment.h:40
ObsTestsFixture.h
test::testLinearity
void testLinearity()
Definition: test/interface/LinearObsOperator.h:59
test::LinearObsOperator::testid
std::string testid() const override
Definition: test/interface/LinearObsOperator.h:323
test::LinearObsOperator::LinearObsOperator
LinearObsOperator()
Definition: test/interface/LinearObsOperator.h:320
test::testTangentLinear
void testTangentLinear()
Definition: test/interface/LinearObsOperator.h:215
oops::ObsDiagnostics
Definition: ObsDiagnostics.h:33
TestEnvironment.h
oops::LinearObsOperator
Definition: oops/interface/LinearObsOperator.h:38
test::testConstructor
void testConstructor()
Tests creation and destruction of ObsErrorCovariances.
Definition: test/base/ObsErrorCovariance.h:32
test::LinearObsOperator
Definition: test/interface/LinearObsOperator.h:317
oops_variables_mod::has
logical function has(this, var)
Definition: variables_mod.F90:140
test::LinearObsOperator::register_tests
void register_tests() const override
Definition: test/interface/LinearObsOperator.h:325
oops::Test
Definition: Test.h:39
oops::Variables
Definition: oops/base/Variables.h:23
test::LinearObsOperator::~LinearObsOperator
virtual ~LinearObsOperator()
Definition: test/interface/LinearObsOperator.h:321
oops::GeoVaLs
Definition: oops/interface/GeoVaLs.h:32
Variables.h