In [2]:

    

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate
%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.
    """
    if sens:
        n = X0.shape[0]
        x0 = np.hstack((X0, np.eye(n).flatten(), np.zeros(len(p) * n)))

        def g(t, X, p):
            var = X[:n]
            var2 = X[: n * 3]
            return np.hstack((f(t, var, p), dfx(t, var2, p), dfp(t, X, p)))
        ivp = scipy.integrate.ode(g)
    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)

    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]])
    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 not ivp.successful():
            result.success = False
            break

    if store_trajectory:
        result.trajectory_t = np.array(times)
        result.trajectory_x = np.array(points)
    return result
'''
def dX_dt(t, X, p):
    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.dot(np.array([alpha - beta * F, - beta * R, gamma * F,\
                            gamma * R -  delta]).reshape((2, 2)),\
                  X[2:].reshape((2, 2))).flatten()

def dfp(t, X, p):
    R = X[0]
    F = X[1]
    alpha, beta, gamma, delta = p
    foo = np.dot(np.array([alpha - beta * F, - beta * R, gamma * F,\
                           gamma * R -  delta]).reshape((2, 2)), \
                 X[6:].reshape((2, 4)))
    return (foo + np.array([[R, - F * R, 0., 0.],[0., 0.,  R * F,\
                                                  -F]])).flatten()


    

In [2]:

    

# initial condition
X0 = np.array([20.,10.])
param = (.2, .01, .001, .1)
result = solve_ivp(dX_dt, np.array([0, 300]), X0, p=param,\
                   integrator='dopri5', store_trajectory=True, sens=True, dfx=dfx, dfp=dfp)


    

In [3]:

    

xs = result.xs[-1, :]
x_ = xs[:2]
gx = xs[2:6].reshape((2, 2))
gp = xs[6:].reshape((2, 4))
print('The result  for the rabbits and foxes:')
print(x_)
print('G_x(t)')
print(gx)
print('G_p(t)')
print(gp)


    

    




The result  for the rabbits and foxes:
[ 12.64880398  19.2015986 ]
G_x(t)
[[ 0.72830829  1.69448833]
 [-2.47569061 -4.09255314]]
G_p(t)
[[ -3.47553958e+01   1.69448833e+03   1.91736179e+03   1.83915592e+02]
 [ -9.91575326e+02  -6.01271300e+03  -4.95138121e+04  -1.95234998e+03]]

In [4]:

    

def g(X, t=0):
    p = param
    var = X[:2]
    var2 = X[: 2 * 3]
    return np.hstack((dX_dt(t, var, p), dfx(t, var2, p), dfp(t, X, p)))


n = 2
X0 = np.hstack((np.array([20., 10.]), np.eye(2).flatten(), np.zeros(8)))
t = np.linspace(0,300,100000)
X, infodict = scipy.integrate.odeint(g, X0, t, full_output=True)
infodict['message']


    

    
Out[4]:





'Integration successful.'

In [5]:

    

xs = X[-1, :]
x__fixed = xs[:2]
gx_fixed = xs[2:6].reshape((2, 2))
gp_fixed = xs[6:].reshape((2, 4))
print('The result  for the rabbits and foxes:')
print(x__fixed)
print('G_x(t)')
print(gx_fixed)
print('G_p(t)')
print(gp_fixed)


    

    




The result  for the rabbits and foxes:
[ 12.64881061  19.20161368]
G_x(t)
[[ 0.72827082  1.69444497]
 [-2.47565671 -4.09253794]]
G_p(t)
[[ -3.47718101e+01   1.69444497e+03   1.91660572e+03   1.83890604e+02]
 [ -9.91559500e+02  -6.01269931e+03  -4.95131342e+04  -1.95233004e+03]]

In [6]:

    

print('\nDelta R, Delta F')
print(x_ - x__fixed)
print('\nDelta G(t)_x')
print(gx - gx_fixed)
print('\nDelta G(t)_p')
print(gp - gp_fixed)


    

    




Delta R, Delta F
[ -6.63744237e-06  -1.50769547e-05]

Delta G(t)_x
[[  3.74715911e-05   4.33563723e-05]
 [ -3.38975057e-05  -1.52020022e-05]]

Delta G(t)_p
[[ 0.01641431  0.04335637  0.75606926  0.02498857]
 [-0.01582585 -0.01369431 -0.67795011 -0.01993601]]

In [61]:

    

import scipy.integrate
import scipy.sparse
import numpy as np
import collections

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


    

In [62]:

    

def dX_dt(t, X, p):
    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]])


# initial condition
X0 = np.array([20.,10.])
param = np.array((.2, .01, .001, .1))
result = solve_ivp(dX_dt, np.array([0, 300]), X0, p=param,\
                   integrator='dopri5', store_trajectory=True, sens=True, dfx=dfx, dfp=dfp)


    

In [75]:

    

result.Gps[1].shape


    

    
Out[75]:





(2, 4)

In [140]:

    

def evaluate_parest_single_shooting(f, dfx, dfp, meas_times, meas_values, s0, p):
    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)
    F = (meas_values - result.xs).reshape((m * nx, 1))
    #F = np.linalg.norm(F)
    J = np.zeros((m * nx, nv))
    #print (J.shape)
    dy = np.eye(nx)
    for i in range(m):
        for j in range(nx):
            #print(2 * m + nx)
            J[2 * i + j, :] = np.hstack((dy[[j], :].dot(result.Gxs[i]), dy[[j], :].dot(result.Gps[i])))
    return F, J

#def gauss_newton(F_J, x0, itmax=100, tol=1e-7, verbose=1):


    

In [129]:

    

np.random.seed(31)
meas_times = np.arange(21) * 5
noise = 10. * np.random.rand(21, 2)
x0 = np.array([20., 10.])
result = solve_ivp(dX_dt, meas_times, X0, p=param)
meas_values = result.xs + noise


    

In [130]:

    

F, J= evaluate_parest_single_shooting(dX_dt, dfx, dfp, meas_times, meas_values, X0, param)


    

    




(42, 6)

In [184]:

    

np.random.seed(31)
meas_times = np.arange(21) * 5
noise = 20. * np.random.rand(21, 2)

x0_init = np.array([20., 10.])
param_init = np.array((.2, .01, .001, .1))
result = solve_ivp(dX_dt, meas_times, x0_init, p=param_init)
meas_values = result.xs + noise


x0 = x0_init * (1 + .3 * np.random.rand(len(x0)))
param = param_init * (1 + .02 * np.random.rand(len(param.shape))) 
tau = .5
for i in range(20):
    F, J= evaluate_parest_single_shooting(dX_dt, dfx, dfp, meas_times, meas_values, x0, param)
    Dx = np.linalg.lstsq(J, F)[0]
    #print (Dx.shape)
    x0 += tau * Dx.flatten()[:2]
    param += tau * Dx.flatten()[2:]
    print(np.linalg.norm(F))


    

    




183.272960207
93.6848438009
60.8908096508
50.7285669842
48.0687713419
47.4020148869
47.2275837077
47.1775962455
47.161169119
47.1547520877
47.1517863019
47.1502358752
47.1493648828
47.1488574209
47.1485566668
47.1483770482
47.1482694187
47.1482048398
47.1481660747
47.1481428041

In [185]:

    

print((param - param_init) / param_init)
print((x0 - x0_init)/ x0_init)
print(param)
print(param_init)


    

    




[ 0.05839425 -0.19220837 -0.14597882 -0.09251752]
[ 0.14365051  0.47141482]
[ 0.21167885  0.00807792  0.00085402  0.09074825]
[ 0.2    0.01   0.001  0.1  ]

In [150]:

    

Dx


    

    
Out[150]:





array([[  1.08952223e-11],
       [ -4.39763683e-11],
       [ -3.87416206e-13],
       [ -5.59538464e-15],
       [  1.07963928e-15],
       [  1.31417969e-13]])

In [70]:

    

meas_values.T.reshape((meas_values.size, 1))


    

    
Out[70]:





array([[  22.86053822],
       [  43.5660351 ],
       [  74.05288828],
       [ 158.77635974],
       [ 285.42187874],
       [ 313.49878537],
       [ 102.20730423],
       [  26.44792153],
       [  18.32924493],
       [  15.97273941],
       [  20.54387204],
       [  34.52449685],
       [  64.87795615],
       [ 127.66342944],
       [ 248.74670075],
       [ 347.09229066],
       [ 145.19726741],
       [  35.05759042],
       [  20.53381241],
       [  13.56541286],
       [  19.15769563],
       [  19.58105567],
       [  16.8045809 ],
       [   6.8121239 ],
       [   6.30084112],
       [  15.34704083],
       [  34.42365375],
       [  54.74501082],
       [  41.44363318],
       [  26.8760131 ],
       [  26.42033959],
       [  13.77331686],
       [  15.81696425],
       [   8.38235889],
       [  15.10177521],
       [  10.436421  ],
       [  26.65913688],
       [  56.13496698],
       [  45.55913752],
       [  29.97999634],
       [  19.83158892],
       [  21.45407851]])

In [68]:

    

meas_values.reshape(21, 2)


    

    
Out[68]:





array([[  22.86053822,   19.58105567],
       [  43.5660351 ,   16.8045809 ],
       [  74.05288828,    6.8121239 ],
       [ 158.77635974,    6.30084112],
       [ 285.42187874,   15.34704083],
       [ 313.49878537,   34.42365375],
       [ 102.20730423,   54.74501082],
       [  26.44792153,   41.44363318],
       [  18.32924493,   26.8760131 ],
       [  15.97273941,   26.42033959],
       [  20.54387204,   13.77331686],
       [  34.52449685,   15.81696425],
       [  64.87795615,    8.38235889],
       [ 127.66342944,   15.10177521],
       [ 248.74670075,   10.436421  ],
       [ 347.09229066,   26.65913688],
       [ 145.19726741,   56.13496698],
       [  35.05759042,   45.55913752],
       [  20.53381241,   29.97999634],
       [  13.56541286,   19.83158892],
       [  19.15769563,   21.45407851]])

In [69]:

    

meas_values


    

    
Out[69]:





array([[  22.86053822,   19.58105567],
       [  43.5660351 ,   16.8045809 ],
       [  74.05288828,    6.8121239 ],
       [ 158.77635974,    6.30084112],
       [ 285.42187874,   15.34704083],
       [ 313.49878537,   34.42365375],
       [ 102.20730423,   54.74501082],
       [  26.44792153,   41.44363318],
       [  18.32924493,   26.8760131 ],
       [  15.97273941,   26.42033959],
       [  20.54387204,   13.77331686],
       [  34.52449685,   15.81696425],
       [  64.87795615,    8.38235889],
       [ 127.66342944,   15.10177521],
       [ 248.74670075,   10.436421  ],
       [ 347.09229066,   26.65913688],
       [ 145.19726741,   56.13496698],
       [  35.05759042,   45.55913752],
       [  20.53381241,   29.97999634],
       [  13.56541286,   19.83158892],
       [  19.15769563,   21.45407851]])

In [71]:

    

v = np.arange(8)


    

In [73]:

    

v = 0


    

In [74]:

    

v


    

    
Out[74]:





0

In [ ]: