From 634cef741ecfcf7b516e07fc770d981e059c8bc7 Mon Sep 17 00:00:00 2001 From: Jrenaud-Desk Date: Tue, 12 Nov 2024 21:54:01 -0500 Subject: [PATCH 1/9] Importing numpy arrays properly --- CyRK/cy/cysolver_api.pyx | 1 + CyRK/cy/pysolver.pyx | 2 -- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CyRK/cy/cysolver_api.pyx b/CyRK/cy/cysolver_api.pyx index 264a273..edf57a4 100644 --- a/CyRK/cy/cysolver_api.pyx +++ b/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) diff --git a/CyRK/cy/pysolver.pyx b/CyRK/cy/pysolver.pyx index 939754b..c6d9d12 100644 --- a/CyRK/cy/pysolver.pyx +++ b/CyRK/cy/pysolver.pyx @@ -7,8 +7,6 @@ from libcpp.cmath cimport fmin, fabs from CyRK.utils.memory cimport make_shared from CyRK.cy.cysolver_api cimport find_expected_size, WrapCySolverResult, INF, EPS_100, Y_LIMIT, DY_LIMIT -cimport numpy as np - import numpy as np np.import_array() From 907fad49ed3c1606f240a0feca30cf414d2ae347 Mon Sep 17 00:00:00 2001 From: Jrenaud-Desk Date: Wed, 13 Nov 2024 00:13:00 -0500 Subject: [PATCH 2/9] Implemented dense output for extra outputs --- CHANGES.md | 11 + CyRK/cy/cysolution.cpp | 8 +- CyRK/cy/cysolution.hpp | 4 +- CyRK/cy/cysolver.cpp | 20 +- CyRK/cy/cysolver_api.pyx | 6 +- CyRK/cy/dense.cpp | 61 +- CyRK/cy/dense.hpp | 24 +- CyRK/cy/rk.cpp | 8 +- Documentation/Extra Output.md | 15 +- Tests/Dense Extra Output Testing.ipynb | 891 +++++++++++++++++++++++++ _build_cyrk.py | 2 +- pyproject.toml | 2 +- 12 files changed, 1024 insertions(+), 28 deletions(-) create mode 100644 Tests/Dense Extra Output Testing.ipynb diff --git a/CHANGES.md b/CHANGES.md index 2db25b0..b1505d6 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -2,6 +2,17 @@ ## 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. + +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: diff --git a/CyRK/cy/cysolution.cpp b/CyRK/cy/cysolution.cpp index da6aa55..603c3cf 100644 --- a/CyRK/cy/cysolution.cpp +++ b/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); } diff --git a/CyRK/cy/cysolution.hpp b/CyRK/cy/cysolution.hpp index dd0c50f..08188e2 100644 --- a/CyRK/cy/cysolution.hpp +++ b/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); }; diff --git a/CyRK/cy/cysolver.cpp b/CyRK/cy/cysolver.cpp index 75ecfe1..80ef826 100644 --- a/CyRK/cy/cysolver.cpp +++ b/CyRK/cy/cysolver.cpp @@ -252,9 +252,6 @@ void CySolverBase::reset() this->reset_called = true; } -#include -#include - void CySolverBase::take_step() { if (!this->reset_called) [[unlikely]] @@ -316,8 +313,15 @@ 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->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); @@ -518,7 +522,13 @@ 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->t_now_ptr, + this->y_now_ptr, + this->dy_now_ptr ); } diff --git a/CyRK/cy/cysolver_api.pyx b/CyRK/cy/cysolver_api.pyx index edf57a4..2f96aab 100644 --- a/CyRK/cy/cysolver_api.pyx +++ b/CyRK/cy/cysolver_api.pyx @@ -28,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] @@ -43,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): diff --git a/CyRK/cy/dense.cpp b/CyRK/cy/dense.cpp index 7587122..b9b2009 100644 --- a/CyRK/cy/dense.cpp +++ b/CyRK/cy/dense.cpp @@ -1,15 +1,30 @@ #include "dense.hpp" +#include + 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 cysolver_diffeq_ptr, + 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), + 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) @@ -119,10 +134,48 @@ void CySolverDense::call(double t_interp, double* y_intepret) [[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++) + 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, 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); + // 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); + 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_intepret 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_intepret[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, num_dy); + memcpy(this->cysolver_dy_now_ptr, dy_tmp_ptr, num_dy); } } \ No newline at end of file diff --git a/CyRK/cy/dense.hpp b/CyRK/cy/dense.hpp index e6de7d6..4cbe0fe 100644 --- a/CyRK/cy/dense.hpp +++ b/CyRK/cy/dense.hpp @@ -4,10 +4,15 @@ * S o by avoiding calls to New we greatly improve performance. */ +#include #include #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 +class CySolverBase; + + class CySolverDense { /* Attributes */ @@ -26,6 +31,14 @@ class CySolverDense // y and t state info unsigned int num_y = 0; + unsigned int num_extra = 0; + + // Pointer to the CySolverBase class + CySolverBase* cysolver_instance_ptr = nullptr; + std::function cysolver_diffeq_ptr = nullptr; + double* cysolver_t_now_ptr = nullptr; + double* cysolver_y_now_ptr = nullptr; + double* cysolver_dy_now_ptr = nullptr; // Time step info double step = 0.0; @@ -52,7 +65,14 @@ class CySolverDense 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 cysolver_diffeq_ptr, + double* cysolver_t_now_ptr, + double* cysolver_y_now_ptr, + double* cysolver_dy_now_ptr + ); - virtual void call(double t_interp, double* y_intepret); + virtual void call(double t_interp, double* y_interp_ptr); }; diff --git a/CyRK/cy/rk.cpp b/CyRK/cy/rk.cpp index 70efddf..fd8934c 100644 --- a/CyRK/cy/rk.cpp +++ b/CyRK/cy/rk.cpp @@ -1021,7 +1021,13 @@ CySolverDense* RKSolver::p_dense_output_heap() this->t_now_ptr[0], this->y_old_ptr, this->num_y, - this->len_Pcols); + this->num_extra, + this->len_Pcols, + this, + this->diffeq, + this->t_now_ptr, + this->y_now_ptr, + this->dy_now_ptr); // Update Q this->p_update_Q(dense_output->Q_ptr); diff --git a/Documentation/Extra Output.md b/Documentation/Extra Output.md index 7c2fa36..c2551ac 100644 --- a/Documentation/Extra Output.md +++ b/Documentation/Extra Output.md @@ -32,8 +32,7 @@ The values of these extra parameters are not analyzed by the solver when determi Current limitations of this feature as of v0.4.0: - Only numerical parameters can be used as extra outputs (no strings or booleans). -- All extra outputs must have the same _type_ as the input `y`s. This means if you are using `y`s which are floats, but you need to capture a complex number, -then you will either need to change the `y`s to complex-valued (with a zero imaginary portion) or, a better option, return the real and imaginary portions of the extra parameter separately. +- All extra outputs must have the same _type_ as the input `y`s. (for `pysolve_ivp` and `cysolve_ivp` extra outputs can only be doubles). ## How to use with `CyRK.nbsolve_ivp` (Numba-based) @@ -127,12 +126,18 @@ extra_parameter_0 = result.y[2, :] extra_parameter_1 = result.y[3, :] ``` +# Interpolating extra outputs -## Additional Considerations When Using `t_eval` (for numba-based method) +## Interpolating for cython-based `cysolve_ivp` and `pysolve_ivp` -_This is applicable to either the numba- or cython-based solver._ +When using either `dense_output` or `t_eval` then interpolations must be performed. For the dependent y variables, CyRK will utilize the user-specified differential equation method's approach to interpolation. However, this only works for _dependent_ variables. +If you are capturing extra output they will not be included in this interpolation. Instead, both `cysolve_ivp` and `pysolve_ivp` will perform said dependent-variable interpolation and then use the results to make additional calls (one per interpolation, i.e., once per value in t_eval) to the differential equation to find the values of the extra outputs at the requested time steps. -By setting the `t_eval` argument for either the `nbsolve_ivp` or `cyrk_ode` solver, an interpolation will occur at the end of integration. +This approach is similar to "option 1" discussed in the next section for the numba-based solver. + +## Interpolating for numba-based `nbsolve_ivp` + +By setting the `t_eval` argument for either the `nbsolve_ivp` solver, an interpolation will occur at the end of integration. This uses the solved `y`s and `time_domain` to find a new reduced `y_reduced` at `t_eval` intervals using a method similar to `numpy.interp` function. Since we are potentially storing extra parameters during integration, we need to tell the solvers how to handle any potential interpolation on these new parameters. diff --git a/Tests/Dense Extra Output Testing.ipynb b/Tests/Dense Extra Output Testing.ipynb new file mode 100644 index 0000000..b049e11 --- /dev/null +++ b/Tests/Dense Extra Output Testing.ipynb @@ -0,0 +1,891 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "f42995c9-51f5-42fd-b59c-8ebc78418f64", + "metadata": {}, + "outputs": [], + "source": [ + "import Cython" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6eacd895-3692-4402-9c99-1081a047108b", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext cython" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "155cdcf8-3aea-4564-8d79-2b5fa18258ee", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Content of stdout:\n", + "_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cpp\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\numpy\\core\\include\\numpy\\npy_1_7_deprecated_api.h(14) : Warning Msg: Using deprecated NumPy API, disable it with #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\dense.cpp(134): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(53): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(68): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(98): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(130): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(341): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(184): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(257): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(267): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(273): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(285): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(545): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(50): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(52): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(57): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(75): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(85): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(160): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(226): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(444): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(510): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(574): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(579): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(603): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(621): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(989): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(111): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(265): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(270): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(297): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cpp(24954): warning C4551: function call missing argument list\n", + " Creating library C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cp311-win_amd64.lib and object C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cp311-win_amd64.exp\n", + "Generating code\n", + "Finished generating code\n", + "Integration Done!\n", + "teval i = 0 ys= 4.526273291353615\t|\t2.130934826873202\t|\t0.25\t|\t0.15\n", + "teval i = 1 ys= 1.8716084470161582\t|\t21.05384255519958\t|\t0.25\t|\t2.0\n", + "teval i = 2 ys= 0.13355624397337607\t|\t10.312017312712717\t|\t0.25\t|\t4.0\n", + "teval i = 3 ys= 1.578715913166198\t|\t48.55169176359066\t|\t0.25\t|\t4.55\n", + "teval i = 4 ys= -1.015908770745836\t|\t57.06848507859018\t|\t0.25\t|\t4.95\n", + "INTERP AT r=0.3: 4.526\t|\t2.131\t|\t0.250\t|\t0.150\n", + "INTERP AT r=4.0: 1.872\t|\t21.054\t|\t0.250\t|\t2.000\n", + "INTERP AT r=8.0: 0.134\t|\t10.312\t|\t0.250\t|\t4.000\n", + "INTERP AT r=9.1: 1.579\t|\t48.552\t|\t0.250\t|\t4.550\n", + "INTERP AT r=9.9: -1.016\t|\t57.068\t|\t0.250\t|\t4.950\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "\n", + " \n", + " Cython: _cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.pyx\n", + " \n", + "\n", + "\n", + "

