Introduce move semantics to adouble, clean-up tape types, clean-up include files

ADOL-C/boost-test/CMakeLists.txt

@@ -16,6 +16,7 @@ endif()
 set(SOURCE_FILES
     adouble.cpp
     main.cpp
+    pdouble.cpp
     traceCompositeTests.cpp
     tracelessCompositeTests.cpp
     tracelessOperatorScalar.cpp
@@ -31,6 +32,10 @@ set(SOURCE_FILES
 file(GLOB UNI5_FOR_FILES "uni5_for/*.cpp")
 list(APPEND SOURCE_FILES ${UNI5_FOR_FILES})
 
+# Add all source files from integration_tests
+file(GLOB INTEGRATION_TESTS_FILES "integration_tests/*.cpp")
+list(APPEND SOURCE_FILES ${INTEGRATION_TESTS_FILES})
+
 add_executable(boost-test-adolc ${SOURCE_FILES})
 target_include_directories(boost-test-adolc PRIVATE "${ADOLC_INCLUDE_DIR}")
 target_link_libraries(boost-test-adolc PRIVATE

ADOL-C/boost-test/integration_tests/checkpointing.cpp

@@ -0,0 +1,185 @@
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+namespace tt = boost::test_tools;
+
+#include "../const.h"
+#include <adolc/adolc.h>
+#include <array>
+#include <numeric>
+#include <vector>
+
+// Euler step (double version)
+int euler_step(size_t n, double *y) {
+  y[0] = y[0] + 0.01 * y[0];
+  y[1] = y[1] + 0.01 * 2 * y[1];
+  return 1;
+}
+
+// Euler step (adouble version)
+int euler_step_act(size_t n, adouble *y) {
+  y[0] = y[0] + 0.01 * y[0];
+  y[1] = y[1] + 0.01 * 2 * y[1];
+  return 1;
+}
+
+BOOST_AUTO_TEST_SUITE(test_checkpoint_example)
+BOOST_AUTO_TEST_CASE(Checkpointing_Gradient_Comparison) {
+  const int16_t tag_full = 1;  // Tag for full taping
+  const int16_t tag_part = 2;  // Tag for partial taping with checkpointing
+  const int16_t tag_check = 3; // Tag for checkpointing
+  const size_t n = 2;          // Number of state variables
+  const int steps = 100;       // Number of time steps
+  const double t0 = 0.0;       // Initial time
+  const double tf = 1.0;       // Final time
+
+  // State variables (double and adouble versions)
+  std::vector<double> y_double(n);
+  ensureContiguousLocations(2 * n);
+  std::vector<adouble> y_adouble_1(n);
+  std::vector<adouble> y_adouble_2(n);
+
+  // Control variables (double and adouble versions)
+  std::vector<double> conp{1.0, 1.0}; // Initial control values
+  std::vector<adouble> con(n);
+
+  // Target value and gradient
+  std::vector<double> out(2);
+  std::vector<double> grad_full(n); // Gradient from full taping
+  std::vector<double> grad_part(n); // Gradient from checkpointing
+
+  // Full taping of the time step loop
+  trace_on(tag_full);
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y_adouble_1[0] = con[0];
+  y_adouble_1[1] = con[1];
+
+  for (int i = 0; i < steps; i++) {
+    euler_step_act(n, y_adouble_1.data());
+  }
+  y_adouble_1[0] + y_adouble_1[1] >>= out[0];
+  trace_off();
+
+  // Compute gradient using full taping
+  gradient(tag_full, n, conp.data(), grad_full.data());
+
+  // Checkpointing setup
+  CP_Context cpc(euler_step_act);    // Checkpointing context
+  cpc.setDoubleFct(euler_step);      // Double version of the time step function
+  cpc.setNumberOfSteps(steps);       // Number of time steps
+  cpc.setNumberOfCheckpoints(5);     // Number of checkpoints
+  cpc.setDimensionXY(n);             // Dimension of input/output
+  cpc.setInput(y_adouble_2.data());  // Input vector
+  cpc.setOutput(y_adouble_2.data()); // Output vector
+  cpc.setTapeNumber(tag_check);      // Tape number for checkpointing
+  cpc.setAlwaysRetaping(false);      // Do not always retape
+
+  // Partial taping with checkpointing
+  trace_on(tag_part);
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y_adouble_2[0] = con[0];
+  y_adouble_2[1] = con[1];
+
+  cpc.checkpointing(); // Perform checkpointing
+
+  y_adouble_2[0] + y_adouble_2[1] >>= out[1];
+  trace_off();
+
+  // test if both taping results are equal
+  BOOST_TEST(out[0] == out[1], tt::tolerance(tol));
+  // Compute gradient using checkpointing
+  gradient(tag_part, n, conp.data(), grad_part.data());
+  // Compare gradients from full taping and checkpointing
+  for (size_t i = 0; i < n; i++) {
+    BOOST_TEST(grad_full[i] == grad_part[i], tt::tolerance(tol));
+  }
+}
+BOOST_AUTO_TEST_CASE(Checkpointing_fov_reverse) {
+  const int16_t tag_full = 1;  // Tag for full taping
+  const int16_t tag_part = 2;  // Tag for partial taping with checkpointing
+  const int16_t tag_check = 3; // Tag for checkpointing
+  const size_t n = 2;          // Number of state variables
+  const int steps = 100;       // Number of time steps
+  const double t0 = 0.0;       // Initial time
+  const double tf = 1.0;       // Final time
+
+  // State variables (double and adouble versions)
+  std::vector<double> y_double(n);
+  ensureContiguousLocations(2 * n);
+  std::vector<adouble> y_adouble_1(n);
+  std::vector<adouble> y_adouble_2(n);
+
+  // Control variables (double and adouble versions)
+  std::vector<double> conp{1.0, 1.0}; // Initial control values
+  std::vector<adouble> con(n);
+
+  // Target value and gradient
+  std::vector<double> out(2);
+  std::vector<double> grad_full(n); // Gradient from full taping
+  std::vector<double> grad_part(n); // Gradient from checkpointing
+
+  // Full taping of the time step loop
+  trace_on(tag_full, 1);
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y_adouble_1[0] = con[0];
+  y_adouble_1[1] = con[1];
+
+  for (int i = 0; i < steps; i++) {
+    euler_step_act(n, y_adouble_1.data());
+  }
+  y_adouble_1[0] + y_adouble_1[1] >>= out[0];
+  trace_off();
+
+  // weights
+  double **U = myalloc2(2, 1);
+  U[0][0] = 1.0;
+  U[1][0] = -1.0;
+
+  // outputs
+  double **Z_full = myalloc2(2, 2);
+  double **Z_part = myalloc2(2, 2);
+
+  // Compute vector-mode reverse
+  fov_reverse(tag_full, 1, 2, 2, U, Z_full);
+
+  // Checkpointing setup
+  CP_Context cpc(euler_step_act);    // Checkpointing context
+  cpc.setDoubleFct(euler_step);      // Double version of the time step function
+  cpc.setNumberOfSteps(steps);       // Number of time steps
+  cpc.setNumberOfCheckpoints(5);     // Number of checkpoints
+  cpc.setDimensionXY(n);             // Dimension of input/output
+  cpc.setInput(y_adouble_2.data());  // Input vector
+  cpc.setOutput(y_adouble_2.data()); // Output vector
+  cpc.setTapeNumber(tag_check);      // Tape number for checkpointing
+  cpc.setAlwaysRetaping(true);       // Do always retape
+
+  // Partial taping with checkpointing
+  trace_on(tag_part, 1);
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y_adouble_2[0] = con[0];
+  y_adouble_2[1] = con[1];
+
+  cpc.checkpointing(); // Perform checkpointing
+
+  y_adouble_2[0] + y_adouble_2[1] >>= out[1];
+  trace_off();
+
+  // test if both taping results are equal
+  BOOST_TEST(out[0] == out[1], tt::tolerance(tol));
+
+  // Compute gradient using checkpointing
+  fov_reverse(tag_part, 1, 2, 2, U, Z_part);
+  // Compare gradients from full taping and checkpointing
+  for (size_t i = 0; i < 2; ++i) {
+    for (size_t j = 0; j < 2; ++j)
+      BOOST_TEST(Z_full[i][j] == Z_part[i][j], tt::tolerance(tol));
+  }
+  myfree2(U);
+  myfree2(Z_full);
+  myfree2(Z_part);
+}
+BOOST_AUTO_TEST_SUITE_END()

ADOL-C/boost-test/integration_tests/external_fct.cpp

@@ -0,0 +1,145 @@
+#include "../const.h"
+#include <adolc/adolc.h>
+#include <boost/test/unit_test.hpp>
+#include <memory>
+#include <vector>
+
+namespace tt = boost::test_tools;
+
+BOOST_AUTO_TEST_SUITE(ExternalFunctionTests)
+
+const double h = 0.01;
+const int steps = 100;
+
+const short tag_full = 1;
+const short tag_part = 2;
+const short tag_ext_fct = 3;
+
+// Control parameters
+std::vector<double> conp = {1.0, 1.0};
+std::vector<double> grad_full(2);
+std::vector<double> grad_ext(2);
+
+ext_diff_fct *edf;
+std::vector<double> yp = {0};
+std::vector<double> ynewp = {0};
+std::vector<double> u = {1.0, 1.0};
+std::vector<double> z = {0};
+
+void euler_step_act(size_t n, adouble *yin, size_t m, adouble *yout) {
+  yout[0] = yin[0] + h * yin[0];
+  yout[1] = yin[1] + h * 2 * yin[1];
+}
+
+int euler_step(size_t n, double *yin, size_t m, double *yout) {
+  yout[0] = yin[0] + h * yin[0];
+  yout[1] = yin[1] + h * 2 * yin[1];
+  return 1;
+}
+
+int zos_for_euler_step(size_t n, double *yin, size_t m, double *yout) {
+  int rc;
+  set_nested_ctx(tag_ext_fct, true);
+  rc = zos_forward(tag_ext_fct, 2, 2, 0, yin, yout);
+  set_nested_ctx(tag_ext_fct, false);
+  return rc;
+}
+
+int fos_rev_euler_step(size_t n, double *u, size_t m, double *z, double *,
+                       double *) {
+  int rc;
+  set_nested_ctx(tag_ext_fct, true);
+  zos_forward(tag_ext_fct, 2, 2, 1, edf->dp_x, edf->dp_y);
+  rc = fos_reverse(tag_ext_fct, 2, 2, u, z);
+  set_nested_ctx(tag_ext_fct, false);
+  return rc;
+}
+
+void setup_full_taping(const short tag_full, const std::vector<double> &conp) {
+  trace_on(tag_full);
+  std::vector<adouble> y(2), ynew(2);
+  std::vector<adouble> con(2);
+
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y[0] = con[0];
+  y[1] = con[1];
+
+  for (int i = 0; i < steps; i++) {
+    euler_step_act(2, y.data(), 2, ynew.data());
+    y[0] = ynew[0];
+    y[1] = ynew[1];
+  }
+
+  adouble f;
+  f = y[0] + y[1];
+  double f_out;
+  f >>= f_out;
+  trace_off();
+}
+
+void setup_external_function(const short tag_ext_fct,
+                             const std::vector<double> &conp) {
+  trace_on(tag_ext_fct);
+  std::vector<adouble> y(2), ynew(2);
+  y[0] <<= conp[0];
+  y[1] <<= conp[1];
+  euler_step_act(2, y.data(), 2, ynew.data());
+
+  double f_out;
+  ynew[0] >>= f_out; // Dummy output
+  ynew[1] >>= f_out; // Dummy output
+  trace_off();
+
+  edf = reg_ext_fct(euler_step);
+  edf->zos_forward = zos_for_euler_step;
+  edf->dp_x = yp.data();
+  edf->dp_y = ynewp.data();
+  edf->fos_reverse = fos_rev_euler_step;
+  edf->dp_U = u.data();
+  edf->dp_Z = z.data();
+}
+
+void setup_external_taping(size_t tag_part, std::vector<double> conp) {
+  trace_on(tag_part);
+  ensureContiguousLocations(4);
+  std::vector<adouble> y(2), ynew(2);
+  std::vector<adouble> con(2);
+
+  con[0] <<= conp[0];
+  con[1] <<= conp[1];
+  y[0] = con[0];
+  y[1] = con[1];
+
+  for (int i = 0; i < steps; i++) {
+    call_ext_fct(edf, 2, y.data(), 2, ynew.data());
+    y[0] = ynew[0];
+    y[1] = ynew[1];
+  }
+
+  adouble f;
+  f = y[0] + y[1];
+  double f_out;
+  f >>= f_out;
+  trace_off();
+}
+
+BOOST_AUTO_TEST_CASE(CompareFullAndExternalGradients) {
+  const double expected0 = exp(h * steps);
+  const double expected1 = exp(2 * h * steps);
+
+  setup_full_taping(tag_full, conp);
+  gradient(tag_full, 2, conp.data(), grad_full.data());
+
+  setup_external_function(tag_ext_fct, conp);
+  setup_external_taping(tag_part, conp);
+
+  gradient(tag_part, 2, conp.data(), grad_ext.data());
+
+  // Verify gradients match
+  for (int i = 0; i < 2; ++i) {
+    BOOST_TEST(grad_full[i] == grad_ext[i], tt::tolerance(tol));
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()