Updated documentation and plots

CyRK/__init__.pxd

@@ -1 +1 @@
-from CyRK.cy.cysolverNew cimport cysolve_ivp, cysolve_ivp_gil, DiffeqFuncType, CySolverResult, WrapCySolverResult, CySolverBase, CySolveOutput, RK23_METHOD_INT, RK45_METHOD_INT, DOP853_METHOD_INT
+from CyRK.cy.cysolverNew cimport cysolve_ivp, cysolve_ivp_gil, DiffeqFuncType, PreEvalFunc, CySolverResult, WrapCySolverResult, CySolverBase, CySolveOutput, RK23_METHOD_INT, RK45_METHOD_INT, DOP853_METHOD_INT

Documentation/Advanced CySolver.md

@@ -1,29 +1,28 @@
-# TODO: This document is very sparse. There are currently many ways to interface with CyRK's C++ backend as well as other features and gotchas. This document will be expanded in the future to highlight all of these. In the mean time, please feel free to post questions as GitHub Issues.
+# TODO: This document needs expansion. There are currently many ways to interface with CyRK's C++ backend as well as other features and gotchas. This document will be expanded in the future to highlight all of these. In the mean time, please feel free to post questions as GitHub Issues.
 
-_Many of the examples below can be found in the interactive "Demos/Advanced CySolver Examples.ipynb"._
+_Many of the examples below can be found in the interactive "Demos/Advanced CySolver Examples.ipynb" jupyter notebook._
 
 # Arbitrary Additional Arguments
 
-`cysolve_ivp` allows users to specify arbitrary additional arguments which are passed to the differential equation at each evaluation. Many problems require additional arguments and they often take the form of numerical parameters (double precision floating point numbers). However, occasionally more advanced information may need to be passed to the diffeq such as boolean flags, complex numbers, and perhaps even strings or whole other classes or structs. Below we outline how to utilize this feature and any limitations.
+`cysolve_ivp` allows users to specify arbitrary additional arguments which are passed to the differential equation function at each evaluation. Many problems require additional arguments that often take the form of numerical parameters (double precision floating point numbers). However, occasionally more advanced information may be required such as: boolean flags, complex numbers, and perhaps even strings or whole other classes and structs. Below we outline how to utilize arbitrary types in CySolver's additional arguments and discuss any limitations.
 
 Before we describe arbitrary arguments, first we will outline basic argument usage assuming a user has a differential equation that requires an array of floating point numbers.
+
+_Note: This assumes you are working in a Jupyter notebook. If you are not then you can exclude the "%%cython --force" at the beginning of each cell._
+
 ```cython
 %%cython --force 
 # distutils: language = c++
-# distutils: extra_compile_args = /std:c++20
 # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
 import numpy as np
 cimport numpy as np
 np.import_array()
 
 from libc.math cimport sin
 
-from CyRK.cy.cysolverNew cimport (
-    cysolve_ivp, find_expected_size, WrapCySolverResult, DiffeqFuncType, MAX_STEP, EPS_100, INF,
-    CySolverResult, CySolveOutput
-    )
+from CyRK.cy.cysolverNew cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType,MAX_STEP, CySolveOutput, PreEvalFunc
 
-cdef void pendulum_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, PreEvalFunc pre_eval_func) noexcept nogil:
     # Arguments are passed in as a void pointer to preallocated and instantiated memory. 
     # This memory could be stack or heap allocated. If it is not valid then you will run into seg faults.
     cdef double* args_dbl_ptr = <double*>args_ptr
@@ -87,12 +86,11 @@ print("\n\nIntegration success =", result_reg.success, "\n\tNumber of adaptive t
 print("Integration message:", result_reg.message)
 ```
 
