Advanced CySolver Examples

This notebook showcases some of the advanced functionality of the CySolver backend distributed with CyRK. In order to run these cells you need to make sure Cython, Jupyter, CyRK, and Numpy are installed.

import matplotlib.pyplot as plt
plt.style.use('dark_background')
# Hack the cython magic so that it can set the correct C++ standard.
%load_ext Cython
import sys
import platform
import numpy as np
import CyRK
print("CyRK version:", CyRK.version)

The Cython extension is already loaded. To reload it, use:
  %reload_ext Cython
CyRK version: 0.17.0

Hack IPython Magic

In order to make this notebook work across platforms without changes we need to hack the cython cells to include some additional header information. This is done using the script jupyter_cyhack.py file contained in the same directory as these demos. Hopefully it works fine for you, but keep in mind that the hack below is required to get these cells to work. So if you are making your own cythonized notebook you will want to see what headers should be included for your operating system.

from Cython.Build import IpythonMagic

from jupyter_cyhack import build_hack

# Get the current cython Ipython parser.
original_cython_parser = IpythonMagic.CythonMagics.cython
# Apply hacks to headers and compile lines.
patched_cython_parser = build_hack(original_cython_parser)
# Patch it in
IpythonMagic.CythonMagics.cython = patched_cython_parser
# Reload
ip = get_ipython()
ip.register_magics(IpythonMagic.CythonMagics)

Parallelizing cysolve_ivp

Below is an example of how you can parallelize over cysolve_ivp using Cython’s prange.

Note that in this example we use cysolve_ivp_noreturn instead of cysolve_ivp. It should work with both. Please report a bug if it does not!

In this example we build a unique CySolveOutput for each job and batch send them to threads to work through the jobs. In the next example we show how you can instead have one CySolveOutput per thread.

%%cython --force
# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
from cython.parallel cimport prange

from libcpp.vector cimport vector
from libcpp.memory cimport unique_ptr, make_unique
from libcpp.utility cimport move

from CyRK.cy.cysolver_api cimport cysolve_ivp_noreturn, CySolverResult, CySolveOutput, ODEMethod
from CyRK.cy.cysolver_test cimport lotkavolterra_diffeq

cdef extern from "<chrono>" namespace "std::chrono" nogil:
    cdef cppclass microseconds:
        long long count()

    # We need a generic duration type to handle the result of subtraction
    cdef cppclass duration "std::chrono::high_resolution_clock::duration":
        pass

    microseconds duration_cast_microseconds "std::chrono::duration_cast<std::chrono::microseconds>"(duration)

cdef extern from "<chrono>" namespace "std::chrono::high_resolution_clock" nogil:
    cdef cppclass time_point:
        duration operator-(time_point)

    time_point now()

# Note the example below uses a `def` vs. `cdef` function but it should work with either

def run_prange(int num_threads = 2):
    """Run a prange test where args are common across all runs."""

    cdef int num_runs = 1000
    cdef double t_start = 0.0
    # Pick a large time so each time step takes a significant amount of time so it can benefit from parallelize.
    cdef double t_end = 2000.0

    cdef vector[vector[double]] y0_vecs = vector[vector[double]]()
    y0_vecs.resize(num_runs)
    cdef size_t i
    for i in range(num_runs):
        # Set initial conditions; all the same for these tests
        # 2 dependent variables
        y0_vecs[i].resize(2)
        y0_vecs[i][0] = 10.0
        y0_vecs[i][1] = 5.0

    # --- CHANGE: Single constant vector instead of vector[vector] ---
    cdef vector[char] constant_args = vector[char]()
    # arg vector
    constant_args.resize(4 * 8)  # Number of double * size of double

    # Recast <char*> args to double pointer for storage of args
    cdef double* arg_dbl_ptr = <double*> constant_args.data()
    arg_dbl_ptr[0] = 1.5
    arg_dbl_ptr[1] = 1.0
    arg_dbl_ptr[2] = 3.0
    arg_dbl_ptr[3] = 1.0

    cdef vector[CySolveOutput] results = vector[CySolveOutput]()
    results.resize(num_runs)
    for i in range(num_runs):
        results[i] = move(make_unique[CySolverResult](ODEMethod.RK45))

    # Main Loop
    cdef time_point start_time, end_time
    cdef Py_ssize_t p_i  # prange has to used a signed integer type

    start_time = now()
    for p_i in prange(num_runs, nogil=True, num_threads=num_threads, use_threads_if=num_threads>1, schedule='static'):
        cysolve_ivp_noreturn(
            results[p_i].get(), # Pass in CySolverResult pointer which is found with CySolveOutput::get()
            lotkavolterra_diffeq,
            t_start,
            t_end,
            y0_vecs[p_i],
            1.0e-5, # rtol
            1.0e-8, # atol
            constant_args # Pass the single shared vector
        )
    end_time = now()

    # Check results
    cdef CySolverResult* reference_result_ptr = results[0].get()
    # Check the reference first.
    assert reference_result_ptr.success
    assert reference_result_ptr.size > 0

    cdef CySolverResult* check_result_ptr
    cdef size_t j
    for i in range(num_runs):

        if i == 0:
            continue
        check_result_ptr = results[i].get()

        assert check_result_ptr.success
        assert check_result_ptr.size > 0
        assert reference_result_ptr.size == check_result_ptr.size

        for j in range(reference_result_ptr.size):
            assert check_result_ptr.time_domain_vec[j] == reference_result_ptr.time_domain_vec[j]
            assert check_result_ptr.solution[2 * j + 0] == reference_result_ptr.solution[2 * j + 0]
            assert check_result_ptr.solution[2 * j + 1] == reference_result_ptr.solution[2 * j + 1]


    cdef double loop_time_ms = <double>(<int>duration_cast_microseconds(end_time - start_time).count()) / 1000.0
    print(f"CyRK prange test `run_prange_constant_args_test` finished in {loop_time_ms:0.1f} ms (If you are seeing this then tests passed too).")
    print(f"\tAvg time: {loop_time_ms/num_runs:0.1f}ms/loop.")

    return loop_time_ms

