Events¶
Integration “events” are points in time or state where the user may want something recorded or have the entire integration to terminate early.
CyRK largely follows the implementation used by SciPy’s solve_ivp.
Event Functions¶
pysolve_ivp Version¶
An event is tracked using a user-provided function. The function must take the form of event(t, y(t), *args). “args”
are the same additional arguments provided to the differential equation function. The event then must return a float. Following SciPy’s
approach: “[CyRK] will find an accurate value of \(t\) at which event(t, y(t), *args) = 0 using a root-finding
algorithm. By default, all zeros will be found. The solver looks for a sign change over each step, so if multiple
zero crossings occur within one step, events may be missed.
Additionally each event function might have the following attributes:”
`event.terminal`: bool or int, optional
If provided as a bool = True, then it will tell the integrator to stop integration if this
event is ever 0.
If an int > 1, then termination occurs after this number of occurrences of this event.
If not defined then integrator will not terminate if this event is hit (implicitly False).
`event.direction`: float, optional
Direction of a zero crossing.
If direction is positive, event will only trigger when going from negative to positive or
from 0 to 0, and vice versa if direction is negative.
If 0, then either direction will trigger event.
Implicitly 0 if not assigned.
The above is largely copied from SciPy’s documentation here.
Example¶
import numpy as np
from CyRK import pysolve_ivp
import matplotlib.pyplot as plt
from numba import njit
# Example with SciPy
@njit
def event_func_1(t, y, a, b, c):
# Check if t greater than or equal to 5.0
if t >= 5.0:
return 0.0
else:
return 1.0
@njit
def event_func_2(t, y, a, b, c):
# 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[0] < -10.0:
return 0.0
elif y[2] > 30.0:
return 0.0
else:
return 1.0
@njit
def event_func_3(t, y, a, b, c):
# We won't actually use the args but lets just check they are correct.
args_correct = False
if a == 10.0 and b == 28.0 and c == 8.0 / 3.0:
args_correct = True
# Then return events if args are correct and t greater than 8
if args_correct and t > 8.0:
return 0.0
else:
return 1.0
@njit
def lorenz_diffeq(dy, t, y, a, b, c):
# Unpack y
y0 = y[0]
y1 = y[1]
y2 = y[2]
dy[0] = a * (y1 - y0)
dy[1] = y0 * (b - y2) - y1
dy[2] = y0 * y1 - c * y2
def run_pysove_with_events(terminate = False):
event_func_1.direction = 0
if terminate:
event_func_1.terminal = 1
event_func_2.direction = 0
event_func_3.direction = 0
time_span = (0.0, 10.0)
y0 = np.asarray([1.0, 0.0, 0.0])
args = (10.0, 28.0, 8.0/3.0)
solution = pysolve_ivp(lorenz_diffeq, time_span, y0, method="RK45",
args=args, rtol=1.0e-6, atol=1.0e-6, t_eval=None, dense_output=False, pass_dy_as_arg=True,
events=(event_func_1, event_func_2, event_func_3))
return solution
cysolve_ivp Implementation¶
Please see the Advanced CySolver demo notebook to see this in action.
%%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
To plot the data we can use a python cell.
# Compare to the solve_ivp and pysolve_ivp examples in "1 - Getting Started" Notebook
# 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()
Differences from SciPy’s Implementation¶
CyRK is capable of capturing extra parameters during integration (see Extra Output for more details).
These will also be captured during event triggers and can be used in an event to determine if an event should trigger.
They are appended to the y array passed to the event function after the dependent variables.
CyRK outputs y values as the first index of y_events compared to SciPy outputting them as the second.
scipy_solution.y_events == cyrk_solution.y_events.T
Additional Arguments to Event Functions¶
pysolve_ivp Example:¶
def event_func(t, y, *args):
y0 = y[0]
y1 = y[1]
y2 = y[2]
# Only 3 dependent ys in this example.
extra0 = y[3]
extra1 = y[4]
extra2 = y[5]
# These extras can now be used to trigger an event!
cysolve_ivp Example:¶
cdef event_func(double t, double* y, char* args):
cdef double y0 = y[0]
cdef double y1 = y[1]
cdef double y2 = y[2]
# Only 3 dependent ys in this example.
cdef double extra0 = y[3]
cdef double extra1 = y[4]
cdef double extra2 = y[5]
# These extras can now be used to trigger an event!
Event Diagnostics¶
CyRK provides several event-related data that can be used to diagnose issues or just learn about the integration results.
solution = pysolve_ivp(... events=(my_event_fuc1, my_event_func2), ...)
# Number of event functions provided, in this case 2.
solution.num_events
# Boolean flag if the integration was terminated by an event.
solution.event_terminated
# If the integration was terminated, then the index of the function the first triggered the termination is provided by the following:
solution.event_terminate_index
# All of this information (plus a lot more!) will be printed if you print the solution diagnostics
solution.print_diagnostics()
Citation¶
Use of CyRK’s integrators with “events” on should cite SciPy since this implementation is based on their work.