Generated by Cython 3.0.11

\n", + "

\n", + " Yellow lines hint at Python interaction.
\n", + " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", + "

\n", + "
+01: # distutils: language = c++
\n", + "
  __pyx_t_7 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 1, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_t_7);\n",
+       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_7) < 0) __PYX_ERR(0, 1, __pyx_L1_error)\n",
+       "  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "
 02: # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
\n", + "
 03:
\n", + "
+04: import numpy as np
\n", + "
  __pyx_t_7 = __Pyx_ImportDottedModule(__pyx_n_s_numpy, NULL); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 4, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_t_7);\n",
+       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_7) < 0) __PYX_ERR(0, 4, __pyx_L1_error)\n",
+       "  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "
 05: cimport numpy as np
\n", + "
 06:
\n", + "
+07: np.import_array()
\n", + "
  __pyx_t_9 = __pyx_f_5numpy_import_array(); if (unlikely(__pyx_t_9 == ((int)-1))) __PYX_ERR(0, 7, __pyx_L1_error)\n",
+       "
 08:
\n", + "
 09: from CyRK cimport PreEvalFunc, cysolve_ivp, CySolveOutput
\n", + "
+10: cdef void diffeq(double* dy, double t, double* y, const void* args, PreEvalFunc NA) noexcept nogil:
\n", + "
static void __pyx_f_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_diffeq(double *__pyx_v_dy, double __pyx_v_t, double *__pyx_v_y, CYTHON_UNUSED void const *__pyx_v_args, CYTHON_UNUSED PreEvalFunc __pyx_v_NA) {\n",
+       "/* … */\n",
+       "  /* function exit code */\n",
+       "}\n",
+       "
 11:     # Real dy
\n", + "
+12:     dy[0] = 3.1 * t - y[1]
\n", + "
  (__pyx_v_dy[0]) = ((3.1 * __pyx_v_t) - (__pyx_v_y[1]));\n",
+       "
+13:     dy[1] = y[0] * (0.3 * t * y[1])
\n", + "
  (__pyx_v_dy[1]) = ((__pyx_v_y[0]) * ((0.3 * __pyx_v_t) * (__pyx_v_y[1])));\n",
+       "
 14:     # Extra
\n", + "
+15:     dy[2] = 0.25
\n", + "
  (__pyx_v_dy[2]) = 0.25;\n",
+       "
+16:     dy[3] = t / 2.
\n", + "
  (__pyx_v_dy[3]) = (__pyx_v_t / 2.);\n",
+       "
 17:
\n", + "
+18: cdef double[2] t_span = [0.0, 10.0]
\n", + "
  __pyx_t_10[0] = 0.0;\n",
+       "  __pyx_t_10[1] = 10.0;\n",
+       "  memcpy(&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_span[0]), __pyx_t_10, sizeof(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_span[0]) * (2));\n",
+       "
+19: cdef double* t_span_ptr = &t_span[0]
\n", + "
  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_span_ptr = (&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_span[0]));\n",
