Return value for 2n-dimensional ODE system

Questions
1

I want to estimate parameters of a 2n-dimensional ODE system, where n is the number of my data points. I post a minimal example.

n=2

def odefunc(y, t, p):
    y1 = y[:n]
    y2 = y[n:]
    dy1 = -p[0]*y1*y2
    dy2 = p[0]*y1*y2 
    return [dy1[0], dy1[1], dy2[0], dy2[1]] # generalize this line

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

The ODE equations are simplified and do not make sense. What I wonder is how should rewrite the return line such that it works for any n. I tried return [dy1, dy2] but this does not work. I think, I have to concatenate dy1 and dy2 somehow. How do I do this?

2

The return shape must be the same as the shape of y.

print(y.shape)
...
result = numpy.concatenate([y1, y2]).T
print(result.shape)
3

Thanks! I’m not sure whether I fully understand. I tried:

n=2

def odefunc(y, t, p):
    y1 = y[:n]
    y2 = y[n:]
    dy1 = -p[0]*y1*y2
    dy2 = p[0]*y1*y2 
    return np.concatenate([dy1, dy2]).T

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

which throws an error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-152-11eea9a67857> in <module>
      8     return np.concatenate([dy1, dy2]).T
      9 
---> 10 pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

~/software/pymc3/pymc3/ode/ode.py in __init__(self, func, times, n_states, n_theta, t0)
     72         # Private
     73         self._augmented_times = np.insert(times, 0, t0).astype(floatX)
---> 74         self._augmented_func = utils.augment_system(func, self.n_states, self.n_theta)
     75         self._sens_ic = utils.make_sens_ic(self.n_states, self.n_theta, floatX)
     76 

~/software/pymc3/pymc3/ode/utils.py in augment_system(ode_func, n_states, n_theta)
     89 
     90     # Get symbolic representation of the ODEs by passing tensors for y, t and theta
---> 91     yhat = ode_func(t_y, t_t, t_p[n_states:])
     92     # Stack the results of the ode_func into a single tensor variable
     93     if not isinstance(yhat, (list, tuple)):

<ipython-input-152-11eea9a67857> in odefunc(y, t, p)
      6     dy1 = -p[0]*y1*y2
      7     dy2 = p[0]*y1*y2
----> 8     return np.concatenate([dy1, dy2]).T
      9 
     10 pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: zero-dimensional arrays cannot be concatenated
4

Looks like there might be an issue with how the sensitivies are computed. Traceback shows there is a problem in pymc3.ode.utils.augment_system which creates a system of ODEs for the gradients of the solution with respect to the parameters.

@michaelosthege Might be worth noting that

n=10

def odefunc(y, t, p):
    y1 = y[:n]
    y2 = y[n:]
    dy1 = -p[0]*y1*y2
    dy2 = p[0]*y1*y2 
    print(dy1.shape)
    print(dy2.shape)
    return np.concatenate((dy1,dy2))

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

Returns

Shape.0
Shape.0

So the shape of dy1 and dy2 are zero, which i don’t think is right. Thoughts? I can get to this a little later this evening.

5

I’m not sure, but is this the same error:

n=2

def odefunc(y, t, p):
    dy = -p[0]*y
    return dy

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=n, n_theta=1) 

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-7-63b1db971a54> in <module>
      5     return dy
      6 
----> 7 pm.ode.DifferentialEquation(odefunc, times=[1], n_states=n, n_theta=1)

~/software/pymc3/pymc3/ode/ode.py in __init__(self, func, times, n_states, n_theta, t0)
     72         # Private
     73         self._augmented_times = np.insert(times, 0, t0).astype(floatX)
---> 74         self._augmented_func = utils.augment_system(func, self.n_states, self.n_theta)
     75         self._sens_ic = utils.make_sens_ic(self.n_states, self.n_theta, floatX)
     76 

~/software/pymc3/pymc3/ode/utils.py in augment_system(ode_func, n_states, n_theta)
     96 
     97     # Now compute gradients
---> 98     J = tt.jacobian(t_yhat, t_y)
     99 
    100     Jdfdy = tt.dot(J, dydp)

~/miniconda3/envs/spols191122/lib/python3.7/site-packages/theano/gradient.py in jacobian(expression, wrt, consider_constant, disconnected_inputs)
   1876         "tensor.jacobian expects a Variable as `expression`"
   1877     assert expression.ndim < 2, \
-> 1878         ("tensor.jacobian expects a 1 dimensional variable as "
   1879          "`expression`. If not use flatten to make it a vector")
   1880 

AssertionError: tensor.jacobian expects a 1 dimensional variable as `expression`. If not use flatten to make it a vector

And I came up with this ugly workaround:

n=3

def odefunc(y, t, p):
    y1 = y[:n]
    y2 = y[n:]
    dy1 = -p[0]*y1*y2
    dy2 = p[0]*y1*y2 
    
    return_string = '[' 
    for i in range(n):
        return_string += f'dy1[{i}], '
    for i in range(n):
        return_string += f'dy2[{i}], '
    return_string += ']'
    
    print(return_string)
    
    return eval(return_string)

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

Output:

[dy1[0], dy1[1], dy1[2], dy2[0], dy2[1], dy2[2], ]
<pymc3.ode.ode.DifferentialEquation at 0x7fe4bb280860>
6

So the solution to this one is pretty straight forward. In the notebook describing how to use DifferentialEquation I write

The argument func needs to be written as if y and p are vectors. So even when your model has one state and/or one parameter, you should explicitly write y[0] and/or p[0] .

So the correct way to write this would be

def odefunc(y, t, p):
    dy = -p[0]*y[0]
    return dy

DifferentialEquation(odefunc, times=[1], n_states=1, n_theta=1)
1 Like
7

The problem seems to be with using numpy.concatenate. Consider the following two examples:

def odefunc_numpy(y, t, p):
    
    dy1 = -p[0]*y[0]
    dy2 = -p[1]*y[1]
    
    return np.concatenate([dy1, dy2])

def odefunc(y, t, p):
    
    dy1 = -p[0]*y[0]
    dy2 = -p[1]*y[1]
    
    return [dy1, dy2]

#No Error
DifferentialEquation(odefunc, times=[1], n_states=2, n_theta=2) 
#Error
DifferentialEquation(odefunc_numpy, times=[1], n_states=2, n_theta=2)

As I said before, I think the problem has to do with how the sensitivities are computed (in this case pymc3.ode.utils.augment_system). It is likely something about how Theano is handling the function. Still investigating.

8

However, in my example y was a 2-dimensional vector:

Hence, I would expect dy to be 2-dimensional.

9

Looks like theano does not support numpy.concatenate.

Either you can do a list comprehension (similar to what you did with the strings), or the more efficient way would be to do:

if isinstance(y, tt.TensorVariable):
   return tt.concatenate(...)
1 Like
10

Oh and Shape.0 does not mean that something is 0-dimensional. It’s the name of the TensorVariable that symbolically represents the shape of a TensorVariable. At the time of the print statement, the shape was not defined because augment_system called the ode function symbolically.
That’s probably also why numpy complained about 0-dim arrays - they were numpy.array(TensorVariable), hence 0-dimensional.

11

Ah, thanks.

12

Thanks! However, I can’t get it working. What am I doing wrong?

n=2

def odefunc(y, t, p):
    y1 = y[:n]
    y2 = y[n:]
    dy1 = -p[0]*y1*y2
    dy2 = p[0]*y1*y2 
    if isinstance(y, tt.TensorVariable):
        return tt.concatenate([dy1, dy2])

pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-337-aae7cc6efcde> in <module>
      9         return tt.concatenate([dy1, dy2])
     10 
---> 11 pm.ode.DifferentialEquation(odefunc, times=[1], n_states=2*n, n_theta=1)

~/software/pymc3/pymc3/ode/ode.py in __init__(self, func, times, n_states, n_theta, t0)
     72         # Private
     73         self._augmented_times = np.insert(times, 0, t0).astype(floatX)
---> 74         self._augmented_func = utils.augment_system(func, self.n_states, self.n_theta)
     75         self._sens_ic = utils.make_sens_ic(self.n_states, self.n_theta, floatX)
     76 

~/software/pymc3/pymc3/ode/utils.py in augment_system(ode_func, n_states, n_theta)
     96 
     97     # Now compute gradients
---> 98     J = tt.jacobian(t_yhat, t_y)
     99 
    100     Jdfdy = tt.dot(J, dydp)

~/miniconda3/envs/spols191122/lib/python3.7/site-packages/theano/gradient.py in jacobian(expression, wrt, consider_constant, disconnected_inputs)
   1876         "tensor.jacobian expects a Variable as `expression`"
   1877     assert expression.ndim < 2, \
-> 1878         ("tensor.jacobian expects a 1 dimensional variable as "
   1879          "`expression`. If not use flatten to make it a vector")
   1880 

AssertionError: tensor.jacobian expects a 1 dimensional variable as `expression`. If not use flatten to make it a vector
13

This is an old post, but the same problem still seems to exist. At least I am running into it with pymc3 3.11.4. I think the issue is caused by the following lines in pymc/ode/utils.py (pymc/utils.py at 71b18e6517f21563a8caad462bc7caea9f06f41d · pymc-devs/pymc · GitHub)

if not isinstance(yhat, (list, tuple)):
    yhat = (yhat,)
t_yhat = at.stack(yhat, axis=0)

When I change them to

if isinstance(yhat, at.TensorVariable):
    t_yhat = yhat
else:
    if not isinstance(yhat, (list, tuple)):
        yhat = (yhat,)
    t_yhat = tt.stack(yhat, axis=0)

I no longer get the error message. There should probably also be an explicit check whether yhat is a numpy array. Maybe someone who knows more about pymc3 and theano can comment on this.

14

Looks like the docstring on ode_func doesn’t do a great job of communicating the expected types of the return value (yhat).

It should be something like this:

  • list/tuple of scalars or 1D ndarray/TensorVariable → elements correspond to the state variables
  • scalar → there’s just one state variable
    The second case is what the if not isinstance() deals with.

The if/else there can definitely be improved, and so can our test coverage of these scenarios.
@mirko-m it would be great if you can submit a PR to fix it on main. We could then cherry-pick the commit also to the v3 branch so the fix can make it into the 3.11.5 release.

15

Sounds good. I submitted the PR: fix: check the output type of ode_func by mirko-m · Pull Request #5414 · pymc-devs/pymc · GitHub