import os
cdef int cpu_count = os.cpu_count()
cdef double result
cdef int threads
cdef double t1_result, last_result
print("\n\nStarting prange Demo!")
for threads in range(1, min(cpu_count, 8)):
    print(f"Working with N={threads}...")
    result = run_prange(threads)
    if threads == 1:
        print(f"Total Time: {result:0.3f}ms.\n\n")
        t1_result = result
    else:
        print(f"Total Time: {result:0.3f}ms ({100.0*(result - last_result)/last_result:0.1f} % change from N={threads-1}; {100.0*(result - t1_result)/t1_result:0.1f} % change from N=1).\n\n")
    last_result = result

Content of stdout:
_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(25554): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(25555): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(25678): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(28820): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(28827): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(29327): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cpp(29328): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_b079f79e5601112d690eb8710dc64b86fc48f37951916b65c17d7897418586c0.cp312-win_amd64.exp
Generating code
Finished generating code

Starting prange Demo!
Working with N=1...
CyRK prange test `run_prange_constant_args_test` finished in 1900.1 ms (If you are seeing this then tests passed too).
        Avg time: 1.9ms/loop.
Total Time: 1900.120ms.


Working with N=2...
CyRK prange test `run_prange_constant_args_test` finished in 984.2 ms (If you are seeing this then tests passed too).
        Avg time: 1.0ms/loop.
Total Time: 984.242ms (-48.2 % change from N=1; -48.2 % change from N=1).


Working with N=3...
CyRK prange test `run_prange_constant_args_test` finished in 641.9 ms (If you are seeing this then tests passed too).
        Avg time: 0.6ms/loop.
Total Time: 641.888ms (-34.8 % change from N=2; -66.2 % change from N=1).


Working with N=4...
CyRK prange test `run_prange_constant_args_test` finished in 491.1 ms (If you are seeing this then tests passed too).
        Avg time: 0.5ms/loop.
Total Time: 491.051ms (-23.5 % change from N=3; -74.2 % change from N=1).


Working with N=5...
CyRK prange test `run_prange_constant_args_test` finished in 389.7 ms (If you are seeing this then tests passed too).
        Avg time: 0.4ms/loop.
Total Time: 389.738ms (-20.6 % change from N=4; -79.5 % change from N=1).


Working with N=6...
CyRK prange test `run_prange_constant_args_test` finished in 338.3 ms (If you are seeing this then tests passed too).
        Avg time: 0.3ms/loop.
Total Time: 338.321ms (-13.2 % change from N=5; -82.2 % change from N=1).


Working with N=7...
CyRK prange test `run_prange_constant_args_test` finished in 298.4 ms (If you are seeing this then tests passed too).
        Avg time: 0.3ms/loop.
Total Time: 298.442ms (-11.8 % change from N=6; -84.3 % change from N=1).

Batching CySolveOutputs Over Multiple Threads

In the parallelized example above we created a CySolveOutput for each of our runs. This is safe because each thread is guaranteed to only work with a unique CySolveOutput so it can write data as it wishes. However, CySolveOutput does carry a little overhead that may cause memory issues for problems with many dependent y variables.

It may be advantageous to make use of “Re-Uses” (see other examples in this demo) with cysolve_ivp_noreturn. Where only one CySolveOutput per thread is built and used. However, to do any useful work we still have to build storage arrays for our y, t, and perhaps other metadata like integration success, for each job. That is a large part of the memory footprint of CySolveOutput in the first place. To emphasize many problems will not see an improvement moving from a global CySolveOutput to a local one with multithreaded workflows (as we see in the example below)!

Examples where this approach may be faster (you should still test both approaches!):

• The ODE system has many dependent y variables so the internal structures of CySolveOutput start to take up a considerable memory footprint.

• You only care about the output at the end of the integration or at certain time steps: then the storage arrays would be much smaller than dealing with the storage arrays inside CySolveOutput.

The example below is just one example of an approach you could take.

%%cython --force
# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False

from cython.parallel cimport prange

from libc.string cimport memcpy

from libcpp cimport bool as cpp_bool
from libcpp.vector cimport vector
from libcpp.memory cimport unique_ptr, make_unique
from libcpp.utility cimport move

from CyRK.cy.cysolver_api cimport cysolve_ivp_noreturn, CySolverResult, ODEMethod
from CyRK.cy.cysolver_test cimport lotkavolterra_diffeq

cdef extern from "<chrono>" namespace "std::chrono" nogil:
    cdef cppclass microseconds:
        long long count()
    cdef cppclass duration "std::chrono::high_resolution_clock::duration":
        pass
    microseconds duration_cast_microseconds "std::chrono::duration_cast<std::chrono::microseconds>"(duration)

