From ecba8eb1fcd36e3f489016e21de2d6e2fdd77b2f Mon Sep 17 00:00:00 2001 From: "Corey R. Randall" Date: Sun, 22 Mar 2026 21:06:44 -0600 Subject: [PATCH 1/4] Catch and propagate python exceptions without storing in AuxData --- src/sksundae/_cy_common.pxd | 6 +++ src/sksundae/_cy_common.pyx | 37 ++++++++++++++++ src/sksundae/_cy_cvode.pyx | 78 ++++++++++++---------------------- src/sksundae/_cy_ida.pyx | 84 +++++++++++++------------------------ src/sksundae/c_cvode.pxd | 14 +++---- src/sksundae/c_ida.pxd | 14 +++---- 6 files changed, 112 insertions(+), 121 deletions(-) diff --git a/src/sksundae/_cy_common.pxd b/src/sksundae/_cy_common.pxd index 133de42..f2dff01 100644 --- a/src/sksundae/_cy_common.pxd +++ b/src/sksundae/_cy_common.pxd @@ -6,6 +6,12 @@ cimport numpy as np # Extern cdef headers from .c_sundials cimport * # Access to C types +# Propagate python exceptions or print SUNDIALS error messages +cdef _pyerr_handler() +cdef void _sunerr_handler( + int line, const char* func, const char* file, const char* msg, int err_code, + void* err_user_data, SUNContext ctx) except * + # Convert between N_Vector and numpy array cdef svec2np(N_Vector nvec, np.ndarray[DTYPE_t, ndim=1] np_array) cdef np2svec(np.ndarray[DTYPE_t, ndim=1] np_array, N_Vector nvec) diff --git a/src/sksundae/_cy_common.pyx b/src/sksundae/_cy_common.pyx index 23252c8..7836a0c 100644 --- a/src/sksundae/_cy_common.pyx +++ b/src/sksundae/_cy_common.pyx @@ -3,6 +3,14 @@ # Dependencies cimport numpy as np +from cpython.exc cimport ( + PyObject, PyErr_Occurred, PyErr_Fetch, PyErr_Restore, + # PyErr_GetRaisedException, PyErr_SetRaisedException, +) + +# PyErr_Fetch and PyErr_Restore are deprecated at 3.12. When support for <3.12 +# is dropped, replace with PyErr_GetRaisedException/PyErr_SetRaisedException. + # Extern cdef headers from .c_sundials cimport * # Access to C types from .c_nvector cimport * # Access to N_Vector functions @@ -37,6 +45,35 @@ elif SUNDIALS_INT_TYPE == "long int": from numpy import int64 as INT_TYPE +cdef _pyerr_handler(): + """Catch and re-raise Python exceptions in Cython code.""" + cdef PyObject *errtype, *errvalue, *errtraceback + + PyErr_Fetch(&errtype, &errvalue, &errtraceback) + PyErr_Restore(errtype, errvalue, errtraceback) + + # Update to the following when support for Python <3.12 is dropped: + # cdef PyObject *exc = PyErr_GetRaisedException() + # PyErr_SetRaisedException(exc) + # raise exc + + raise errvalue + + +cdef void _sunerr_handler(int line, const char* func, const char* file, + const char* msg, int err_code, void* err_user_data, + SUNContext ctx) except *: + """Custom error handler for shorter messages (no line or file).""" + + if PyErr_Occurred(): + pass + + else: + decoded_func = func.decode("utf-8") + decoded_msg = msg.decode("utf-8").replace(", ,", ",").strip() + print(f"\n[{decoded_func}, Error: {err_code}] {decoded_msg}\n") + + cdef svec2np(N_Vector nvec, np.ndarray[DTYPE_t, ndim=1] np_array): """Fill a numpy array with values from an N_Vector.""" cdef sunrealtype* nv_ptr diff --git a/src/sksundae/_cy_cvode.pyx b/src/sksundae/_cy_cvode.pyx index 8cded6a..f78ef87 100644 --- a/src/sksundae/_cy_cvode.pyx +++ b/src/sksundae/_cy_cvode.pyx @@ -17,13 +17,7 @@ cimport numpy as np from scipy import sparse as sp from scipy.optimize._numdiff import group_columns -from cpython.exc cimport ( - PyErr_Fetch, PyErr_NormalizeException, - PyObject, PyErr_CheckSignals, PyErr_Occurred, # PyErr_GetRaisedException, -) - -# PyErr_Fetch and PyErr_NormalizeException are deprecated at 3.12. When support -# for <3.12 is dropped, replace with PyErr_GetRaisedException. +from cpython.exc cimport PyErr_Occurred # Extern cdef headers from .c_cvode cimport * @@ -96,7 +90,7 @@ LSMESSAGES = { cdef int _rhsfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, - void* data) except? -1: + void* data) except -1: """Wraps 'rhsfn' by converting between N_Vector and ndarray types.""" aux = data @@ -114,7 +108,7 @@ cdef int _rhsfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, cdef int _eventsfn_wrapper(sunrealtype t, N_Vector yy, sunrealtype* ee, - void* data) except? -1: + void* data) except -1: """Wraps 'eventsfn' by converting between N_Vector and ndarray types.""" aux = data @@ -133,7 +127,7 @@ cdef int _eventsfn_wrapper(sunrealtype t, N_Vector yy, sunrealtype* ee, cdef int _jacfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, SUNMatrix JJ, void* data, N_Vector tmp1, N_Vector tmp2, - N_Vector tmp3) except? -1: + N_Vector tmp3) except -1: """Wraps 'jacfn' by converting between N_Vector and ndarray types.""" aux = data @@ -153,7 +147,7 @@ cdef int _jacfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, SUNMatrix JJ, cdef int _psetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, sunbooleantype jok, sunbooleantype* jcurPtr, - sunrealtype gamma, void* data) except? -1: + sunrealtype gamma, void* data) except -1: """Wraps 'psetup' by converting between N_Vector and ndarray types.""" aux = data @@ -176,7 +170,7 @@ cdef int _psetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, cdef int _psolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rv, N_Vector zv, sunrealtype gamma, sunrealtype delta, - int lr, void* data) except? -1: + int lr, void* data) except -1: """Wraps 'psolve' by converting between N_Vector and ndarray types.""" aux = data @@ -199,7 +193,7 @@ cdef int _psolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rv, cdef int _jvsetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, - void* data) except? -1: + void* data) except -1: """Wraps 'jvsetup' by converting between N_Vector and ndarray types.""" aux = data @@ -217,7 +211,7 @@ cdef int _jvsetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, cdef int _jvsolve_wrapper(N_Vector vv, N_Vector Jv, sunrealtype t, N_Vector yy, - N_Vector yp, void* data, N_Vector tmp) except? -1: + N_Vector yp, void* data, N_Vector tmp) except -1: """Wraps 'jvsolve' by converting between N_Vector and ndarray types.""" aux = data @@ -237,27 +231,6 @@ cdef int _jvsolve_wrapper(N_Vector vv, N_Vector Jv, sunrealtype t, N_Vector yy, return 0 -cdef void _err_handler(int line, const char* func, const char* file, - const char* msg, int err_code, void* err_user_data, - SUNContext ctx) except *: - """Custom error handler for shorter messages (no line or file).""" - cdef PyObject *errtype, *errvalue, *errtraceback - - if PyErr_Occurred(): - aux = err_user_data - # aux.pyerr = PyErr_GetRaisedException() - - PyErr_Fetch(&errtype, &errvalue, &errtraceback) - PyErr_NormalizeException(&errtype, &errvalue, &errtraceback) - - aux.pyerr = errvalue - - else: - decoded_func = func.decode("utf-8") - decoded_msg = msg.decode("utf-8").replace(", ,", ",").strip() - print(f"\n[{decoded_func}, Error: {err_code}] {decoded_msg}\n") - - cdef class AuxData: """ Auxiliary data. @@ -278,7 +251,6 @@ cdef class AuxData: cdef bint with_userdata cdef bint is_constrained - cdef object pyerr # Exception cdef object rhsfn # Callable cdef object userdata # Any cdef object eventsfn # Callable @@ -289,7 +261,6 @@ cdef class AuxData: cdef object jactimes # CVODEJacTimes def __cinit__(self, sunindextype NEQ, object options): - self.pyerr = None self.np_yy = np.empty(NEQ, DTYPE) self.np_yp = np.empty(NEQ, DTYPE) @@ -454,6 +425,9 @@ cdef class _cvLSSparseDQJac: """ self.mem = mem + def __dealloc__(self): + self.mem = NULL + class CVODEResult(RichResult): _order_keys = ["message", "success", "status", "t", "y", "i_events", @@ -756,7 +730,7 @@ cdef class CVODE: # 16) Set optional inputs SUNContext_ClearErrHandlers(self.ctx) - SUNContext_PushErrHandler(self.ctx, _err_handler, self.aux) + SUNContext_PushErrHandler(self.ctx, _sunerr_handler, NULL) cdef sunrealtype first_step = self._options["first_step"] flag = CVodeSetInitStep(self.mem, first_step) @@ -884,9 +858,11 @@ cdef class CVODE: return result - cdef _normal_solve(self, np.ndarray[DTYPE_t, ndim=1] tspan, - np.ndarray[DTYPE_t, ndim=1] y0, - ): + cdef _normal_solve( + self, + np.ndarray[DTYPE_t, ndim=1] tspan, + np.ndarray[DTYPE_t, ndim=1] y0, + ): cdef int ind cdef int flag @@ -937,10 +913,8 @@ cdef class CVODE: ind += 1 - if self.aux.pyerr is not None: - raise self.aux.pyerr - elif PyErr_CheckSignals() == -1: - return + if PyErr_Occurred(): + _pyerr_handler() elif stop: break @@ -963,9 +937,11 @@ cdef class CVODE: return result - cdef _onestep_solve(self, np.ndarray[DTYPE_t, ndim=1] tspan, - np.ndarray[DTYPE_t, ndim=1] y0, - ): + cdef _onestep_solve( + self, + np.ndarray[DTYPE_t, ndim=1] tspan, + np.ndarray[DTYPE_t, ndim=1] y0, + ): cdef int ind cdef int flag @@ -1022,10 +998,8 @@ cdef class CVODE: ind += 1 - if self.aux.pyerr is not None: - raise self.aux.pyerr - elif PyErr_CheckSignals() == -1: - return + if PyErr_Occurred(): + _pyerr_handler() elif stop: break diff --git a/src/sksundae/_cy_ida.pyx b/src/sksundae/_cy_ida.pyx index bbcec58..458572b 100644 --- a/src/sksundae/_cy_ida.pyx +++ b/src/sksundae/_cy_ida.pyx @@ -17,13 +17,7 @@ cimport numpy as np from scipy import sparse as sp from scipy.optimize._numdiff import group_columns -from cpython.exc cimport ( - PyErr_Fetch, PyErr_NormalizeException, - PyObject, PyErr_CheckSignals, PyErr_Occurred, # PyErr_GetRaisedException, -) - -# PyErr_Fetch and PyErr_NormalizeException are deprecated at 3.12. When support -# for <3.12 is dropped, replace with PyErr_GetRaisedException. +from cpython.exc cimport PyErr_Occurred # Extern cdef headers from .c_ida cimport * @@ -94,7 +88,7 @@ LSMESSAGES = { cdef int _resfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, - void* data) except? -1: + void* data) except -1: """Wraps 'resfn' by converting between N_Vector and ndarray types.""" aux = data @@ -113,7 +107,7 @@ cdef int _resfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, cdef int _eventsfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, - sunrealtype* ee, void* data) except? -1: + sunrealtype* ee, void* data) except -1: """Wraps 'eventsfn' by converting between N_Vector and ndarray types.""" aux = data @@ -133,7 +127,7 @@ cdef int _eventsfn_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, cdef int _jacfn_wrapper(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector rr, SUNMatrix JJ, void* data, N_Vector tmp1, - N_Vector tmp2, N_Vector tmp3) except? -1: + N_Vector tmp2, N_Vector tmp3) except -1: """Wraps 'jacfn' by converting between N_Vector and ndarray types.""" aux = data @@ -154,7 +148,7 @@ cdef int _jacfn_wrapper(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, cdef int _psetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, - sunrealtype cj, void* data) except? -1: + sunrealtype cj, void* data) except -1: """Wraps 'psetup' by converting between N_Vector and ndarray types.""" aux = data @@ -174,7 +168,7 @@ cdef int _psetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, cdef int _psolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector rv, N_Vector zv, sunrealtype cj, - sunrealtype delta, void* data) except? -1: + sunrealtype delta, void* data) except -1: """Wraps 'psolve' by converting between N_Vector and ndarray types.""" aux = data @@ -198,7 +192,7 @@ cdef int _psolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, cdef int _jvsetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, - sunrealtype cj, void* data) except? -1: + sunrealtype cj, void* data) except -1: """Wraps 'jvsolve' by converting between N_Vector and ndarray types.""" aux = data @@ -218,7 +212,7 @@ cdef int _jvsetup_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, cdef int _jvsolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector vv, N_Vector Jv, sunrealtype cj, void* data, - N_Vector tmp1, N_Vector tmp2) except? -1: + N_Vector tmp1, N_Vector tmp2) except -1: """Wraps 'jvsolve' by converting between N_Vector and ndarray types.""" aux = data @@ -241,27 +235,6 @@ cdef int _jvsolve_wrapper(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, return 0 -cdef void _err_handler(int line, const char* func, const char* file, - const char* msg, int err_code, void* err_user_data, - SUNContext ctx) except *: - """Custom error handler for shorter messages (no line or file).""" - cdef PyObject *errtype, *errvalue, *errtraceback - - if PyErr_Occurred(): - aux = err_user_data - # aux.pyerr = PyErr_GetRaisedException() - - PyErr_Fetch(&errtype, &errvalue, &errtraceback) - PyErr_NormalizeException(&errtype, &errvalue, &errtraceback) - - aux.pyerr = errvalue - - else: - decoded_func = func.decode("utf-8") - decoded_msg = msg.decode("utf-8").replace(", ,", ",").strip() - print(f"\n[{decoded_func}, Error: {err_code}] {decoded_msg}\n") - - cdef class AuxData: """ Auxiliary data. @@ -283,7 +256,6 @@ cdef class AuxData: cdef bint with_userdata cdef bint is_constrained - cdef object pyerr # Exception cdef object resfn # Callable cdef object userdata # Any cdef object eventsfn # Callable @@ -294,7 +266,6 @@ cdef class AuxData: cdef object jactimes # IDAJacTimes def __cinit__(self, sunindextype NEQ, object options): - self.pyerr = None self.np_yy = np.empty(NEQ, DTYPE) self.np_yp = np.empty(NEQ, DTYPE) self.np_rr = np.empty(NEQ, DTYPE) @@ -468,6 +439,9 @@ cdef class _idaLSSparseDQJac: """ self.mem = mem + + def __dealloc__(self): + self.mem = NULL class IDAResult(RichResult): @@ -780,9 +754,9 @@ cdef class IDA: # 15) Set optional inputs SUNContext_ClearErrHandlers(self.ctx) - SUNContext_PushErrHandler(self.ctx, _err_handler, self.aux) + SUNContext_PushErrHandler(self.ctx, _sunerr_handler, NULL) - # Set algebraic variable indices + # Set algebraic variable indices np_algidx = np.ones(self.NEQ, DTYPE) if self._options["algebraic_idx"] is not None: @@ -951,10 +925,12 @@ cdef class IDA: return result - cdef _normal_solve(self, np.ndarray[DTYPE_t, ndim=1] tspan, - np.ndarray[DTYPE_t, ndim=1] y0, - np.ndarray[DTYPE_t, ndim=1] yp0 - ): + cdef _normal_solve( + self, + np.ndarray[DTYPE_t, ndim=1] tspan, + np.ndarray[DTYPE_t, ndim=1] y0, + np.ndarray[DTYPE_t, ndim=1] yp0, + ): cdef int ind cdef int flag @@ -1010,10 +986,8 @@ cdef class IDA: ind += 1 - if self.aux.pyerr is not None: - raise self.aux.pyerr - elif PyErr_CheckSignals() == -1: - return + if PyErr_Occurred(): + _pyerr_handler() elif stop: break @@ -1037,10 +1011,12 @@ cdef class IDA: return result - cdef _onestep_solve(self, np.ndarray[DTYPE_t, ndim=1] tspan, - np.ndarray[DTYPE_t, ndim=1] y0, - np.ndarray[DTYPE_t, ndim=1] yp0 - ): + cdef _onestep_solve( + self, + np.ndarray[DTYPE_t, ndim=1] tspan, + np.ndarray[DTYPE_t, ndim=1] y0, + np.ndarray[DTYPE_t, ndim=1] yp0, + ): cdef int ind cdef int flag @@ -1103,10 +1079,8 @@ cdef class IDA: ind += 1 - if self.aux.pyerr is not None: - raise self.aux.pyerr - elif PyErr_CheckSignals() == -1: - return + if PyErr_Occurred(): + _pyerr_handler() elif stop: break diff --git a/src/sksundae/c_cvode.pxd b/src/sksundae/c_cvode.pxd index 23782dd..56039be 100644 --- a/src/sksundae/c_cvode.pxd +++ b/src/sksundae/c_cvode.pxd @@ -6,8 +6,8 @@ from .c_sundials cimport * # Access to types cdef extern from "cvode/cvode.h": # user-supplied functions - ctypedef int (*CVRhsFn)(sunrealtype t, N_Vector yy, N_Vector yp, void* data) except? -1 - ctypedef int (*CVRootFn)(sunrealtype t, N_Vector yy, sunrealtype* ee, void* data) except? -1 + ctypedef int (*CVRhsFn)(sunrealtype t, N_Vector yy, N_Vector yp, void* data) except -1 + ctypedef int (*CVRootFn)(sunrealtype t, N_Vector yy, sunrealtype* ee, void* data) except -1 # imethod int CV_ADAMS @@ -69,22 +69,22 @@ cdef extern from "cvode/cvode_ls.h": # user-supplied functions ctypedef int (*CVLsJacFn)( sunrealtype tt, N_Vector yy, N_Vector yp, SUNMatrix JJ, void* data, - N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1 + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except -1 ctypedef int (*CVLsPrecSetupFn)( sunrealtype tt, N_Vector yy, N_Vector yp, sunbooleantype jok, sunbooleantype* jcurPtr, - sunrealtype gamma, void* data) except? -1 + sunrealtype gamma, void* data) except -1 ctypedef int (*CVLsPrecSolveFn)( sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rv, N_Vector zv, sunrealtype gamma, - sunrealtype delta, int lr, void* data) except? -1 + sunrealtype delta, int lr, void* data) except -1 ctypedef int (*CVLsJacTimesSetupFn)( - sunrealtype tt, N_Vector yy, N_Vector yp, void* data) except? -1 + sunrealtype tt, N_Vector yy, N_Vector yp, void* data) except -1 ctypedef int (*CVLsJacTimesVecFn)( N_Vector vv, N_Vector Jv, sunrealtype tt, N_Vector yy, N_Vector yp, void* data, - N_Vector tmp) except? -1 + N_Vector tmp) except -1 # exported functions int CVodeSetLinearSolver(void* mem, SUNLinearSolver LS, SUNMatrix A) diff --git a/src/sksundae/c_ida.pxd b/src/sksundae/c_ida.pxd index b80e59e..83276cd 100644 --- a/src/sksundae/c_ida.pxd +++ b/src/sksundae/c_ida.pxd @@ -6,8 +6,8 @@ from .c_sundials cimport * cdef extern from "ida/ida.h": # user-supplied functions - ctypedef int (*IDAResFn)(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void* data) except? -1 - ctypedef int (*IDARootFn)(sunrealtype t, N_Vector yy, N_Vector yp, sunrealtype* ee, void* data) except? -1 + ctypedef int (*IDAResFn)(sunrealtype t, N_Vector yy, N_Vector yp, N_Vector rr, void* data) except -1 + ctypedef int (*IDARootFn)(sunrealtype t, N_Vector yy, N_Vector yp, sunrealtype* ee, void* data) except -1 # itask int IDA_NORMAL @@ -75,23 +75,23 @@ cdef extern from "ida/ida_ls.h": # user-supplied functions ctypedef int (*IDALsJacFn)( sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector rr, - SUNMatrix JJ, void* data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except? -1 + SUNMatrix JJ, void* data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) except -1 ctypedef int (*IDALsPrecSetupFn)( sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, sunrealtype cj, - void* data) except? -1 + void* data) except -1 ctypedef int (*IDALsPrecSolveFn)( sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector rv, - N_Vector zv, sunrealtype cj, sunrealtype delta, void* data) except? -1 + N_Vector zv, sunrealtype cj, sunrealtype delta, void* data) except -1 ctypedef int (*IDALsJacTimesSetupFn)( sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, sunrealtype cj, - void* data) except? -1 + void* data) except -1 ctypedef int (*IDALsJacTimesVecFn)( sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector vv, N_Vector Jv, - sunrealtype cj, void* data, N_Vector tmp1, N_Vector tmp2) except? -1 + sunrealtype cj, void* data, N_Vector tmp1, N_Vector tmp2) except -1 # exported functions int IDASetLinearSolver(void* mem, SUNLinearSolver LS, SUNMatrix A) From dc4bc669c3ea71a5e111bcb1b113b8072e01b7b8 Mon Sep 17 00:00:00 2001 From: "Corey R. Randall" Date: Sun, 22 Mar 2026 22:00:57 -0600 Subject: [PATCH 2/4] Add _pyerr_handler() to init_step() and step() methods --- src/sksundae/_cy_cvode.pyx | 17 +++++++++++------ src/sksundae/_cy_ida.pyx | 25 +++++++++++++++++-------- 2 files changed, 28 insertions(+), 14 deletions(-) diff --git a/src/sksundae/_cy_cvode.pyx b/src/sksundae/_cy_cvode.pyx index f78ef87..1797288 100644 --- a/src/sksundae/_cy_cvode.pyx +++ b/src/sksundae/_cy_cvode.pyx @@ -834,6 +834,9 @@ cdef class CVODE: # 17) Advance solution in time flag = CVode(self.mem, tt, self.yy, &tout, itask) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) if flag == CV_ROOT_RETURN: @@ -894,6 +897,9 @@ cdef class CVODE: flag = CVode(self.mem, tend, self.yy, &tt, CV_NORMAL) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) if flag == CV_ROOT_RETURN: @@ -913,9 +919,7 @@ cdef class CVODE: ind += 1 - if PyErr_Occurred(): - _pyerr_handler() - elif stop: + if stop: break if self.aux.eventsfn: @@ -977,6 +981,9 @@ cdef class CVODE: while True: flag = CVode(self.mem, tend, self.yy, &tt, CV_ONE_STEP) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) if flag == CV_ROOT_RETURN: @@ -998,9 +1005,7 @@ cdef class CVODE: ind += 1 - if PyErr_Occurred(): - _pyerr_handler() - elif stop: + if stop: break if self.aux.eventsfn: diff --git a/src/sksundae/_cy_ida.pyx b/src/sksundae/_cy_ida.pyx index 458572b..4f370f2 100644 --- a/src/sksundae/_cy_ida.pyx +++ b/src/sksundae/_cy_ida.pyx @@ -854,11 +854,15 @@ cdef class IDA: if calc_initcond: flag = IDACalcIC(self.mem, ic_opt, ic_t0) - if flag < 0: + if PyErr_Occurred(): + _pyerr_handler() + elif flag < 0: raise RuntimeError("IDACalcIC - " + IDAMESSAGES[flag]) flag = IDAGetConsistentIC(self.mem, self.yy, self.yp) - if flag < 0: + if PyErr_Occurred(): + _pyerr_handler() + elif flag < 0: raise RuntimeError("IDAGetConsistentIC - " + IDAMESSAGES[flag]) self._initialized = True @@ -899,6 +903,9 @@ cdef class IDA: # 17) Advance solution in time flag = IDASolve(self.mem, tt, &tout, self.yy, self.yp, itask) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) svec2np(self.yp, yp_tmp) @@ -965,6 +972,9 @@ cdef class IDA: flag = IDASolve(self.mem, tend, &tt, self.yy, self.yp, IDA_NORMAL) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) svec2np(self.yp, yp_tmp) @@ -986,9 +996,7 @@ cdef class IDA: ind += 1 - if PyErr_Occurred(): - _pyerr_handler() - elif stop: + if stop: break if self.aux.eventsfn: @@ -1055,6 +1063,9 @@ cdef class IDA: while True: flag = IDASolve(self.mem, tend, &tt, self.yy, self.yp, IDA_ONE_STEP) + if PyErr_Occurred(): + _pyerr_handler() + svec2np(self.yy, yy_tmp) svec2np(self.yp, yp_tmp) @@ -1079,9 +1090,7 @@ cdef class IDA: ind += 1 - if PyErr_Occurred(): - _pyerr_handler() - elif stop: + if stop: break if self.aux.eventsfn: From b0bd0317091db49cca1e27a970d8dca64eb90cad Mon Sep 17 00:00:00 2001 From: "Corey R. Randall" Date: Mon, 23 Mar 2026 17:52:30 -0600 Subject: [PATCH 3/4] Remove unnecessary sparsity attr from *LSSparseDQJac classes --- docs/source/user_guide/linear_solvers.rst | 2 +- noxfile.py | 2 +- src/sksundae/_cy_cvode.pyx | 5 +---- src/sksundae/_cy_ida.pyx | 5 +---- src/sksundae/cvode/_solver.py | 4 ++-- src/sksundae/ida/_solver.py | 4 ++-- 6 files changed, 8 insertions(+), 14 deletions(-) diff --git a/docs/source/user_guide/linear_solvers.rst b/docs/source/user_guide/linear_solvers.rst index 78eb1a6..1b12c33 100644 --- a/docs/source/user_guide/linear_solvers.rst +++ b/docs/source/user_guide/linear_solvers.rst @@ -62,7 +62,7 @@ Switching between SUNDIALS direct linear solvers in scikit-SUNDAE is straightfor U = 0 # upper bandwidth solver = CVODE(rhsfn, linsolver='band', lband=L, uband=U) -Ensure you provide both the lower bandwidth `lband` and upper bandwidth `uband` when using the `band` linear solver. Each bandwidth defines the LARGEST distance between a non-zero element and the main diagonal, on either side, as shown in the figure below. Forgetting to set either bandwidth will raise an error. If `lband + uband` matches the dimension `N` of the matrix, the performance of the `band` and `dense` linear solvers will be approximately the same. +Ensure you provide both the lower bandwidth `lband` and upper bandwidth `uband` when using the `band` linear solver. Each bandwidth defines the LARGEST distance between a non-zero element and the main diagonal, on either side, as shown in the figure below. Forgetting to set either bandwidth will raise an error. If `lband + uband + 1` matches the dimension `N` of the matrix, the performance of the `band` and `dense` linear solvers will be approximately the same. .. figure:: figures/banded_jacobian.png :width: 50% diff --git a/noxfile.py b/noxfile.py index 123150d..5cc1b66 100644 --- a/noxfile.py +++ b/noxfile.py @@ -21,7 +21,7 @@ def run_cleanup(_) -> None: @nox.session(name='linter', python=False) def run_ruff(session: nox.Session) -> None: """ - Run ruff to check for linting errors. + Run ruff to check for linting errors Use the optional 'format' or 'format-unsafe' arguments to run ruff with the --fix or --unsafe-fixes option prior to the linter. You can also use 'stats' diff --git a/src/sksundae/_cy_cvode.pyx b/src/sksundae/_cy_cvode.pyx index 1797288..5e4533b 100644 --- a/src/sksundae/_cy_cvode.pyx +++ b/src/sksundae/_cy_cvode.pyx @@ -335,9 +335,7 @@ cdef class _cvLSSparseDQJac: """ cdef void* mem cdef AuxData aux - cdef object groups # dict[int, np.ndarray[int]] - cdef object sparsity # sparse.csc_matrix, shape(NEQ, NEQ) def __cinit__(self, AuxData aux): @@ -351,7 +349,6 @@ cdef class _cvLSSparseDQJac: self.aux = aux self.groups = groups - self.sparsity = aux.sparsity def __call__( self, @@ -368,7 +365,7 @@ cdef class _cvLSSparseDQJac: cdef np.ndarray[DTYPE_t, ndim=1] diff, inc, inc_inv, ytemp, yptemp aux = self.aux - sparsity = self.sparsity + sparsity = aux.sparsity ytemp = y.copy() yptemp = yp.copy() diff --git a/src/sksundae/_cy_ida.pyx b/src/sksundae/_cy_ida.pyx index 4f370f2..9764c1c 100644 --- a/src/sksundae/_cy_ida.pyx +++ b/src/sksundae/_cy_ida.pyx @@ -341,9 +341,7 @@ cdef class _idaLSSparseDQJac: """ cdef void* mem cdef AuxData aux - cdef object groups # dict[int, np.ndarray[int]] - cdef object sparsity # csc_matrix, shape(NEQ, NEQ) def __cinit__(self, AuxData aux): @@ -357,7 +355,6 @@ cdef class _idaLSSparseDQJac: self.aux = aux self.groups = groups - self.sparsity = aux.sparsity def __call__( self, @@ -377,7 +374,7 @@ cdef class _idaLSSparseDQJac: cdef np.ndarray[DTYPE_t, ndim=1] ytemp, yptemp, rtemp aux = self.aux - sparsity = self.sparsity + sparsity = aux.sparsity ytemp = y.copy() yptemp = yp.copy() diff --git a/src/sksundae/cvode/_solver.py b/src/sksundae/cvode/_solver.py index f2f6322..5a0e174 100644 --- a/src/sksundae/cvode/_solver.py +++ b/src/sksundae/cvode/_solver.py @@ -99,7 +99,7 @@ def __init__(self, rhsfn: Callable, **options) -> None: If 'constraints_idx' is not None, then this option must include an array of equal length specifying the types of constraints to apply. Values should be in `{-2, -1, 1, 2}` which apply `y[i] < 0`, - `y[i] <= 0`, `y[i] >=0`, and `y[i] > 0`, respectively. The + `y[i] <= 0`, `y[i] >= 0`, and `y[i] > 0`, respectively. The default is None. eventsfn : Callable or None, optional Events function with signature `g(t, y, events[, userdata])`. @@ -124,7 +124,7 @@ def __init__(self, rhsfn: Callable, **options) -> None: Number of events to track. The default is 0. jacfn : Callable or None, optional Jacobian function like `J(t, y, yp, JJ[, userdata])`. Fills the - pre-allocated memory for 'JJ' with values defined by the Jacobian + pre-allocated memory for `JJ` with values defined by the Jacobian `JJ[i,j] = dyp_i/dy_j`. An internal finite difference method is applied when None (default). Note that the template for `JJ` is determined by the linear solver choice: a 1D array for 'sparse' or diff --git a/src/sksundae/ida/_solver.py b/src/sksundae/ida/_solver.py index 054f1c1..7d87657 100644 --- a/src/sksundae/ida/_solver.py +++ b/src/sksundae/ida/_solver.py @@ -105,7 +105,7 @@ def __init__(self, resfn: Callable, **options) -> None: If 'constraints_idx' is not None, then this option must include an array of equal length specifying the types of constraints to apply. Values should be in `{-2, -1, 1, 2}` which apply `y[i] < 0`, - `y[i] <= 0`, `y[i] >=0`, and `y[i] > 0`, respectively. The + `y[i] <= 0`, `y[i] >= 0`, and `y[i] > 0`, respectively. The default is None. eventsfn : Callable or None, optional Events function with signature `g(t, y, yp, events[, userdata])`. @@ -130,7 +130,7 @@ def __init__(self, resfn: Callable, **options) -> None: Number of events to track. The default is 0. jacfn : Callable or None, optional Jacobian function like `J(t, y, yp, res, cj, JJ[, userdata])`. - The function should fill the pre-allocated memory for 'JJ' with the + The function should fill the pre-allocated memory for `JJ` with the values defined by `JJ[i,j] = dres_i/dy_j + cj*dres_i/dyp_j`. An internal finite difference method is applied when None (default). Note that the template for `JJ` is determined by the linear solver From 255b71c20ea4eaa944dcf57120ceb13ef473af57 Mon Sep 17 00:00:00 2001 From: "Corey R. Randall" Date: Mon, 23 Mar 2026 17:59:11 -0600 Subject: [PATCH 4/4] Update CHANGELOG for PR#51 --- .github/ISSUE_TEMPLATE/bug_report.yml | 2 +- .github/ISSUE_TEMPLATE/feature_request.yml | 1 + CHANGELOG.md | 3 +++ 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/.github/ISSUE_TEMPLATE/bug_report.yml b/.github/ISSUE_TEMPLATE/bug_report.yml index a8e0df5..38db148 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.yml +++ b/.github/ISSUE_TEMPLATE/bug_report.yml @@ -1,6 +1,6 @@ name: Bug Report description: File a bug report -title: "[Bug]: " +title: "[bug]: " labels: ["bug"] body: diff --git a/.github/ISSUE_TEMPLATE/feature_request.yml b/.github/ISSUE_TEMPLATE/feature_request.yml index f90d9ea..c787e66 100644 --- a/.github/ISSUE_TEMPLATE/feature_request.yml +++ b/.github/ISSUE_TEMPLATE/feature_request.yml @@ -1,5 +1,6 @@ name: Feature Request description: Suggest an enhancement or new feature +title: "[feature]: " labels: ["feature"] body: diff --git a/CHANGELOG.md b/CHANGELOG.md index 71dadfa..f29b461 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,12 +7,15 @@ - Custom `__reduce__` methods, allowing solvers to be serialized ([#38](https://github.com/NatLabRockies/scikit-sundae/pull/38)) ### Optimizations +None. ### Bug Fixes +- Address memory leak with raised exceptions caused by persisting solver/data objects ([#51](https://github.com/NatLabRockies/scikit-sundae/pull/51)) - Fixes issue where `jacfn` is ignored when using `sparse` linear solver ([#48](https://github.com/NatLabRockies/scikit-sundae/pull/48)) - Ensures exception propagations work correctly with numpy 2.4 release ([#41](https://github.com/NatLabRockies/scikit-sundae/pull/41)) ### Breaking Changes +None. ### Chores - Move to using `ruff` for linting and start tracking `Cython` addition issue ([#45](https://github.com/NatLabRockies/scikit-sundae/pull/45))