Return value for 2n-dimensional ODE system

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?

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

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

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

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.

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>

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)

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.

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

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

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

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.

Ah, thanks.

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