cdef extern from "<chrono>" namespace "std::chrono::high_resolution_clock" nogil:
    cdef cppclass time_point:
        duration operator-(time_point)
    time_point now()

cdef void single_thread_run(
        bint* success_storage_ptr,
        vector[double]* t_storage_vec_ptr,
        vector[double]* y_storage_vec_ptr,
        vector[double]* y0_vec_ptr,
        size_t num_dy,
        vector[char]& args_vec,
        double t_start,
        double t_end,
        size_t num_jobs
        ) noexcept nogil:

    # Thread-Local Allocation
    # This solver exists only for this thread and is reused for every job in the batch.
    # For this example we will make it easy and just hard code the integration method.
    cdef unique_ptr[CySolverResult] solver_uptr = make_unique[CySolverResult](ODEMethod.RK45)
    cdef CySolverResult* solver_ptr = solver_uptr.get()

    cdef size_t solution_size

    # Have this thread work on each job sequentially
    cdef size_t job_idx
    for job_idx in range(num_jobs):
        # All arrays will be provided starting at this particular thread's first job.
        # It does not need to know any details about what the other threads are doing.
        # So we start from 0 and end when we are done with our work.

        # Reuse the solver pointer
        cysolve_ivp_noreturn(
            solver_ptr,
            lotkavolterra_diffeq,  # Again this is hardcoded, as are many of the other arguments
            t_start,
            t_end,
            y0_vec_ptr[job_idx],
            1.0e-5, # rtol
            1.0e-8, # atol
            args_vec
        )

        # D. Extract Results to Shared Arrays
        # Store success status
        success_storage_ptr[job_idx] = solver_ptr.success

        if solver_ptr.success:
            solution_size = solver_ptr.time_domain_vec.size()
            # Copy the time_domain vector from the solution into our storage array.
            t_storage_vec_ptr[job_idx].resize(solution_size)
            memcpy(t_storage_vec_ptr[job_idx].data(), solver_ptr.time_domain_vec.data(), solution_size * sizeof(double))

            # Repeat for y's
            # There are multiple y's stored in a flattened solution array in the solver_ptr. But we don't actually
            #  care about their shape right now because we want to copy the exact shape over to our final array.
            # All we care about is the total size so we can use memcpy.
            # We use num_dy not num_y because the user might be capturing extra outputs
            y_storage_vec_ptr[job_idx].resize(solution_size * num_dy)
            memcpy(y_storage_vec_ptr[job_idx].data(), solver_ptr.solution.data(), solution_size * num_dy * sizeof(double))


def run_prange_batched(int num_threads=2):

    if num_threads < 1:
        num_threads = 1

    cdef int num_runs = 1000
    cdef double t_start = 0.0
    cdef double t_end = 2000.0

    # 1. Setup Inputs
    cdef vector[vector[double]] y0_vecs = vector[vector[double]]()
    y0_vecs.resize(num_runs)
    cdef size_t num_dy = 2
    cdef size_t num_y = 2
    cdef size_t i
    for i in range(num_runs):
        y0_vecs[i].resize(2)
        y0_vecs[i][0] = 10.0
        y0_vecs[i][1] = 5.0

    # Setup Constant Args
    cdef vector[char] constant_args = vector[char]()
    constant_args.resize(4 * 8)
    cdef double* arg_dbl_ptr = <double*>constant_args.data()
    arg_dbl_ptr[0] = 1.5
    arg_dbl_ptr[1] = 1.0
    arg_dbl_ptr[2] = 3.0
    arg_dbl_ptr[3] = 1.0

    # Setup Output storage for our full solution domain (Lightweight vectors instead of full CySolverResult objects)
    # We pre-allocate these so threads can write to them without locking.
    # Since each thread writes to a unique index 'job_idx', this is thread-safe.
    cdef vector[vector[double]] final_time_domains = vector[vector[double]](num_runs)
    cdef vector[vector[double]] final_ys           = vector[vector[double]](num_runs)

    # Important! We use cpp_bool a lot throughout CyRK but c++ bool is just a bit and can't have a memory address.
    # Since we need to take the pointer addresses of elements of this vector we need to use another struct. bint should work.
    cdef vector[bint] success_flags = vector[bint](num_runs)

    # Determine each thread's starting index and number of jobs
    cdef Py_ssize_t thread_idx
    cdef vector[size_t] thread_jobs        = vector[size_t](num_threads)
    cdef vector[size_t] thread_start_index = vector[size_t](num_threads)

    cdef size_t start_idx
    cdef size_t end_idx
    for thread_idx in range(num_threads):
        # Integer division chunking: (i*N)//T to ((i+1)*N)//T
        start_idx = (thread_idx * num_runs) // num_threads
        end_idx   = ((thread_idx + 1) * num_runs) // num_threads

        thread_start_index[thread_idx] = start_idx
        thread_jobs[thread_idx]        = end_idx - start_idx

    # Loop over the number of threads, NOT the number of runs.
    print(f"Starting Batched Run ({num_threads} threads; jobs per thread = {num_runs // num_threads})...")
    cdef time_point start_time, end_time
    start_time = now()
    for thread_idx in prange(num_threads, nogil=True, schedule='static'):
        single_thread_run(
            &success_flags[thread_start_index[thread_idx]],
            &final_time_domains[thread_start_index[thread_idx]],
            &final_ys[thread_start_index[thread_idx]],
            &y0_vecs[thread_start_index[thread_idx]],
            num_dy,
            constant_args,
            t_start,
            t_end,
            thread_jobs[thread_idx]
            )
    end_time = now()

    # 3. Validation
    for i in range(num_runs):
        assert success_flags[i] == True

    cdef double loop_time_ms = <double>(<int>duration_cast_microseconds(end_time - start_time).count()) / 1000.0
    print(f"Batched test finished in {loop_time_ms:0.1f} ms.")
    print(f"\tAvg time: {loop_time_ms/num_runs:0.3f}ms/loop.")

    return loop_time_ms