+       "
 20:
\n", + "
+21: cdef double[2] y0 = [5., 2.]
\n", + "
  __pyx_t_11[0] = 5.;\n",
+       "  __pyx_t_11[1] = 2.;\n",
+       "  memcpy(&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y0[0]), __pyx_t_11, sizeof(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y0[0]) * (2));\n",
+       "
+22: cdef double* y0_ptr = &y0[0]
\n", + "
  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y0_ptr = (&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y0[0]));\n",
+       "
+23: cdef double[5] t_eval = [0.3, 4.0, 8., 9.10, 9.9]
\n", + "
  __pyx_t_12[0] = 0.3;\n",
+       "  __pyx_t_12[1] = 4.0;\n",
+       "  __pyx_t_12[2] = 8.;\n",
+       "  __pyx_t_12[3] = 9.10;\n",
+       "  __pyx_t_12[4] = 9.9;\n",
+       "  memcpy(&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_eval[0]), __pyx_t_12, sizeof(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_eval[0]) * (5));\n",
+       "
+24: cdef double* t_eval_ptr = &t_eval[0]
\n", + "
  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_eval_ptr = (&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_eval[0]));\n",
+       "
+25: cdef int num_t_eval = 5
\n", + "
  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_num_t_eval = 5;\n",
+       "
 26:
