Extra ouput interpolation

CHANGES.md

@@ -2,6 +2,18 @@
 
 ## 2024
 
+#### v0.11.2 (2024-11-12)
+
+New:
+* Using extra output (via `num_extra`) with `cysolve_ivp` and `pysolve_ivp` now works when `dense_output` is set to true. CySolverSolution will now make additional calls to the differential equation to determine correct values for extra outputs which are provided alongside the interpolated y values.
+  * Added relevant tests.
+
+Other:
+* Refactored some misspellings in the cysolver c++ backend.
+
+Fix:
+* Fixed missing `np.import_array` in cysolver backend.
+
 #### v0.11.1 (2024-11-11)
 
 Fixes:

CyRK/array/interp.pyx

@@ -6,8 +6,6 @@ import numpy as np
 
 from libc.math cimport isnan, floor
 
-from libc.stdio cimport printf
-
 # Get machine precision.
 cdef double EPS
 EPS = np.finfo(dtype=np.float64).eps
@@ -504,17 +502,14 @@ cdef void interp_array_ptr(
         x_slope = 1.
 
     cdef double desired_x
-    printf("\n\n!!!DEBUG:: WORKING ON PRANGE\n")
     for index in range(desired_len):
         desired_x = desired_x_array[index]
 
         # Perform binary search with guess
         guess = <Py_ssize_t>floor(x_slope * <double>index)
-        printf("!!!DEBUG:: INDEX = %d; desired_x=%f; guess = %d; len_x = %d\n", index, desired_x, guess, len_x)
         j = c_binary_search_with_guess(desired_x, x_domain, len_x, guess)
 
         # Run interpolation
-        printf("!!!DEBUG:: j= %d\n", j)
         result = interp_ptr(desired_x, x_domain, dependent_values, len_x, provided_j=j)
 
         # Store result

CyRK/cy/cysolution.cpp

@@ -336,7 +336,7 @@ void CySolverResult::update_message(const char* const new_message_ptr)
     std::strcpy(this->message_ptr, new_message_ptr);
 }
 
-void CySolverResult::call(const double t, double* y_interp)
+void CySolverResult::call(const double t, double* y_interp_ptr)
 {
     if (!this->capture_dense_output) [[unlikely]]
     {
@@ -401,18 +401,18 @@ void CySolverResult::call(const double t, double* y_interp)
         }
 
         // Call interpolant to update y
-        this->dense_vector[closest_index]->call(t, y_interp);
+        this->dense_vector[closest_index]->call(t, y_interp_ptr);
     }
 }
 
-void CySolverResult::call_vectorize(const double* t_array_ptr, size_t len_t, double* y_interp)
+void CySolverResult::call_vectorize(const double* t_array_ptr, size_t len_t, double* y_interp_ptr)
 {
     double* y_sub_ptr;
 
     for (size_t i = 0; i < len_t; i++)
     {
         // Assume y is passed as a y0_0, y1_0, y2_0, ... y0_1, y1_1, y2_1, ...
-        y_sub_ptr = &y_interp[this->num_y * i];
+        y_sub_ptr = &y_interp_ptr[this->num_dy * i];
 
         this->call(t_array_ptr[i], y_sub_ptr);
     }

CyRK/cy/cysolution.hpp

@@ -110,6 +110,6 @@ class CySolverResult {
     void finalize();
     void reset();
     void update_message(const char* const new_message_ptr);
-    void call(const double t, double* y_interp);
-    void call_vectorize(const double* t_array_ptr, size_t len_t, double* y_interp);
+    void call(const double t, double* y_interp_ptr);
+    void call_vectorize(const double* t_array_ptr, size_t len_t, double* y_interp_ptr);
 };

CyRK/cy/cysolver.cpp

@@ -147,7 +147,7 @@ CySolverBase::CySolverBase(
 CySolverBase::~CySolverBase()
 {
     this->storage_ptr = nullptr;
-    if (this->use_pysolver)
+    if (this->deconstruct_python)
     {
         // Decrease reference count on the cython extension class instance
         Py_XDECREF(this->cython_extension_class_instance);
@@ -252,9 +252,6 @@ void CySolverBase::reset()
     this->reset_called = true;
 }
 
-#include <cstdio>
-#include <iostream>
-
 void CySolverBase::take_step()
 {    
     if (!this->reset_called) [[unlikely]]
@@ -316,8 +313,16 @@ void CySolverBase::take_step()
                     this->t_now_ptr[0],
                     this->y_old_ptr,
                     this->num_y,
-                    0 // Fake Q order just for consistent constructor call
+                    this->num_extra,
+                    0, // Fake Q order just for consistent constructor call
+                    this,
+                    this->diffeq,
+                    this->cython_extension_class_instance,
+                    this->t_now_ptr,
+                    this->y_now_ptr,
+                    this->dy_now_ptr
                     );
+
                 // Update the dense output class with integrator-specific data
                 this->p_dense_output_stack(dense_output);
 
@@ -327,7 +332,7 @@ void CySolverBase::take_step()
                 // Check if there are any t_eval steps between this new index and the last index.
                 // Get lowest and highest indices
                 auto lower_i = std::lower_bound(this->t_eval_vec.begin(), this->t_eval_vec.end(), this->t_now_ptr[0]) - this->t_eval_vec.begin();
-                auto upper_i  = std::upper_bound(this->t_eval_vec.begin(), this->t_eval_vec.end(), this->t_now_ptr[0]) - this->t_eval_vec.begin();
+                auto upper_i = std::upper_bound(this->t_eval_vec.begin(), this->t_eval_vec.end(), this->t_now_ptr[0]) - this->t_eval_vec.begin();
 
                 size_t t_eval_index_new;
                 if (lower_i == upper_i)
@@ -518,7 +523,14 @@ CySolverDense* CySolverBase::p_dense_output_heap()
         this->t_now_ptr[0],
         this->y_old_ptr,
         this->num_y,
-        0 // Fake Q order just for consistent constructor call
+        this->num_extra,
+        0, // Fake Q order just for consistent constructor call
+        this,
+        this->diffeq,
+        this->cython_extension_class_instance,
+        this->t_now_ptr,
+        this->y_now_ptr,
+        this->dy_now_ptr
         );
 }
 
@@ -551,8 +563,8 @@ void CySolverBase::set_cython_extension_instance(PyObject* cython_extension_clas
         }
         else
         {
-            // TODO: Do we need to decref this at some point? During CySolver's deconstruction?
             Py_XINCREF(this->cython_extension_class_instance);
+            this->deconstruct_python = true;
         }
     }
 }

