pymc-devs / sunode Goto Github PK
View Code? Open in Web Editor NEWSolve ODEs fast, with support for PyMC
Home Page: https://sunode.readthedocs.io
License: MIT License
Solve ODEs fast, with support for PyMC
Home Page: https://sunode.readthedocs.io
License: MIT License
There are hard-codes dependencies that force to use conda. For example, here: build_cvodes.py
Hi
I think the docs are missing an example of how we can use sunode to compute gradients wrt to model parameters using the adjoint solver.
I'm not even sure what gradients are computed - of the solution ? Sundials talks about a "derived function", which I take to be , for example, a loss function. But not sure how this applies to sunode.
I understand I should use solve_backward; but the function is not documents. I'm especially confused about what the "grads" argument is, and what the output variables are - "grad_out" vs "lamda_out".
If you explain this I would be happy to share a notebook that uses sunode to fit an ODE - this could be useful for the documentation.
Thanks!
any plan to support m1 mac ?
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14574:14: error: implicit declaration of function 'SUNDenseLinearSolver' is invalid in C99 [-Werror,-Wimplicit-function-declaration]
{ result = SUNDenseLinearSolver(x0, x1); }
^
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14574:12: warning: incompatible integer to pointer conversion assigning to 'struct _generic_SUNLinearSolver *' from 'int' [-Wint-conversion]
{ result = SUNDenseLinearSolver(x0, x1); }
^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14589:31: error: too few arguments to function call, expected 3, have 2
return SUNDenseMatrix(x0, x1);
~~~~~~~~~~~~~~ ^
/opt/homebrew/include/sunmatrix/sunmatrix_dense.h:79:27: note: 'SUNDenseMatrix' declared here
SUNDIALS_EXPORT SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N, SUNContext sunctx);
^
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14615:35: error: too few arguments to function call, expected 3, have 2
{ result = SUNDenseMatrix(x0, x1); }
~~~~~~~~~~~~~~ ^
/opt/homebrew/include/sunmatrix/sunmatrix_dense.h:79:27: note: 'SUNDenseMatrix' declared here
SUNDIALS_EXPORT SUNMatrix SUNDenseMatrix(sunindextype M, sunindextype N, SUNContext sunctx);
^
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14918:10: error: implicit declaration of function 'SUNLapackDense' is invalid in C99 [-Werror,-Wimplicit-function-declaration]
return SUNLapackDense(x0, x1);
^
build/temp.macosx-11.0-arm64-cpython-310/_sundials_cvodes.c:14918:10: warning: incompatible integer to pointer conversion returning 'int' from a function with result type 'struct _generic_SUNLinearSolver *' [-Wint-conversion]
return SUNLapackDense(x0, x1);
^~~~~~~~~~~~~~~~~~~~~~
fatal error: too many errors emitted, stopping now [-ferror-limit=]
17 warnings and 20 errors generated.
error: command '/usr/bin/clang' failed with exit code 1
[end of output]
note: This error originates from a subprocess, and is likely not a problem with pip.
Hi, I would like to know if it is possible for sunode to be installed and used on linux-aarch64, as it is not mentioned in the conda forge environments?
Hello,
I tried to install sunode but was stuck with solving environment step. Could you have a look please?
To reproduce:
'''
conda config --add channels conda-forge
conda config --set channel_priority strict
conda install sunode
'''
Then I saw this in my terminal. It was stuck at solving environment and didn't manage to succeed.
Thanks!
Converting the solution of an ODE with multidimensional states to an xarray dataset fails if I choose unstack_state=True
:
import numpy as np
import matplotlib.pyplot as plt
import sunode
def multi_dim_ode(t, y, p):
return {"y": -p.a * y.y}
params = {
"a": (),
}
states = {
"y": ("ydim",),
}
problem = sunode.SympyProblem(
params=params,
states=states,
rhs_sympy=multi_dim_ode,
derivative_params=(),
coords={"ydim": np.array([0, 2])},
)
solver = sunode.solver.Solver(problem, compute_sens=False, solver="BDF")
y0 = np.zeros((), dtype=problem.state_dtype)
y0["y"] = np.array([1, 3])
tvals = np.linspace(0, 10, num=5)
solver.set_params_dict(
{
"a": 0.1,
}
)
output = solver.make_output_buffers(tvals)
solver.solve(t0=0, tvals=tvals, y0=y0, y_out=output)
# Works
solver.as_xarray(tvals, output, unstack_state=False)
# Does not work
solver.as_xarray(tvals, output, unstack_state=True)
It seems that the additional dimensions are not passed to xarray. The error I get is:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/variable.py in as_variable(obj, name)
101 try:
--> 102 obj = Variable(*obj)
103 except (TypeError, ValueError) as error:
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/variable.py in __init__(self, dims, data, attrs, encoding, fastpath)
303 self._data = as_compatible_data(data, fastpath=fastpath)
--> 304 self._dims = self._parse_dimensions(dims)
305 self._attrs = None
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/variable.py in _parse_dimensions(self, dims)
495 if len(dims) != self.ndim:
--> 496 raise ValueError(
497 "dimensions %s must have the same length as the "
ValueError: dimensions ('time',) must have the same length as the number of data dimensions, ndim=2
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
<ipython-input-13-f84c5315b76f> in <module>
40
41 # Does not work
---> 42 solver.as_xarray(tvals, output, unstack_state=True)
~/git_repos/sunode/sunode/solver.py in as_xarray(self, tvals, out, sens_out, unstack_state, unstack_params)
324
325 def as_xarray(self, tvals, out, sens_out=None, unstack_state=True, unstack_params=True):
--> 326 return self._problem.solution_to_xarray(
327 tvals, out, self._user_data,
328 sensitivity=sens_out,
~/git_repos/sunode/sunode/problem.py in solution_to_xarray(self, tvals, solution, user_data, sensitivity, unstack_state, unstack_params)
127 for name in state:
128 assert name not in data
--> 129 data[name] = ('time', state[name])
130 else:
131 data['solution'] = ('time', solution)
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/dataset.py in __setitem__(self, key, value)
1272 )
1273
-> 1274 self.update({key: value})
1275
1276 def __delitem__(self, key: Hashable) -> None:
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/dataset.py in update(self, other, inplace)
3573 """
3574 _check_inplace(inplace)
-> 3575 merge_result = dataset_update_method(self, other)
3576 return self._replace(inplace=True, **merge_result._asdict())
3577
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/merge.py in dataset_update_method(dataset, other)
862 other[key] = value.drop_vars(coord_names)
863
--> 864 return merge_core([dataset, other], priority_arg=1, indexes=dataset.indexes)
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/merge.py in merge_core(objects, compat, join, priority_arg, explicit_coords, indexes, fill_value)
550 coerced, join=join, copy=False, indexes=indexes, fill_value=fill_value
551 )
--> 552 collected = collect_variables_and_indexes(aligned)
553
554 prioritized = _get_priority_vars_and_indexes(aligned, priority_arg, compat=compat)
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/merge.py in collect_variables_and_indexes(list_of_mappings)
275 append_all(coords, indexes)
276
--> 277 variable = as_variable(variable, name=name)
278 if variable.dims == (name,):
279 variable = variable.to_index_variable()
~/miniconda3/envs/sunode/lib/python3.8/site-packages/xarray/core/variable.py in as_variable(obj, name)
103 except (TypeError, ValueError) as error:
104 # use .format() instead of % because it handles tuples consistently
--> 105 raise error.__class__(
106 "Could not convert tuple of form "
107 "(dims, data[, attrs, encoding]): "
ValueError: Could not convert tuple of form (dims, data[, attrs, encoding]): ('time', array([[1. , 3. ],
[0.77880078, 2.33640235],
[0.60653066, 1.81959198],
[0.47236655, 1.41709966],
[0.36787944, 1.10363832]])) to Variable.
While using coords I have noticed that the time to create SympyProblem
increases enormously with larger number of coords (see example below). Up to 40 it does not take too long (<2sec). But after 50, it takes up to minutes, which is unacceptable for me as I would like to use at least 100 coords. I am not sure if this is a bug related to sunode or putting to much stress on sympy...
Nonetheless, thank you in advance for looking in to it.
import numpy as np
import sunode
import sympy as sym
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
def I_func(I0_e, μ, dz, c_list):
I_list = sym.symarray('I', len(c_list))
for n,c in enumerate(c_list):
if n == 0:
I_list[n] = I0_e * sym.exp(-μ * dz * c)
else:
I_list[n] = I_list[n-1] * sym.exp(-μ * dz * c)
return I_list
def R_func(Φ,I,μ,c,dz):
return - (Φ * I * (1 - sym.exp(-μ * dz * c)) / dz)
R_func_v = np.vectorize(R_func)
# Define the right hand side of the ode
def rhs_sympy(t, c, p):
I_list = I_func(p.I0_e,p.μ,p.dz,c.A)
R_list = R_func_v(p.Φ,I_list,p.μ,c.A,p.dz)
return {
'A' : R_list
}
dz = 1e-3
z_list = np.arange(0,0.01,dz)
coords = {
'z_idx': pd.Index(np.arange(len(z_list)))
}
params = {
'Φ': (),
'μ': (),
'I0_e': (),
'dz': (),
}
c0 = {
'A' : 'z_idx'
}
print('Number of coords is:', len(z_list))
Number of coords is: 10
%%time
ode = sunode.symode.SympyProblem(
params,
c0,
rhs_sympy,
derivative_params=[],
coords=coords
)
CPU times: user 17.3 ms, sys: 0 ns, total: 17.3 ms
Wall time: 16.7 ms
solver = sunode.solver.Solver(ode)
solver.set_params_dict({
'Φ': 0.65,
'μ': 7675,
'I0_e': 0.0038,
'dz': dz,
})
tvals = np.linspace(0, 1., 1001)
c0 = np.full(ode.n_states,0.07)
out = solver.make_output_buffers(tvals)
%%time
solver.solve(0, tvals, c0, out)
CPU times: user 4.28 ms, sys: 0 ns, total: 4.28 ms
Wall time: 2.95 ms
plt.figure()
ax = plt.axes()
plt.plot(tvals,out[:,:]);
ax.set_xlabel('t [s]')
ax.set_ylabel('C$_{A}$')
ax.set_xscale('log')
for n in [2,5,10,20,30,40,50,75]:
df = 0.01
dz = df/n
z_list = np.arange(0,df,dz)
coords = {'z_idx': pd.Index(np.arange(len(z_list)))}
print('Number of coords is:', len(z_list))
%time ode = sunode.symode.SympyProblem(params, c0, rhs_sympy, derivative_params=[], coords=coords)
Number of coords is: 2
CPU times: user 14.3 ms, sys: 0 ns, total: 14.3 ms
Wall time: 13.3 ms
Number of coords is: 5
CPU times: user 30.2 ms, sys: 0 ns, total: 30.2 ms
Wall time: 28.5 ms
Number of coords is: 10
CPU times: user 71.6 ms, sys: 1.88 ms, total: 73.4 ms
Wall time: 73 ms
Number of coords is: 20
CPU times: user 279 ms, sys: 0 ns, total: 279 ms
Wall time: 278 ms
Number of coords is: 30
CPU times: user 670 ms, sys: 0 ns, total: 670 ms
Wall time: 670 ms
Number of coords is: 40
CPU times: user 1.81 s, sys: 0 ns, total: 1.81 s
Wall time: 1.81 s
Number of coords is: 50
CPU times: user 40.2 s, sys: 17 ms, total: 40.2 s
Wall time: 40.2 s
Number of coords is: 75
CPU times: user 5min 4s, sys: 0 ns, total: 5min 4s
Wall time: 5min 4s
from sinfo import sinfo
sinfo()
-----
matplotlib 3.5.1
numpy 1.21.5
pandas 1.3.5
sinfo 0.3.1
sunode 0.2.1
sympy 1.9
-----
IPython 7.30.1
jupyter_client 7.1.0
jupyter_core 4.9.1
jupyterlab 3.2.5
notebook 6.4.6
-----
Python 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:53) [GCC 9.4.0]
Linux-5.10.74.3-microsoft-standard-WSL2-x86_64-with-debian-bullseye-sid
12 logical CPU cores, x86_64
-----
Session information updated at 2022-01-05 15:13
The license file (MIT) and the license specified in meta.yml
(Apache 2) don't match.
@aseyboldt which one is correct?
Hello,
Is there a way to specify the external input to the RHS of the ODE? I.e., my ODE has an external force term x_dot = f(x, u, param).
Hi,
I've been using sunode for over 2 months now. It's a great project and I've really benefited from the speed.
I'm interested in using sunode with pyMC4. I ran the test script in the readMe and recieved a strange error that I'm not sure how to fix.
Traceback (most recent call last):
idata = pm.sample(tune=1000, draws=1000, chains=6, cores=6)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/pymc/sampling.py", line 524, in sample
step = assign_step_methods(model, step, methods=pm.STEP_METHODS, step_kwargs=kwargs)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/pymc/sampling.py", line 215, in assign_step_methods
tg.grad(model_logp, var)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 635, in grad
rval = _populate_grad_dict(var_to_app_to_idx, grad_dict, wrt, cost_name)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1440, in _populate_grad_dict
rval = [access_grad_cache(elem) for elem in wrt]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1440, in <listcomp>
rval = [access_grad_cache(elem) for elem in wrt]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1393, in access_grad_cache
term = access_term_cache(node)[idx]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1393, in access_grad_cache
term = access_term_cache(node)[idx]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1393, in access_grad_cache
term = access_term_cache(node)[idx]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1393, in access_grad_cache
term = access_term_cache(node)[idx]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in access_term_cache
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1064, in <listcomp>
output_grads = [access_grad_cache(var) for var in node.outputs]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1393, in access_grad_cache
term = access_term_cache(node)[idx]
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/gradient.py", line 1219, in access_term_cache
input_grads = node.op.L_op(inputs, node.outputs, new_output_grads)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/graph/op.py", line 371, in L_op
return self.grad(inputs, output_grads)
File "/home/aarcher/research/pdo_pathway_pymc/sunode/sunode/wrappers/as_aesara.py", line 261, in grad
lamda, gradient = backward(y0, params, params_fixed, g)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/graph/op.py", line 296, in __call__
node = self.make_node(*inputs, **kwargs)
File "/home/aarcher/anaconda3/envs/pymc_env/lib/python3.10/site-packages/aesara/graph/op.py", line 240, in make_node
raise TypeError(
TypeError: Invalid input types for Op SolveODEAdjointBackward{_solver_id=140214287550256, _t0=1900, _tvals_id=140214298068912}:
Input 3/4: Expected TensorType(float64, (None,)), got TensorType(float64, (1,))
Your help would be much appreciated. Thank you in advance!
Hi and thanks for this package!
I want to solve a set of ODEs where one of the parameters is a step function with prescribed times at which it can jump. What would be the recommended way of implementing something like that with sunode?
Excuse me, I want to solve a single DOF oscillator with a nonlinear spring, whose nonlinear stiffness is modelled as Piecewise Function. How I define ODE?
I try use sympy.Piecewise. While there were no errors and warnings reported, the results appeared to be wrong.
I appended a picture of the correct result calculated by scipy.solve_ivp. The ODE is
p = [1, 0.1, 100, 1000, 0]
def impact_1dof(t, x, p):
return [x[1],
1 / p[0] * ( - p[1] * x[1] - p[2] * x[0] - p[3] * (x[0] > p[4]) * (x[0] - p[4]))]
The reslut is correct as follow.
Next is code using sunode, and its reslut.
I need your help, if you understand what's going on. THANKS.
def SDOF(t, y, p):
# nonlinear_damper = pt.switch(y.x > p.delta1, p.cnl * y.v, 0.0)
return {
'x': y.v,
'v': -1/p.m*(p.k*y.x+p.c*y.v +p.knl* sp.Piecewise((0, y.x <= 0), (y.x, y.x > 0)))
}
problem = SympyProblem(
params={
'm': (),
'k': (),
'c': (),
'knl':(),
},
states={
'x': (),
'v': (),
},
rhs_sympy=SDOF,
derivative_params=[
('m',),
('c',),
('k',),
('knl',),
],
)
solver = Solver(problem, sens_mode=None, solver='BDF')
tvals = np.linspace(0, 5,10000)
y0 = np.zeros((), dtype=problem.state_dtype)
y0['x'] = -0.1
y0['v'] = 0
# We can also specify the parameters by name:
solver.set_params_dict({
'm': 1 ,
'k': 100,
'c': 0.1,
'knl': 100,
})
output = solver.make_output_buffers(tvals)
solver.solve(t0=0, tvals=tvals, y0=y0, y_out=output)
# We can convert the solution to an xarray Dataset
solx= solver.as_xarray(tvals, output).solution_x
solv = solver.as_xarray(tvals, output).solution_v
plt.figure()
plt.plot(solx,solv)
Hi!
What is the recommended way to install Sunode on Apple Silicon via Conda, please?
It seems sundoe does not have a conda-forge package for Apple Silicon: https://anaconda.org/conda-forge/sunode does not show any osx-arm64
label, unlike e.g. pymc:
Thanks for any help!
Hi everyone, I'm trying to run the example code in "Quickstart with PyMC3 — sunode documentation" but met the error "NotImplementedError: Cannot convert hares_start ~ HalfNormal to a tensor variable." But I can use "DifferentialEquation" in pymc3.ode to solve ODEs with pymc successfully.
I have pymc3.11.4, sunode 0.2.2, theano-py 1.1.2 (rather than theano), aesara 2.2.6 installed on macOS Big Sur11.6. I'm using python 3.9.7. I installed pymc using the command in their official documentation https://github.com/pymc-devs/pymc/wiki/Installation-Guide-(MacOS):
conda create -c conda-forge -n pymc3_env python=3.9 pymc3 theano-pymc mkl mkl-service
conda activate pymc3_en
and then installed sunode in the pymc3_env using
conda config --set channel_priority strict
conda install sunode
I also tried to uninstall aesara and run again, but then I got another error :name "aesara" is not defined.
Could you please provide some suggestions? Thank you a lot!
Currently the rhs function is defined using sympy and then compiled using numba.
Now that aesara has a numba backend we can use aesara to compile this function instead, which should speed up the compilation of larger models a lot and also makes it easier to define the ODE, because we follow numpy more closely.
With sympy 1.8 and sunode 0.1.1 (both installed with conda from conda-forge) I get an error when importing sunode:
import sunode
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
<ipython-input-11-d67c482f6196> in <module>
----> 1 import sunode
~/miniconda3/envs/sp1-rosa/lib/python3.8/site-packages/sunode/__init__.py in <module>
3 from sunode.vector import empty_vector, from_numpy
4 from sunode.matrix import empty_matrix
----> 5 from sunode.symode import SympyProblem
6 import sunode.solver
7
~/miniconda3/envs/sp1-rosa/lib/python3.8/site-packages/sunode/symode/__init__.py in <module>
----> 1 from sunode.symode.problem import SympyProblem
2
3 __all__ = ['SympyProblem']
~/miniconda3/envs/sp1-rosa/lib/python3.8/site-packages/sunode/symode/problem.py in <module>
12
13 from sunode import problem, basic, dtypesubset
---> 14 from sunode.symode.lambdify import lambdify_consts
15
16
~/miniconda3/envs/sp1-rosa/lib/python3.8/site-packages/sunode/symode/lambdify.py in <module>
7 from functools import partial
8
----> 9 from sympy.printing.pycode import SciPyPrinter
10 import sympy
11 import numpy as np
ImportError: cannot import name 'SciPyPrinter' from 'sympy.printing.pycode' (/home/anna/miniconda3/envs/sp1-rosa/lib/python3.8/site-packages/sympy/printing/pycode.py)
The SciPyPrinter
has been moved to sympy.printing.numpy.SciPyPrinter
in sympy/sympy@8406514 so updating the import should probably fix this.
y0 = {'A': (u ~ Uniform, ()), 'B': array(1.), 'C': array(1.)}
params = {'alpha': array(1.), 'beta': array(1.), 'extra': array([0.])}
solution, flat_solution, problem, sol, y0_flat, params_subs_flat, flat_sens, wrapper = sunode.wrappers.as_theano.solve_ivp(
y0=y0,
params=params,
rhs=dydt_dict,
tvals=t,
t0=t[0],
derivatives="forward",
solver_kwargs=dict(sens_mode="simultaneous", compute_sens=True)
)
Traceback:
[...]
solution, flat_solution, problem, sol, y0_flat, params_subs_flat, flat_sens, wrapper = sunode.wrappers.as_theano.solve_ivp(
/opt/conda/lib/python3.8/site-packages/sunode/wrappers/as_theano.py:123: in solve_ivp
sol = solver.Solver(problem, **solver_kwargs)
/opt/conda/lib/python3.8/site-packages/sunode/solver.py:252: in __init__
sens_rhs = self._problem.make_sundials_sensitivity_rhs()
/opt/conda/lib/python3.8/site-packages/sunode/problem.py:262: in make_sundials_sensitivity_rhs
sens_rhs = self.make_sensitivity_rhs()
/opt/conda/lib/python3.8/site-packages/sunode/symode/problem.py:473: in make_sensitivity_rhs
dydp = lambdify_consts(
/opt/conda/lib/python3.8/site-packages/sunode/symode/lambdify.py:265: in lambdify_consts
spec.loader.exec_module(module)
<frozen importlib._bootstrap_external>:779: in exec_module
???
/opt/conda/lib/python3.8/site-packages/sunode/symode/lambdify.py:197: in get_code
return self.source_to_code(self._asts[fullname])
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
data = <_ast.Module object at 0x7f2eb9bb3190>, path = '<string>'
@staticmethod
def source_to_code(data, path='<string>'):
"""Compile 'data' into a code object.
The 'data' argument can be anything that compile() can handle. The'path'
argument should be where the data was retrieved (when applicable)."""
> return compile(data, path, 'exec', dont_inherit=True)
E ValueError: empty body on FunctionDef
/opt/conda/lib/python3.8/importlib/abc.py:244: ValueError
Hi. I have suffered the same problem (infinite loops on CVODE) as you in this post, (https://githubmemory.com/repo/LLNL/sundials/issues/44)
(The interesting point to add is that, sometimes it gets certain values after taking plenty of times...)
Is there any update on this problem?
I have used the version:
pymc3 3.9.3,
sunode 0.1.1, and
sundials 5.3.0
If there is an update on this problem, I plan to update the newest version.
Very simple question - can I use a non symbolic version of the rhs function and if so, how?
Specifically, a case of ode(t ,y, p) where t, y and p are lists or numpy arrays, and the function returns dYdt (a list or numpy array where each element is a different state equation)
When running the tutorial code (https://sunode.readthedocs.io/en/latest/without_pymc.html) I get the following error from this line:
solver = sunode.solver.Solver(problem, compute_sens=False, solver='BDF')
> __init__() got an unexpected keyword argument 'compute_sens'
And also the following error:
plt.plot(output.view(tvals, problem.state_dtype)['hares'])
> ValueError: Type must be a sub-type of ndarray type
I use:
sundials 5.8.0 h7ddfbfd_0 conda-forge
sunode 0.3.0 py38h5e27569_0 conda-forge
When computing the gradient with respect to initial conditions, I find that the adjoint and fwd gradient computations differ slightly. Following equation 14 of CVODES manual, the discrepancy is in the adjoint equation and can be corrected by adding the constant, -np.matmul(sens0, lambda_out - grads[0, :]), to grad_out from solve_backward. sens0 is the initial sensitivities, lambda_out is the adjoint variables at time 0 and grads[0, :]) is the derivative of the likelihood with respect to the state variables at time 0.
I have to adjust lambda at 0 by -grads[0, :] because line 691 of sunode/solver.py appears to not loop over the initial time point of grads.
Does SolveODEAdjointBackward in sunode/sunode/wrappers/as_aesara.py implement a similar correction constant when computing the gradient wrt initial conditions?
I attached some code below as an example. It is based on the example script of the readMe. It computes the gradient of the sum of Hares and Lynx at time 0, 5 and 10 wrt alpha, beta and hares0 where hares0 is log10(Hares(0)).
import numpy as np
import sunode
import sunode.wrappers.as_aesara
import pymc as pm
import matplotlib.pyplot as plt
lib = sunode._cvodes.lib
def lotka_volterra(t, y, p):
"""Right hand side of Lotka-Volterra equation.
All inputs are dataclasses of sympy variables, or in the case
of non-scalar variables numpy arrays of sympy variables.
"""
return {
'hares': p.alpha * y.hares - p.beta * y.lynx * y.hares,
'lynx': p.delta * y.hares * y.lynx - p.gamma * y.lynx,
}
# initialize problem
problem = sunode.symode.SympyProblem(
params={
# We need to specify the shape of each parameter.
# Any empty tuple corresponds to a scalar value.
'alpha': (),
'beta': (),
'gamma': (),
'delta': (),
'hares0': ()
},
states={
# The same for all state variables
'hares': (),
'lynx': (),
},
rhs_sympy=lotka_volterra,
derivative_params=[
# We need to specify with respect to which variables
# gradients should be computed.
('alpha',),
('beta',),
('hares0',),
],
)
tvals = np.linspace(0, 10, 3)
y0 = np.zeros((), dtype=problem.state_dtype)
y0['hares'] = 1e0
y0['lynx'] = 0.1
params_dict = {
'alpha': 0.1,
'beta': 0.2,
'gamma': 0.3,
'delta': 0.4,
'hares0': 1e0
}
sens0 = np.zeros((3, 2))
sens0[2,0] = np.log(10)*1e0
solver = sunode.solver.Solver(problem, solver='BDF', sens_mode='simultaneous')
yout, sens_out = solver.make_output_buffers(tvals)
# gradient via fwd senstivity
solver.set_params_dict(params_dict)
output = solver.make_output_buffers(tvals)
solver.solve(t0=0, tvals=tvals, y0=y0, y_out=yout, sens0=sens0, sens_out=sens_out)
grad_out_fwd = [ sens_out[:,j,:].sum() for j in range(3)]
print(grad_out_fwd)
# gradient via adj senstivity
solver = sunode.solver.AdjointSolver(problem, solver='BDF')
solver.set_params_dict({
'alpha': 0.1,
'beta': 0.2,
'gamma': 0.3,
'delta': 0.4,
'hares0': 1e0
})
tvals_expanded = np.linspace(0, 10, 21)
yout, grad_out, lambda_out = solver.make_output_buffers(tvals_expanded)
lib.CVodeSetMaxNumSteps(solver._ode, 10000)
solver.solve_forward(t0=0, tvals=tvals, y0=y0, y_out=yout)
grads = np.zeros_like(yout)
grads[::10,:] = 1
solver.solve_backward(t0=tvals_expanded[-1], tend=tvals_expanded[0], tvals=tvals_expanded[1:-1],
grads=grads, grad_out=grad_out, lamda_out=lambda_out)
grad_out_adj = -np.matmul(sens0, lambda_out -grads[0, :]) + grad_out
print(grad_out_adj)
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.