\n", + "
 27: cdef CySolveOutput cyresult = \\
\n", + "
+28:     cysolve_ivp(
\n", + "
  __pyx_t_14.__pyx_n = 16;\n",
+       "  __pyx_t_14.method = 1;\n",
+       "  __pyx_t_14.rtol = 1.0e-3;\n",
+       "  __pyx_t_14.atol = 1.0e-6;\n",
+       "  __pyx_t_14.args_ptr = NULL;\n",
+       "  __pyx_t_14.num_extra = 2;\n",
+       "  __pyx_t_14.max_num_steps = 0;\n",
+       "  __pyx_t_14.max_ram_MB = 0x7D0;\n",
+       "  __pyx_t_14.dense_output = 1;\n",
+       "  __pyx_t_14.t_eval = __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_eval_ptr;\n",
+       "  __pyx_t_14.len_t_eval = __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_num_t_eval;\n",
+       "  __pyx_t_14.pre_eval_func = NULL;\n",
+       "  __pyx_t_14.rtols_ptr = NULL;\n",
+       "  __pyx_t_14.atols_ptr = NULL;\n",
+       "  __pyx_t_14.max_step = 10000.;\n",
+       "  __pyx_t_14.first_step = 0.0;\n",
+       "  __pyx_t_14.expected_size = 0;\n",
+       "  __pyx_t_13 = __pyx_f_4CyRK_2cy_12cysolver_api_cysolve_ivp(__pyx_f_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_diffeq, __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_t_span_ptr, __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y0_ptr, 2, &__pyx_t_14); \n",
+       "  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult = __PYX_STD_MOVE_IF_SUPPORTED(__pyx_t_13);\n",
+       "
 29:         diffeq,
\n", + "
 30:         t_span_ptr,
\n", + "
 31:         y0_ptr,
\n", + "
 32:         2,
\n", + "
 33:         1,
\n", + "
 34:         1.0e-3,
\n", + "
 35:         1.0e-6,
\n", + "
 36:         NULL,
\n", + "
 37:         2, # num extra
\n", + "
 38:         0,
\n", + "
 39:         2000,
\n", + "
 40:         True,
\n", + "
 41:         t_eval_ptr,
\n", + "
 42:         num_t_eval,
\n", + "
 43:         NULL,
\n", + "
 44:         NULL,
\n", + "
 45:         NULL,
\n", + "
 46:         10000.,
\n", + "
 47:         0.0,
\n", + "
 48:         0)
\n", + "
 49:
\n", + "
+50: print('')
\n", + "
  __pyx_t_7 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_tuple__23, NULL); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 50, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_t_7);\n",
+       "  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "/* … */\n",
+       "  __pyx_tuple__23 = PyTuple_Pack(1, __pyx_kp_u__22); if (unlikely(!__pyx_tuple__23)) __PYX_ERR(0, 50, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_tuple__23);\n",
+       "  __Pyx_GIVEREF(__pyx_tuple__23);\n",
+       "
+51: print('Integration Done!')
\n", + "
  __pyx_t_7 = __Pyx_PyObject_Call(__pyx_builtin_print, __pyx_tuple__24, NULL); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 51, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_t_7);\n",
+       "  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "/* … */\n",
+       "  __pyx_tuple__24 = PyTuple_Pack(1, __pyx_kp_u_Integration_Done); if (unlikely(!__pyx_tuple__24)) __PYX_ERR(0, 51, __pyx_L1_error)\n",
+       "  __Pyx_GOTREF(__pyx_tuple__24);\n",
+       "  __Pyx_GIVEREF(__pyx_tuple__24);\n",
+       "
 52:
\n", + "
+53: cdef double[4] y_arr = [-999, -999, -999, -999]
\n", + "
  __pyx_t_15[0] = -999.0;\n",
+       "  __pyx_t_15[1] = -999.0;\n",
+       "  __pyx_t_15[2] = -999.0;\n",
+       "  __pyx_t_15[3] = -999.0;\n",
+       "  memcpy(&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr[0]), __pyx_t_15, sizeof(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr[0]) * (4));\n",
+       "
+54: cdef double* y_arr_ptr = &y_arr[0]
\n", + "
  __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr = (&(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr[0]));\n",
+       "
 55:
\n", + "
 56: cdef int i
\n", + "
+57: for i in range(5):
\n", + "
  for (__pyx_t_9 = 0; __pyx_t_9 < 5; __pyx_t_9+=1) {\n",
+       "    __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i = __pyx_t_9;\n",
+       "
+58:     print(f"teval i = {i} ys= {cyresult.get().solution[4*i]}\\t|\\t{cyresult.get().solution[4*i+1]}\\t|\\t{cyresult.get().solution[4*i+2]}\\t|\\t{cyresult.get().solution[4*i+3]}")
\n", + "
    __pyx_t_7 = PyTuple_New(10); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_7);\n",
+       "    __pyx_t_16 = 0;\n",
+       "    __pyx_t_17 = 127;\n",
+       "    __Pyx_INCREF(__pyx_kp_u_teval_i);\n",
+       "    __pyx_t_16 += 10;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u_teval_i);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_kp_u_teval_i);\n",
+       "    __pyx_t_4 = __Pyx_PyUnicode_From_int(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i, 0, ' ', 'd'); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __pyx_t_16 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_4);\n",
+       "    __Pyx_GIVEREF(__pyx_t_4);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_4);\n",
+       "    __pyx_t_4 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u_ys);\n",
+       "    __pyx_t_16 += 5;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u_ys);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 2, __pyx_kp_u_ys);\n",
+       "    __pyx_t_4 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult.get()->solution[(4 * __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i)])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __pyx_t_5 = __Pyx_PyObject_FormatSimple(__pyx_t_4, __pyx_empty_unicode); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) : __pyx_t_17;\n",
+       "    __pyx_t_16 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_5);\n",
+       "    __Pyx_GIVEREF(__pyx_t_5);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 3, __pyx_t_5);\n",
+       "    __pyx_t_5 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_16 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 4, __pyx_kp_u__25);\n",
+       "    __pyx_t_5 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult.get()->solution[((4 * __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i) + 1)])); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __pyx_t_4 = __Pyx_PyObject_FormatSimple(__pyx_t_5, __pyx_empty_unicode); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_4) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_4) : __pyx_t_17;\n",
+       "    __pyx_t_16 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_4);\n",
+       "    __Pyx_GIVEREF(__pyx_t_4);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 5, __pyx_t_4);\n",
+       "    __pyx_t_4 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_16 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 6, __pyx_kp_u__25);\n",
+       "    __pyx_t_4 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult.get()->solution[((4 * __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i) + 2)])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __pyx_t_5 = __Pyx_PyObject_FormatSimple(__pyx_t_4, __pyx_empty_unicode); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) : __pyx_t_17;\n",
+       "    __pyx_t_16 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_5);\n",
+       "    __Pyx_GIVEREF(__pyx_t_5);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 7, __pyx_t_5);\n",
+       "    __pyx_t_5 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_16 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 8, __pyx_kp_u__25);\n",
+       "    __pyx_t_5 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult.get()->solution[((4 * __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_i) + 3)])); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __pyx_t_4 = __Pyx_PyObject_FormatSimple(__pyx_t_5, __pyx_empty_unicode); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_4) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_4) : __pyx_t_17;\n",
+       "    __pyx_t_16 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_4);\n",
+       "    __Pyx_GIVEREF(__pyx_t_4);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_7, 9, __pyx_t_4);\n",
+       "    __pyx_t_4 = 0;\n",
+       "    __pyx_t_4 = __Pyx_PyUnicode_Join(__pyx_t_7, 10, __pyx_t_16, __pyx_t_17); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "    __pyx_t_7 = __Pyx_PyObject_CallOneArg(__pyx_builtin_print, __pyx_t_4); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 58, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_7);\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "  }\n",
+       "
 59:
\n", + "
 60: cdef double rr
\n", + "
+61: for r in [0.3, 4.0, 8., 9.10, 9.9]:
\n", + "
  __pyx_t_7 = __pyx_tuple__26; __Pyx_INCREF(__pyx_t_7);\n",
