Skip to content

Commit

Permalink
merged upstream changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@ayush9pandey
ayush9pandey committed Mar 23, 2020
2 parents 8c5abc8 + e238583 commit 9d9d879
Show file tree
Hide file tree
Showing 31 changed files with 40,075 additions and 894 deletions.
98 changes: 56 additions & 42 deletions auto_reduce/local_sensitivity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import numpy as np
# from .ode import ODE
from scipy.integrate import solve_ivp, odeint
from sympy import lambdify

class SSM(System):
'''
Expand All @@ -13,16 +14,12 @@ class SSM(System):
order central difference as given in the paper.
'''
def __init__(self, x, f, params = None, C = None, g = None, h = None, u = None,
params_values = [], x_init = [], timepoints = None, mode = None):
params_values = [], x_init = [], timepoints = None):
super().__init__(x, f, params, C, g, h, u, params_values, x_init)
if timepoints is None:
timepoints = []
else:
self.timepoints = timepoints
if mode is None:
self.mode = 'approximate'
else:
self.mode = mode

return

Expand All @@ -32,8 +29,10 @@ def compute_Zj(self, x, j, **kwargs):
Returns a vector of size n x 1.
Use mode = 'accurate' for this object attribute to use accurate computations using numdifftools.
'''
if self.mode == 'accurate':
return self.sensitivity_to_parameter(**kwargs)
# if 'mode' in kwargs:
# if kwargs.get('mode') == 'accurate':
# del kwargs['mode']
# return self.sensitivity_to_parameter(x, j, **kwargs)
# initialize Z
Z = np.zeros(self.n)
P_holder = self.params_values
Expand Down Expand Up @@ -71,37 +70,53 @@ def compute_J(self, x, **kwargs):
Returns a matrix of size n x n.
Use mode = 'accurate' for this object attribute to use accurate computations using numdifftools.
'''
if self.mode == 'accurate':
import numdifftools as nd
ode = kwargs.get('ode')
params = kwargs.get('params')
ode_func = lambda t, xs: ode(t, xs, params)
return nd.jacobian(lambda x: ode_func(0, x))
if 'fun' in kwargs:
fun = kwargs.get('fun')
else:
fun = self.f
if 'var' in kwargs:
var = kwargs.get('var')
else:
var = x
# initialize J
J = np.zeros( (self.n, self.n) )
J = np.zeros( (self.n, len(var)) )
P = self.params_values
# store x
X = x
u = self.u
if 'mode' in kwargs:
if kwargs.get('mode') == 'accurate':
del kwargs['mode']
try:
import numdifftools as nd
except:
raise ValueError('The package numdifftools is not installed for this method to work.')
fun_l = lambdify((self.x, self.params), fun)
def fun_ode(t, x, params):
y = fun_l(x, params)
return np.array(y)
jfun = nd.Jacobian(lambda x: fun_ode(0, x, P), **kwargs)
return jfun(x)
# store the variable with respect to which we approximate the differentiation (df/dvar)
X = var
for i in range(self.n):
for j in range(self.n):
for j in range(len(var)):
F = np.zeros( (4,1) )
h = X[j]*0.01
# Gets O(4) central difference on dfi/dxj
# Gets O(4) central difference on dfi/dvarj
if h != 0:
x = X
x[j] = X[j] + 2*h
f = self.evaluate(self.f, x, P)
var = X
var[j] = X[j] + 2*h
f = self.evaluate(fun, var, P, u)
F[0] = f[i]
x[j] = X[j] + h
f = self.evaluate(self.f, x, P)
var[j] = X[j] + h
f = self.evaluate(fun, var, P, u)
F[1] = f[i]
x[j] = X[j] - h
f = self.evaluate(self.f, x, P)
var[j] = X[j] - h
f = self.evaluate(fun, var, P, u)
F[2] = f[i]
x[j] = X[j] - 2*h
f = self.evaluate(self.f, x, P)
var[j] = X[j] - 2*h
f = self.evaluate(fun, var, P, u)
F[3] = f[i]
#Store approx. dfi/dxj into J
#Store approvar. dfi/dvarj into J
J[i,j]= (-F[0] + 8*F[1] - 8*F[2] + F[3])/(12*h)
# print(J[i,j])
# if J[i,j] == np.Inf:
Expand All @@ -117,8 +132,9 @@ def compute_SSM(self, normalize = False, **kwargs):
If normalize argument is true, the coefficients are normalized by the nominal value of each paramneter.
Use mode = 'accurate' for this object attribute to use accurate computations using numdifftools.
'''
if self.mode == 'accurate':
return self.solve_extended_ode(**kwargs)
if 'mode' in kwargs:
if kwargs.get('mode') == 'accurate_SSM':
return self.solve_extended_ode(**kwargs)

def sens_func(t, x, J, Z):
# forms ODE to solve for sensitivity coefficient S
Expand All @@ -140,13 +156,13 @@ def sens_func(t, x, J, Z):
if len(timepoints) == 1:
continue
# get the jacobian matrix
J = self.compute_J(xs[k,:])
J = self.compute_J(xs[k,:], **kwargs)
#Solve for S = dx/dp for all x and all P (or theta, the parameters) at time point k
for j in range(len(P)):
utils.printProgressBar(int(j + k*len(P)), len(self.timepoints)*len(P) - 1, prefix = 'SSM Progress:', suffix = 'Complete', length = 50)
# print('for parameter',P[j])
# get the pmatrix
Zj = self.compute_Zj(xs[k,:], j)
Zj = self.compute_Zj(xs[k,:], j, **kwargs)
# solve for S
sens_func_ode = lambda t, x : sens_func(t, x, J, Zj)
sol = odeint(sens_func_ode, S0, timepoints, tfirst = True)
Expand Down Expand Up @@ -179,13 +195,12 @@ def get_system(self):

############## Sam Clamons ###########

def sensitivity_to_parameter(self, ode_sol = None, ode_jac = None, ode = None, params = None, n_vars = None, p = None, t_min = None,
t_max = None):
def sensitivity_to_parameter(self, x, j, **kwargs):
'''
Calculates the response of each derivative (defined by ode) to changes
in a single parameter (the jth one) over time.
in a single parameter (the jth one) at point x.
params:
keyword argument options?:
ode_sol - An OdeSolution object holding a continuously-interpolated
solution for ode.
ode_jac - Jacobian of the ode, as calculated by numdifftools.Jacobian.
Expand All @@ -205,14 +220,13 @@ def dS_dt(t, s):
xs = ode_sol(t)
# Wrapper to let numdifftools calculate df/dp.
def ode_as_parameter_call(param):
call_params = copy.deepcopy(params)
call_params[p] = param
call_params = copy.deepcopy(self.params_values)
call_params[j] = self.params_values
return ode(t, xs, call_params)
df_dp = lambda xs: nd.Jacobian(ode_as_parameter_call)(xs).transpose()[:,0]
return df_dp(params[p]) + np.matmul(ode_jac(xs), s)
sol = solve_ivp(dS_dt, (t_min, t_max), np.zeros(n_vars),
dense_output = True)
return sol.sol
return df_dp(params[j]) + np.matmul(ode_jac(xs), s)
sol = odeint(dS_dt, np.zeros(n_vars), self.timepoints, **kwargs)
return sol


############## Sam Clamons ###########
Expand Down

0 comments on commit 9d9d879

Please sign in to comment.