ODE with non-analytical time-dependent parameters in PyMC3

Question

I'm working on solving the following ODE with PyMC3:

def production( y, t, p ):
    return p[0]*getBeam( t ) - p[1]*y[0]

The getBeam( t ) is my time dependent coefficient. Those coefficients are given by an array of data which is accessed by the time index as follows:

def getBeam( t ):
    nBeam = I[int(t/10)]*pow( 10, -6 )/q_e
    return nBeam

I have successfully implemented it by using the scipy.integrate.odeint, but I have hard time to do it with pymc3.ode. In fact, by using the following:

ode_model = DifferentialEquation(func=production, times=x, n_states=1, n_theta=3, t0=0)
with pm.Model() as model:
    a = pm.Uniform( "S-Factor", lower=0.01, upper=100 )
    ode_solution = ode_model(y0=[0], theta=[a, Yield, lambd])

I obviously get the error TypeError: __trunc__ returned non-Integral (type TensorVariable), as the t is a TensorVariable, thus can not be used to access the array in which the coefficients are stored.

Is there a way to overcome this difficulty? I thought about using the theano.function but I can not get it working since, unfortunately, the coefficients can not be expressed by any analytical function: they are just stored inside the array I which index represents the time variable t.

Thanks

Why do you not give the numerical constant `pow( 10, -6 )` directly as `1e-6`? — Lutz Lehmann, Apr 21 '21 at 17:59
Yes, you're right, I can write it in that way. Thank you, but the problem still persists. — Jakub Skowroński, Apr 21 '21 at 18:11

score 0 · Answer 1 · answered May 18 '21 at 20:29

Since you already have a working implementation with scipy.integrate.odeint, you could use theano.compile.ops.as_op, though it comes with some inconveniences (see how to fit a method belonging to an instance with pymc3? and How to write a custom Deterministic or Stochastic in pymc3 with theano.op?)

Using your exact definitions for production and getBeam, the following code seems to work for me:

from scipy.integrate import odeint
from theano.compile.ops import as_op
import theano.tensor as tt
import pymc3 as pm

def ScipySolveODE(a):
    return odeint(production, y0=[0], t=x, args=([a, Yield, lambd],)).flatten()

@as_op(itypes=[tt.dscalar], otypes=[tt.dvector])
def TheanoSolveODE(a):
    return ScipySolveODE(a)

with pm.Model() as model:
    a = pm.Uniform( "S-Factor", lower=0.01, upper=100 )
    ode_solution = TheanoSolveODE(a)

Sorry I know this is more of a workaround than an actual solution...

ODE with non-analytical time-dependent parameters in PyMC3

1 Answers1