Используйте solve_ivp с различными параметрами - PullRequest
Новая
       50

Используйте solve_ivp с различными параметрами

0 голосов
/ 10 мая 2018

Я пытаюсь решить следующую простую линейную систему уравнений:

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

где x - это n-матричный вектор, а A_eps (t) - матрица, которая меняется со временем.

Вот что я пробовал до сих пор после поста:

scipy ode update set_f_params внутри функции, установленной как set_solout

Сначала я определил правую часть как функцию: 

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)

Затем я просто настроил некоторые параметры для 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)

Для настройки решателя я сделал: 

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

Я хотел бы написать что-то вроде: 

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)

где next_enviroment (current_env, Q): 

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)

- это функция, которая принимает текущее состояние и выбирает случайное новое состояние в соответствии с заданным правилом. Моя проблема сложена вдвое:

  • Во-первых, я не знаю, как читать текущую среду в решателе.
  • Во-вторых, учитывая текущую среду, как передать ее в функцию.

Спасибо за помощь.

1 Ответ

0 голосов
/ 10 мая 2018

Я нашел ответ. Вот код: 

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)

Вот график моделирования:

Эволюция двух фенотипов в случайной случайной среде с двумя состояниями. Штат j благоприятствует фенотипу j

...