CyRK/cy/cysolver.hpp

@@ -123,6 +123,7 @@ class CySolverBase {
     bool reset_called   = false;
     bool capture_extra  = false;
     bool user_provided_max_num_steps = false;
+    bool deconstruct_python = false;
 
     // Dense (Interpolation) Attributes
     bool use_dense_output = false;

CyRK/cy/cysolver_api.pyx

@@ -2,6 +2,7 @@
 # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
 
 import numpy as np
+np.import_array()
 
 # =====================================================================================================================
 # Import CySolverResult (container for integration results)
@@ -27,7 +28,7 @@ cdef class WrapCySolverResult:
         if not self.cyresult_ptr.capture_dense_output:
             raise AttributeError("Can not call WrapCySolverResult when dense_output set to False.")
 
-        y_interp_array = np.empty(self.cyresult_ptr.num_y, dtype=np.float64, order='C')
+        y_interp_array = np.empty(self.cyresult_ptr.num_dy, dtype=np.float64, order='C')
         cdef double[::1] y_interp_view = y_interp_array
         cdef double* y_interp_ptr      = &y_interp_view[0]
 
@@ -42,13 +43,13 @@ cdef class WrapCySolverResult:
 
         cdef size_t len_t = len(t_view)
 
-        y_interp_array = np.empty(self.cyresult_ptr.num_y * len_t, dtype=np.float64, order='C')
+        y_interp_array = np.empty(self.cyresult_ptr.num_dy * len_t, dtype=np.float64, order='C')
         cdef double[::1] y_interp_view = y_interp_array
         cdef double* y_interp_ptr = &y_interp_view[0]
         cdef double* t_array_ptr  = &t_view[0]
 
         self.cyresult_ptr.call_vectorize(t_array_ptr, len_t, y_interp_ptr)
-        return y_interp_array.reshape(len_t, self.cyresult_ptr.num_y).T
+        return y_interp_array.reshape(len_t, self.cyresult_ptr.num_dy).T
 
     @property
     def success(self):
@@ -87,7 +88,7 @@ cdef class WrapCySolverResult:
         if type(t) == np.ndarray:
             return self.call_vectorize(t)
         else:
-            return self.call(t).reshape(self.cyresult_ptr.num_y, 1)
+            return self.call(t).reshape(self.cyresult_ptr.num_dy, 1)
 
 
 # =====================================================================================================================

CyRK/cy/dense.cpp

@@ -1,15 +1,31 @@
 
 #include "dense.hpp"
 
+// Constructors
 CySolverDense::CySolverDense(
         int integrator_int,
         double t_old,
         double t_now,
         double* y_in_ptr,
         unsigned int num_y,
-        unsigned int Q_order) :
+        unsigned int num_extra,
+        unsigned int Q_order,
+        CySolverBase* cysolver_instance_ptr,
+        std::function<void (CySolverBase *)> cysolver_diffeq_ptr,
+        PyObject* cython_extension_class_instance,
+        double* cysolver_t_now_ptr,
+        double* cysolver_y_now_ptr,
+        double* cysolver_dy_now_ptr
+        ) :
             integrator_int(integrator_int),
             num_y(num_y),
