working on pre-eval

jrenaud90 · Jul 22, 2024 · e8d6ce9 · e8d6ce9
1 parent 965e928
commit e8d6ce9
Show file tree

Hide file tree

Showing 16 changed files with 335 additions and 84 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -7,10 +7,13 @@
 C++ Backend:
 * Changed optional args from double pointer to void pointers to allow for arbitrary objects to be passed in.
   * Added description of this feature to "Documentation/Advanced CySolver.md" documentation and "Demos/Advanced CySolver Examples.ipynb" jupyter notebook.
+* Allow users to specify a "Pre-Eval" function that can be passed to the differential equation. This function should take in time, y, and args and update an output pointer which can then be used by the diffeq to solve for dydt.
 
 `cysolve_ivp`:
+* Change call signature to accept new `pre_eval_func` function.
 * Added more differential equations to tests.
-* Added tests to check new void arg system.
+* Added tests to check new void arg feature.
+* Added tests to check new pre-eval function feature.
 
 ### v0.10.0 (2024-07-17)
 

diff --git a/CyRK/cy/common.hpp b/CyRK/cy/common.hpp
@@ -39,7 +39,11 @@ static const double DYNAMIC_GROWTH_RATE = 1.618;
 static constexpr double SIZE_MAX_DBL = 0.99 * SIZE_MAX;
 
 
-typedef void (*DiffeqFuncType)(double*, double, double*, const void*);
+
+typedef void (*PreEvalFunc)(void*, double, double*, const void*);
+
+typedef void (*DiffeqFuncType)(double*, double, double*, const void*, PreEvalFunc);
+
 
 
 struct MaxNumStepsOutput

diff --git a/CyRK/cy/cysolve.cpp b/CyRK/cy/cysolve.cpp
@@ -16,6 +16,7 @@ std::shared_ptr<CySolverResult> baseline_cysolve_ivp(
     const bool dense_output,
     const double* t_eval,
     const size_t len_t_eval,
+    PreEvalFunc pre_eval_func,
     // rk optional arguments
     const double rtol,
     const double atol,
@@ -85,7 +86,7 @@ std::shared_ptr<CySolverResult> baseline_cysolve_ivp(
         solver = new RK23(
             // Common Inputs
             diffeq_ptr, solution_ptr, t_start, t_end, y0_ptr, num_y, num_extra, args_ptr, max_num_steps, max_ram_MB,
-            dense_output, t_eval, len_t_eval,
+            dense_output, t_eval, len_t_eval, pre_eval_func,
             // RK Inputs
             rtol, atol, rtols_ptr, atols_ptr, max_step_size, first_step_size
         );
@@ -95,7 +96,7 @@ std::shared_ptr<CySolverResult> baseline_cysolve_ivp(
         solver = new RK45(
             // Common Inputs
             diffeq_ptr, solution_ptr, t_start, t_end, y0_ptr, num_y, num_extra, args_ptr, max_num_steps, max_ram_MB,
-            dense_output, t_eval, len_t_eval,
+            dense_output, t_eval, len_t_eval, pre_eval_func,
             // RK Inputs
             rtol, atol, rtols_ptr, atols_ptr, max_step_size, first_step_size
         );
@@ -105,7 +106,7 @@ std::shared_ptr<CySolverResult> baseline_cysolve_ivp(
         solver = new DOP853(
             // Common Inputs
             diffeq_ptr, solution_ptr, t_start, t_end, y0_ptr, num_y, num_extra, args_ptr, max_num_steps, max_ram_MB,
-            dense_output, t_eval, len_t_eval,
+            dense_output, t_eval, len_t_eval, pre_eval_func,
             // RK Inputs
             rtol, atol, rtols_ptr, atols_ptr, max_step_size, first_step_size
         );
@@ -181,6 +182,9 @@ PySolver::PySolver(
     // We need to pass a fake diffeq pointer (diffeq ptr is unused in python-based solver)
     DiffeqFuncType diffeq_ptr = nullptr;
 
+    // We also need to pass a fake pre-eval function
+    PreEvalFunc pre_eval_func = nullptr;
+
     // Build the solver class. This must be heap allocated to take advantage of polymorphism.
     switch (this->integration_method)
     {
@@ -201,6 +205,7 @@ PySolver::PySolver(
             dense_output,
             t_eval,
             len_t_eval,
+            pre_eval_func,
             // rk optional arguments
             rtol,
             atol,
@@ -226,6 +231,7 @@ PySolver::PySolver(
             dense_output,
             t_eval,
             len_t_eval,
+            pre_eval_func,
             // rk optional arguments
             rtol,
             atol,
@@ -251,6 +257,7 @@ PySolver::PySolver(
             dense_output,
             t_eval,
             len_t_eval,
+            pre_eval_func,
             // rk optional arguments
             rtol,
             atol,

diff --git a/CyRK/cy/cysolve.hpp b/CyRK/cy/cysolve.hpp
@@ -22,6 +22,7 @@ std::shared_ptr<CySolverResult> baseline_cysolve_ivp(
     const bool dense_output = false,
     const double* t_eval = nullptr,
     const size_t len_t_eval = 0,
+    PreEvalFunc pre_eval_func = nullptr,
     // rk optional arguments
     const double rtol = 1.0e-3,
     const double atol = 1.0e-6,

diff --git a/CyRK/cy/cysolver.cpp b/CyRK/cy/cysolver.cpp
@@ -38,7 +38,8 @@ CySolverBase::CySolverBase(
         const size_t max_ram_MB,
         const bool use_dense_output,
         const double* t_eval,
-        const size_t len_t_eval) :
+        const size_t len_t_eval,
+        PreEvalFunc pre_eval_func) :
             status(0),
             num_y(num_y),
             num_extra(num_extra),
@@ -48,7 +49,8 @@ CySolverBase::CySolverBase(
             diffeq_ptr(diffeq_ptr),
             args_ptr(args_ptr),
             use_dense_output(use_dense_output),
-            len_t_eval(len_t_eval)
+            len_t_eval(len_t_eval),
+            pre_eval_func(pre_eval_func)
 {
     // Parse inputs
     this->capture_extra = num_extra > 0;
@@ -194,7 +196,7 @@ bool CySolverBase::check_status() const
 void CySolverBase::cy_diffeq() noexcept
 {
     // Call c function
-    this->diffeq_ptr(this->dy_now_ptr, this->t_now_ptr[0], this->y_now_ptr, this->args_ptr);
+    this->diffeq_ptr(this->dy_now_ptr, this->t_now_ptr[0], this->y_now_ptr, this->args_ptr, this->pre_eval_func);
 }
 
 void CySolverBase::reset()

diff --git a/CyRK/cy/cysolver.hpp b/CyRK/cy/cysolver.hpp
@@ -52,7 +52,8 @@ class CySolverBase {
         const size_t max_ram_MB = 2000,
         const bool use_dense_output = false,
         const double* t_eval = nullptr,
-        const size_t len_t_eval = 0
+        const size_t len_t_eval = 0,
+        PreEvalFunc pre_eval_func = nullptr
     );
 
     void change_storage(std::shared_ptr<CySolverResult> new_storage_ptr, bool auto_reset = true);
@@ -114,6 +115,9 @@ class CySolverBase {
     // Information on capturing extra information during integration.
     int num_extra = 0;
 
+    // Function to send to diffeq which is called before dy is calculated
+    PreEvalFunc pre_eval_func = nullptr;
+
     // Keep bools together to reduce size
     bool direction_flag = false;
     bool reset_called   = false;

diff --git a/CyRK/cy/cysolverNew.pxd b/CyRK/cy/cysolverNew.pxd
@@ -20,7 +20,8 @@ cdef extern from "common.cpp" nogil:
     const unsigned int DY_LIMIT
     const double MAX_STEP
 
-    ctypedef void (*DiffeqFuncType)(double*, double, double*, const void*)
+    ctypedef void (*PreEvalFunc)(void*, double, double*, const void*)
+    ctypedef void (*DiffeqFuncType)(double*, double, double*, const void*, PreEvalFunc)
 
     cdef size_t find_expected_size(        
         int num_y,
@@ -121,7 +122,8 @@ cdef extern from "cysolver.cpp" nogil:
             const size_t max_ram_MB,
             const cpp_bool use_dense_output,
             const double* t_eval,
-            const size_t len_t_eval
+            const size_t len_t_eval,
+            PreEvalFunc pre_eval_func
         )
 
         shared_ptr[CySolverResult] storage_ptr
@@ -190,6 +192,7 @@ cdef extern from "rk.cpp" nogil:
             const cpp_bool use_dense_output,
             const double* t_eval,
             const size_t len_t_eval,
+            PreEvalFunc pre_eval_func,
             const double rtol,
             const double atol,
             const double* rtols_ptr,
@@ -220,6 +223,7 @@ cdef extern from "rk.cpp" nogil:
             const cpp_bool use_dense_output,
             const double* t_eval,
             const size_t len_t_eval,
+            PreEvalFunc pre_eval_func,
             const double rtol,
             const double atol,
             const double* rtols_ptr,
@@ -246,6 +250,7 @@ cdef extern from "rk.cpp" nogil:
             const cpp_bool use_dense_output,
             const double* t_eval,
             const size_t len_t_eval,
+            PreEvalFunc pre_eval_func,
             const double rtol,
             const double atol,
             const double* rtols_ptr,
@@ -272,6 +277,7 @@ cdef extern from "rk.cpp" nogil:
             const cpp_bool use_dense_output,
             const double* t_eval,
             const size_t len_t_eval,
+            PreEvalFunc pre_eval_func,
             const double rtol,
             const double atol,
             const double* rtols_ptr,
@@ -301,6 +307,7 @@ cdef extern from "cysolve.cpp" nogil:
             const cpp_bool dense_output,
             const double* t_eval,
             const size_t len_t_eval,
+            PreEvalFunc pre_eval_func,
             const double rtol,
             const double atol,
             const double* rtols_ptr,
@@ -361,6 +368,7 @@ cdef CySolveOutput cysolve_ivp(
     bint dense_output = *,
     double* t_eval = *,
     size_t len_t_eval = *,
+    PreEvalFunc pre_eval_func = *,
     double* rtols_ptr = *,
     double* atols_ptr = *,
     double max_step = *,
@@ -383,6 +391,7 @@ cdef CySolveOutput cysolve_ivp_gil(
     bint dense_output = *,
     double* t_eval = *,
     size_t len_t_eval = *,
+    PreEvalFunc pre_eval_func = *,
     double* rtols_ptr = *,
     double* atols_ptr = *,
     double max_step = *,

diff --git a/CyRK/cy/cysolverNew.pyx b/CyRK/cy/cysolverNew.pyx
@@ -103,6 +103,7 @@ cdef CySolveOutput cysolve_ivp(
             bint dense_output = False,
             double* t_eval = NULL,
             size_t len_t_eval = 0,
+            PreEvalFunc pre_eval_func = NULL,
             double* rtols_ptr = NULL,
             double* atols_ptr = NULL,
             double max_step = MAX_STEP,
@@ -124,6 +125,7 @@ cdef CySolveOutput cysolve_ivp(
         dense_output,
         t_eval,
         len_t_eval,
+        pre_eval_func,
         rtol,
         atol,
         rtols_ptr,
@@ -149,6 +151,7 @@ cdef CySolveOutput cysolve_ivp_gil(
             bint dense_output = False,
             double* t_eval = NULL,
             size_t len_t_eval = 0,
+            PreEvalFunc pre_eval_func = NULL,
             double* rtols_ptr = NULL,
             double* atols_ptr = NULL,
             double max_step = MAX_STEP,
@@ -170,6 +173,7 @@ cdef CySolveOutput cysolve_ivp_gil(
         dense_output,
         t_eval,
         len_t_eval,
+        pre_eval_func,
         rtol,
         atol,
         rtols_ptr,

diff --git a/CyRK/cy/cysolverNew_test.pxd b/CyRK/cy/cysolverNew_test.pxd
@@ -1,9 +1,13 @@
 from libcpp cimport bool as cpp_bool
 
-cdef void baseline_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void accuracy_test_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void extraoutput_test_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void lorenz_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void lorenz_extraoutput_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void lotkavolterra_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
-cdef void pendulum_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr) noexcept nogil
+from libc.math cimport sin, cos, fabs, fmin, fmax
+
+from CyRK.cy.cysolverNew cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType,MAX_STEP, CySolveOutput, PreEvalFunc
+
+cdef void baseline_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void accuracy_test_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void extraoutput_test_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void lorenz_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void lorenz_extraoutput_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void lotkavolterra_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil
+cdef void pendulum_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil