Add Python interfaces based on nanobind #796
Conversation
Steven-Roberts
left a comment
Steven-Roberts
left a comment
There was a problem hiding this comment.
Finished first pass through the PR
gardner48
left a comment
gardner48
left a comment
There was a problem hiding this comment.
Finished sundials directory
| headers: | ||
| - ../../include/sundials/sundials_errors.h | ||
| enum_exclude_by_name__regex: | ||
| - "SUNErrCode_" # not needed in Python |
There was a problem hiding this comment.
Do users still have access to the SUN_ERR_* constants?
There was a problem hiding this comment.
No. Unfortunately, the x macro we use in sundials_error.h makes it to where the generator is unable to parse them.
| { | ||
| if (!sunctx->python) | ||
| { | ||
| sunctx->python = SUNContextFunctionTable_Alloc(); |
There was a problem hiding this comment.
Should we explicitly write the Create wrapper and always allocate the function table?
| [](SUNContext sunctx, | ||
| std::function<std::remove_pointer_t<SUNErrHandlerFn>> err_fn) | ||
| { | ||
| if (!sunctx->python) |
There was a problem hiding this comment.
Check for NULL function
| if (!sunctx->python) | |
| if (!err_fn) { throw sundials4py::illegal_value("err_fn was null"); } | |
| if (!sunctx->python) |
| { | ||
| auto fn_table = static_cast<SUNLinearSolverFunctionTable*>( | ||
| std::malloc(sizeof(SUNLinearSolverFunctionTable))); | ||
| std::memset(fn_table, 0, sizeof(SUNLinearSolverFunctionTable)); |
| m.def( | ||
| "SUNLogger_QueueMsg", | ||
| [](SUNLogger logger, SUNLogLevel lvl, const char* scope, const char* label, | ||
| const char* msg_txt) -> SUNErrCode |
There was a problem hiding this comment.
Does this handle the variadic inputs?
gardner48
left a comment
gardner48
left a comment
There was a problem hiding this comment.
Finished sundials classes starting on test directory
| # Programmer(s): Cody J. Balos @ LLNL | ||
| # ----------------------------------------------------------------- | ||
| # SUNDIALS Copyright Start | ||
| # Copyright (c) 2002-2025, Lawrence Livermore National Security |
| # Programmer(s): Cody J. Balos @ LLNL | ||
| # ----------------------------------------------------------------- | ||
| # SUNDIALS Copyright Start | ||
| # Copyright (c) 2002-2025, Lawrence Livermore National Security |
| # Programmer(s): Cody J. Balos @ LLNL | ||
| # ----------------------------------------------------------------- | ||
| # SUNDIALS Copyright Start | ||
| # Copyright (c) 2002-2025, Lawrence Livermore National Security |
| - "^SUNBandMatrix_Data$" | ||
| # We don't interface these function in Python. Instead users can index the numpy array returned by Data. | ||
| - "^SUNBandMatrix_Cols$" | ||
| - "^SUNBandMatrix_Column$" |
There was a problem hiding this comment.
If it's not already we should note in the user guide which functions are not available from python
| auto owner = nb::find(A); | ||
| auto ptr = SUNDenseMatrix_Data(A); | ||
| // SUNDenseMatrix_Data returns a column-major ordered array (i.e., Fortran style) | ||
| return nb::ndarray<sunrealtype, nb::numpy, nb::ndim<2>, |
There was a problem hiding this comment.
Returning this as a 2D array contradicts the comment here about matching the C function. However, reshaping is more natural in the dense case than the banded.
| rows, cols, nnz = 3, 3, 4 | ||
| A = SUNSparseMatrix(rows, cols, nnz, SUN_CSR_MAT, sunctx) | ||
| ret = SUNMatZero(A) | ||
| assert ret == 0 |
There was a problem hiding this comment.
| assert ret == 0 | |
| assert ret == SUN_SUCCESS | |
| dataA = SUNSparseMatrix_Data(A) | |
| np.testing.assert_array_equal(dataA, np.zeros_like(dataA)) |
| assert np.allclose(dataB, dataA) | ||
| assert np.allclose(idx_valsB, idx_vals) | ||
| assert np.allclose(idx_ptrsB, idx_ptrs) |
There was a problem hiding this comment.
| assert np.allclose(dataB, dataA) | |
| assert np.allclose(idx_valsB, idx_vals) | |
| assert np.allclose(idx_ptrsB, idx_ptrs) | |
| np.testing.assert_allclose(dataB, dataA) | |
| np.testing.assert_allclose(idx_valsB, idx_vals) | |
| np.testing.assert_allclose(idx_ptrsB, idx_ptrs) |
| ret = SUNMatScaleAdd(3.0, A, B) | ||
| assert ret == 0 | ||
| # 3*A + B = [3+2, 3+2, 3+2, 6+4] = [5, 5, 5, 10] | ||
| assert np.allclose(dataA, [5.0, 5.0, 5.0, 10.0]) |
There was a problem hiding this comment.
| assert np.allclose(dataA, [5.0, 5.0, 5.0, 10.0]) | |
| np.testing.assert_allclose(dataA, [5.0, 5.0, 5.0, 10.0]) |
| # Diagonal elements should be 3*2+1=7, off-diagonal unchanged | ||
| # So dataA = [7.0, 7.0, 7.0, 2.0] (assuming diagonal at idx_vals[0:3]) | ||
| # But since idx_vals = [0,1,2,0], diagonal is at positions 0,1,2 | ||
| assert np.allclose(dataA[:3], [7.0, 7.0, 7.0]) |
There was a problem hiding this comment.
Guessing this is AI not quite getting it right
| # Diagonal elements should be 3*2+1=7, off-diagonal unchanged | |
| # So dataA = [7.0, 7.0, 7.0, 2.0] (assuming diagonal at idx_vals[0:3]) | |
| # But since idx_vals = [0,1,2,0], diagonal is at positions 0,1,2 | |
| assert np.allclose(dataA[:3], [7.0, 7.0, 7.0]) | |
| # Diagonal elements should be 3*2+1=7, off-diagonal should be 3*2=6 | |
| np.testing.assert_allclose(dataA, [7.0, 7.0, 7.0, 6.0]) |
| # y[0] = 1.0*x[0] = 1.0 | ||
| # y[1] = 1.0*x[1] = 1.0 | ||
| # y[2] = 1.0*x[2] + 2.0*x[0] = 1.0 + 2.0 = 3.0 | ||
| assert np.allclose(N_VGetArrayPointer(y), [1.0, 1.0, 3.0]) |
There was a problem hiding this comment.
| assert np.allclose(N_VGetArrayPointer(y), [1.0, 1.0, 3.0]) | |
| np.testing.assert_allclose(N_VGetArrayPointer(y), [1.0, 1.0, 3.0]) |
gardner48
left a comment
gardner48
left a comment
There was a problem hiding this comment.
Finished all of the sundials4py directory. Only the doc directory to go but I probably won't get to that today, so feel free to start pushing updates.
| sol = N_VClone(y) | ||
| ode_problem.solution(y, sol, tret) | ||
|
|
||
| assert np.allclose(N_VGetArrayPointer(sol), N_VGetArrayPointer(y), atol=1e-2) |
There was a problem hiding this comment.
Find and replace is probably easier than suggesting all the places this occurs. My understanding is using np.testing.assert_allclose has helpful error messages when the check fails.
|
|
||
|
|
||
| @pytest.mark.skipif( | ||
| sunrealtype == np.float32, reason="Test not supported for sunrealtype=np.float32" |
There was a problem hiding this comment.
Looks like the tolerances are precision-dependent, what's the issue testing with single precision?
| status = ARKodeSStolerances(ark.get(), SUNREALTYPE_RTOL, SUNREALTYPE_ATOL) | ||
| assert status == ARK_SUCCESS | ||
|
|
||
| status = ARKodeSetLinearSolver(ark.get(), ls, None) |
There was a problem hiding this comment.
Might be useful to have a direct solver version of this test as well to check attaching a matrix (and maybe a Jacobian function)
| sol = N_VClone(y) | ||
| ode_problem.solution(y, sol, tret) | ||
|
|
||
| assert np.allclose(N_VGetArrayPointer(sol), N_VGetArrayPointer(y), atol=1e-2) |
There was a problem hiding this comment.
Should the comparison tolerance be tied to the precision?
|
|
||
| [project] | ||
| name = "sundials4py" | ||
| version = "7.6.0" |
There was a problem hiding this comment.
Update scripts/updateVersion.sh to update this file
gardner48
left a comment
gardner48
left a comment
There was a problem hiding this comment.
Finished all of the sundials4py directory. Only the doc directory to go but I probably won't get to that today, so feel free to start pushing updates.
This PR adds sundials4py, which is a Python interface to SUNDIALS.
The generator code can be found here: https://github.com/sundials-codes/sundials4py-generate.
Rendered docs: https://sundials--796.org.readthedocs.build/en/796/.
Review Requests
@drreynolds:
doc/bindings/sundials4py/arkodebindings/sundials4py/idasbindings/sundials4py/kinsolbindings/sundials4py/sundomeigestbindings/sundials4py/sunlinsolbindings/sundials4py/sunmatrixbindings/sundials4py/examples@gardner48:
Everything but can skip:
_generated.hppfiles (Steven and I will review these)bindings/sundials4py/examples(Steven and Dan are both reviewing these)bindings/sundials4py/arkode(Steven and Dan both reviewing these)bindings/sundials4py/idas(Steven and Dan both reviewing these)@Steven-Roberts
_generated.hppfiles inbindings/sundials4pybindings/sundials4py/arkodebindings/sundials4py/idasbindings/sundials4py/includebindings/sundials4py/sundialsbindings/sundials4py/nvectorbindings/sundials4py/sunadaptcontrollerbindings/sundials4py/sunadjointcheckpointschemebindings/sundials4py/sunnonlinsolbindings/sundials4py/sunmemorybindings/sundials4py/examples