In [19]:
import scipy.integrate
import scipy.sparse
import numpy as np
import collections
import matplotlib.pyplot as plt
%matplotlib inline
class IVPResult:
pass
def solve_ivp(f, ts, x0, p=None, integrator='dopri5', store_trajectory=False,
sens=False, dfx=None, dfp=None):
"""
Solve initial value problem
d/dt x = f(t, x, p); x(t0) = x0.
Evaluate Solution at time points specified in ts, with ts[0] = t0.
Compute sensitivity matrices evaluated at time points in ts if sens=True.
"""
if sens:
def f_vode(t, xGxGp, p):
x = xGxGp[:x0.shape[0]]
Gx = xGxGp[x0.shape[0]:x0.shape[0]+x0.shape[0]*x0.shape[0]]
Gx = Gx.reshape([x0.shape[0], x0.shape[0]])
Gp = xGxGp[x0.shape[0]+x0.shape[0]*x0.shape[0]:]
Gp = Gp.reshape([x0.shape[0], p.shape[0]])
dx = f(t, x, p)
dfxev = dfx(t, x, p)
dfpev = dfp(t, x, p)
dGx = dfxev.dot(Gx)
dGp = dfxev.dot(Gp) + dfpev
return np.concatenate(
[dx,
dGx.reshape(-1),
dGp.reshape(-1)])
ivp = scipy.integrate.ode(f_vode)
else:
ivp = scipy.integrate.ode(f)
ivp.set_integrator(integrator)
if store_trajectory:
times = []
points = []
def solout(t, x):
if len(times) == 0 or t != times[-1]:
times.append(t)
points.append(np.copy(x[:x0.shape[0]]))
ivp.set_solout(solout)
if sens:
ivp.set_initial_value(np.concatenate(
[x0,
np.eye(x0.shape[0]).reshape(-1),
np.zeros([x0.shape[0], p.shape[0]]).reshape(-1)
]), ts[0])
else:
ivp.set_initial_value(x0, ts[0])
ivp.set_f_params(p)
result = IVPResult()
result.ts = ts
result.xs = np.zeros([ts.shape[0], x0.shape[0]])
if sens:
result.Gxs = np.zeros([ts.shape[0], x0.shape[0], x0.shape[0]])
result.Gps = np.zeros([ts.shape[0], x0.shape[0], p.shape[0]])
result.success = True
result.xs[0,:] = x0
for ii in range(1,ts.shape[0]):
ivp.integrate(ts[ii])
result.xs[ii,:] = ivp.y[:x0.shape[0]]
if sens:
result.Gxs[ii,:,:] = ivp.y[x0.shape[0]:x0.shape[0]+x0.shape[0]*x0.shape[0]].reshape([x0.shape[0], x0.shape[0]])
result.Gps[ii,:,:] = ivp.y[x0.shape[0]+x0.shape[0]*x0.shape[0]:].reshape([x0.shape[0], p.shape[0]])
if not ivp.successful():
result.success = False
break
if store_trajectory:
result.trajectory_t = np.array(times)
result.trajectory_x = np.array(points)
return result
class MSResult:
pass
# delta x should be a vector of dimension equal to the number
# of variables, which is true only if we do not change the parameters.
def evaluate_parest_single_shooting(f, dfx, dfp, meas_times, meas_values, s, p):
m, nx = s.shape
x_traject = []
t_traject = []
for l in range(m - 1):
interval = np.array((meas_times[l], meas_times[l + 1]))
result = solve_ivp(f, interval, s[l], p=p, sens=True, dfx=dfx, dfp=dfp, store_trajectory=True)
print(result.Gps.shape)
x_traject.append(result.trajectory_x)
t_traject.append(result.trajectory_t)
multi_result = MSResult()
multi_result.F1 = (meas_values - s).reshape((m * nx, 1))
multi_result.t_traject = t_traject
multi_result.x_traject = x_traject
multi_result.J1 = np.hstack((np.tile(np.eye(nx), m), np.zeros((nx, len(p)))))
return multi_result
'''
result = solve_ivp(f, meas_times, s0, p=p, sens=True, dfx=dfx, dfp=dfp)
m = len(meas_times)
nx = len(s0)
nv = nx + len(p)
# make F a column vector with [y11, y12, y21, y22, y31, y32, ... ym1, ym2]
F = (meas_values - result.xs).reshape((m * nx, 1))
# each row of J will be [dFij / dx01, dFij / dx02, dFij / dp1, ..., dFij / dp6 ]
J = np.zeros((m * nx, nv))
dy = np.eye(nx)
for i in range(m):
for j in range(nx):
# dhi / dy = dyi / dy = all zeros except for the ith element
# that will be 1
# For dFij / dp I dhi / dp = dyi(tj, yi, p) / dp is equal to 0
J[2 * i + j, :] = np.hstack((dy[[j], :].dot(result.Gxs[i]), dy[[j], :].dot(result.Gps[i])))
return F, J
'''
def dX_dt(t, X, p): #predator-prey DEs
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
dRdt = alpha * R - beta * F * R
dFdt = gamma * R * F - delta * F
return np.array([dRdt,dFdt])
def dfx(t, X, p):
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
return np.array([alpha - beta * F, - beta * R, gamma * F,\
gamma * R - delta]).reshape((2, 2))
def dfp(t, X, p):
R = X[0]
F = X[1]
alpha, beta, gamma, delta = p
return np.array([[R, - F * R, 0., 0.],[0., 0., R * F, -F]])
In [20]:
np.random.seed(32)
meas_times = np.arange(21) * 5
noise = 10. * np.random.rand(21, 2)
s0_init = np.array([20., 10.])
param_init = np.array((.2, .01, .001, .1))
result = solve_ivp(dX_dt, meas_times, s0_init, p=param_init)
meas_values = result.xs + noise
In [21]:
multi_result = evaluate_parest_single_shooting(dX_dt, dfx, dfp, meas_times, meas_values, meas_values, param_init)
In [22]:
t = multi_result.t_traject
x = multi_result.x_traject
plt.plot(np.hstack(t), np.vstack(x))
Out[22]:
In [22]:
multi_result.J1
Out[22]:
In [30]:
t[0].shape
Out[30]:
In [15]:
meas_values.shape
Out[15]:
In [5]:
np.kron(np.eye(2), np.ones(3))
Out[5]:
In [11]:
np.tile(np.eye(2),21),
Out[11]:
In [12]:
result.Gp
In [ ]: