How to read a system of differential equations from a text file to solve the system with scipy.odeint?

Question

I have a large (>2000 equations) system of ODE's that I want to solve with python scipy's odeint.

I have three problems that I want to solve (maybe I will have to ask 3 different questions?). For simplicity, I will explain them here with a toy model, but please keep in mind that my system is large. Suppose I have the following system of ODE's:

dS/dt = -beta*S
dI/dt = beta*S - gamma*I
dR/dt = gamma*I

with beta = cpI

where c, p and gamma are parameters that I want to pass to odeint.

odeint is expecting a file like this:

def myODEs(y, t, params):
    c,p, gamma = params
    beta = c*p
    S = y[0]
    I = y[1]
    R = y[2]
    dydt = [-beta*S*I,
           beta*S*I - gamma*I,
           - gamma*I]  
    return dydt

that then can be passed to odeint like this:

myoutput = odeint(myODEs, [1000, 1, 0], np.linspace(0, 100, 50), args = ([c,p,gamma], ))

I generated a text file in Mathematica, say myOdes.txt, where each line of the file corresponds to the RHS of my system of ODE's, so it looks like this

#myODEs.txt

-beta*S*I
beta*S*I - gamma*I
- gamma*I

My text file looks similar to what odeint is expecting, but I am not quite there yet. I have three main problems:

1. How can I pass my text file so that odeint understands that this is the RHS of my system?
2. How can I define my variables in a smart way, that is, in a systematic way? Since there are >2000 of them, I cannot manually define them. Ideally I would define them in a separate file and read that as well.
3. How can I pass the parameters (there are a lot of them) as a text file too?

I read this question that is close to my problems 1 and 2 and tried to copy it (I directly put values for the parameters so that I didn't have to worry about my point 3 above):

    systemOfEquations = []
    with open("myODEs.txt", "r") as fp :
        for line in fp :
            systemOfEquations.append(line)

    def dX_dt(X, t):
        vals = dict(S=X[0], I=X[1], R=X[2], t=t)
        return [eq for eq in systemOfEquations]

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 1, 5))

but I got the error:

    odepack.error: Result from function call is not a proper array of          floats.
    ValueError: could not convert string to float: -((12*0.01/1000)*I*S),

Edit: I modified my code to:

    systemOfEquations = []
    with open("SIREquationsMathematica2.txt", "r") as fp :
        for line in fp :
               pattern = regex.compile(r'.+?\s+=\s+(.+?)$')
               expressionString = regex.search(pattern, line) 
               systemOfEquations.append( sympy.sympify( expressionString) )


    def dX_dt(X, t):
        vals = dict(S=X[0], I=X[1], R=X[2], t=t)
        return [eq for eq in systemOfEquations]

    out = odeint(dX_dt, [1000,1,0], np.linspace(0, 100, 50), )

and this works (I don't quite get what the first two lines of the for loop are doing). However, I would like to do the process of defining the variables more automatic, and I still don't know how to use this solution and pass parameters in a text file. Along the same lines, how can I define parameters (that will depend on the variables) inside the dX_dt function?

Thanks in advance!

Answer 1

This isn't a full answer, but rather some observations/questions, but they are too long for comments.

dX_dt is called many times by odeint with a 1d array y and tuple t. You provide t via the args parameter. y is generated by odeint and varies with each step. dX_dt should be streamlined so it runs fast.

Usually an expresion like [eq for eq in systemOfEquations] can be simplified to systemOfEquations. [eq for eq...] doesn't do anything meaningful. But there may be something about systemOfEquations that requires it.

I'd suggest you print out systemOfEquations (for this small 3 line case), both for your benefit and ours. You are using sympy to translated the strings from the file into equations. We need to see what it produces.

Note that myODEs is a function, not a file. It may be imported from a module, which of course is a file.

The point to vals = dict(S=X[0], I=X[1], R=X[2], t=t) is to produce a dictionary that the sympy expressions can work with. A more direct (and I think faster) dX_dt function would look like:

def myODEs(y, t, params):
    c,p, gamma = params
    beta = c*p
    dydt = [-beta*y[0]*y[1],
           beta*y[0]*y[1] - gamma*y[1],
           - gamma*y[1]]  
    return dydt

I suspect that the dX_dt that runs sympy generated expressions will be a lot slower than a 'hardcoded' one like this.

I'm going add sympy tag, because, as written, that is the key to translating your text file into a function that odeint can use.

I'd be inclined to put the equation variability in the t parameters, rather a list of sympy expressions.

That is replace:

    dydt = [-beta*y[0]*y[1],
           beta*y[0]*y[1] - gamma*y[1],
           - gamma*y[1]]  

with something like

    arg12=np.array([-beta, beta, 0])
    arg1 = np.array([0, -gamma, -gamma])
    arg0 = np.array([0,0,0])
    dydt = arg12*y[0]*y[1] + arg1*y[1] + arg0*y[0]

Once this is right, then the argxx definitions can be move outside dX_dt, and passed via args. Now dX_dt is just a simple, and fast, calculation.

This whole sympy approach may work fine, but I'm afraid that in practice it will be slow. But someone with more sympy experience may have other insights.