-Now assume that we have a new diffeq that is similar to before but it allows for atmospheric drag in the system. For the purposes of demonstration, assume the user wants to turn on or off that drag by passing in a boolean flag. The previous example will not work because the arg array is an array of doubles and can not contain other types of data such as bools. Instead the user needs to build a custom data structure and pass the address of that structure to cysolve_ivp. Finally, inside the diffeq function the arg pointer must be cast back to the known data type in order to access the data.
+Now assume that we have a new diffeq that is similar to before but it allows for atmospheric drag in the system. For the purposes of demonstration, assume the user wants to turn on or off that drag by passing in a boolean flag. The previous example will not work because the `args_arr` is an array of doubles and can not contain other types of data such as booleans. Instead the user will need to build a custom data structure and pass the address of that structure to `cysolve_ivp`. Finally, inside the diffeq function the arg pointer must be cast back to that known data type in order to access the data.
 
 ```cython
 %%cython --force 
 # distutils: language = c++
-# distutils: extra_compile_args = /std:c++20
 # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
 import numpy as np
 cimport numpy as np
@@ -101,10 +99,7 @@ np.import_array()
 from libc.math cimport sin
 from libcpp cimport bool as cpp_bool
 
-from CyRK.cy.cysolverNew cimport (
-    cysolve_ivp, find_expected_size, WrapCySolverResult, DiffeqFuncType, MAX_STEP, EPS_100, INF,
-    CySolverResult, CySolveOutput
-    )
+from CyRK.cy.cysolverNew cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType,MAX_STEP, CySolveOutput, PreEvalFunc
 
 cdef struct PendulumArgs:
     # Structure that contains heterogeneous types
@@ -113,7 +108,7 @@ cdef struct PendulumArgs:
     double l
     double m
 
-cdef void pendulum_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, PreEvalFunc pre_eval_func) noexcept nogil:
     # Arg pointer still must be listed as a void pointer or it will not work with cysolve_ivp.
     # But now the user can recast that void pointer to the structure they wish.
     cdef PendulumArgs* pendulum_args_ptr = <PendulumArgs*>args_ptr
@@ -183,8 +178,131 @@ print("Integration message:", result_drag.message)
 
 # Pre-Evaluation Function
 
-It is occasionally advantageous for users to define differential equation functions that utilize a "pre-evaluation" function that take in the current state and provide parameters which are then used to find dydt. This allows different models to be defined as modified versions of this pre-eval function without changing the rest of the differential equation.
+It is occasionally advantageous for users to define differential equation functions that utilize a "pre-evaluation" function that will use the current state to perform calculations that are then used by the diffeq function to find dydt. While this functionality could be hard coded into the diffeq, having a pre-eval function allows for different models to be defined without changing the rest of the differential equation. A large set of different pre-eval functions could then be passed in to subsequent runs of the solver to compare various models. 
+
+Consider the previous example of a pendulum that is experiencing drag versus one that is not. The equations of motion for the pendulum are only slightly modified between the two cases. Rather than writing a different differential equation, or build an additional argument structure as we used above, we can instead write one diffeq and two different pre-evaluation functions. 
+
+```cython
+%%cython --force 
+# distutils: language = c++
+# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
+import numpy as np
+cimport numpy as np
+np.import_array()
+
+from libc.math cimport sin
+
+from CyRK.cy.cysolverNew cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType,MAX_STEP, CySolveOutput, PreEvalFunc
+
+cdef void pendulum_preeval_nodrag(void* output_ptr, double time, double* y_ptr, const void* args_ptr) noexcept nogil:
+    # Unpack args (in this example we do not need these but they follow the same rules as the Arbitrary Args section discussed above)
+    cdef double* args_dbl_ptr = <double*>args_ptr
+
+    # External torque
+    cdef double torque = 0.1 * sin(time)
 
