Problem simulating a 2 dimensional system

[fccoelho]

Hi, I am trying to use nsim to simulate the following model:

"""
SIR Epidemic Model

S --> I --> R

X(t) = [S(t), I(t)]

Equations:

dS = -r0*S*I dt -sqrt(r0*S*I) dW1
dI = r0*S*I dt + sqrt(r0*S*I) dW1 -sqrt(I) dW2

"""

import nsim
import numpy as np

class SIR(nsim.ItoModel):
    r0 = 1.5
    y0 = np.array([0.9, 0.1]).T

    def f(self, y, t):
        A = np.array([[-self.r0*y[0]*y[1]],
                       [self.r0*y[0]*y[1] - y[1]]])
        return A

    def G(self, y, t):
        B = np.array([[-np.sqrt(self.r0*y[0]*y[1]), 0],
                       [np.sqrt(self.r0*y[0]*y[1]), - y[1]]])
        return B

To simulate it I tried the following:

s = nsim.Simulation(SIR, T=520.0)

It seems to work but when I try to plot the results, I get this error:

Traceback (most recent call last):
  File "SDE_SIR.py", line 38, in <module>
    print(s.output, dir(s))
  File "/usr/local/lib/python3.5/dist-packages/nsim/nsim.py", line 1330, in output
    self.compute()
  File "/usr/local/lib/python3.5/dist-packages/nsim/nsim.py", line 1313, in compute
    ar = self.system.integrate(tspan)
  File "/usr/local/lib/python3.5/dist-packages/nsim/nsim.py", line 870, in integrate
    ar = self.integrator[0](self.f, self.G, self.y0, tspan)
  File "/usr/local/lib/python3.5/dist-packages/sdeint/integrate.py", line 125, in itoint
    return chosenAlgorithm(f, G, y0, tspan)
  File "/usr/local/lib/python3.5/dist-packages/sdeint/integrate.py", line 305, in itoSRI2
    return _Roessler2010_SRK2(f, G, y0, tspan, Imethod, dW, I)
  File "/usr/local/lib/python3.5/dist-packages/sdeint/integrate.py", line 439, in _Roessler2010_SRK2
    H20b = np.reshape(H20, (d, 1))
  File "/usr/local/lib/python3.5/dist-packages/numpy/core/fromnumeric.py", line 225, in reshape
    return reshape(newshape, order=order)
ValueError: total size of new array must be unchanged

Is this a bug, or maybe a mistake in my model?

[mattja]

For a 2-dimensional system nsim expects the state y to be an array of shape (2,)
Your function f() returns an array of shape (2,1).

If you remove the extra brackets in f then it works:

    def f(self, y, t):
        A = np.array([-self.r0*y[0]*y[1],
                       self.r0*y[0]*y[1] - y[1]])
        return A

[mattja]

The fact that you got an error message from numpy instead of a more helpful message means that nsim could use some better input validation. So this is a bug in nsim.

It's reasonable for people to give a (2,1) array as a column vector. I'll probably fix this bug by enhancing nsim so that the state y can be an array with arbitrary number of dimensions (using numpy.ravel).
Then your code would work as-is.

[fccoelho]

Thanks! It is now working.