import os
cdef int cpu_count = os.cpu_count()
cdef double result
cdef int threads
cdef double t1_result, last_result
print("\n\nStarting Batched prange Demo!")
for threads in range(1, min(cpu_count, 8)):
    print(f"Working with N={threads}...")
    result = run_prange_batched(threads)
    if threads == 1:
        print(f"Total Time: {result:0.3f}ms.\n\n")
        t1_result = result
    else:
        print(f"Total Time: {result:0.3f}ms ({100.0*(result - last_result)/last_result:0.1f} % change from N={threads-1}; {100.0*(result - t1_result)/t1_result:0.1f} % change from N=1).\n\n")
    last_result = result

Content of stdout:
_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(25738): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(25739): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(25862): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(29004): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(29011): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(29511): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cpp(29512): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_161526d290fbaed164bb95774fbc88f57d73a0433fdc7e9e0c59e7bd722866c8.cp312-win_amd64.exp
Generating code
Finished generating code

Starting Batched prange Demo!
Working with N=1...
Starting Batched Run (1 threads; jobs per thread = 1000)...
Batched test finished in 1897.6 ms.
        Avg time: 1.898ms/loop.
Total Time: 1897.574ms.


Working with N=2...
Starting Batched Run (2 threads; jobs per thread = 500)...
Batched test finished in 979.7 ms.
        Avg time: 0.980ms/loop.
Total Time: 979.657ms (-48.4 % change from N=1; -48.4 % change from N=1).


Working with N=3...
Starting Batched Run (3 threads; jobs per thread = 333)...
Batched test finished in 666.8 ms.
        Avg time: 0.667ms/loop.
Total Time: 666.806ms (-31.9 % change from N=2; -64.9 % change from N=1).


Working with N=4...
Starting Batched Run (4 threads; jobs per thread = 250)...
Batched test finished in 505.2 ms.
        Avg time: 0.505ms/loop.
Total Time: 505.152ms (-24.2 % change from N=3; -73.4 % change from N=1).


Working with N=5...
Starting Batched Run (5 threads; jobs per thread = 200)...
Batched test finished in 414.2 ms.
        Avg time: 0.414ms/loop.
Total Time: 414.221ms (-18.0 % change from N=4; -78.2 % change from N=1).


Working with N=6...
Starting Batched Run (6 threads; jobs per thread = 166)...
Batched test finished in 357.3 ms.
        Avg time: 0.357ms/loop.
Total Time: 357.297ms (-13.7 % change from N=5; -81.2 % change from N=1).


Working with N=7...
Starting Batched Run (7 threads; jobs per thread = 142)...
Batched test finished in 326.5 ms.
        Avg time: 0.327ms/loop.
Total Time: 326.510ms (-8.6 % change from N=6; -82.8 % change from N=1).

Arbitrary Additional Arguments

Additional Args = Array of Doubles

See discussion in “Advanced CySolver.md”

%%cython --force
# 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 libcpp.utility cimport move
from libcpp.vector cimport vector

from CyRK cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType, MAX_STEP, CySolveOutput, PreEvalFunc, ODEMethod

cdef void pendulum_diffeq(double* dy_ptr, double t, double* y_ptr, char* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil:
    # Arguments are passed in as char pointers.
    # These point to preallocated and instantiated memory in the CySolver class that is filled with user-provided data.

    # The user can then recast the generic char pointers back to original format of args_ptr. In this case, an array of doubles.
    # Seg faults will occur if this function recasts to the incorrect type.
    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))

    # Unpack y
    cdef double y0, y1, torque
    y0 = y_ptr[0]
    y1 = y_ptr[1]

    # External torque
    torque = 0.1 * sin(t)

    dy_ptr[0] = y1
    dy_ptr[1] = coeff_1 * sin(y0) + coeff_2 * torque


# Now we define a function that builds our inputs and calls `cysolve_ivp` to solve the problem
def run():
    cdef DiffeqFuncType diffeq = pendulum_diffeq

    # Define time domain
    cdef double t_start = 0.0
    cdef double t_end = 10.0

    # Define initial conditions
    cdef vector[double] y0_vec = vector[double](2)
    y0_vec[0] = 0.01
    y0_vec[1] = 0.0

    # Define our arguments in a char vector
    # The size of the vector is the size of the underlying object(s) x number of them.
    cdef vector[char] args_vec = vector[char](sizeof(double) * 3)
    # To make it easier to populate the args_vec we recast it as a double pointer.
    cdef double* args_dbl_ptr = <double*>args_vec.data()
    args_dbl_ptr[0] = 1.0
    args_dbl_ptr[1] = 1.0
    args_dbl_ptr[2] = 9.81

    cdef CySolveOutput result = cysolve_ivp(
        diffeq,
        t_start,
        t_end,
        y0_vec,
        method = ODEMethod.RK45, # Integration method
        rtol = 1.0e-6,
        atol = 1.0e-8,
        args_vec = args_vec
        )

    # 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(move(result))

    return pysafe_result