+       "  __pyx_t_16 = 0;\n",
+       "  for (;;) {\n",
+       "    if (__pyx_t_16 >= 5) break;\n",
+       "    #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS\n",
+       "    __pyx_t_4 = PyTuple_GET_ITEM(__pyx_t_7, __pyx_t_16); __Pyx_INCREF(__pyx_t_4); __pyx_t_16++; if (unlikely((0 < 0))) __PYX_ERR(0, 61, __pyx_L1_error)\n",
+       "    #else\n",
+       "    __pyx_t_4 = __Pyx_PySequence_ITEM(__pyx_t_7, __pyx_t_16); __pyx_t_16++; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 61, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    #endif\n",
+       "    if (PyDict_SetItem(__pyx_d, __pyx_n_s_r, __pyx_t_4) < 0) __PYX_ERR(0, 61, __pyx_L1_error)\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "/* … */\n",
+       "  }\n",
+       "  __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;\n",
+       "
+62:     rr = r
\n", + "
    __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_r); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 62, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __pyx_t_18 = __pyx_PyFloat_AsDouble(__pyx_t_4); if (unlikely((__pyx_t_18 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 62, __pyx_L1_error)\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "    __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_rr = __pyx_t_18;\n",
+       "
+63:     cyresult.get().call(rr, y_arr_ptr)
\n", + "
    __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_cyresult.get()->call(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_rr, __pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr);\n",
+       "
+64:     print(f"INTERP AT r={rr}: {y_arr_ptr[0]:0.3f}\\t|\\t{y_arr_ptr[1]:0.3f}\\t|\\t{y_arr_ptr[2]:0.3f}\\t|\\t{y_arr_ptr[3]:0.3f}")
\n", + "
    __pyx_t_4 = PyTuple_New(10); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __pyx_t_19 = 0;\n",
+       "    __pyx_t_17 = 127;\n",
+       "    __Pyx_INCREF(__pyx_kp_u_INTERP_AT_r);\n",
+       "    __pyx_t_19 += 12;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u_INTERP_AT_r);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_kp_u_INTERP_AT_r);\n",
+       "    __pyx_t_5 = PyFloat_FromDouble(__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_rr); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __pyx_t_20 = __Pyx_PyObject_FormatSimple(__pyx_t_5, __pyx_empty_unicode); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) : __pyx_t_17;\n",
+       "    __pyx_t_19 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_20);\n",
+       "    __Pyx_GIVEREF(__pyx_t_20);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_20);\n",
+       "    __pyx_t_20 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u_);\n",
+       "    __pyx_t_19 += 2;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u_);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 2, __pyx_kp_u_);\n",
+       "    __pyx_t_20 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr[0])); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __pyx_t_5 = __Pyx_PyObject_Format(__pyx_t_20, __pyx_kp_u_0_3f); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __Pyx_DECREF(__pyx_t_20); __pyx_t_20 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) : __pyx_t_17;\n",
+       "    __pyx_t_19 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_5);\n",
+       "    __Pyx_GIVEREF(__pyx_t_5);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 3, __pyx_t_5);\n",
+       "    __pyx_t_5 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_19 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 4, __pyx_kp_u__25);\n",
+       "    __pyx_t_5 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr[1])); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __pyx_t_20 = __Pyx_PyObject_Format(__pyx_t_5, __pyx_kp_u_0_3f); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) : __pyx_t_17;\n",
+       "    __pyx_t_19 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_20);\n",
+       "    __Pyx_GIVEREF(__pyx_t_20);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 5, __pyx_t_20);\n",
+       "    __pyx_t_20 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_19 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 6, __pyx_kp_u__25);\n",
+       "    __pyx_t_20 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr[2])); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __pyx_t_5 = __Pyx_PyObject_Format(__pyx_t_20, __pyx_kp_u_0_3f); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __Pyx_DECREF(__pyx_t_20); __pyx_t_20 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_5) : __pyx_t_17;\n",
+       "    __pyx_t_19 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_5);\n",
+       "    __Pyx_GIVEREF(__pyx_t_5);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 7, __pyx_t_5);\n",
+       "    __pyx_t_5 = 0;\n",
+       "    __Pyx_INCREF(__pyx_kp_u__25);\n",
+       "    __pyx_t_19 += 3;\n",
+       "    __Pyx_GIVEREF(__pyx_kp_u__25);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 8, __pyx_kp_u__25);\n",
+       "    __pyx_t_5 = PyFloat_FromDouble((__pyx_v_54_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a_y_arr_ptr[3])); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_5);\n",
+       "    __pyx_t_20 = __Pyx_PyObject_Format(__pyx_t_5, __pyx_kp_u_0_3f); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;\n",
+       "    __pyx_t_17 = (__Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) > __pyx_t_17) ? __Pyx_PyUnicode_MAX_CHAR_VALUE(__pyx_t_20) : __pyx_t_17;\n",
+       "    __pyx_t_19 += __Pyx_PyUnicode_GET_LENGTH(__pyx_t_20);\n",
+       "    __Pyx_GIVEREF(__pyx_t_20);\n",
+       "    PyTuple_SET_ITEM(__pyx_t_4, 9, __pyx_t_20);\n",
+       "    __pyx_t_20 = 0;\n",
+       "    __pyx_t_20 = __Pyx_PyUnicode_Join(__pyx_t_4, 10, __pyx_t_19, __pyx_t_17); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_20);\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "    __pyx_t_4 = __Pyx_PyObject_CallOneArg(__pyx_builtin_print, __pyx_t_20); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 64, __pyx_L1_error)\n",
+       "    __Pyx_GOTREF(__pyx_t_4);\n",
+       "    __Pyx_DECREF(__pyx_t_20); __pyx_t_20 = 0;\n",
+       "    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;\n",
+       "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%cython -a -f\n", + "# distutils: language = c++\n", + "# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False\n", + "\n", + "import numpy as np\n", + "cimport numpy as np\n", + "\n", + "np.import_array()\n", + "\n", + "from CyRK cimport PreEvalFunc, cysolve_ivp, CySolveOutput\n", + "cdef void diffeq(double* dy, double t, double* y, const void* args, PreEvalFunc NA) noexcept nogil:\n", + " # Real dy\n", + " dy[0] = 3.1 * t - y[1]\n", + " dy[1] = y[0] * (0.3 * t * y[1])\n", + " # Extra\n", + " dy[2] = 0.25\n", + " dy[3] = t / 2.\n", + "\n", + "cdef double[2] t_span = [0.0, 10.0]\n", + "cdef double* t_span_ptr = &t_span[0]\n", + "\n", + "cdef double[2] y0 = [5., 2.]\n", + "cdef double* y0_ptr = &y0[0]\n", + "cdef double[5] t_eval = [0.3, 4.0, 8., 9.10, 9.9]\n", + "cdef double* t_eval_ptr = &t_eval[0]\n", + "cdef int num_t_eval = 5\n", + "\n", + "cdef CySolveOutput cyresult = \\\n", + " cysolve_ivp(\n", + " diffeq,\n", + " t_span_ptr,\n", + " y0_ptr,\n", + " 2,\n", + " 1,\n", + " 1.0e-3,\n", + " 1.0e-6,\n", + " NULL,\n", + " 2, # num extra\n", + " 0,\n", + " 2000,\n", + " True,\n", + " t_eval_ptr,\n", + " num_t_eval,\n", + " NULL,\n", + " NULL,\n", + " NULL,\n", + " 10000.,\n", + " 0.0,\n", + " 0)\n", + "\n", + "print('')\n", + "print('Integration Done!')\n", + "\n", + "cdef double[4] y_arr = [-999, -999, -999, -999]\n", + "cdef double* y_arr_ptr = &y_arr[0]\n", + "\n", + "cdef int i\n", + "for i in range(5):\n", + " print(f\"teval i = {i} ys= {cyresult.get().solution[4*i]}\\t|\\t{cyresult.get().solution[4*i+1]}\\t|\\t{cyresult.get().solution[4*i+2]}\\t|\\t{cyresult.get().solution[4*i+3]}\")\n", + "\n", + "cdef double rr\n", + "for r in [0.3, 4.0, 8., 9.10, 9.9]:\n", + " rr = r\n", + " cyresult.get().call(rr, y_arr_ptr)\n", + " print(f\"INTERP AT r={rr}: {y_arr_ptr[0]:0.3f}\\t|\\t{y_arr_ptr[1]:0.3f}\\t|\\t{y_arr_ptr[2]:0.3f}\\t|\\t{y_arr_ptr[3]:0.3f}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c4a6148-51a4-4013-97bd-2062746107bb", + "metadata": {}, + "outputs": [], + "source": [ + "# 0.11.1 result:\n", + "Finished generating code\n", + "Integration Done!\n", + "INTERP AT r=0.3: 4.526\t|\t2.131\t|\t-999.000\t|\t-999.000\n", + "INTERP AT r=4.0: 1.872\t|\t21.054\t|\t-999.000\t|\t-999.000\n", + "INTERP AT r=8.22: 2.500\t|\t23.547\t|\t-999.000\t|\t-999.000\n", + "INTERP AT r=9.1: 1.579\t|\t48.552\t|\t-999.000\t|\t-999.000" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/_build_cyrk.py b/_build_cyrk.py index 84c2050..70771f6 100644 --- a/_build_cyrk.py +++ b/_build_cyrk.py @@ -11,7 +11,7 @@ import numpy as np num_procs = os.cpu_count() -num_threads = int(math.floor(num_procs * 0.75)) +num_threads = max(1, num_procs - 1) install_platform = platform.system() diff --git a/pyproject.toml b/pyproject.toml index 329e9ab..dfb889e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name='CyRK' -version = '0.11.2a0.dev0' +version = '0.11.2a1.dev3' description='Runge-Kutta ODE Integrator Implemented in Cython and Numba.' authors= [ {name = 'Joe P. Renaud', email = 'joe.p.renaud@gmail.com'} From 30005a16ab0b63bb67080573722ebf8e6f5e6c4e Mon Sep 17 00:00:00 2001 From: Jrenaud-Desk Date: Wed, 13 Nov 2024 07:56:30 -0500 Subject: [PATCH 3/9] Fixed issue where python may deconstruct the diffeq while dense output still had reference and could call it --- CyRK/cy/cysolver.cpp | 2 + CyRK/cy/cysolver_api.pyx | 2 +- CyRK/cy/dense.cpp | 36 +- CyRK/cy/dense.hpp | 3 + CyRK/cy/pysolver_cyhook.pyx | 2 + CyRK/cy/rk.cpp | 1 + Tests/Dense Extra Output Testing.ipynb | 620 ++++++++++++++++--------- pyproject.toml | 2 +- 8 files changed, 433 insertions(+), 235 deletions(-) diff --git a/CyRK/cy/cysolver.cpp b/CyRK/cy/cysolver.cpp index 80ef826..ffe4f02 100644 --- a/CyRK/cy/cysolver.cpp +++ b/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 diff --git a/CyRK/cy/cysolver_api.pyx b/CyRK/cy/cysolver_api.pyx index 2f96aab..610ea60 100644 --- a/CyRK/cy/cysolver_api.pyx +++ b/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) # ===================================================================================================================== diff --git a/CyRK/cy/dense.cpp b/CyRK/cy/dense.cpp index b9b2009..f0865aa 100644 --- a/CyRK/cy/dense.cpp +++ b/CyRK/cy/dense.cpp @@ -13,6 +13,7 @@ CySolverDense::CySolverDense( unsigned int Q_order, CySolverBase* cysolver_instance_ptr, std::function 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); } } \ No newline at end of file diff --git a/CyRK/cy/dense.hpp b/CyRK/cy/dense.hpp index 4cbe0fe..1833da1 100644 --- a/CyRK/cy/dense.hpp +++ b/CyRK/cy/dense.hpp @@ -7,6 +7,7 @@ #include #include +#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 cysolver_diffeq_ptr, + PyObject* cython_extension_class_instance, double* cysolver_t_now_ptr, double* cysolver_y_now_ptr, double* cysolver_dy_now_ptr diff --git a/CyRK/cy/pysolver_cyhook.pyx b/CyRK/cy/pysolver_cyhook.pyx index 9658bd5..c70da68 100644 --- a/CyRK/cy/pysolver_cyhook.pyx +++ b/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) diff --git a/CyRK/cy/rk.cpp b/CyRK/cy/rk.cpp index fd8934c..4137ee6 100644 --- a/CyRK/cy/rk.cpp +++ b/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); diff --git a/Tests/Dense Extra Output Testing.ipynb b/Tests/Dense Extra Output Testing.ipynb index b049e11..598c71c 100644 --- a/Tests/Dense Extra Output Testing.ipynb +++ b/Tests/Dense Extra Output Testing.ipynb @@ -22,18 +22,20 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 3, "id": "155cdcf8-3aea-4564-8d79-2b5fa18258ee", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Content of stdout:\n", - "_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cpp\n", + "_cython_magic_8ec2014aa0d4675425dd699e02c0b83159c42fcf.cpp\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\numpy\\core\\include\\numpy\\npy_1_7_deprecated_api.h(14) : Warning Msg: Using deprecated NumPy API, disable it with #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION\n", - "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\dense.cpp(134): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\dense.cpp(144): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(53): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(68): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolution.cpp(98): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", @@ -44,7 +46,7 @@ "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(267): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(273): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(285): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", - "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(545): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", + "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolver.cpp(547): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(50): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(52): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\rk.cpp(57): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", @@ -63,8 +65,8 @@ "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(265): warning C5051: attribute [[unlikely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(270): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", "C:\\Users\\joepr\\anaconda3\\envs\\cyrk11py311\\Lib\\site-packages\\CyRK\\cy\\cysolve.cpp(297): warning C5051: attribute [[likely]] requires at least '/std:c++20'; ignored\n", - "C:\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cpp(24954): warning C4551: function call missing argument list\n", - " Creating library C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cp311-win_amd64.lib and object C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.cp311-win_amd64.exp\n", + "C:\\Users\\joepr\\.ipython\\cython\\_cython_magic_8ec2014aa0d4675425dd699e02c0b83159c42fcf.cpp(25001): warning C4551: function call missing argument list\n", + " Creating library C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_8ec2014aa0d4675425dd699e02c0b83159c42fcf.cp311-win_amd64.lib and object C:\\Users\\joepr\\.ipython\\cython\\Users\\joepr\\.ipython\\cython\\_cython_magic_8ec2014aa0d4675425dd699e02c0b83159c42fcf.cp311-win_amd64.exp\n", "Generating code\n", "Finished generating code\n", "Integration Done!\n", @@ -88,7 +90,7 @@ "\n", "\n", " \n", - " Cython: _cython_magic_9308ae100906791caefb8b06e931293b7a570f1a.pyx\n", + " Cython: _cython_magic_8ec2014aa0d4675425dd699e02c0b83159c42fcf.pyx\n", "