jrenaud90 / cyrk Goto Github PK

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.

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:

• Have the cython integrator dynamically determine the type of y0 and use that type throughout the calculation
• Have a wrapper that is called before cyrk_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.
  • There would need to be two different cyrk_odes: 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.

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

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:

• python=3.9.12
• numba=0.55.1
• numpy=1.21.5
• cython=0.29.32
• cyrk=0.1.1

Everything but cyrk is installed from conda standard channel.

Add unit tests and add link to GitHub workflow

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.

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:

• Tests never failed on Windows.
• More tests would fail on MacOS vs. Ubuntu.
• Setting 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.

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.

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:

• cyrk`

Can someone please let me know why this package cannot be found? Thanks!

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:

• Have the cython integrator dynamically determine the type of y0 and use that type throughout the calculation
• Have a wrapper that is called before cyrk_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.
  • There would need to be two different cyrk_odes: 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.

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.

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".

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.

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.

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.

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 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.

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

• Implement intermediate saving for cyrk_ode
  • Implement basic saving
  • Implement interpolation for t_eval
  • Create tests
• Implement intermediate saving for nbrk_ode
  • Implement basic saving
  • Implement interpolation for t_eval
  • Create tests

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).