result_reg = run()
print("\n\nIntegration success =", result_reg.success, "\n\tNumber of adaptive time steps required:", result_reg.size)
print("Integration message:", result_reg.message)

Content of stdout:
_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(25021): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(25022): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(25145): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(28009): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(28010): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(28630): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cpp(28637): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_31382ce44063f0dd07dd67fcf2e5fc782bfe4eca24cc5b848ceff7d9914b376a.cp312-win_amd64.exp
Generating code
Finished generating code

Integration success = True
        Number of adaptive time steps required: 152
Integration message: Integration completed without issue.

fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(result_reg.t, result_reg.y[0], c='C0')
ax2.plot(result_reg.t, result_reg.y[1], c='C3')
ax.set_ylabel("Angle [Rad]", c='C0')
ax.set_xlabel("Time [s]")
ax2.set_ylabel("Angular Velocity [Rad s-1]", c='C3')
plt.show()

Additional Args = Cython C-Struct

See discussion in “Advanced CySolver.md”

%%cython --force
# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
import numpy as np
cimport numpy as np
np.import_array()

from libc.string cimport memcpy
from libc.math cimport sin
from libcpp cimport bool as cpp_bool
from libcpp.utility cimport move
from libcpp.vector cimport vector

from CyRK cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType,MAX_STEP, CySolveOutput, PreEvalFunc, ODEMethod

cdef struct PendulumArgs:
    # Structure that contains heterogeneous types
    cpp_bool use_drag
    double g
    double l
    double m

cdef void pendulum_diffeq(double* dy_ptr, double t, double* y_ptr, char* args_ptr, PreEvalFunc pre_eval_func) noexcept nogil:
    # Arg pointer still must be listed as a char pointer or it will not work with cysolve_ivp.
    # But now the user can recast that char pointer to the structure they wish.
    cdef PendulumArgs* pendulum_args_ptr = <PendulumArgs*>args_ptr
    # And easily access its members which can be many heterogeneous types.
    cdef double l = pendulum_args_ptr.l
    cdef double m = pendulum_args_ptr.m
    cdef double g = pendulum_args_ptr.g
    cdef double use_drag = pendulum_args_ptr.use_drag

    cdef double coeff_1 = (-3. * g / (2. * l))
    cdef double coeff_2 = (3. / (m * l**2))

    # Unpack y
    cdef double y0, y1, torque
    y0 = y_ptr[0]
    y1 = y_ptr[1]

    # External torque
    torque = 0.1 * sin(t)

    dy_ptr[0] = y1
    dy_ptr[1] = coeff_1 * sin(y0) + coeff_2 * torque

    if use_drag:
        dy_ptr[1] -= 1.5 * y1

# Now we define a function that builds our inputs and calls `cysolve_ivp` to solve the problem
def run():
    cdef DiffeqFuncType diffeq = pendulum_diffeq

    # Define time domain
    cdef double t_start = 0.0
    cdef double t_end = 10.0

    # Define initial conditions
    cdef vector[double] y0_vec = vector[double](2)
    y0_vec[0] = 0.01
    y0_vec[1] = 0.0

    # Define our arguments.
    # We now have a a structure that we need to allocate memory for.
    # For this example, let's do it on the stack.
    cdef PendulumArgs pendulum_args = PendulumArgs(True, 9.81, 1.0, 1.0)
    # We need to pass in a char vector to cysolve_ivp, so let's create a vector of the correct size
    cdef vector[char] args_vec = vector[char](sizeof(PendulumArgs))
    # Then we need to copy over the contents to this vector.
    memcpy(args_vec.data(), &pendulum_args, sizeof(PendulumArgs))

    cdef CySolveOutput result = cysolve_ivp(
        diffeq,
        t_start,
        t_end,
        y0_vec,
        method = ODEMethod.RK45,
        rtol = 1.0e-6,
        atol = 1.0e-8,
        args_vec = args_vec
        )

    # 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(move(result))

    return pysafe_result

result_drag = run()
print("\n\nIntegration success =", result_drag.success, "\n\tNumber of adaptive time steps required:", result_drag.size)
print("Integration message:", result_drag.message)

Content of stdout:
_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(25074): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(25075): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(25198): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(28062): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(28063): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(28683): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cpp(28690): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_d4bc01dc0f3accbc4586830a5c49ad06d41619e57ccc498e020a9d049f406bdf.cp312-win_amd64.exp
Generating code
Finished generating code

Integration success = True
        Number of adaptive time steps required: 112
Integration message: Integration completed without issue.

fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(result_drag.t, result_drag.y[0], c='C0')
ax2.plot(result_drag.t, result_drag.y[1], c='C3')
ax.set_ylabel("Angle [Rad]", c='C0')
ax.set_xlabel("Time [s]")
ax2.set_ylabel("Angular Velocity [Rad s-1]", c='C3')
plt.show()

Pre-Evaluation Functions

See discussion in “Advanced CySolver.md”

%%cython --force
# 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 libcpp.utility cimport move
from libcpp.vector cimport vector

