CODEWITHSUNDEEP
pythonscipyscientific-computing
Zee
Zee

Reputation: 41

Use solve_ivp with varying parameters

I am trying to solve the following simple linear system of equations:

x'(t) = A_eps(t) x(t)

where x is an n-diemnsional vector and A_eps(t) is a matrix that changes over time.

Here is what I have tried so far following the post:

scipy ode update set_f_params inside function set as set_solout

First I defined the right hand side as a function: 

from functools import partial
from scipy.integrate import solve_ivp

def Aeps_x(t, x, enviroment, Q, F, H):

    '''
    inputs: 
        - t = time.
        - x[j] = # of cell with phenotype j in an enviroment. 
        - enviroment = number of enviroment. There are E in total.  
        - Q = ExE transition matrix of a Markov Chain.  
        - F = ExP matrix with growth rate of phenotypes for each enviroment. 
        - H = PxP matrix with the switching rates from phenotype i to j.
    '''
    E = Q.shape[0]
    P = F.shape[1]

    A = np.diag(F[enviroment]) - H

    dx = A.dot(x)

    return(dx)

Then, I just set up some parameters for the r.h.s.: 

EMat = np.array([[0, 1, 1, 1], 
             [1, 0, 1, 1],
             [1, 1, 0, 1],
             [1, 1, 1, 0]])
E0 = EMat.shape[0]

row_sums = EMat.sum(axis=1).reshape(E0, 1)
Q = EMat / row_sums

F = linalg.toeplitz([1, 0.1, 0.1, 0.1]) # only one strong phenotype for each 
enviroment
F = F[0:E0, ] # only 4 enviroments

P0 = F.shape[1]
H = np.diag([0.5]*P0)

To set up the solver, I did: 

sol_param = solve_ivp(
    partial(Aeps_x, enviroment = 2, Q=Q, F=F, H=H), (0, 4, 20), x0, dense_output=True)

I would like to write something like: 

sol_param = solve_ivp(
    partial(Aeps_x, enviroment = next_enviroment(current_env, Q), 
    Q=Q, F=F, H=H), (0, 4, 20), x0, dense_output=True)

where next_enviroment(current_env, Q) is: 

def next_enviroment(current_env, Q):
    '''
    Given an current state, computes the next state in a markov chain with transition matrix Q. 
    '''
    sample = np.random.multinomial(1, Q[intitial_env], size = 1)
    next_env = np.argmax(sample)
    return(next_env)

is a function that takes the current state and choose a random new state according to a given rule. My problem is two folded:

Thanks for the help.

Upvotes: 3

Views: 2216

Answers (1)

Zee
Zee

Reputation: 41

I found an answer. Here's the code: 

EMat = np.array([[10, 0.1], 
                 [0.1, 10]])

# number of enviroments and phenotypes
E = 2
P = 2

row_sums = EMat.sum(axis=1).reshape(E, 1)
Q = EMat / row_sums

H = np.array([[0, 0.05],
              [1e-6, 0]])

F = np.array([[2, -0.05], 
              [-7, -0.05]]) 


import scipy

N = 1000
tfinal = 25
t = np.linspace(start=0, stop=tfinal, num=N)
t0 = 0
x0 = [20, 20]
e0 = 0


solver = scipy.integrate.ode(Aeps_x).set_integrator('dopri5', nsteps=100)
solver.set_initial_value(x0, t0).set_f_params(e0, Q, F, H)

sol = np.zeros((N, E))
sol[0] = x0
Enviroments = np.zeros(N)

k = 1
while solver.successful() and solver.t < tfinal:
    solver.integrate(t[k])
    sol[k] = solver.y
    k += 1

    Enviroments[k-1] = e0

    e0 =  next_enviroment(e0, Q=Q)
    solver.set_f_params(e0, Q, F, H)

Here's a pplot of the simulation:

Evolution of two phenotypes in an stochastic random enviroment with two states. State j favors phenotype j

Upvotes: 1

Related Questions