Fixed issue where python may deconstruct the diffeq while dense output still had reference and could call it

Author: jrenaud90

CyRK/cy/cysolver.cpp

@@ -317,6 +317,7 @@ void CySolverBase::take_step()
                     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
@@ -526,6 +527,7 @@ CySolverDense* CySolverBase::p_dense_output_heap()
         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

CyRK/cy/cysolver_api.pyx

@@ -88,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

@@ -13,6 +13,7 @@ CySolverDense::CySolverDense(
         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
@@ -22,6 +23,7 @@ CySolverDense::CySolverDense(
             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),
@@ -33,9 +35,17 @@ 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)
+    {
+        // TODO: Do we need to decref this at some point? During CySolver's deconstruction?
+        Py_XINCREF(this->cython_extension_class_instance);
+    }
+    
 }
 
-void CySolverDense::call(double t_interp, double* y_intepret)
+void CySolverDense::call(double t_interp, double* y_interp_ptr)
 {
     double step_factor = (t_interp - this->t_old) / this->step;
 
@@ -65,7 +75,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;
 
@@ -90,7 +100,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;
 
@@ -127,13 +137,13 @@ 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);
+        std::memcpy(y_interp_ptr, this->y_stored_ptr, sizeof(double) * this->num_y);
         break;
     }
 
@@ -151,31 +161,33 @@ void CySolverDense::call(double t_interp, double* y_intepret)
         // y array
         double y_tmp[Y_LIMIT];
         double* y_tmp_ptr = &y_tmp[0];
-        memcpy(y_tmp_ptr, this->cysolver_y_now_ptr, this->num_y);
+        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, num_dy);
+        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_intepret, this->num_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
+        printf("DEBUG!! About to call diffeq from dense...\n");
         this->cysolver_diffeq_ptr(this->cysolver_instance_ptr);
+        printf("DEBUG!! After diffeq call\n");
 
-        // Capture extra output and add to the y_intepret array
+        // 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] = this->cysolver_dy_now_ptr[i];
+            y_interp_ptr[i] = this->cysolver_dy_now_ptr[i];
         }
 
         // Reset CySolver state to what it was before
         cysolver_t_now_ptr[0] = t_tmp;
-        memcpy(this->cysolver_y_now_ptr, y_tmp_ptr, num_dy);
-        memcpy(this->cysolver_dy_now_ptr, dy_tmp_ptr, num_dy);
+        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);
     }
 }

CyRK/cy/dense.hpp

@@ -7,6 +7,7 @@
 #include <functional>
 #include <cstring>
 
+#include "Python.h"
 #include "common.hpp"
 
 // We need a pointer to the CySolverBase class. But that file includes this one. So we need to do a forward declaration
@@ -39,6 +40,7 @@ class CySolverDense
     double* cysolver_t_now_ptr  = nullptr;
     double* cysolver_y_now_ptr  = nullptr;
     double* cysolver_dy_now_ptr = nullptr;
+    PyObject* cython_extension_class_instance = nullptr;
 
     // Time step info
     double step = 0.0;
@@ -69,6 +71,7 @@ class CySolverDense
         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

CyRK/cy/pysolver_cyhook.pyx

@@ -1,6 +1,7 @@
 # distutils: language = c++
 # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
 
+from libc.stdio cimport printf
 
 cdef public api void call_diffeq_from_cython(object py_instance, DiffeqMethod diffeq):
     """Callback function used by the C++ model.
@@ -9,4 +10,5 @@ cdef public api void call_diffeq_from_cython(object py_instance, DiffeqMethod di
     """
 
     # Call the python diffeq.
+    printf("DEBUG!!! \t\t Calling python diffeq from api'd func.\n")
     diffeq(py_instance)

CyRK/cy/rk.cpp

@@ -1025,6 +1025,7 @@ CySolverDense* RKSolver::p_dense_output_heap()
         this->len_Pcols,
         this,
         this->diffeq,
+        this->cython_extension_class_instance,
         this->t_now_ptr,
         this->y_now_ptr,
         this->dy_now_ptr);