from CyRK cimport cysolve_ivp, WrapCySolverResult, DiffeqFuncType, MAX_STEP, CySolveOutput, PreEvalFunc, ODEMethod

cdef void pendulum_preeval_nodrag(char* output_ptr, double time, double* y_ptr, char* 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] = 0.0  # No Drag

cdef void pendulum_preeval_withdrag(char* output_ptr, double time, double* y_ptr, char* 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, char* 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 char so we can call function
    cdef char* pre_eval_storage_char_ptr = <char*>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_char_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 char* here.
    # Even though pre_eval_func was passed the char* 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 t_start = 0.0
    cdef double t_end   = 10.0

    # Define initial conditions
    cdef vector[double] y0_vec = vector[double](2)
    y0_vec[0] = 0.01
    y0_vec[1] = 0.0

    # Define our arguments.
    cdef vector[char] args_vec = vector[char](sizeof(double) * 3)
    cdef double* args_dbl_ptr = <double*>args_vec.data()
    args_dbl_ptr[0] = 1.0
    args_dbl_ptr[1] = 1.0
    args_dbl_ptr[2] = 9.81

    cdef CySolveOutput result = cysolve_ivp(
        diffeq,
        t_start,
        t_end,
        y0_vec,
        method = ODEMethod.RK45, # Integration method
        rtol = 1.0e-6,
        atol = 1.0e-8,
        args_vec = args_vec,
        num_extra = 0,
        max_num_steps = 100_000_000,
        max_ram_MB = 2000,
        dense_output = False,
        t_eval_vec = vector[double](), # C++ / Cython require you to provide args in order even if they unused like t_eval_vec here.
        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(move(result))

    return pysafe_result

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)

Content of stdout:
_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(25254): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(25255): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(25378): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(28242): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(28243): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(28863): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cpp(28870): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_9138999a69b1ea35f98f8614c03d533204edb47e27667122523e452dcfd94b6b.cp312-win_amd64.exp
Generating code
Finished generating code

Integration success = True
        Number of adaptive time steps required: 112
Integration message: Integration completed without issue.

fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(result_preeval.t, result_preeval.y[0], c='C0')
ax2.plot(result_preeval.t, result_preeval.y[1], c='C3')
ax.set_ylabel("Angle [Rad]", c='C0')
ax.set_xlabel("Time [s]")
ax2.set_ylabel("Angular Velocity [Rad s-1]", c='C3')
plt.show()

PySolver and CySolver Solution Reuses

See discussion in “CySolver & PySolver Reuses.md”

PySolver

# Note if using this format, `dy` must be the first argument. Additionally, a special flag must be set to True when calling pysolve_ivp, see below.
def cy_diffeq(dy, t, y):
    dy[0] = (1. - 0.01 * y[1]) * y[0]
    dy[1] = (0.02 * y[0] - 1.) * y[1]

# Since this is pure python we can njit it safely
from numba import njit
cy_diffeq = njit(cy_diffeq)

import numpy as np
from CyRK import pysolve_ivp

initial_conds = np.asarray((20., 20.), dtype=np.float64, order='C')
time_span = (0., 50.)
rtol = 1.0e-7
atol = 1.0e-8

result = \
    pysolve_ivp(cy_diffeq, time_span, initial_conds,
                method="RK45", rtol=rtol, atol=atol, pass_dy_as_arg=True)

# Use the result's memory allocation to re-call pysolve_ivp
result = \
    pysolve_ivp(cy_diffeq, time_span, initial_conds,
                method="RK45", rtol=rtol, atol=atol, pass_dy_as_arg=True,
                solution_reuse=result) # Set the reuse argument.

print("Was Integration was successful?", result.success)
print(result.message)
print("Size of solution: ", result.size)
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(result.t, result.y[0], c='C0')
ax.plot(result.t, result.y[1], c='C3')
plt.show()

# We can change parameters too
new_initial_conds = np.asarray((3., 0.1), dtype=np.float64, order='C')
result = \
    pysolve_ivp(cy_diffeq, time_span, new_initial_conds,
                method="RK45", rtol=rtol, atol=atol, pass_dy_as_arg=True,
                solution_reuse=result) # Set the reuse argument.

print("Was Integration was successful?", result.success)
print(result.message)
print("Size of solution: ", result.size)
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(result.t, result.y[0], c='C0')
ax.plot(result.t, result.y[1], c='C3')
plt.show()

Was Integration was successful? True
Integration completed without issue.
Size of solution:  360

Was Integration was successful? True
Integration completed without issue.
Size of solution:  601

CySolver

%%cython --force
# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False
from libcpp.utility cimport move
from libcpp.vector cimport vector

import numpy as np
cimport numpy as np
np.import_array()

from CyRK cimport cysolve_ivp, cysolve_ivp_noreturn, DiffeqFuncType, WrapCySolverResult, CySolveOutput, PreEvalFunc

cdef void cython_diffeq(double* dy, double t, double* y, char* args, PreEvalFunc pre_eval_func) noexcept nogil:

    # Build Coeffs
    cdef double coeff_1 = (1. - 0.01 * y[1])
    cdef double coeff_2 = (0.02 * y[0] - 1.)

    # Store results
    dy[0] = coeff_1 * y[0]
    dy[1] = coeff_2 * y[1]

# Import the required functions from CyRK
from CyRK cimport cysolve_ivp, DiffeqFuncType, WrapCySolverResult, CySolveOutput

# Let's get the integration number for the RK45 method
from CyRK cimport ODEMethod

def run_reuse_cysolver(tuple t_span, double[::1] y0):

    # Cast our diffeq to the accepted format
    cdef DiffeqFuncType diffeq = cython_diffeq

    # Convert the python user input to pure C types
    cdef size_t num_y          = len(y0)
    cdef double t_start        = t_span[0]
    cdef double t_end          = t_span[1]
    cdef vector[double] y0_vec = vector[double](num_y)
    cdef size_t yi
    for yi in range(num_y):
        y0_vec[yi] = y0[yi]

    # Run the integrator!
    cdef CySolveOutput result = cysolve_ivp(
        diffeq,
        t_start,
        t_end,
        y0_vec,
        method = ODEMethod.RK45,
        rtol = 1.0e-7,
        atol = 1.0e-8
    )

    # Run again using the noreturn version so we can reuse the integration storage.
    cysolve_ivp_noreturn(
        result.get(),  # result is a unique pointer; we need to pass the raw pointer to the baseline_cysolve_ivp
        diffeq,
        t_start,
        t_end,
        y0_vec,
        rtol = 1.0e-7,
        atol = 1.0e-8
    )

    # Like with PySolver, we can change the inputs
    cdef double t_end_2 = 2.0 * t_end
    cysolve_ivp_noreturn(
        result.get(),  # result is a unique pointer; we need to pass the raw pointer to the baseline_cysolve_ivp
        diffeq,
        t_start,
        t_end_2,
        y0_vec,
        rtol = 1.0e-7,
        atol = 1.0e-8
    )

    # The CySolveOutput is not accesible via Python. We need to wrap it first
    cdef WrapCySolverResult pysafe_result = WrapCySolverResult()
    pysafe_result.set_cyresult_pointer(move(result))

    return pysafe_result

Content of stdout:
_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(24996): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(24997): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(25120): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(28761): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(28768): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(29126): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cpp(29127): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_e8dc2852613e7617e71ec8816dba3ba288779843c465bd29263e43cbadf90870.cp312-win_amd64.exp
Generating code
Finished generating code

# Assume we are working in a Jupyter notebook so we don't need to import `run_cysolver` if it was defined in an earlier cell.
# from my_cython_code import run_cysolver

import numpy as np
initial_conds = np.asarray((20., 20.), dtype=np.float64, order='C')
time_span = (0., 50.)

result = run_reuse_cysolver(time_span, initial_conds)
result.print_diagnostics()

print("Was Integration was successful?", result.success)
print(result.message)
print("Size of solution: ", result.size)
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(result.t, result.y[0], c='C0')
ax.plot(result.t, result.y[1], c='C3')
plt.show()

----------------------------------------------------
CyRK v0.17.0 - WrapCySolverResult Diagnostic.
----------------------------------------------------
# of y:             2.
# of dy:            2.
# of events:        0.
Success:            True.
Error Code:         1.
Status:             Integration completed without issue..
Size:               705.
Steps Taken:        704.
Event Termination:  False.
Integrator Message:
        Integration completed without issue.

----------------- CySolverResult -------------------
Capture Extra:          False.
Capture Dense Output:   False.
Integration Direction:  Forward.
Integration Method:     Explicit Runge-Kutta method of order 5(4).
# of Interpolates:      0.

---- Additional Argument Info ----
args size (bytes):   0.
args size (doubles): 0.
Args Pointer is Null.
End of Additional Argument Info.

------------------ CySolverBase --------------------
Integration Method: 4.
# of y:             2.
# of dy:            2.
PySolver:           False.
---- Current State Info ----
t_now: 100.0.
y_now:
        y0  = 1.40214e+02.
        y1  = 3.76947e+01.
dy_now:
        dy0 = 8.73608e+01.
        dy1 = 6.80120e+01.
End of Current State Info.

-------------- Diagnostic Complete -----------------

Was Integration was successful? True
Integration completed without issue.
Size of solution:  705

Using Events with CySolver

%%cython --force
# cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True, initializedcheck=False

from libcpp cimport bool as cpp_bool
from libcpp.vector cimport vector
from libcpp.limits cimport numeric_limits

import numpy as np
cimport numpy as cnp
cnp.import_array()

from CyRK cimport CyrkErrorCodes, PreEvalFunc, MAX_STEP, DiffeqFuncType, Event, EventFunc, cysolve_ivp_noreturn, WrapCySolverResult, ODEMethod, CySolverResult
from CyRK.cy.cysolver_test cimport lorenz_diffeq

cdef double lorenz_event_func_1(double t, double* y_ptr, char* unused_args) noexcept nogil:

    # Check if t greater than or equal to 5.0
    if t >= 5.0:
        return 0.0
    else:
        return 1.0

cdef double lorenz_event_func_2(double t, double* y_ptr, char* unused_args) noexcept nogil:

    # Check y values.
    # In the time span [0,10]:
    #    y_0  starts at 1, spikes then goes below zero and oscillates with a min below -10. Have this return if y_0 < -10
    if y_ptr[0] < -10.0:
        return 0.0
    elif y_ptr[2] > 30.0:
        return 0.0
    else:
        return 1.0

cdef double lorenz_event_func_3(double t, double* y_ptr, char* args) noexcept nogil:

    # Use arguments to set threshold
    cdef double* args_dbl_ptr = <double*>args

    # We won't actually use the args but lets just check they are correct.
    cdef cpp_bool args_correct = False
    if args_dbl_ptr[0] == 10.0 and args_dbl_ptr[1] == 28.0 and args_dbl_ptr[2] == 8.0 / 3.0:
        args_correct = True

    # Then return events if args are correct and t greater than 4
    if args_correct and t >= 4.0:
        if t <= 4.1:
            return 0.0
        elif t >= 4.5 and t <= 4.6:
            return 0.0
        elif t >= 5.0 and t <= 5.1:
            return 0.0
        elif t >= 5.5 and t <= 5.6:
            return 0.0
        elif t >= 6.0 and t <= 6.1:
            return 0.0
        else:
            return 1.0
    else:
        return 1.0

def run_cysolver_with_events(bint use_termination = False):

    # Inputs for lorenz diffeq
    cdef double t_start = 0.0
    cdef double t_end = 10.0

    cdef size_t num_y = 3
    cdef vector[double] y0_vec = vector[double](num_y)
    y0_vec[0] = 1.0
    y0_vec[1] = 0.0
    y0_vec[2] = 0.0

    cdef vector[double] t_eval_vec = vector[double]()
    cdef vector[char] args_vec = vector[char](3 * sizeof(double))
    args_dbl_ptr = <double*>args_vec.data()
    args_dbl_ptr[0] = 10.0
    args_dbl_ptr[1] = 28.0
    args_dbl_ptr[2] = 8.0 / 3.0

    cdef PreEvalFunc pre_eval_func = NULL

    # Build events
    cdef vector[Event] events_vec = vector[Event]()
    # Set the maximum number of events allowed before termination to the max size of size_t (effectively infinite)
    cdef size_t max_allowed = numeric_limits[size_t].max()
    cdef int direction = 0
    if use_termination:
        events_vec.emplace_back(lorenz_event_func_1, 1, 0)
    else:
        events_vec.emplace_back(lorenz_event_func_1, max_allowed, direction)
    events_vec.emplace_back(lorenz_event_func_2, max_allowed, direction)
    events_vec.emplace_back(lorenz_event_func_3, max_allowed, direction)

    cdef ODEMethod integration_method = ODEMethod.RK45
    cdef WrapCySolverResult solution = WrapCySolverResult()
    solution.build_cyresult(integration_method)
    cdef CySolverResult* solution_ptr = solution.cyresult_uptr.get()

    cdef DiffeqFuncType diffeq = lorenz_diffeq
    cdef size_t num_extra = 0
    cdef cpp_bool use_dense = False

     # Run Solver
    cysolve_ivp_noreturn(
        solution_ptr,
        diffeq,
        t_start,
        t_end,
        y0_vec,
        rtol = 1.0e-6,
        atol = 1.0e-6,
        args_vec = args_vec,
        num_extra = num_extra,
        max_num_steps = 100000,
        max_ram_MB = 2000,
        dense_output = use_dense,
        t_eval_vec = t_eval_vec,
        pre_eval_func = pre_eval_func,
        events_vec = events_vec,
        rtols_vec = vector[double](),
        atols_vec = vector[double](),
        max_step = MAX_STEP,
        first_step = 0.0,
        expected_size = 128
        )

    solution.finalize()
    return solution

Content of stdout:
_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(25643): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(25644): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(25767): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(28631): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(28632): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(29252): warning C4551: function call missing argument list
C:\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cpp(29259): warning C4551: function call missing argument list
   Creating library C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cp312-win_amd64.lib and object C:\Users\joepr\.ipython\cython\Users\joepr\.ipython\cython\_cython_magic_be0781bd329588d8007d63393b0fe6608ce0b7554b2f6b534e40a773d1238870.cp312-win_amd64.exp
Generating code
Finished generating code

# Compare to the solve_ivp and pysolve_ivp examples in "1 - Getting Started" Notebook

# CyRK v0.16.0
#    Terminate off: 61.1us; 62.0us; 62.6us
#    Terminate on:  34.4us; 33.7us; 34.4us
# CyRK v0.17.0
#    Terminate off: 57.9us; 56.9us; 58.5us
#    Terminate on:  29.4us; 30.6us; 29.6us
terminate = True

solution = run_cysolver_with_events(terminate)
print("Solution Succeeded:", solution.success)

%timeit run_cysolver_with_events(terminate)

fig, ax = plt.subplots()
ax.plot(solution.t, solution.y[0], label='y0', c='C0')
ax.plot(solution.t, solution.y[1], label='y1', c='C3')
ax.plot(solution.t, solution.y[2], label='y2', c='C6')
# Event 1 will be dots.
# Event 2 will be X's
# Event 3 will be +'s
ax.scatter(solution.t_events[1], solution.y_events[1][0, :], c='C0', marker='o')
ax.scatter(solution.t_events[0], solution.y_events[0][1, :], c='C3', marker='x', s=120)
ax.scatter(solution.t_events[1], solution.y_events[1][2, :], c='C6', marker='o')
ax.scatter(solution.t_events[2], solution.y_events[2][0, :], c='C0', marker='+', s=120)
ax.set(xlabel='Time', ylabel='y')
ax.legend(loc="upper right")
plt.show()

Solution Succeeded: True
27.3 μs ± 376 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)