+            num_extra(num_extra),
+            cysolver_instance_ptr(cysolver_instance_ptr),
+            cysolver_diffeq_ptr(cysolver_diffeq_ptr),
+            cython_extension_class_instance(cython_extension_class_instance),
+            cysolver_t_now_ptr(cysolver_t_now_ptr),
+            cysolver_y_now_ptr(cysolver_y_now_ptr),
+            cysolver_dy_now_ptr(cysolver_dy_now_ptr),
             t_old(t_old),
             t_now(t_now),
             Q_order(Q_order)
@@ -18,9 +34,27 @@ CySolverDense::CySolverDense(
     std::memcpy(this->y_stored_ptr, y_in_ptr, sizeof(double) * this->num_y);
     // Calculate step
     this->step = this->t_now - this->t_old;
+
+    // Make a strong reference to the python class (if this dense output was built using the python hooks).
+    if (cython_extension_class_instance)
+    {
+        Py_XINCREF(this->cython_extension_class_instance);
+        this->deconstruct_python = true;
+    }
+
 }
 
-void CySolverDense::call(double t_interp, double* y_intepret)
+// Destructors
+CySolverDense::~CySolverDense()
+{
+    if (this->deconstruct_python)
+    {
+        // Decrease reference count on the cython extension class instance
+        Py_XDECREF(this->cython_extension_class_instance);
+    }
+}
+
+void CySolverDense::call(double t_interp, double* y_interp_ptr)
 {
     double step_factor = (t_interp - this->t_old) / this->step;
 
@@ -50,7 +84,7 @@ void CySolverDense::call(double t_interp, double* y_intepret)
             // Finally multiply by step
             temp_double *= this->step;
 
-            y_intepret[y_i] = this->y_stored_ptr[y_i] + temp_double;
+            y_interp_ptr[y_i] = this->y_stored_ptr[y_i] + temp_double;
         }
         break;
 
@@ -75,7 +109,7 @@ void CySolverDense::call(double t_interp, double* y_intepret)
             // Finally multiply by step
             temp_double *= this->step;
 
-            y_intepret[y_i] = this->y_stored_ptr[y_i] + temp_double;
+            y_interp_ptr[y_i] = this->y_stored_ptr[y_i] + temp_double;
         }
         break;
 
@@ -112,17 +146,55 @@ void CySolverDense::call(double t_interp, double* y_intepret)
             temp_double += this->Q_ptr[Q_stride + 6];
             temp_double *= step_factor;
 
-            y_intepret[y_i] = this->y_stored_ptr[y_i] + temp_double;
+            y_interp_ptr[y_i] = this->y_stored_ptr[y_i] + temp_double;
         }
         break;
 
     [[unlikely]] default:
         // Don't know the model. Just return the input.
-        std::memcpy(y_intepret, this->y_stored_ptr, sizeof(double) * this->num_y);
-        for (size_t i = 0; i < this->num_y; i++)
+        std::memcpy(y_interp_ptr, this->y_stored_ptr, sizeof(double) * this->num_y);
+        break;
+    }
+
+    if (this->num_extra > 0)
+    {
+        // We have interpolated the dependent y-values but have not handled any extra outputs
+        // We can not use the RK (or any other integration method's) fancy interpolation because extra outputs are
+        // not included in the, for example, Q matrix building process.
+        // TODO: Perhaps we could include them in that? 
+        // For now, we will make an additional call to the diffeq using the y0 we just found above and t_interp.
+
+        size_t num_dy = this->num_y + this->num_extra;
+
+        // Store a copy of dy_now, t_now, and y_now into old vectors so we can make the call non destructively.
+        // y array
+        double y_tmp[Y_LIMIT];
+        double* y_tmp_ptr = &y_tmp[0];
+        memcpy(y_tmp_ptr, this->cysolver_y_now_ptr, sizeof(double) * this->num_y);
+        // dy array
+        double dy_tmp[DY_LIMIT];
+        double* dy_tmp_ptr = &dy_tmp[0];
+        memcpy(dy_tmp_ptr, this->cysolver_dy_now_ptr, sizeof(double) * num_dy);
+        // t
+        double t_tmp = cysolver_t_now_ptr[0];
+
+        // Load new values into t and y
+        memcpy(this->cysolver_y_now_ptr, y_interp_ptr, sizeof(double) * this->num_y);
+        cysolver_t_now_ptr[0] = t_interp;
+
+        // Call diffeq to update dy_now pointer
+        this->cysolver_diffeq_ptr(this->cysolver_instance_ptr);
+
+        // Capture extra output and add to the y_interp_ptr array
+        // We already have y interpolated from above so start at num_y
+        for (size_t i = this->num_y; i < num_dy; i++)
         {
-            y_intepret[i] = 0.75;
+            y_interp_ptr[i] = this->cysolver_dy_now_ptr[i];
         }
-        break;
+
+        // Reset CySolver state to what it was before
+        cysolver_t_now_ptr[0] = t_tmp;
+        memcpy(this->cysolver_y_now_ptr, y_tmp_ptr, sizeof(double) * num_y);
+        memcpy(this->cysolver_dy_now_ptr, dy_tmp_ptr, sizeof(double) * num_dy);
     }
 }