jrenaud90 / cyrk Goto Github PK
View Code? Open in Web Editor NEWRunge-Kutta ODE Integrator Implemented in Cython and Numba
License: Other
Runge-Kutta ODE Integrator Implemented in Cython and Numba
License: Other
Had to set cache_njit=False in the test suite for the helper.py conversion functions (see this commit) due to very long run times for Windows when running Python 3.11. This behavior was new with the development of CyRK v0.8.3 and PR #39.
Very hard to debug the issue as it was not occurring on my local machine (also Win 10 & Py 3.11). Perhaps is is a different numba version or something like that. Did not look into the specific versions of dependencies. #a675494
For systems of ODEs, each y could have a different absolute tolerance compared to one another. Right now only a constant atol is allowed and is applied equally to all y-variables.
I'm excited to try out the cython compiled ODE solver but having difficulty with what appears to be binary compatibility with numpy.
Error:
ValueError Traceback (most recent call last)
Input In [111], in <cell line: 2>()
1 from numba import njit
----> 2 from CyRK import cyrk_ode
File ~\miniconda3\envs\arc\lib\site-packages\CyRK\__init__.py:9, in <module>
6 from .nb.nbrk import nbrk_ode
8 # Import cython solver
----> 9 from CyRK.cy.cyrk import cyrk_ode
File ~\miniconda3\envs\arc\lib\site-packages\CyRK\cy\_cyrk.pyx:1, in init CyRK.cy.cyrk()
ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 96 from C header, got 88 from PyObject
Relevant installed package versions:
Everything but cyrk is installed from conda standard channel.
Create a cython-based cdef class integrator to improve performance by allowing users to define cython-based diffeq.
The form could look like:
def BaseClass:
# ... Integration functionality.
def diffeq(self):
# ... Placeholder diffeq.
def ProblemSolver(BaseClass):
def diffeq(self):
# ... Specific diffeq.
CyRK's tests are currently failing on ubuntu systems. Digging into it, it appears only some of the tests are failing. For example if you run a diffeq using the DOP853 (rk_method=2
) these appear to work fine. The test using RK23 or RK45 fails unless you bring the timespan down to (0, 10) or lower. Otherwise the integrator just sits there and hangs. It does not appear to be crashing but the test will eventually timeout and fail.
CyRK's CySolver only allows np.float64
y0 and dydt. This is due to limitations between fused types (used in cyrk_ode
) and cython cdef classes.
A workaround is to convert your N-dimensional complex system of ODEs into a 2N-dimensional floating point system of ODEs. Pass this onto CySolver, and then convert the output back into complex values.
As of v0.2.0, CyRK's cython integrator only allows for complex numbers. Non-complex differential equations have to be passed into the integrator by setting all of the imaginary portions to 0. This works fine but roughly doubles the memory usage and computation steps for various parts of the calculation. To improve performance, the integrator should allow for pure floats if complex calculations are not required.
This could either be done by:
y0
and use that type throughout the calculationcyrk_ode
is called. This wrapper (which could be pure python or numba code) would then check the type of y0
and then pass the arguments the appropriate cyrk_ode.
cyrk_ode_complex
and cyrk_ode_float
.The first method would be more ideal as the second would require maintaining two very similar functions. It would also introduce a (small?) overhead in that initial python wrapper function call.
See details (and hopefully a future solution) about this issue here.
Current workaround is to not build or upload wheels during CI. Users will have to recompile the cython code locally as of CyRK v0.2.0.
Hey. I am trying to install CyRK on an HPC environment. I have pip version 21.3.1.
I tried
pip3 install --user CyRK
and got
Collecting CyRK Could not find a version that satisfies the requirement CyRK (from versions: ) No matching distribution found for CyRK
Then I tried to use Conda to install the package on my local machine, and got
`PackagesNotFoundError: The following packages are not available from current channels:
Can someone please let me know why this package cannot be found? Thanks!
As of version 0.5.0, the cyrk_ode
function utilizes constants defined in a separate module CyRK.rk
. These same constants are used in the nbrk_ode
solver (however, the numba version expects np.ndarrays and this file defines things in terms of carrays). It would be good to have nbrk_ode
import these Runge-Kutta constants so that they are not defined in two different locations.
Currently the cyrk_ode
function uses python lists to store values found during integration. However, python lists lead to a big performance hit.
It was found during the creation of the CySolver
class that if numpy arrays are allocated at some large size then data can be stored in them instead. If the amount of data exceeds the size of these arrays then they can be expanded with a much lower performance hit than dealing with python lists.
This has already been coded up in the CySolver
class. It just needs to be ported over to cyrk_ode
and tested. I expect there will be a good amount of performance gained. Good key words to search for in CySolver
to see how it was implemented are "expected_size" and "time_domain_array_view".
Prior to v0.5.0, numba's fastmath
variable was turned on for nbrk_ode
. With the fix to issue #20 in PR #23 we noticed that the standard set of tests would fail on MacOS and Ubuntu systems. It was very random, even rerunning a test without a change in code could cause different results.
Some consistent findings:
fastmath=False
(in CyRK.nb.nbrk at the njit decorator) allowed tests to pass.PR #23 changed the way minimum step size was calculated. Prior to this change it was set to:
10 * np.finfo(np.float64).eps
now it is set, correctly, to
10. * abs(np.nextafter(t, direction * np.inf) - t)
For small t
this value is quite small (5e-323 vs. a EPS of ~1e-16). I believe there is something not playing well between this small value floating around (pun) while fastmath is on. At least on some operating systems.
It is more important that the step size is set correctly so, as of 0.5.0, fastmath is turned off for nbrk_ode
there is a slight performance hit with this change.
See the pytest.skip
in the accuracy test in "test_a_pysolve_ivp.py" test file. Most (all?) of these are failing at the np.allclose check. I have manually checked some other cases with backward integration and it worked fine. Not sure what the trigger is.
This is likely to do with the backend's messy indexing when it comes to forward and backward integration. Code could really use a second pair of eyes to help clean it up and squash those bugs.
Will need to check that any index fixes still work for backward integration with t_eval and dense_output are set.
During the original implementation of CyRK I misunderstood how the minimum step size check was supposed to work. Scipy defines the minimum step size based on the actual value of time,
import numpy as np
min_step = 10 * np.abs(np.nextafter(t, direction * np.inf) - t)
whereas CyRK originally had a constant step size based off of the floating point error,
defines the minimum step size based on the actual value of time,
import numpy as np
min_step = 10 * np.finfo(np.float64).eps
This latter approach works for when t is close to zero, but is not accurate when t is far from zero since the number of floats between numbers decreases at higher values.
Since min_step
is used to set the step_size from becoming too small during integration - having an incorrect value could lead to incorrect or unexpected behavior.
This bug affects both the cyrk_ode
and nbrk_ode
.
As can be seen in the README.md figure, when the time span is large (in that example, when the end step is >~ 1e4) the numba solver becomes equal to if not worse than the scipy solver. This also occurs when the relative integration tolerance is very small. So it appears to get worse with longer integration times in general.
This is not an issue for the cython solver.
The macos workflow for building wheels is defaulting to whatever the macos-latest
runner architecture is (which changed recently from x64-86 to arm64).
I think you could easily distribute both architectures by running the macos workflow on a macos-latest
and a macos-13
runner (they have arm64 and x64 architectures natively).
Right now there is no way (for either the Numba or Cython solver) to capture other variables besides the y-variables and time domain. However, there may be secondary variables that are calculated which the user may want to record without having to rerun functions.
Solution Idea:
If a differential equation is designed as the following where z
is a container for independent variables (much like dy),
def diffeq(t, y, dy, z, arg_0, arg_1, ...):
id_0 = ...
id_1 = ...
dy[0] = id_0 ...
dy[1] = id_1 ....
z[0] = id_0
z[1] = id_1
Then the solvers (at least the cython one) should be able to pass in that z argument and track it. The trick will be how to handle the output of the solver to keep it consistent for when the user does not want to track any independent variables.
Add unit tests and add link to GitHub workflow
As of v0.2.0, CyRK's cython integrator only allows for complex numbers. Non-complex differential equations have to be passed into the integrator by setting all of the imaginary portions to 0. This works fine but roughly doubles the memory usage and computation steps for various parts of the calculation. To improve performance, the integrator should allow for pure floats if complex calculations are not required.
This could either be done by:
y0
and use that type throughout the calculationcyrk_ode
is called. This wrapper (which could be pure python or numba code) would then check the type of y0
and then pass the arguments the appropriate cyrk_ode.
cyrk_ode_complex
and cyrk_ode_float
.The first method would be more ideal as the second would require maintaining two very similar functions. It would also introduce a (small?) overhead in that initial python wrapper function call.
I just figured out that t_eval
isn't behaving in the same way as scipy.integrate.solve_ivp
:
solve_ivp
solves the ode for the points in t_eval
Times at which to store the computed solution, must be sorted and lie within t_span. If None (default), use points selected by the solver.
CyRK
interpolates the solution t_eval
Both solvers uses an adaptive time stepping protocol based on the recent error at each step. This results in a final non-uniform time domain of variable size. If the user would like the results at specific time steps then they can provide a np.ndarray array at the desired steps via t_eval. The solver will then interpolate the results to fit this array.
In my kind of special case CyRK
is actually slower than solve_ivp
.
For performance and simplicity reasons, it was decided that the optional arguments (args
) passed to both the CySolver
and cyrk_ode
solvers must be given as floating point numbers (no complex values, strings, booleans, etc). It would be nice to offer the ability to pass other types of values via the args
parameter.
For the same reasons, the optional capture of additional outputs for both solvers is also limited to certain types (floats for CySolver; floats or complex for cyrk_ode).
In the current implementation of CyRK the only variables that are stored during integration are the dependent-y variables. For example,
# y = dy(t, y0)
t, y, success, message = nbrk_ode(dy, time_span, y0)
However, it is often useful to a user for additional parameters to be stored. For instance, a
and b
below.
def dy(t, y):
dy_ = np.empty(2, dtype=np.complex128)
a = g(t, y) # Some function of time and y
b = h(t, y) # Some function of time and y
dy_[0] = a * y[0]
dy_[1] = b * y[1]
return dy_
The values of a
and b
are lost after the integration is completed. The user would then need to make separate calls to the g
and h
function to retrieve their values. This is duplicating calculations and forcing the user to write additional code.
Propose CyRK stores these intermediate results during integration and returns them to the user.
One way this could be done is by having the user define the differential equation like,
def dy(t, y):
dy_ = np.empty(2, dtype=np.complex128)
a = g(t, y) # Some function of time and y
b = h(t, y) # Some function of time and y
dy_[0] = a * y[0]
dy_[1] = b * y[1]
return dy_, (a, b)
Then both cyrk_ode and nbrk_ode would store the result (probably best as a np.array) and return the result to the user.
To-Do
cyrk_ode
t_eval
nbrk_ode
t_eval
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.