finite_difference_tests.cpp
Go to the documentation of this file.
1 // Copyright 2019, University of Maryland and the MANGO development team.
2 //
3 // This file is part of MANGO.
4 //
5 // MANGO is free software: you can redistribute it and/or modify it
6 // under the terms of the GNU Lesser General Public License as
7 // published by the Free Software Foundation, either version 3 of the
8 // License, or (at your option) any later version.
9 //
10 // MANGO is distributed in the hope that it will be useful, but
11 // WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with MANGO. If not, see
17 // <https://www.gnu.org/licenses/>.
18 
19 #include "catch.hpp"
20 #include "mango.hpp"
21 #include "Solver.hpp"
22 #include "Least_squares_solver.hpp"
23 
24 #include <cassert>
25 #include <cmath>
26 #include <iostream>
27 #include <iomanip>
28 
29 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
30 // Test finite-difference gradient, for a non-least-squares problem.
31 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
32 
33 void objective_function_1(int* N_parameters, const double* x, double* f, int* failed_int, mango::Problem* problem, void* user_data) {
34  assert(*N_parameters == 3);
35  *f = exp(x[0] * x[0] - exp(x[1]) + sin(x[2])); // A 3-D model function, for testing.
36  *failed_int = false;
37 }
38 
39 
40 TEST_CASE_METHOD(mango::Solver, "Solver::finite_difference_gradient()","[Solver][finite difference]") {
41  // The Catch2 macros automatically call the mango::Solver() constructor (the version with no arguments).
42  N_parameters = 3;
43  best_state_vector = new double[N_parameters];
44  state_vector = new double[N_parameters];
45  double* gradient = new double[N_parameters];
46  objective_function = &objective_function_1;
47  function_evaluations = 0;
48  verbose = 0;
49  double base_case_objective_function;
50 
51  // Set up MPI:
52  mpi_partition = new mango::MPI_Partition();
53  auto N_worker_groups_requested = GENERATE(range(1,5)); // Scan over N_worker_groups
54  mpi_partition->set_N_worker_groups(N_worker_groups_requested);
55  mpi_partition->init(MPI_COMM_WORLD);
56 
57  // Set the point about which we will compute the gradient:
58  state_vector[0] = 1.2;
59  state_vector[1] = 0.9;
60  state_vector[2] = -0.4;
61 
62  finite_difference_step_size = 1.0e-7;
63  double correct_objective_function = 2.443823056453063e-01;
64 
65  SECTION("1-sided differences") {
66  centered_differences = false;
67 
68  if (mpi_partition->get_proc0_world()) {
69  // Case of proc0_world
70  finite_difference_gradient(state_vector, &base_case_objective_function, gradient);
71  // Tell group leaders to exit.
72  int data = -1;
73  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
74  } else {
75  // Case for group leaders:
76  if (mpi_partition->get_proc0_worker_groups()) {
77  group_leaders_loop();
78  } else {
79  // Everybody else, i.e. workers. Nothing to do here.
80  }
81  }
82 
83  if (mpi_partition->get_proc0_world()) {
84  //std::cout << "base_case_objective_function: " << std::scientific << std::setw(24) << std::setprecision(15) << base_case_objective_function << ", 1-sided gradient: " << gradient[0] << ", " << gradient[1] << ", " << gradient[2] << std::endl;
85  // Finally, see if the results are correct:
86  CHECK( function_evaluations == 4);
87  CHECK(base_case_objective_function == Approx(correct_objective_function).epsilon(1e-14));
88  CHECK( gradient[0] == Approx( 5.865176283537110e-01).epsilon(1e-13));
89  CHECK( gradient[1] == Approx(-6.010834349701177e-01).epsilon(1e-13));
90  CHECK( gradient[2] == Approx( 2.250910244305793e-01).epsilon(1e-13));
91  // Also check best_state_vector, best_objective_function and best_function_evaluation?
92  }
93  }
94  SECTION("centered differences") {
95  centered_differences = true;
96 
97  if (mpi_partition->get_proc0_world()) {
98  // Case of proc0_world
99  finite_difference_gradient(state_vector, &base_case_objective_function, gradient);
100  // Tell group leaders to exit.
101  int data = -1;
102  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
103  } else {
104  // Case for group leaders:
105  if (mpi_partition->get_proc0_worker_groups()) {
106  group_leaders_loop();
107  } else {
108  // Everybody else, i.e. workers. Nothing to do here.
109  }
110  }
111 
112  if (mpi_partition->get_proc0_world()) {
113  //std::cout << "base_case_objective_function: " << std::scientific << std::setw(24) << std::setprecision(15) << base_case_objective_function << ", centered gradient: " << gradient[0] << ", " << gradient[1] << ", " << gradient[2] << std::endl;
114  // Finally, see if the results are correct:
115  CHECK( function_evaluations == 7);
116  CHECK(base_case_objective_function == Approx(correct_objective_function).epsilon(1e-14));
117  CHECK( gradient[0] == Approx( 5.865175337071982e-01).epsilon(1e-13));
118  CHECK( gradient[1] == Approx(-6.010834789627051e-01).epsilon(1e-13));
119  CHECK( gradient[2] == Approx( 2.250910093037906e-01).epsilon(1e-13));
120  }
121  }
122 
123  delete[] state_vector;
124  delete[] gradient;
125 }
126 
127 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
128 // Now consider a least-squares problem, and test both finite_difference_Jacobian() and
129 // finite_difference_gradient().
130 ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
131 
132 void residual_function_1(int* N_parameters, const double* x, int* N_terms, double* f, int* failed_int, mango::Problem* problem, void* user_data) {
133  assert(*N_parameters == 2);
134  assert(*N_terms == 4);
135  for (int j = 0; j < *N_terms; j++) {
136  f[j] = exp(j + x[0] * x[0] - exp(x[1]));
137  }
138  *failed_int = false;
139 }
140 
141 
142 TEST_CASE_METHOD(mango::Least_squares_solver, "Least_squares_solver::finite_difference_Jacobian() and problem::finite_difference_gradient()","[problem][finite difference]") {
143  // The Catch2 macros automatically call the mango::problem() constructor (the version with no arguments).
144  N_parameters = 2;
145  N_terms = 4;
146  best_state_vector = new double[N_parameters];
147  residuals = new double[N_terms]; // We must allocate this variable since the destructor will delete it.
148  double* base_case_residuals = new double[N_terms];
149  state_vector = new double[N_parameters];
150  double* Jacobian = new double[N_parameters * N_terms];
151  targets = new double[N_terms];
152  sigmas = new double[N_terms];
153  best_residual_function = new double[N_terms];
154  double* gradient = new double[N_parameters];
155  double base_case_objective_function;
156  mpi_partition = new mango::MPI_Partition();
157  residual_function = &residual_function_1;
158  function_evaluations = 0;
159  verbose = 0;
160 
161  // Set up MPI:
162  auto N_worker_groups_requested = GENERATE(range(1,5)); // Scan over N_worker_groups
163  //int N_worker_groups_requested = 1;
164  mpi_partition->set_N_worker_groups(N_worker_groups_requested);
165  mpi_partition->init(MPI_COMM_WORLD);
166 
167  // Set the point about which we will compute the derivatives:
168  state_vector[0] = 1.2;
169  state_vector[1] = 0.9;
170 
171  // Initialize targets and sigmas
172  for (int j=0; j<N_terms; j++) {
173  targets[j] = 1.5 + 2 * j;
174  sigmas[j] = 0.8 + 1.3 * j;
175  }
176 
177  finite_difference_step_size = 1.0e-7;
178 
179  // Here are all the reference values to compare against:
180  double correct_residuals[] = {3.607380846860443e-01, 9.805877804351946e-01, 2.665513944765978e+00, 7.245618119561541e+00};
181  double correct_d_residuals_d_x0_1sided[] = { 8.657715439008840e-01, 2.353411054922816e+00, 6.397234502131255e+00, 1.738948636642590e+01};
182  double correct_d_residuals_d_x0_centered[] = { 8.657714037352271e-01, 2.353410674116319e+00, 6.397233469623842e+00, 1.738948351093228e+01};
183  double correct_d_residuals_d_x1_1sided[] = {-8.872724499564555e-01, -2.411856577788640e+00, -6.556105911492693e+00, -1.782134355643450e+01};
184  double correct_d_residuals_d_x1_centered[] = {-8.872725151820582e-01, -2.411856754314101e+00, -6.556106388888594e+00, -1.782134486205678e+01};
185 
186  // Based on the above data, compute what the objective function and gradient should be:
187  double temp;
188  double correct_objective_function = 0;
189  double correct_gradient_1sided[] = {0.0, 0.0};
190  double correct_gradient_centered[] = {0.0, 0.0};
191  for (int j=0; j<N_terms; j++) {
192  temp = (correct_residuals[j] - targets[j]) / sigmas[j];
193  correct_objective_function += temp * temp;
194  correct_gradient_1sided[0] += 2 * (correct_residuals[j] - targets[j]) / (sigmas[j] * sigmas[j]) * correct_d_residuals_d_x0_1sided[j];
195  correct_gradient_1sided[1] += 2 * (correct_residuals[j] - targets[j]) / (sigmas[j] * sigmas[j]) * correct_d_residuals_d_x1_1sided[j];
196  correct_gradient_centered[0] += 2 * (correct_residuals[j] - targets[j]) / (sigmas[j] * sigmas[j]) * correct_d_residuals_d_x0_centered[j];
197  correct_gradient_centered[1] += 2 * (correct_residuals[j] - targets[j]) / (sigmas[j] * sigmas[j]) * correct_d_residuals_d_x1_centered[j];
198  }
199 
200  SECTION("1-sided differences, Jacobian") { // This section tests finite_difference_Jacobian()
201  centered_differences = false;
202 
203  if (mpi_partition->get_proc0_world()) {
204  // Case of proc0_world
205  finite_difference_Jacobian(state_vector, base_case_residuals, Jacobian);
206  // Tell group leaders to exit.
207  int data = -1;
208  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
209  } else {
210  // Case for group leaders:
211  if (mpi_partition->get_proc0_worker_groups()) {
212  group_leaders_loop();
213  } else {
214  // Everybody else, i.e. workers. Nothing to do here.
215  }
216  }
217 
218  if (mpi_partition->get_proc0_world()) {
219  // Finally, see if the results are correct:
220 
221  //std::cout << "residuals:" << std::scientific << std::setprecision(15);
222  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << base_case_residuals[k];
223  //std::cout << std::endl;
224  //std::cout << "d(residuals)/d(x[0]):";
225  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << Jacobian[k];
226  //std::cout << std::endl;
227  //std::cout << "d(residuals)/d(x[1]):";
228  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << Jacobian[k+N_terms];
229  //std::cout << std::endl;
230 
231  CHECK(function_evaluations == 3);
232 
233  //for (int k=0; k<N_terms; k++) CHECK(residuals[k] == Approx(correct_residuals[k]).epsilon(1e-14));
234  CHECK(base_case_residuals[0] == Approx(correct_residuals[0]).epsilon(1e-14));
235  CHECK(base_case_residuals[1] == Approx(correct_residuals[1]).epsilon(1e-14));
236  CHECK(base_case_residuals[2] == Approx(correct_residuals[2]).epsilon(1e-14));
237  CHECK(base_case_residuals[3] == Approx(correct_residuals[3]).epsilon(1e-14));
238 
239  // Order of Jacobian entries is [j_parameter*N_terms + j_term], so term is least significant, parameter is most significant.
240  CHECK(Jacobian[0] == Approx(correct_d_residuals_d_x0_1sided[0]).epsilon(1e-13));
241  CHECK(Jacobian[1] == Approx(correct_d_residuals_d_x0_1sided[1]).epsilon(1e-13));
242  CHECK(Jacobian[2] == Approx(correct_d_residuals_d_x0_1sided[2]).epsilon(1e-13));
243  CHECK(Jacobian[3] == Approx(correct_d_residuals_d_x0_1sided[3]).epsilon(1e-13));
244 
245  CHECK(Jacobian[4] == Approx(correct_d_residuals_d_x1_1sided[0]).epsilon(1e-13));
246  CHECK(Jacobian[5] == Approx(correct_d_residuals_d_x1_1sided[1]).epsilon(1e-13));
247  CHECK(Jacobian[6] == Approx(correct_d_residuals_d_x1_1sided[2]).epsilon(1e-13));
248  CHECK(Jacobian[7] == Approx(correct_d_residuals_d_x1_1sided[3]).epsilon(1e-13));
249 
250  }
251  }
252  SECTION("Centered differences, Jacobian") { // This section tests finite_difference_Jacobian()
253  centered_differences = true;
254 
255  if (mpi_partition->get_proc0_world()) {
256  // Case of proc0_world
257  finite_difference_Jacobian(state_vector, base_case_residuals, Jacobian);
258  // Tell group leaders to exit.
259  int data = -1;
260  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
261  } else {
262  // Case for group leaders:
263  if (mpi_partition->get_proc0_worker_groups()) {
264  group_leaders_loop();
265  } else {
266  // Everybody else, i.e. workers. Nothing to do here.
267  }
268  }
269 
270  if (mpi_partition->get_proc0_world()) {
271  // Finally, see if the results are correct:
272 
273  //std::cout << "residuals:" << std::scientific << std::setprecision(15);
274  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << base_case_residuals[k];
275  //std::cout << std::endl;
276  //std::cout << "d(residuals)/d(x[0]):";
277  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << Jacobian[k];
278  //std::cout << std::endl;
279  //std::cout << "d(residuals)/d(x[1]):";
280  //for (int k=0; k<N_terms; k++) std::cout << " " << std::setw(24) << Jacobian[k+N_terms];
281  //std::cout << std::endl;
282 
283  CHECK(function_evaluations == 5);
284 
285  //for (int k=0; k<N_terms; k++) CHECK(residuals[k] == Approx(correct_residuals[k]).epsilon(1e-14));
286  CHECK(base_case_residuals[0] == Approx(correct_residuals[0]).epsilon(1e-14));
287  CHECK(base_case_residuals[1] == Approx(correct_residuals[1]).epsilon(1e-14));
288  CHECK(base_case_residuals[2] == Approx(correct_residuals[2]).epsilon(1e-14));
289  CHECK(base_case_residuals[3] == Approx(correct_residuals[3]).epsilon(1e-14));
290 
291  // Order of Jacobian entries is [j_parameter*N_terms + j_term], so term is least significant, parameter is most significant.
292  CHECK(Jacobian[0] == Approx(correct_d_residuals_d_x0_centered[0]).epsilon(1e-13));
293  CHECK(Jacobian[1] == Approx(correct_d_residuals_d_x0_centered[1]).epsilon(1e-13));
294  CHECK(Jacobian[2] == Approx(correct_d_residuals_d_x0_centered[2]).epsilon(1e-13));
295  CHECK(Jacobian[3] == Approx(correct_d_residuals_d_x0_centered[3]).epsilon(1e-13));
296 
297  CHECK(Jacobian[4] == Approx(correct_d_residuals_d_x1_centered[0]).epsilon(1e-13));
298  CHECK(Jacobian[5] == Approx(correct_d_residuals_d_x1_centered[1]).epsilon(1e-13));
299  CHECK(Jacobian[6] == Approx(correct_d_residuals_d_x1_centered[2]).epsilon(1e-13));
300  CHECK(Jacobian[7] == Approx(correct_d_residuals_d_x1_centered[3]).epsilon(1e-13));
301 
302  }
303  }
304  SECTION("1-sided differences, gradient") { // This section tests finite_difference_gradient()
305  centered_differences = false;
306 
307  if (mpi_partition->get_proc0_world()) {
308  // Case of proc0_world
309  finite_difference_gradient(state_vector, &base_case_objective_function, gradient);
310  // Tell group leaders to exit.
311  int data = -1;
312  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
313  } else {
314  // Case for group leaders:
315  if (mpi_partition->get_proc0_worker_groups()) {
316  group_leaders_loop();
317  } else {
318  // Everybody else, i.e. workers. Nothing to do here.
319  }
320  }
321 
322  if (mpi_partition->get_proc0_world()) {
323  // Finally, see if the results are correct:
324  CHECK(function_evaluations == 3);
325 
326  CHECK(base_case_objective_function == Approx(correct_objective_function).epsilon(1e-14));
327 
328  CHECK(gradient[0] == Approx(correct_gradient_1sided[0]).epsilon(1e-13));
329  CHECK(gradient[1] == Approx(correct_gradient_1sided[1]).epsilon(1e-13));
330  }
331  }
332  SECTION("Centered differences, gradient") { // This section tests finite_difference_gradient()
333  centered_differences = true;
334 
335  if (mpi_partition->get_proc0_world()) {
336  // Case of proc0_world
337  finite_difference_gradient(state_vector, &base_case_objective_function, gradient);
338  // Tell group leaders to exit.
339  int data = -1;
340  MPI_Bcast(&data,1,MPI_INT,0,mpi_partition->get_comm_group_leaders());
341  } else {
342  // Case for group leaders:
343  if (mpi_partition->get_proc0_worker_groups()) {
344  group_leaders_loop();
345  } else {
346  // Everybody else, i.e. workers. Nothing to do here.
347  }
348  }
349 
350  if (mpi_partition->get_proc0_world()) {
351  // Finally, see if the results are correct:
352  CHECK(function_evaluations == 5);
353 
354  CHECK(base_case_objective_function == Approx(correct_objective_function).epsilon(1e-14));
355 
356  CHECK(gradient[0] == Approx(correct_gradient_centered[0]).epsilon(1e-13));
357  CHECK(gradient[1] == Approx(correct_gradient_centered[1]).epsilon(1e-13));
358  }
359  }
360 
361  delete[] state_vector;
362  delete[] Jacobian;
363  delete[] targets;
364  delete[] sigmas;
365  delete[] best_residual_function;
366  delete[] base_case_residuals;
367  delete[] gradient;
368  // best_state_vector and residuals will be deleted by destructor.
369 }
370 
371 
Least_squares_solver.hpp
residual_function_1
void residual_function_1(int *N_parameters, const double *x, int *N_terms, double *f, int *failed_int, mango::Problem *problem, void *user_data)
Definition: finite_difference_tests.cpp:132
objective_function_1
void objective_function_1(int *N_parameters, const double *x, double *f, int *failed_int, mango::Problem *problem, void *user_data)
Definition: finite_difference_tests.cpp:33
mango::Least_squares_solver
Definition: Least_squares_solver.hpp:28
Solver.hpp
mango.hpp
TEST_CASE_METHOD
TEST_CASE_METHOD(mango::Solver, "Solver::finite_difference_gradient()","[Solver][finite difference]")
Definition: finite_difference_tests.cpp:40
mango::Solver
Definition: Solver.hpp:31
mango::MPI_Partition
A class for dividing the set of MPI processes into worker groups.
Definition: mango.hpp:195
mango::Problem
Definition: mango.hpp:442