ReFRACtor
nlls_solver_gsl_sm_helical_valley.cc
Go to the documentation of this file.
1 #include <bard_nlls_problem.h>
2 #include <brown_nlls_problem.h>
3 #include <freudenstein_roth_nlls_problem.h>
4 #include <helical_valley_nlls_problem.h>
5 #include <jennrich_sampson_nlls_problem.h>
6 #include <meyer_nlls_problem.h>
7 #include <powell_nlls_problem.h>
8 #include <powell_singular_nlls_problem.h>
9 #include <rosenbrock2_nlls_problem.h>
10 #include <unit_test_support.h>
11 #include <fp_exception.h>
12 #include <nlls_solver_gsl_sm.h>
13 
14 #include <bard_ml_problem.h>
15 #include <meyer_ml_problem.h>
16 #include <nlls_max_likelihood.h>
17 
18 
19 using namespace FullPhysics;
20 using namespace blitz;
21 
22 BOOST_FIXTURE_TEST_SUITE(nlls_solver_gsl_sm_freudenstein_roth_b, GlobalFixture)
23 
24 /* convergence check thresholds */
25 double x_tol=1.0e-6, g_tol=6.0555e-06, f_tol=1.0e-6;
26 bool verbose=false;
27 
28 
29 
30 BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__default)
31 {
32  // Create the problem
33  //
34  Array<double, 1> x0(2); x0 = 6.0, 7.0;
35  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
36 
37  // Set the problem at the initial guess.
38  //
39  nlls->parameters(x0);
40 
41  // Create the solver with almost all default setting,
42  // and solve the NLLS problem.
43  //
44  NLLSSolverGSLSM solver(nlls, 100);
45  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::UNTRIED);
46  solver.solve();
47 
48  int n_f_calls = nlls->num_residual_evaluations();
49  int n_j_calls = nlls->num_jacobian_evaluations();
50  double cst = nlls->cost();
51 
52  std::cout
53  << "Testing NLLSSolverGSLSM with Freudenstein/Roth function:" << std::endl
54  << " Number of residual function evaluations = " << n_f_calls << std::endl
55  << " Number of jacobian function evaluations = " << n_j_calls << std::endl
56  << " Final solver status = " << solver.status_str() << std::endl
57  << " Final problem status (point) = " << nlls->parameters() << std::endl
58  << " Final problem status (cost value) = " << cst << std::endl
59  << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
60  for( int i=0; i<=solver.num_accepted_steps(); i++ )
61  std::cout
62  << " ========================================" << std::endl
63  << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
64  << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
65  << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
66 
67  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
68  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
69  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
70 
71  int iLast = solver.num_accepted_steps();
72  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
73  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
74  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
75 
76  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::SUCCESS);
77  BOOST_CHECK(n_f_calls < 10);
78  BOOST_CHECK(n_j_calls <= n_f_calls);
79  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
80  if( abs(cst-0.0) < 0.0000001 ) {
81  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
82  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
83  } else {
84  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
85  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
86  }
87 }
88 
89 
90 BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__lm_more_svd)
91 {
92  // Create the problem
93  //
94  Array<double, 1> x0(2); x0 = 6.0, 7.0;
95  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
96 
97  // Set the problem at the initial guess.
98  //
99  nlls->parameters(x0);
100 
101  // Create the solver and solve the NLLS problem.
102  // Solver is set to the Lev/Mar trust region method
103  // with More's scaling, and it will use SVD for solving
104  // linear systems.
105  //
106  gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
107  fdf_params.trs = gsl_multifit_nlinear_trs_lm;
108  fdf_params.scale = gsl_multifit_nlinear_scale_more;
109  fdf_params.solver = gsl_multifit_nlinear_solver_svd;
110  //
111  NLLSSolverGSLSM solver(nlls, 100, fdf_params);
112  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::UNTRIED);
113  solver.solve();
114 
115  int n_f_calls = nlls->num_residual_evaluations();
116  int n_j_calls = nlls->num_jacobian_evaluations();
117  double cst = nlls->cost();
118 
119  std::cout
120  << "Testing NLLSSolverGSLSM with Freudenstein/Roth function:" << std::endl
121  << " Number of residual function evaluations = " << n_f_calls << std::endl
122  << " Number of jacobian function evaluations = " << n_j_calls << std::endl
123  << " Final solver status = " << solver.status_str() << std::endl
124  << " Final problem status (point) = " << nlls->parameters() << std::endl
125  << " Final problem status (cost value) = " << cst << std::endl
126  << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
127  for( int i=0; i<=solver.num_accepted_steps(); i++ )
128  std::cout
129  << " ========================================" << std::endl
130  << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
131  << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
132  << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
133 
134  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
135  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
136  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
137 
138  int iLast = solver.num_accepted_steps();
139  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
140  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
141  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
142 
143  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::SUCCESS);
144  BOOST_CHECK(n_f_calls < 10);
145  BOOST_CHECK(n_j_calls <= n_f_calls);
146  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
147  if( abs(cst-0.0) < 0.0000001 ) {
148  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
149  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
150  } else {
151  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
152  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
153  }
154 }
155 
156 
157 BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__lmaccel_more_svd)
158 {
159  // Create the problem
160  //
161  Array<double, 1> x0(2); x0 = 6.0, 7.0;
162  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
163 
164  // Set the problem at the initial guess.
165  //
166  nlls->parameters(x0);
167 
168  // Create the solver and solve the NLLS problem. Solver
169  // is set to the Lev/Mar accelerated trust region method
170  // with More's scaling, and it will use SVD for solving
171  // linear systems.
172  //
173  gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
174  fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
175  fdf_params.scale = gsl_multifit_nlinear_scale_more;
176  fdf_params.solver = gsl_multifit_nlinear_solver_svd;
177  fdf_params.fdtype = GSL_MULTIFIT_NLINEAR_CTRDIFF;
178  //
179  NLLSSolverGSLSM solver(nlls, 100, fdf_params);
180  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::UNTRIED);
181  solver.solve();
182 
183  int n_f_calls = nlls->num_residual_evaluations();
184  int n_j_calls = nlls->num_jacobian_evaluations();
185  double cst = nlls->cost();
186 
187  std::cout
188  << "Testing NLLSSolverGSLSM with Freudenstein/Roth function:" << std::endl
189  << " Number of residual function evaluations = " << n_f_calls << std::endl
190  << " Number of jacobian function evaluations = " << n_j_calls << std::endl
191  << " Final solver status = " << solver.status_str() << std::endl
192  << " Final problem status (point) = " << nlls->parameters() << std::endl
193  << " Final problem status (cost value) = " << cst << std::endl
194  << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
195  for( int i=0; i<=solver.num_accepted_steps(); i++ )
196  std::cout
197  << " ========================================" << std::endl
198  << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
199  << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
200  << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
201 
202  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
203  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
204  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
205 
206  int iLast = solver.num_accepted_steps();
207  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
208  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
209  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
210 
211  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::SUCCESS);
212  BOOST_CHECK(n_f_calls < 25);
213  BOOST_CHECK(n_j_calls <= n_f_calls);
214  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
215  if( abs(cst-0.0) < 0.0000001 ) {
216  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
217  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
218  } else {
219  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
220  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
221  }
222 }
223 
224 
225 BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__subspace2D_more_svd)
226 {
227  // Create the problem
228  //
229  Array<double, 1> x0(2); x0 = 6.0, 7.0;
230  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
231 
232  // Set the problem at the initial guess.
233  //
234  nlls->parameters(x0);
235 
236  // Create the solver and solve the NLLS problem.
237  // Solver is set to the Subspace2D trust region method
238  // with More's scaling, and it will use SVD for solving
239  // linear systems.
240  //
241  gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
242  fdf_params.trs = gsl_multifit_nlinear_trs_subspace2D;
243  fdf_params.scale = gsl_multifit_nlinear_scale_more;
244  fdf_params.solver = gsl_multifit_nlinear_solver_svd;
245  //
246  NLLSSolverGSLSM solver(nlls, 100, fdf_params);
247  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::UNTRIED);
248  solver.solve();
249 
250  int n_f_calls = nlls->num_residual_evaluations();
251  int n_j_calls = nlls->num_jacobian_evaluations();
252  double cst = nlls->cost();
253 
254  std::cout
255  << "Testing NLLSSolverGSLSM with Freudenstein/Roth function:" << std::endl
256  << " Number of residual function evaluations = " << n_f_calls << std::endl
257  << " Number of jacobian function evaluations = " << n_j_calls << std::endl
258  << " Final solver status = " << solver.status_str() << std::endl
259  << " Final problem status (point) = " << nlls->parameters() << std::endl
260  << " Final problem status (cost value) = " << cst << std::endl
261  << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
262  for( int i=0; i<=solver.num_accepted_steps(); i++ )
263  std::cout
264  << " ========================================" << std::endl
265  << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
266  << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
267  << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
268 
269  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
270  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
271  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
272 
273  int iLast = solver.num_accepted_steps();
274  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
275  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
276  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
277 
278  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::SUCCESS);
279  BOOST_CHECK(n_f_calls < 10);
280  BOOST_CHECK(n_j_calls <= n_f_calls);
281  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
282  if( abs(cst-0.0) < 0.0000001 ) {
283  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
284  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
285  } else {
286  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
287  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
288  }
289 }
290 
291 
292 BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__ddogleg_more_svd)
293 {
294  // Create the problem
295  //
296  Array<double, 1> x0(2); x0 = 6.0, 7.0;
297  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
298 
299  // Set the problem at the initial guess.
300  //
301  nlls->parameters(x0);
302 
303  // Create the solver and solve the NLLS problem.
304  // Solver is set to the DDogLeg trust region method
305  // with More's scaling, and it will use SVD for solving
306  // linear systems.
307  //
308  gsl_multifit_nlinear_parameters fdf_params=gsl_multifit_nlinear_default_parameters();
309  fdf_params.trs = gsl_multifit_nlinear_trs_ddogleg;
310  fdf_params.scale = gsl_multifit_nlinear_scale_more;
311  fdf_params.solver = gsl_multifit_nlinear_solver_svd;
312  //
313  NLLSSolverGSLSM solver(nlls, 100, fdf_params);
314  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::UNTRIED);
315  solver.solve();
316 
317  int n_f_calls = nlls->num_residual_evaluations();
318  int n_j_calls = nlls->num_jacobian_evaluations();
319  double cst = nlls->cost();
320 
321  std::cout
322  << "Testing NLLSSolverGSLSM with Freudenstein/Roth function:" << std::endl
323  << " Number of residual function evaluations = " << n_f_calls << std::endl
324  << " Number of jacobian function evaluations = " << n_j_calls << std::endl
325  << " Final solver status = " << solver.status_str() << std::endl
326  << " Final problem status (point) = " << nlls->parameters() << std::endl
327  << " Final problem status (cost value) = " << cst << std::endl
328  << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
329  for( int i=0; i<=solver.num_accepted_steps(); i++ )
330  std::cout
331  << " ========================================" << std::endl
332  << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
333  << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
334  << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
335 
336  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
337  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
338  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
339 
340  int iLast = solver.num_accepted_steps();
341  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
342  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
343  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
344 
345  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSLSM::SUCCESS);
346  BOOST_CHECK(n_f_calls < 10);
347  BOOST_CHECK(n_j_calls <= n_f_calls);
348  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
349  if( abs(cst-0.0) < 0.0000001 ) {
350  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
351  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
352  } else {
353  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
354  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
355  }
356 }
357 
358 
359 /*
360 BOOST_AUTO_TEST_CASE(freudenstein_roth__b)
361 {
362  Array<double, 1> x0(2); x0 = 6.0, 7.0;
363  boost::shared_ptr<FreudensteinRothNLLSProblem> nlls(new FreudensteinRothNLLSProblem);
364  nlls->parameters(x0);
365  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
366  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
367  solver.solve();
368 
369  int n_f_calls = nlls->num_residual_evaluations();
370  int n_j_calls = nlls->num_jacobian_evaluations();
371  double cst = nlls->cost();
372 
373 // std::cout
374 // << "Testing NLLSSolverGSL with Freudenstein/Roth function:" << std::endl
375 // << " Number of residual function evaluations = " << n_f_calls << std::endl
376 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
377 // << " Final solver status = " << solver.status_str() << std::endl
378 // << " Final problem status (point) = " << nlls->parameters() << std::endl
379 // << " Final problem status (cost value) = " << cst << std::endl
380 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
381 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
382 // std::cout
383 // << " ========================================" << std::endl
384 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
385 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
386 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
387 
388  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
389  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
390  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
391 
392  int iLast = solver.num_accepted_steps();
393  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
394  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
395  BOOST_CHECK_CLOSE(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast])), 0.0, 1e-12);
396 
397  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
398  BOOST_CHECK(n_f_calls < 10);
399  BOOST_CHECK(n_j_calls <= n_f_calls);
400  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) || (abs(cst-24.4921) < 0.001) );
401  if( abs(cst-0.0) < 0.0000001 ) {
402  BOOST_CHECK_CLOSE(nlls->parameters()(0), 5.0, 0.01);
403  BOOST_CHECK_CLOSE(nlls->parameters()(1), 4.0, 0.01);
404  } else {
405  BOOST_CHECK_CLOSE(nlls->parameters()(0), 11.413, 0.01);
406  BOOST_CHECK_CLOSE(nlls->parameters()(1), -0.89681, 0.01);
407  }
408 }
409 
410 
411 BOOST_AUTO_TEST_CASE(helical_valley)
412 {
413  Array<double, 1> x0(3); x0 = -1.0, 0.0, 0.0;
414  boost::shared_ptr<HelicalValleyNLLSProblem> nlls(new HelicalValleyNLLSProblem);
415  nlls->parameters(x0);
416  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
417  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
418  solver.solve();
419 
420  int n_f_calls = nlls->num_residual_evaluations();
421  int n_j_calls = nlls->num_jacobian_evaluations();
422  double cst = nlls->cost();
423 
424 // std::cout
425 // << "Testing NLLSSolverGSL with Helical Valley function:" << std::endl
426 // << " Number of residual function evaluations = " << n_f_calls << std::endl
427 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
428 // << " Final solver status = " << solver.status_str() << std::endl
429 // << " Final problem status (point) = " << nlls->parameters() << std::endl
430 // << " Final problem status (cost value) = " << cst << std::endl
431 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
432 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
433 // std::cout
434 // << " ========================================" << std::endl
435 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
436 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
437 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
438 
439  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
440  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
441  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
442 
443  int iLast = solver.num_accepted_steps();
444  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
445  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
446  BOOST_CHECK_CLOSE(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast])), 0.0, 1e-12);
447 
448  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
449  BOOST_CHECK(n_f_calls < 15);
450  BOOST_CHECK(n_j_calls <= n_f_calls);
451  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) );
452  BOOST_CHECK_CLOSE(nlls->parameters()(0), 1.0, 0.01);
453  BOOST_CHECK(abs(nlls->parameters()(1)-0.0) < 0.0000001);
454  BOOST_CHECK(abs(nlls->parameters()(2)-0.0) < 0.0000001);
455 }
456 
457 
458 BOOST_AUTO_TEST_CASE(jennrich_sampson)
459 {
460  Array<double, 1> x0(2); x0 = 0.3, 0.4;
461  boost::shared_ptr<JennrichSampsonNLLSProblem> nlls(new JennrichSampsonNLLSProblem);
462  nlls->parameters(x0);
463  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
464  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
465  solver.solve();
466 
467  int n_f_calls = nlls->num_residual_evaluations();
468  int n_j_calls = nlls->num_jacobian_evaluations();
469  double cst = nlls->cost();
470 
471 // std::cout
472 // << "Testing NLLSSolverGSL with Jennrich/Sampson function:" << std::endl
473 // << " Number of residual function evaluations = " << n_f_calls << std::endl
474 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
475 // << " Final solver status = " << solver.status_str() << std::endl
476 // << " Final problem status (point) = " << nlls->parameters() << std::endl
477 // << " Final problem status (cost value) = " << cst << std::endl
478 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
479 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
480 // std::cout
481 // << " ========================================" << std::endl
482 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
483 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
484 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
485 
486  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
487  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
488  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
489 
490  int iLast = solver.num_accepted_steps();
491  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
492  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
493  BOOST_CHECK(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast])) < 1e-12);
494 
495  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
496  BOOST_CHECK(n_f_calls < 25);
497  BOOST_CHECK(n_j_calls <= n_f_calls);
498  BOOST_CHECK_CLOSE(cst, 62.1813, 0.01);
499  BOOST_CHECK_CLOSE(nlls->parameters()(0), 0.2578, 0.01);
500  BOOST_CHECK_CLOSE(nlls->parameters()(1), 0.2578, 0.01);
501 }
502 
503 
504 BOOST_AUTO_TEST_CASE(meyer)
505 {
506  Array<double, 1> x0(3); x0 = 0.02, 4000.0, 250.0;
507  boost::shared_ptr<MeyerNLLSProblem> nlls(new MeyerNLLSProblem);
508  nlls->parameters(x0);
509  NLLSSolverGSL solver(nlls, 200, dx_tol_abs, dx_tol_rel, g_tol, verbose);
510  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
511  solver.solve();
512 
513  int n_f_calls = nlls->num_residual_evaluations();
514  int n_j_calls = nlls->num_jacobian_evaluations();
515  double cst = nlls->cost();
516 
517 // std::cout
518 // << "Testing NLLSSolverGSL with Meyer function:" << std::endl
519 // << " Number of residual function evaluations = " << n_f_calls << std::endl
520 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
521 // << " Final solver status = " << solver.status_str() << std::endl
522 // << " Final problem status (point) = " << nlls->parameters() << std::endl
523 // << " Final problem status (cost value) = " << cst << std::endl
524 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
525 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
526 // std::cout
527 // << " ========================================" << std::endl
528 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
529 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
530 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
531 
532  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
533  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
534  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
535 
536  int iLast = solver.num_accepted_steps();
537  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
538  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
539  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-8);
540 
541  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
542  BOOST_CHECK(n_f_calls < 140);
543  BOOST_CHECK(n_j_calls <= n_f_calls);
544  BOOST_CHECK_CLOSE(cst, 43.9729, 0.01);
545  BOOST_CHECK_CLOSE(nlls->parameters()(0), 0.0056096, 0.01);
546  BOOST_CHECK_CLOSE(nlls->parameters()(1), 6181.35, 0.01);
547  BOOST_CHECK_CLOSE(nlls->parameters()(2), 345.224, 0.01);
548 }
549 
550 
551 BOOST_AUTO_TEST_CASE(powell)
552 {
553  Array<double, 1> x0(2); x0 = 0.0, 1.0;
554  boost::shared_ptr<PowellNLLSProblem> nlls(new PowellNLLSProblem);
555  nlls->parameters(x0);
556  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
557  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
558  solver.solve();
559 
560  int n_f_calls = nlls->num_residual_evaluations();
561  int n_j_calls = nlls->num_jacobian_evaluations();
562  double cst = nlls->cost();
563 
564 // std::cout
565 // << "Testing NLLSSolverGSL with Powell function:" << std::endl
566 // << " Number of residual function evaluations = " << n_f_calls << std::endl
567 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
568 // << " Final solver status = " << solver.status_str() << std::endl
569 // << " Final problem status (point) = " << nlls->parameters() << std::endl
570 // << " Final problem status (cost value) = " << cst << std::endl
571 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
572 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
573 // std::cout
574 // << " ========================================" << std::endl
575 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
576 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
577 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
578 
579  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
580  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
581  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
582 
583  int iLast = solver.num_accepted_steps();
584  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
585  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
586  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
587 
588  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
589  BOOST_CHECK(n_f_calls < 20);
590  BOOST_CHECK(n_j_calls <= n_f_calls);
591  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) );
592  BOOST_CHECK_CLOSE(nlls->parameters()(0), 1.09816e-5, 0.01);
593  BOOST_CHECK_CLOSE(nlls->parameters()(1), 9.10615, 0.01);
594 }
595 
596 
597 BOOST_AUTO_TEST_CASE(powell_singular)
598 {
599  Array<double, 1> x0(4); x0 = 3.0, -1.0, 0.0, 1.0;
600  boost::shared_ptr<PowellSingularNLLSProblem> nlls(new PowellSingularNLLSProblem);
601  nlls->parameters(x0);
602  NLLSSolverGSL solver(nlls, 100, 1e-21, 1e-21, 1e-21);
603  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
604  solver.solve();
605 
606  int n_f_calls = nlls->num_residual_evaluations();
607  int n_j_calls = nlls->num_jacobian_evaluations();
608  double cst = nlls->cost();
609 
610 // std::cout
611 // << "Testing NLLSSolverGSL with Powell Singular function:" << std::endl
612 // << " Number of residual function evaluations = " << n_f_calls << std::endl
613 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
614 // << " Final solver status = " << solver.status_str() << std::endl
615 // << " Final problem status (point) = " << nlls->parameters() << std::endl
616 // << " Final problem status (cost value) = " << cst << std::endl
617 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
618 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
619 // std::cout
620 // << " ========================================" << std::endl
621 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
622 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
623 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
624 
625  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
626  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
627  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
628 
629  int iLast = solver.num_accepted_steps();
630  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
631  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
632  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
633 
634  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
635  BOOST_CHECK(n_f_calls < 30);
636  BOOST_CHECK(n_j_calls <= n_f_calls);
637  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) );
638  BOOST_CHECK(abs(nlls->parameters()(0)-0.0) < 0.0000001);
639  BOOST_CHECK(abs(nlls->parameters()(1)-0.0) < 0.0000001);
640  BOOST_CHECK(abs(nlls->parameters()(2)-0.0) < 0.0000001);
641  BOOST_CHECK(abs(nlls->parameters()(3)-0.0) < 0.0000001);
642 }
643 
644 
645 BOOST_AUTO_TEST_CASE(rosenbrock2)
646 {
647  Array<double, 1> x0(2); x0 = -1.2, 1.0;
648  boost::shared_ptr<Rosenbrock2NLLSProblem> nlls(new Rosenbrock2NLLSProblem);
649  nlls->parameters(x0);
650  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
651  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
652  solver.solve();
653 
654  int n_f_calls = nlls->num_residual_evaluations();
655  int n_j_calls = nlls->num_jacobian_evaluations();
656  double cst = nlls->cost();
657 
658 // std::cout
659 // << "Testing NLLSSolverGSL with Rosenbrock function:" << std::endl
660 // << " Number of residual function evaluations = " << n_f_calls << std::endl
661 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
662 // << " Final solver status = " << solver.status_str() << std::endl
663 // << " Final problem status (point) = " << nlls->parameters() << std::endl
664 // << " Final problem status (cost value) = " << cst << std::endl
665 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
666 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
667 // std::cout
668 // << " ========================================" << std::endl
669 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
670 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
671 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
672 
673  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
674  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
675  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
676 
677  int iLast = solver.num_accepted_steps();
678  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
679  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
680  BOOST_CHECK_CLOSE(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast])), 0.0, 1e-12);
681 
682  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
683  BOOST_CHECK(n_f_calls < 25);
684  BOOST_CHECK(n_j_calls <= n_f_calls);
685  BOOST_CHECK( (abs(cst-0.0) < 0.0000001) );
686  BOOST_CHECK_CLOSE(nlls->parameters()(0), 1.0, 0.01);
687  BOOST_CHECK_CLOSE(nlls->parameters()(1), 1.0, 0.01);
688 }
689 
690 
691 BOOST_AUTO_TEST_CASE(powell_singular_multiple)
692 {
693  Array<double, 1> x0(4); x0 = 3.0, -1.0, 0.0, 1.0;
694  boost::shared_ptr<PowellSingularNLLSProblem> nlls(new PowellSingularNLLSProblem);
695  nlls->parameters(x0);
696  NLLSSolverGSL solver1(nlls, 100, 1e-5, 1e-5, 1e-5);
697  NLLSSolverGSL solver2(nlls, 100, 1e-21, 1e-21, 1e-21);
698  BOOST_CHECK_EQUAL((int)solver1.status(), (int)NLLSSolverGSL::UNTRIED);
699  BOOST_CHECK_EQUAL((int)solver2.status(), (int)NLLSSolverGSL::UNTRIED);
700 
701 
702 
703 // std::cout
704 // << std::endl
705 // << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< nlls / solver1 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
706 // << std::endl;
707 
708  solver1.solve();
709 
710  int n_f_calls1 = nlls->num_residual_evaluations();
711  int n_j_calls1 = nlls->num_jacobian_evaluations();
712  double cst1 = nlls->cost();
713 
714 // std::cout
715 // << "Testing NLLSSolverGSL with Powell Singular function:" << std::endl
716 // << " Number of residual function evaluations = " << n_f_calls1 << std::endl
717 // << " Number of jacobian function evaluations = " << n_j_calls1 << std::endl
718 // << " Final solver1 status = " << solver1.status_str() << std::endl
719 // << " Final problem status (point) = " << nlls->parameters() << std::endl
720 // << " Final problem status (cost value) = " << cst1 << std::endl
721 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
722 // for( int i=0; i<=solver1.num_accepted_steps(); i++ )
723 // std::cout
724 // << " ========================================" << std::endl
725 // << " At point["<<i<<"] " << solver1.accepted_points()[i] << std::endl
726 // << " Cost["<<i<<"] = " << solver1.cost_at_accepted_points()[i] << std::endl
727 // << " Grad["<<i<<"] = " << solver1.gradient_at_accepted_points()[i] << std::endl;
728 
729  BOOST_CHECK_EQUAL((int) solver1.accepted_points().size(), solver1.num_accepted_steps()+1);
730  BOOST_CHECK_EQUAL((int) solver1.cost_at_accepted_points().size(), solver1.num_accepted_steps()+1);
731  BOOST_CHECK_EQUAL((int) solver1.gradient_at_accepted_points().size(), solver1.num_accepted_steps()+1);
732 
733  int iLast1 = solver1.num_accepted_steps();
734  BOOST_CHECK_EQUAL(cst1, solver1.cost_at_accepted_points()[iLast1]);
735  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver1.accepted_points()[iLast1])), 0.0, 1e-12);
736  BOOST_CHECK_CLOSE(sum(abs(nlls->gradient()-solver1.gradient_at_accepted_points()[iLast1])), 0.0, 1e-12);
737 
738  BOOST_CHECK_EQUAL((int)solver1.status(), (int)NLLSSolverGSL::SUCCESS);
739  BOOST_CHECK(n_f_calls1 < 15);
740  BOOST_CHECK(n_j_calls1 <= n_f_calls1);
741 
742 
743 
744 // std::cout
745 // << std::endl
746 // << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< nlls / solver2 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
747 // << std::endl;
748 
749  solver2.solve();
750 
751  int n_f_calls2 = nlls->num_residual_evaluations();
752  int n_j_calls2 = nlls->num_jacobian_evaluations();
753  double cst2 = nlls->cost();
754 
755 // std::cout
756 // << "Testing NLLSSolverGSL with Powell Singular function:" << std::endl
757 // << " Number of residual function evaluations = " << n_f_calls2 << std::endl
758 // << " Number of jacobian function evaluations = " << n_j_calls2 << std::endl
759 // << " Final solver2 status = " << solver2.status_str() << std::endl
760 // << " Final problem status (point) = " << nlls->parameters() << std::endl
761 // << " Final problem status (cost value) = " << cst2 << std::endl
762 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
763 // for( int i=0; i<=solver2.num_accepted_steps(); i++ )
764 // std::cout
765 // << " ========================================" << std::endl
766 // << " At point["<<i<<"] " << solver2.accepted_points()[i] << std::endl
767 // << " Cost["<<i<<"] = " << solver2.cost_at_accepted_points()[i] << std::endl
768 // << " Grad["<<i<<"] = " << solver2.gradient_at_accepted_points()[i] << std::endl;
769 
770  BOOST_CHECK_EQUAL((int) solver2.accepted_points().size(), solver2.num_accepted_steps()+1);
771  BOOST_CHECK_EQUAL((int) solver2.cost_at_accepted_points().size(), solver2.num_accepted_steps()+1);
772  BOOST_CHECK_EQUAL((int) solver2.gradient_at_accepted_points().size(), solver2.num_accepted_steps()+1);
773 
774  int iLast2 = solver2.num_accepted_steps();
775  BOOST_CHECK_EQUAL(cst2, solver2.cost_at_accepted_points()[iLast2]);
776  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver2.accepted_points()[iLast2])), 0.0, 1e-12);
777  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver2.gradient_at_accepted_points()[iLast2]))) < 1e-12);
778 
779  BOOST_CHECK_EQUAL((int)solver2.status(), (int)NLLSSolverGSL::SUCCESS);
780  BOOST_CHECK(n_f_calls2 < 30);
781  BOOST_CHECK(n_j_calls2 <= n_f_calls2);
782  BOOST_CHECK( (abs(cst2-0.0) < 0.0000001) );
783  BOOST_CHECK(abs(nlls->parameters()(0)-0.0) < 0.0000001);
784  BOOST_CHECK(abs(nlls->parameters()(1)-0.0) < 0.0000001);
785  BOOST_CHECK(abs(nlls->parameters()(2)-0.0) < 0.0000001);
786  BOOST_CHECK(abs(nlls->parameters()(3)-0.0) < 0.0000001);
787 
788 
789  // Last point recorded by the first solver must be the
790  // same as the first point recorded by the second solver.
791  //
792  BOOST_CHECK_EQUAL(solver1.cost_at_accepted_points()[iLast1], solver2.cost_at_accepted_points()[0]);
793  BOOST_CHECK_CLOSE(sum(abs(solver1.accepted_points()[iLast1]-solver2.accepted_points()[0])), 0.0, 1e-12);
794  BOOST_CHECK_CLOSE(sum(abs(solver1.gradient_at_accepted_points()[iLast1]-solver2.gradient_at_accepted_points()[0])), 0.0, 1e-12);
795 
796 }
797 
798 
799 
800 
801 
802 
803 
804 
805 // The case tests interfacing of
806 // maximum-likelihood/nlls-problem/nlls-solver
807 BOOST_AUTO_TEST_CASE(bard_ml)
808 {
809  Array<double, 1> x0(3); x0 = 1.0, 1.0, 1.0;
810  Array<double, 1> measurement(15);
811  Array<double, 1> measurement_error_cov(15);
812  measurement = 0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.10, 4.39;
813  measurement_error_cov = 1.0;
814 
815  boost::shared_ptr<BardMLProblem> ml(new BardMLProblem(measurement, measurement_error_cov));
816  boost::shared_ptr<NLLSMaxLikelihood> nlls(new NLLSMaxLikelihood(ml));
817  nlls->parameters(x0);
818  NLLSSolverGSL solver(nlls, 100, dx_tol_abs, dx_tol_rel, g_tol, verbose);
819  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
820  solver.solve();
821 
822  int n_f_calls = nlls->num_residual_evaluations();
823  int n_j_calls = nlls->num_jacobian_evaluations();
824  double cst = nlls->cost();
825 
826 // std::cout
827 // << "Testing NLLSSolverGSL with Bard/ML function:" << std::endl
828 // << " Number of residual function evaluations = " << n_f_calls << std::endl
829 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
830 // << " Final solver status = " << solver.status_str() << std::endl
831 // << " Final problem status (point) = " << nlls->parameters() << std::endl
832 // << " Final problem status (cost value) = " << cst << std::endl
833 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
834 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
835 // std::cout
836 // << " ========================================" << std::endl
837 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
838 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
839 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
840 
841  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
842  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
843  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
844 
845  int iLast = solver.num_accepted_steps();
846  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
847  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
848  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-12);
849 
850  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
851  BOOST_CHECK(n_f_calls < 10);
852  BOOST_CHECK(n_j_calls <= n_f_calls);
853  BOOST_CHECK( (abs(cst-0.0041074) < 0.000001) || (abs(cst-8.71431) < 0.001) );
854  if( abs(cst-0.0041074) < 0.000001 ) {
855  BOOST_CHECK_CLOSE(nlls->parameters()(0), 0.082411, 0.01);
856  BOOST_CHECK_CLOSE(nlls->parameters()(1), 1.1330, 0.01);
857  BOOST_CHECK_CLOSE(nlls->parameters()(2), 2.3437, 0.01);
858  }
859 }
860 
861 
862 // The case tests interfacing of
863 // maximum-likelihood/nlls-problem/nlls-solver
864 BOOST_AUTO_TEST_CASE(meyer_ml)
865 {
866  Array<double, 1> x0(3); x0 = 0.02, 4000.0, 250.0;
867  Array<double, 1> measurement(16);
868  Array<double, 1> measurement_error_cov(16);
869  measurement = 34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0, 8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0;
870  measurement_error_cov = 1.0;
871 
872  boost::shared_ptr<MeyerMLProblem> ml(new MeyerMLProblem(measurement, measurement_error_cov));
873  boost::shared_ptr<NLLSMaxLikelihood> nlls(new NLLSMaxLikelihood(ml));
874  nlls->parameters(x0);
875  NLLSSolverGSL solver(nlls, 200, dx_tol_abs, dx_tol_rel, g_tol, verbose);
876  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::UNTRIED);
877  solver.solve();
878 
879  int n_f_calls = nlls->num_residual_evaluations();
880  int n_j_calls = nlls->num_jacobian_evaluations();
881  double cst = nlls->cost();
882 
883 // std::cout
884 // << "Testing NLLSSolverGSL with Meyer/ML function:" << std::endl
885 // << " Number of residual function evaluations = " << n_f_calls << std::endl
886 // << " Number of jacobian function evaluations = " << n_j_calls << std::endl
887 // << " Final solver status = " << solver.status_str() << std::endl
888 // << " Final problem status (point) = " << nlls->parameters() << std::endl
889 // << " Final problem status (cost value) = " << cst << std::endl
890 // << " Final problem status (gradient) = " << nlls->gradient() << std::endl;
891 // for( int i=0; i<=solver.num_accepted_steps(); i++ )
892 // std::cout
893 // << " ========================================" << std::endl
894 // << " At point["<<i<<"] " << solver.accepted_points()[i] << std::endl
895 // << " Cost["<<i<<"] = " << solver.cost_at_accepted_points()[i] << std::endl
896 // << " Grad["<<i<<"] = " << solver.gradient_at_accepted_points()[i] << std::endl;
897 
898  BOOST_CHECK_EQUAL((int) solver.accepted_points().size(), solver.num_accepted_steps()+1);
899  BOOST_CHECK_EQUAL((int) solver.cost_at_accepted_points().size(), solver.num_accepted_steps()+1);
900  BOOST_CHECK_EQUAL((int) solver.gradient_at_accepted_points().size(), solver.num_accepted_steps()+1);
901 
902  int iLast = solver.num_accepted_steps();
903  BOOST_CHECK_EQUAL(cst, solver.cost_at_accepted_points()[iLast]);
904  BOOST_CHECK_CLOSE(sum(abs(nlls->parameters()-solver.accepted_points()[iLast])), 0.0, 1e-12);
905  BOOST_CHECK(fabs(sum(abs(nlls->gradient()-solver.gradient_at_accepted_points()[iLast]))) < 1e-8);
906 
907  BOOST_CHECK_EQUAL((int)solver.status(), (int)NLLSSolverGSL::SUCCESS);
908  BOOST_CHECK(n_f_calls < 140);
909  BOOST_CHECK(n_j_calls <= n_f_calls);
910  BOOST_CHECK_CLOSE(cst, 43.9729, 0.01);
911  BOOST_CHECK_CLOSE(nlls->parameters()(0), 0.0056096, 0.01);
912  BOOST_CHECK_CLOSE(nlls->parameters()(1), 6181.35, 0.01);
913  BOOST_CHECK_CLOSE(nlls->parameters()(2), 345.224, 0.01);
914 }
915 
916 */
917 
918 BOOST_AUTO_TEST_SUITE_END()
FullPhysics::IterativeSolver::status_str
virtual const char *const status_str() const
Returns the string version of the solver status.
Definition: iterative_solver.cc:14
FullPhysics::GlobalFixture
This is a global fixture that is available to all unit tests.
Definition: global_fixture.h:10
FullPhysics::IterativeSolverDer::gradient_at_accepted_points
virtual std::vector< blitz::Array< double, 1 > > gradient_at_accepted_points() const
Returns a vector (std) of gradients evaluated at accepted points.
Definition: iterative_solver_der.h:62
FullPhysics::IterativeSolver::SUCCESS
solve method called and a solution found
Definition: iterative_solver.h:51
FullPhysics::NLLSSolverGSLSM::solve
virtual void solve()
The method that solves the optimization problem.
Definition: nlls_solver_gsl_sm.cc:28
FullPhysics::IterativeSolver::num_accepted_steps
virtual int num_accepted_steps() const
Returns the number of the accepted steps.
Definition: iterative_solver.h:96
blitz
Apply value function to a blitz array.
Definition: auto_derivative.h:40
FullPhysics::IterativeSolver::UNTRIED
solve method not called yet
Definition: iterative_solver.h:55
BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(freudenstein_roth_b_gsl_sm__default)
Definition: nlls_solver_gsl_sm_helical_valley.cc:30
FullPhysics::IterativeSolver::status
virtual status_t status() const
Returns a value of IterativeSolver::status_t type.
Definition: iterative_solver.h:190
FullPhysics::IterativeSolver::accepted_points
virtual std::vector< blitz::Array< double, 1 > > accepted_points() const
Returns a vector (std) of accepted points.
Definition: iterative_solver.h:133
FullPhysics::IterativeSolver::cost_at_accepted_points
virtual std::vector< double > cost_at_accepted_points() const
Returns a vector (std) of cost function values at accepted points.
Definition: iterative_solver.h:157
FullPhysics
Contains classes to abstract away details in various Spurr Radiative Transfer software.
Definition: doxygen_python.h:1
Copyright © 2017, California Institute of Technology.
ALL RIGHTS RESERVED.
U.S. Government Sponsorship acknowledged.
Generated Fri Aug 24 2018 15:44:10