In [1]:
#simulation to perform many runs of the healthy vs. unhealthy dynamics
#this version will start with an initial condition in the diseased state space with an
#initial covariance that is nonzero.
import numpy as np
import matplotlib.pyplot as plt
import math
import sys
from threechoice_dynamics import *
#tf.set_random_seed(101)
pi = math.pi
#np.random.seed(101)
#define the gradients and parameters of the simulation
xdim = 2 #state dimension
udim = 2 #control dimension
T = 10000 #number of steps
#define the noise for the system
dt = 1.0e-4
gamma = 1.0e-4
sigsstate = (1./dt)*(1e-9)#fix these
sigsobs = 1.0e-6
Q = sigsstate*np.eye(xdim)
R = sigsobs*np.eye(xdim)
#the non-tensorflow anonymous functions, for generalizations
true_nontf = lambda x,c: grad_threechoice(x[0],x[1],dt,c)
target_nontf = lambda x: grad_threechoice_healthy(x[0],x[1],dt)
In [8]:
ns = 500 #number of samples
#------- the main simulation loop. I might make this a function if i use it a lot
#make these numpy version
statenoise = np.random.normal(0,sigsstate**0.5,[xdim,T,ns])
obsnoise = np.random.normal(0,sigsobs**0.5,[xdim,T,ns])
G = dt**(0.5)*np.eye(xdim) #system noise matrix, for covariance prediction
xvec = np.zeros((xdim,T,ns))
yvec = np.zeros((xdim,T,ns))
loss = np.zeros((T,ns))
lag = 0
ci = [0.,0.]
for m in range(ns):
x_init = np.random.uniform(0.,0.5,(2,))
x_k = x_init
xvec[:,0,m] = x_init
#print(m)
#go ahead and propagate lag-steps ahead before starting state estimation and such
for k in range(1,T):
#update actual dynamics
grad_cont = true_nontf(xvec[:,k-1,m],ci)
xvec[:,k,m] = xvec[:,k-1,m] + grad_cont + statenoise[:,k,m]
yvec[:,k,m] = xvec[:,k,m] + obsnoise[:,k,m]
loss[k,m] = np.linalg.norm(
true_nontf(xvec[:,k,m],ci)-
target_nontf(xvec[:,k,m]))
In [9]:
#------------ save
import pickle
fname = 'threechoice_initruns_nocont_T10000.p'
index = ['xvec','yvec','loss']
alldata = [index,xvec,yvec,loss]
pickle.dump( alldata, open( fname, "wb" ) )
In [15]:
#set up and run EKF in tensorflow
import tensorflow as tf
from tf_funs import *
X_est = tf.placeholder(shape=xdim,dtype=tf.float32,name='X_est') #the state estimate
PI_est = tf.placeholder(shape = (xdim,xdim),dtype=tf.float32, name = 'PI_est') #estimated covariance
Y_tp1 = tf.placeholder(shape= xdim,dtype=tf.float32, name = 'Y_tp1') #the most recent observation
#now run EKF for last 100 steps and get the state estimates for the final
Q_tf = tf.constant(Q,dtype=tf.float32, name = 'Q') #state noise covariance
R_tf = tf.constant(R,dtype=tf.float32, name = 'R') #observation noise covariance
Control = tf.placeholder(shape = udim, dtype=tf.float32, name='Control')
#graphs for updating state and observation
true_model_est = grad_threechoice_tf(X_est[0],X_est[1],dt,Control)
true_model_est_null = grad_threechoice_tf(X_est[0],X_est[1],dt,[0.,0.]) #state est. gradient null control
target_model_est = grad_threechoice_healthy_tf(X_est[0],X_est[1],dt) #state est. target dynamics
#the non-tensorflow anonymous functions, for generalizations
true_nontf = lambda x,c: grad_threechoice(x[0],x[1],dt,c)
target_nontf = lambda x: grad_threechoice_healthy(x[0],x[1],dt)
X_plus,PI_plus = EKF(X_est,Y_tp1,PI_est,true_model_est,true_model_est_null,Q_tf,R_tf,xdim,dt)
In [19]:
session_conf = tf.ConfigProto(
intra_op_parallelism_threads=4,
inter_op_parallelism_threads=4)
sess = tf.Session(config=session_conf)
init = tf.global_variables_initializer()
sess.run(init)
t_start = T-100
x_estvec = np.zeros((xdim,T-t_start,ns))
PI_estvec = np.zeros((xdim,xdim,T-t_start,ns))
xvec_init = xvec[:,t_start:T,:] #values of x that we are running EKF on
x_inits = xvec[:,t_start,:]
PI_init = [[1.0e-4,0.],[0.,1.0e-4]] #initial covariance
x_estvec[:,0,:] = x_inits
for k in range(ns):
PI_estvec[:,:,0,k] = PI_init
print('starting estimation')
for m in range(ns):
for k in range(1,T-t_start):
test = sess.run([X_plus,PI_plus],
{X_est:x_estvec[:,k-1,m], PI_est:PI_estvec[:,:,k-1,m],
Control:ci, Y_tp1:yvec[:,t_start + k,m]})
x_estvec[:,k,m] = test[0]
PI_estvec[:,:,k,m] = test[1]
#print(m)
In [23]:
import pickle
xvec_init = xvec[:,t_start:T,:] #values of x that we are running EKF on
fname = 'threechoice_initruns_nocont_ekf.p'
index = ['x_estvec','PI_estvec','xvec_init']
alldata = [index,x_estvec,PI_estvec,xvec_init]
pickle.dump( alldata, open( fname, "wb" ) )
In [26]:
#run the whole last 1000 steps to get state esitmate and loss
loss_tf = loss_full(X_est,PI_est,true_model_est,target_model_est)
session_conf = tf.ConfigProto(
intra_op_parallelism_threads=4,
inter_op_parallelism_threads=4)
sess = tf.Session(config=session_conf)
init = tf.global_variables_initializer()
sess.run(init)
t_start = T-1000
x_estvec = np.zeros((xdim,T-t_start,ns))
PI_estvec = np.zeros((xdim,xdim,T-t_start,ns))
xvec_init = xvec[:,t_start:T,:] #values of x that we are running EKF on
loss = np.zeros((4,T-t_start,ns))
x_inits = xvec[:,t_start,:]
PI_init = [[1.0e-4,0.],[0.,1.0e-4]] #initial covariance
x_estvec[:,0,:] = x_inits
for k in range(ns):
PI_estvec[:,:,0,k] = PI_init
print('starting estimation')
for m in range(ns):
print(m)
for k in range(1,T-t_start):
test = sess.run([X_plus,PI_plus],
{X_est:x_estvec[:,k-1,m], PI_est:PI_estvec[:,:,k-1,m],
Control:ci, Y_tp1:yvec[:,t_start + k,m]})
x_estvec[:,k,m] = test[0]
PI_estvec[:,:,k,m] = test[1]
ltest = sess.run(loss_tf,{X_est:x_estvec[:,k-1,m],
PI_est:PI_estvec[:,:,k-1,m],
Control:ci
})
loss[:,k-lag,m] = ltest
#print(m)
In [27]:
#save
xvec_init = xvec[:,t_start:T,:] #values of x that we are running EKF on
fname = 'threechoice_initruns_nocont_ekf_T1000.p'
index = ['x_estvec','PI_estvec','xvec_init','loss']
alldata = [index,x_estvec,PI_estvec,xvec_init,loss]
pickle.dump( alldata, open( fname, "wb" ) )