-Consider the example of a pendulum that is experiencing drag and one that is not. The equations of motion for the pendulum are only slightly modified between the two cases. Rather than writing a different differential equation we can instead write one diffeq and two different pre-evaluation functions. 
+    # Convert output pointer to double pointer so we can store data
+    cdef double* output_dbl_ptr = <double*>output_ptr
+    output_dbl_ptr[0] = torque
+    output_dbl_ptr[1] = 0.0  # No Drag
+
+cdef void pendulum_preeval_withdrag(void* output_ptr, double time, double* y_ptr, const void* args_ptr) noexcept nogil:
+    # Unpack args (in this example we do not need these but they follow the same rules as the Arbitrary Args section discussed above)
+    cdef double* args_dbl_ptr = <double*>args_ptr
+
+    # External torque
+    cdef double torque = 0.1 * sin(time)
+
+    # Convert output pointer to double pointer so we can store data
+    cdef double* output_dbl_ptr = <double*>output_ptr
+    output_dbl_ptr[0] = torque
+    output_dbl_ptr[1] = -1.5 * y_ptr[1]  # With Drag
+
+
+cdef void pendulum_preeval_diffeq(double* dy_ptr, double t, double* y_ptr, const void* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil:
+
+    # Build other parameters that do not depend on the pre-eval func
+    cdef double* args_dbl_ptr = <double*>args_ptr
+    cdef double l = args_dbl_ptr[0]
+    cdef double m = args_dbl_ptr[1]
+    cdef double g = args_dbl_ptr[2]
+
+    cdef double coeff_1 = (-3. * g / (2. * l))
+    cdef double coeff_2 = (3. / (m * l**2))
+
+    # Make stack allocated storage for pre eval output
+    cdef double[4] pre_eval_storage
+    cdef double* pre_eval_storage_ptr = &pre_eval_storage[0]
+
+    # Cast storage to void so we can call function
+    cdef void* pre_eval_storage_void_ptr = <void*>pre_eval_storage_ptr
+
+    # Call Pre-Eval Function
+    # Note that even though CyRK calls this function a "pre-eval" function, it can be placed anywhere inside the diffeq function. 
+    pre_eval_func(pre_eval_storage_void_ptr, t, y_ptr, args_ptr)
+
+    cdef double y0 = y_ptr[0]
+    cdef double y1 = y_ptr[1]
+
+    # Use results of pre-eval function to update dy. Note that we are using the double* not the void* here.
+    # Even though pre_eval_func was passed the void* it updated the memory that the double* pointed to so we can use it below.
+    dy_ptr[0] = y1
+    dy_ptr[1] = coeff_1 * sin(y0) + coeff_2 * pre_eval_storage_ptr[0] + pre_eval_storage_ptr[1]
+
+
+# Now we define a function that builds our inputs and calls `cysolve_ivp` to solve the problem
+def run():
+    cdef DiffeqFuncType diffeq = pendulum_preeval_diffeq
+    
+    # Setup pointer to pre-eval function
+    cdef PreEvalFunc pre_eval_func
+    
+    cdef bint use_drag = True
+    if use_drag:
+        pre_eval_func = pendulum_preeval_withdrag
+    else:
+        pre_eval_func = pendulum_preeval_nodrag
+    
+    # Define time domain
+    cdef double[2] time_span_arr = [0.0, 10.0]
+    cdef double* t_span_ptr = &time_span_arr[0]
+    
+    # Define initial conditions
+    cdef double[2] y0_arr = [0.01, 0.0]
+    cdef double* y0_ptr = &y0_arr[0]
+    
+    # Define our arguments.
+    cdef double[3] args_arr = [1.0, 1.0, 9.81]
+    cdef double* args_dbl_ptr = &args_arr[0]
+
+    # To work with cysolve_ivp, we must cast the args ptr to a void pointer
+    cdef void* args_ptr = <void*>args_dbl_ptr
+
+    cdef CySolveOutput result = cysolve_ivp(
+        diffeq,
+        t_span_ptr,
+        y0_ptr,
+        2,  # Number of dependent variables
+        method = 1, # Integration method
+        rtol = 1.0e-6,
+        atol = 1.0e-8,
+        args_ptr = args_ptr,
+        num_extra = 0,
+        max_num_steps = 100_000_000,
+        max_ram_MB = 2000,
+        dense_output = False,
+        t_eval = NULL,
+        len_t_eval = 0,
+        pre_eval_func = pre_eval_func
+        )
+
+    # If we want to pass the solution back to python we need to wrap it in a CyRK `WrapCySolverResult` class.
+    cdef WrapCySolverResult pysafe_result = WrapCySolverResult()
+    pysafe_result.set_cyresult_pointer(result)
+
+    return pysafe_result
 
-# TODO: Left off: change demo example to have two pre-eval funcs that match the above description. finish this section. 
+result_preeval = run()
+print("\n\nIntegration success =", result_preeval.success, "\n\tNumber of adaptive time steps required:", result_preeval.size)
+print("Integration message:", result_preeval.message)
+```

README.md

@@ -11,7 +11,7 @@
 
 ---
 
-<a href="https://github.com/jrenaud90/CyRK/releases"><img src="https://img.shields.io/badge/CyRK-0.10.0 Alpha-orange" alt="CyRK Version 0.10.0 Alpha" /></a>
+<a href="https://github.com/jrenaud90/CyRK/releases"><img src="https://img.shields.io/badge/CyRK-0.10.1 Alpha-orange" alt="CyRK Version 0.10.1 Alpha" /></a>
 
 
 **Runge-Kutta ODE Integrator Implemented in Cython and Numba**
@@ -28,7 +28,7 @@ The [cython-based](https://cython.org/) `cysolver_ivp` function that works with
 An additional benefit of the two cython implementations is that they are pre-compiled. This avoids most of the start-up performance hit experienced by just-in-time compilers like numba.
 
 
-<img style="text-align: center" src="https://github.com/jrenaud90/CyRK/blob/main/Benchmarks/CyRK_SciPy_Compare_predprey_v0-10-0.png" alt="CyRK Performance Graphic" />
+<img style="text-align: center" src="https://github.com/jrenaud90/CyRK/blob/main/Benchmarks/CyRK_SciPy_Compare_predprey_v0-10-1.png" alt="CyRK Performance Graphic" />
 
 ## Installation
 
@@ -225,11 +225,20 @@ np.import_array()
 # Currently, CyRK only allows additional arguments to be passed in as a double array pointer (they all must be of type double). Mixed type args will be explored in the future if there is demand for it (make a GitHub issue if you'd like to see this feature!).
 # The "noexcept nogil" tells cython that the Python Global Interpretor Lock is not required, and that no exceptions should be raised by the code within this function (both improve performance).
 # If you do need the gil for your differential equation then you must use the `cysolve_ivp_gil` function instead of `cysolve_ivp`
-cdef void cython_diffeq(double* dy, double t, double* y, const double* args) noexcept nogil:
+
+# Import the required functions from CyRK
+from CyRK cimport cysolve_ivp, DiffeqFuncType, WrapCySolverResult, CySolveOutput, PreEvalFunc
+
+# Note that currently you must provide the "const void* args, PreEvalFunc pre_eval_func" as inputs even if they are unused.
+# See "Advanced CySolver.md" in the documentation for information about these parameters.
+cdef void cython_diffeq(double* dy, double t, double* y, const void* args, PreEvalFunc pre_eval_func) noexcept nogil:
     
     # Unpack args
-    cdef double a = args[0]
-    cdef double b = args[1]
+    # CySolver assumes an arbitrary data type for additional arguments. So we must cast them to the array of 
+    # doubles that we want to use for this equation
+    cdef double* args_as_dbls = <double*>args
+    cdef double a = args_as_dbls[0]
+    cdef double b = args_as_dbls[1]
     
     # Build Coeffs
     cdef double coeff_1 = (1. - a * y[1])
@@ -263,7 +272,9 @@ def run_cysolver(tuple t_span, double[::1] y0):
 
     # Assume constant args
     cdef double[2] args   = [0.01, 0.02]
-    cdef double* args_ptr = &args[0]
+    cdef double* args_dbl_ptr = &args[0]
+    # Need to cast the arg double pointer to void
+    cdef void* args_ptr = <void*>args_dbl_ptr
 
     # Run the integrator!
     cdef CySolveOutput result = cysolve_ivp(

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name='CyRK'
-version = '0.10.1a0.dev2'
+version = '0.10.1'
 description='Runge-Kutta ODE Integrator Implemented in Cython and Numba.'
 authors= [
     {name = 'Joe P. Renaud', email = '[email protected]'}