In [1]:
# the tensorflow version of the wong/wang dynamics that only uses a TF version of the controller
#the rest is in numpy
import tensorflow as tf
import numpy as np
import matplotlib.pyplot as plt
import math
from tf_funs import * #EKF, and helper in tensorflow
from wang_dynamics import * #the gradients for wong/wang dynamics
import sys
tf.set_random_seed(101)
pi = math.pi
np.random.seed(101)
sess = tf.Session()
In [49]:
#build graphs
xdim = 2 #state dimension
udim = 2 #control dimension
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
#Q = tf.placeholder(dtype=tf.float32)
#R = tf.placeholder(dtype=tf.float32)
Control = tf.placeholder(shape = udim, dtype=tf.float32, name='Control')
#params for experiment
#wong dynamics
T = 1000 #number of steps
mu0 = 30. #stimulus strength
coh = -6.0 #coherence
gamma = 1e-4 #regularization term
#koulakov dynamics
#set up analogous stimulus to koulakov line
ifact = 7.5e-4 #scaling factor, found by inspection
I_k = [0,0]
I_kbase = ifact*(1+abs(coh)/100.0);
if coh < 0:
I_k = [-I_kbase,I_kbase]
elif coh > 0:
I_k = [I_kbase,-I_kbase]
#define the noise for the system
dt = 1.0e-3
sigsstate = (1./dt)*9e-9 #not sure how strong this should be yet
sigsobs = 1.0e-9 #good around 0.5
Q = sigsstate*np.eye(xdim)
Q_tf = tf.constant(Q,dtype=tf.float32, name = 'Q') #state noise covariance
R = sigsobs*np.eye(xdim)
R_tf = tf.constant(R,dtype=tf.float32, name = 'R') #observation noise covariance
#graphs for updating state and observation
true_model_est = wonggrad(X_est,Control,mu0,coh) #state est. gradient, full myopic
true_model_est_null = wonggrad(X_est,[0.,0.],mu0,coh)#state est. gradient null control
target_model_est = koulakov_line(X_est,I_k) #state est. target dynamics
#the non-tensorflow anonymous functions, for generalizations
true_nontf = lambda x,c: wonggrad_nontf(x,c,mu0,coh)
target_nontf = lambda x: koulakov_line_nontf(x,I_k)
X_plus,PI_plus = EKF(X_est,Y_tp1,PI_est,true_model_est,true_model_est_null,Q_tf,R_tf,xdim,dt)
#myopic controller
#graphs for controller
useMO = 0 #handle to use mean only. If ==1, then use mean-only. otherwise full
if useMO ==1:
print('using mean-only control')
Cnew = myopicController_meanonly(
X_est,PI_est,Control,gamma,true_model_est,
true_model_est_null,target_model_est,xdim,udim)
else:
print('using full myopic control')
Cnew = myopicController(
X_est,PI_est,Control,gamma,true_model_est,
true_model_est_null,target_model_est,xdim,udim)
#covariance prediction update graph
Ak = dynamics_linearized(X_est,true_model_est_null,xdim)
In [50]:
ns = 500 #number of samples
#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
x_estvec = np.zeros((xdim,T,ns))
xvec = np.zeros((xdim,T,ns))
yvec = np.zeros((xdim,T,ns))
x_targvec = np.zeros((xdim,T,ns))
PI_estvec = np.zeros((xdim,xdim,T,ns))
contall = np.zeros((udim,T,ns))
loss = np.zeros((T,ns))
loss_nocont = np.zeros((T,ns))
lag = 0 #how many steps in the past will we receive observations
init = tf.global_variables_initializer()
sess.run(init)
#failure happened at iteration 76
x_init = np.array([ 0.79408495, 0.20955779])
#load the noise instances that caused the crash
import pickle
fname = 'noise_bad_lag10_wang'
alldata = pickle.load( open( fname, "rb" ) )
statenoise = alldata[1]
obsnoise = alldata[2]
for m in [76]:
#x_init = [0.9,0.76] #initial state
#x_init = np.random.uniform(0.,1.,(2,))
PI_init = [[0.,0.],[0.,0.]] #initial covariance
c_init = [0.,0.]
xest_k = x_init
pi_k = PI_init
c_k = c_init
x_k = x_init
x_targ_k = x_init
ykp1 = np.array(x_init)
x_estvec[:,0,m] = x_init
xvec[:,0,m] = x_init
x_targvec[:,0,m] = x_init
PI_estvec[:,:,0,m] = PI_init
print(m)
#go ahead and propagate lag-steps ahead before starting state estimation and such
for k in range(1,lag):
#update actual dynamics
grad_cont = true_nontf(xvec[:,k-1,m],c_init)
grad_targ = target_nontf(x_targvec[:,k-1,m])
xvec[:,k,m] = xvec[:,k-1,m] + grad_cont + statenoise[:,k,m]
x_targvec[:,k,m] = x_targvec[:,k-1,m] + grad_targ + statenoise[:,k,m]
yvec[:,k,m] = xvec[:,k,m] + obsnoise[:,k,m]
#set estimates in beginning lags to initial state
x_estvec[:,k,m] = x_init
PI_estvec[:,:,k,m] = PI_init
for k in range(max(1,lag),T):
#update actual dynamics
grad_cont = true_nontf(xvec[:,k-1,m],contall[:,k-1,m])
grad_targ = target_nontf(x_targvec[:,k-1,m])
xvec[:,k,m] = xvec[:,k-1,m] + grad_cont + statenoise[:,k,m]
x_targvec[:,k,m] = x_targvec[:,k-1,m] + grad_targ + statenoise[:,k,m]
yvec[:,k,m] = xvec[:,k,m] + obsnoise[:,k,m]
#run state estimator to update estimate of state k-lag
#if k==0:
# print([k,x_estvec[:,k-lag-1,m], PI_estvec[:,:,k-lag-1,m],
# contall[:,k-lag-1,m], yvec[:,k-lag,m]])
test = sess.run([X_plus,PI_plus],
{X_est:x_estvec[:,k-lag-1,m], PI_est:PI_estvec[:,:,k-lag-1,m],
Control:contall[:,k-lag-1,m], Y_tp1:yvec[:,k-lag,m]})
x_estvec[:,k-lag,m] = test[0]
PI_estvec[:,:,k-lag,m] = test[1]
#predit lag states in the future to calculate control
x_est_n = x_estvec[:,k-lag,m]
PI_est_n = PI_estvec[:,:,k-lag,m]
for n in range(lag):
#state prediction step
grad_cont = true_nontf(x_est_n,contall[:,k-lag-1+n,m])
#covariance prediction step. calculate jacobian
Ak_n= sess.run(Ak,
{X_est: x_est_n, PI_est: PI_est_n,
Control: contall[:,0,m], Y_tp1:yvec[:,0,m]})
x_est_n = x_est_n + grad_cont
PI_est_n = np.matmul(Ak_n,PI_est_n) + np.matmul(PI_est_n,np.transpose(Ak_n)) + np.matmul(
np.matmul(G,Q),np.transpose(G))
#run myopic controller using predicted state estimated. cov, doesnt matter
#find control for time k
c_k = sess.run(Cnew,{X_est:x_est_n, PI_est:PI_est_n,
Control:contall[:,k-1,m], Y_tp1:yvec[:,k,m]})
#do quick comparison for jumps in control due to a singularity in dynamics
if abs(np.linalg.norm(c_k)) > 1000.:
contall[:,k,m] = contall[:,k-1,m]
print('dynamics got singular. hold tight')
print(k)
else:
contall[:,k,m] = c_k
loss[k-lag,m] = np.linalg.norm(true_nontf(x_estvec[:,k-lag,m],contall[:,k-lag,m])-
target_nontf(x_estvec[:,k-lag,m]))
loss_nocont[k-lag,m] = np.linalg.norm(
true_nontf(x_estvec[:,k-lag,m],[0.,0.])-
target_nontf(x_estvec[:,k-lag,m]))
#set final lag estimate values to esimate
for k in range(lag):
x_targvec[:,T-lag+k,m] = x_targvec[:,T-lag-1,m]
x_estvec[:,T-lag+k,m] = x_estvec[:,T-lag-1,m]
PI_estvec[:,:,T-lag+k,m] = PI_estvec[:,:,T-lag-1,m]
In [48]:
#check some values to see what made the control go nuts
ind = 362
m = 76
#Ak = dynamics_linearized(X_est,true_model_est_null,xdim)
print(sess.run(Ak,{X_est:x_estvec[:,ind,m], PI_est:PI_estvec[:,:,ind,m],
Control:contall[:,ind,m], Y_tp1:yvec[:,ind,m]}))
test = sess.run([X_plus,PI_plus],
{X_est:x_estvec[:,ind,m], PI_est:PI_estvec[:,:,ind,m],
Control:contall[:,ind,m], Y_tp1:yvec[:,ind,m]})
c_k = sess.run(Cnew,{X_est:x_est_n, PI_est:PI_est_n,
Control:contall[:,ind,m], Y_tp1:yvec[:,ind,m]})
print(test[0])
print(test[1])
print(c_k)
In [36]:
contall[:,ind,m]
Out[36]:
In [52]:
#true state
tind = 76
plt.plot(x_estvec[0,:,tind],x_estvec[1,:,tind],label='controlled dynamics')
plt.plot(x_targvec[0,:,tind],x_targvec[1,:,tind],label = 'target dynamics')
plt.legend()
plt.xlabel('S_1')
plt.ylabel('S_2')
plt.title('controlled trajecotory')
plt.show()
plt.plot(x_targvec[0,:,tind],label='x0')
plt.plot(x_targvec[1,:,tind],label = 'x1')
plt.legend()
plt.show()
#plt.plot(test[2][:,0]+ (test[3][:,0,0])**0.5)
#plt.plot(test[2][:,0]- (test[3][:,0,0])**0.5)
#plt.show()
plt.plot(contall[0,1:380,tind],label='control on S_1, full myopic')
plt.plot(contall[1,1:380,tind],label='control on S_2, full myopic')
plt.ylabel('C(t)')
plt.title('control signals')
plt.legend()
plt.show()
#expected state
fullloss = plt.plot(np.cumsum(loss[:,tind]),label='lagged, full myopic control')
nocontloss = plt.plot(np.cumsum(loss_nocont[:,tind]),label='no control')
#nolagloss = plt.plot(np.cumsum(loss_nolag[:,tind]),label='no lag,full myopic')
#loss_mo = plt.plot(loss[:,tind],label='lagged, mean-only myopic')
plt.xlabel('time')
plt.ylabel('loss')
plt.title('cumulative loss over time, sum_i(loss_i)')
plt.legend()
plt.show()
#standard deviation
#plot the histogram of losses for the entire set of simulations
# the histogram of the data
plt.hist(loss[:,tind], 10, normed=1, facecolor='green', alpha=0.75)
plt.show()
In [ ]:
In [ ]: