Compare the predictive performance of using different inference techniques in MRD. The performance is evaluated over unseen T-shirts and postures that are not part of the training dataset.
The metrics for evaluation are RMS error, normalized RMS error and pearson correlation. Statistical significance is evaluated using wilcoxon signed rank sum test.
In this notebook, test inference is performed using 3 different inference techniques and the results are stored along with the evaluation metrics.
In [ ]:
# import the modules
import os
import GPy
import csv
import sys
import time
import filterpy
import numpy as np
import cPickle as pickle
from matplotlib import cm
import scipy.stats as stats
import sklearn.metrics as metrics
import GPy.plotting.Tango as Tango
from matplotlib import pyplot as plt
from filterpy.kalman import KalmanFilter
from sklearn.neighbors import NearestNeighbors
from filterpy.kalman import UnscentedKalmanFilter as UKF
%matplotlib notebook
In [ ]:
# NN Search based inference function
def nnsearchFunc(mrdModel, testData):
# model variables
kKey = 'Cloud'
mKey = 'TopCoord'
# initialize dimension variables
qDim = mrdModel.X.mean.shape[1]
nDimIn = testData[kKey].shape[1]
nDimOut = testData[mKey].shape[1]
nSamples = testData[mKey].shape[0]
latentVals = np.zeros((nSamples,qDim))
predictVals = np.zeros((nSamples,nDimOut))
# obtain the training data latent positions
latentPositions = mrdModel.X.mean.values
# build KD tree to search for nearest neighbors
nn = NearestNeighbors(n_neighbors=5,algorithm='kd_tree').fit(mrdModel.Y0.Y)
# loop over dataset
startTime = time.time()
for n in range(nSamples):
# yIn is y*, yTrueOut is z*
yIn = testData[kKey][n,:]
yTrueOut = testData[mKey][n,:]
# approximate x* with xtrain
_, indices = nn.kneighbors(np.atleast_2d(yIn))
latentVal = latentPositions[indices[0],:].mean(axis=0)
# obtain z* from x* through GP mapping
yOut = mrdModel.predict(np.atleast_2d(latentVal), Yindex=1)
latentVals[n,:] = latentVal
predictVals[n,:] = yOut[0]
# print output
sys.stdout.write('.')
sys.stdout.flush()
stopTime = time.time()
print '\nFinished Strategy NN Search'
nrmse = np.divide(np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values')),
testData[mKey].max(axis=0) - testData[mKey].min(axis=0))
rmse = np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values'))
corr = np.zeros((1,nDimOut))
for d in range(nDimOut):
corr[0,d],_ = stats.pearsonr(testData[mKey][:,d],predictVals[:,d])
results = {}
results['corr'] = corr
results['rmse'] = rmse
results['nrmse'] = nrmse
results['pred'] = predictVals
results['latent'] = latentVals
results['time'] = nSamples/(stopTime - startTime)
return results
In [ ]:
def optimizeFunc(mrdModel, testData):
# model variables
kKey = 'Cloud'
mKey = 'TopCoord'
qDim = mrdModel.X.mean.shape[1]
nDimIn = testData[kKey].shape[1]
nDimOut = testData[mKey].shape[1]
nSamples = testData[mKey].shape[0]
latentVals = np.zeros((nSamples,qDim))
predictVals = np.zeros((nSamples,nDimOut))
# obtain the training data latent positions
latentPositions = mrdModel.X.mean.values
nn = NearestNeighbors(n_neighbors=5,algorithm='kd_tree').fit(mrdModel.Y0.Y)
startTime = time.time()
for n in range(nSamples):
yIn = testData[kKey][n,:]
yTrueOut = testData[mKey][n,:]
# infer x* through variational inference
[xPredict, infX] = mrdModel.Y0.infer_newX(yIn[None,:], optimize=True)
latentVal = xPredict.mean
# obtain z* from x* using GP mapping
yOut = mrdModel.predict(np.atleast_2d(latentVal), Yindex=1)
latentVals[n,:] = latentVal
predictVals[n,:] = yOut[0]
# print output
sys.stdout.write('.')
sys.stdout.flush()
stopTime = time.time()
print '\nFinished Strategy Optimize'
nrmse = np.divide(np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values')),
testData[mKey].max(axis=0) - testData[mKey].min(axis=0))
rmse = np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values'))
corr = np.zeros((1,nDimOut))
for d in range(nDimOut):
corr[0,d],_ = stats.pearsonr(testData[mKey][:,d],predictVals[:,d])
results = {}
results['corr'] = corr
results['rmse'] = rmse
results['nrmse'] = nrmse
results['pred'] = predictVals
results['latent'] = latentVals
results['time'] = nSamples/(stopTime - startTime)
return results
In [ ]:
def filterFunc(mrdModel, testData):
qDim = mrdModel.X.mean.values.shape[1]
scales1 = mrdModel.Y0.kern.input_sensitivity(summarize=False)
scales2 = mrdModel.Y1.kern.input_sensitivity(summarize=False)
scales1 = scales1/scales1.max()
scales2 = scales2/scales2.max()
# get the number of dimensions
yThresh = 0.05
indices = np.asarray(range(qDim))
active1 = indices[scales1 >= yThresh]
active2 = indices[scales2 >= yThresh]
sharedDims = np.intersect1d(active2,active2)
nShared = len(sharedDims)
# get init latent state from optimization
hybridFPS = 15.0
deltaT = 1.0/30.0
# state transition matrix
def f_cv(x, dt):
nShared = len(x)/2
F = np.eye(2*nShared)
F[:nShared,nShared:] = dt*np.eye(nShared)
return np.dot(F, x)
def h_cv(x):
nShared = len(x)/2
return x[:nShared]
# create kalman filter
sigmas = filterpy.kalman.MerweScaledSigmaPoints(n=2*nShared, alpha=0.1, beta=2.0, kappa=1.0)
kf = UKF(dim_x=2*nShared, dim_z=nShared, fx=f_cv, hx=h_cv, dt=deltaT, points=sigmas)
# init state
yIn = testData['Cloud'][0,:]
[xPredict, infX] = mrdModel.Y0.infer_newX(yIn[None,:], optimize=True)
xPredict = xPredict.mean
kf.x = np.zeros((2*nShared))
kf.x[:nShared] = xPredict[0,sharedDims]
# init covariance
kf.P *= 1e-4
# process and measurement noise
kf.Q *= 1e-5
kf.R *= 1e-3
# model variables
kKey = 'Cloud'
mKey = 'TopCoord'
qDim = mrdModel.X.mean.shape[1]
nDimIn = testData[kKey].shape[1]
nDimOut = testData[mKey].shape[1]
nSamples = testData[mKey].shape[0]
latentVals = np.zeros((nSamples,qDim))
predictVals = np.zeros((nSamples,nDimOut))
# obtain the training data latent positions
latentPositions = mrdModel.X.mean
nn = NearestNeighbors(n_neighbors=5,algorithm='kd_tree').fit(mrdModel.Y0.Y)
# main inference loop
startTime = time.time()
for n in range(nSamples):
# yIn is y* and yTrueOut is z*
yIn = testData[kKey][n,:]
yTrueOut = testData[mKey][n,:]
# alternate between nearest neighbor inference and optimization
kf.predict()
if n%hybridFPS == 0:
[xPredict, infX] = mrdModel.Y0.infer_newX(yIn[None,:], optimize=True)
xPredict = xPredict.mean
kf.update(xPredict[0,sharedDims], R=1e-6*np.eye(nShared))
else:
_,indices = nn.kneighbors(np.atleast_2d(yIn))
xPredict = latentPositions[indices[0],:].mean(axis=0)
kf.update(xPredict[sharedDims])
# update the state using unscented kalman filter
latentVal = np.atleast_2d(xPredict)
latentVal[0,sharedDims] = kf.x[:nShared]
# obtain z* from x* using GP mapping
yOut = mrdModel.predict(latentVal, Yindex=1)
latentVals[n,:] = latentVal
predictVals[n,:] = yOut[0]
# print output
sys.stdout.write('.')
sys.stdout.flush()
stopTime = time.time()
print '\nFinished Strategy Hybrid'
nrmse = np.divide(np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values')),
testData[mKey].max(axis=0) - testData[mKey].min(axis=0))
rmse = np.sqrt(metrics.mean_squared_error(testData[mKey],predictVals,multioutput='raw_values'))
corr = np.zeros((1,nDimOut))
for d in range(nDimOut):
corr[0,d],_ = stats.pearsonr(testData[mKey][:,d],predictVals[:,d])
results = {}
results['corr'] = corr
results['rmse'] = rmse
results['nrmse'] = nrmse
results['pred'] = predictVals
results['latent'] = latentVals
results['time'] = nSamples/(stopTime - startTime)
return results
In [ ]:
def plotLatent(inX, options):
s = 100
marker = 'o'
resolution = 50
Tango.reset()
fig = plt.figure()
ax = fig.add_subplot(111)
input1, input2 = options['indices']
for i in range(len(inX)):
if inX[i].shape[0] > options['maxPoints'][i]:
subsample = np.random.choice(inX[i].shape[0], size=options['maxPoints'][i], replace=False)
inX[i] = inX[i][subsample]
xmin, ymin = inX[0][:, [input1, input2]].min(0)
xmax, ymax = inX[0][:, [input1, input2]].max(0)
x_r, y_r = xmax-xmin, ymax-ymin
xmin -= .1*x_r
xmax += .1*x_r
ymin -= .1*y_r
ymax += .1*y_r
print xmin, xmax, ymin, ymax
if options['plotVariance']:
def plotFunction(x):
Xtest_full = np.zeros((x.shape[0], qDim))
Xtest_full[:, [input1, input2]] = x
_, var = model.predict(np.atleast_2d(Xtest_full))
var = var[:, :1]
return -np.log(var)
model = options['model']
qDim = model.X.mean.shape[1]
x, y = np.mgrid[xmin:xmax:1j*resolution, ymin:ymax:1j*resolution]
gridData = np.hstack((x.flatten()[:, None], y.flatten()[:, None]))
gridVariance = (plotFunction(gridData)).reshape((resolution, resolution))
varianceHandle = plt.imshow(gridVariance.T, interpolation='bilinear', origin='lower', cmap=cm.gray,
extent=(xmin, xmax, ymin, ymax))
handles = [None]*len(inX)
for i in range(len(inX)):
handles[i] = ax.scatter(inX[i][:, input1], inX[i][:, input2], marker=marker, s=s, c=options['color'][i], linewidth=.2, edgecolor='k', alpha=0.9)
ax.grid(b=False)
ax.set_aspect('auto')
ax.tick_params(axis='both', labelsize=20)
ax.set_title(options['title'], fontsize=25)
ax.legend(handles,options['legend'],loc=2,fontsize=20)
ax.set_xlabel('Latent Dimension %i' % (input1+1), fontsize=20)
ax.set_ylabel('Latent Dimension %i' % (input2+1), fontsize=20)
ax.set_xlim((xmin, xmax))
ax.set_ylim((ymin, ymax))
fig.canvas.draw()
fig.tight_layout()
fig.canvas.draw()
return ax
In [ ]:
nShr = 4
nPos = 6
names = []
dims = [1,7500,8]
keys = ['Time','Cloud','TopCoord']
for nS in range(nShr):
for nP in range(nPos):
names.append('K1S%dP%dT1' % (nS+1,nP+1))
# create directory for results
dName = '../Results/Exp3'
if not os.path.exists(dName):
os.makedirs(dName)
Data = pickle.load(open('../Data/Data.p','rb'))
In [ ]:
experiments = ['Pos', 'Shr']
functions = ['filter','nnsearch','optimize']
for nS in range(nShr):
for nP in range(nPos):
print 'Cycle:%d,%d' % (nS+1,nP+1)
testInd = {'Pos':nS*nPos+nP, 'Shr':((nS+1)%nShr)*nPos+nP}
# load mrd model and obtain plot indices
mrdModel = pickle.load(open('../Models/Model%d%d.p' % (nS+1,nP+1), 'rb'))
scales = 1./np.asarray(mrdModel.Y0.kern.lengthscale)
indices = scales.argsort()[-2:][::-1]
for exp in ['Pos']:
testData = {}
for key,dim in zip(keys,dims):
testData[key] = Data[names[testInd[exp]]][key]
results = {}
results['filter'] = filterFunc(mrdModel, testData)
results['nnsearch'] = nnsearchFunc(mrdModel, testData)
results['optimize'] = optimizeFunc(mrdModel, testData)
pickle.dump(results,open('../Results/Exp3/Res%s%d%d.p' % (exp,nS+1,nP+1), 'wb'))
# latent space plot
mrdX = []
for func in functions:
mrdX.append(results[func]['latent'])
maxPoints = [300,300,300]
options = {'title':'',
'indices':indices,
'plotVariance':True,
'model':mrdModel.Y0,
'maxPoints':maxPoints,
'legend':['hybrid','nnsearch','optimize'],
'color':[Tango.colorsHex['mediumRed'],Tango.colorsHex['darkGreen'],Tango.colorsHex['mediumBlue']]}
plotLatent(mrdX, options)
plt.savefig('../Results/Exp3/Latent%s%d%d.pdf' % (exp,nS+1,nP+1), format='pdf')
# inference plot
mKey = 'TopCoord'
times = testData['Time'][:,1]
nSamplesTest = testData[mKey].shape[0]
legendInds = [0]
plotInds = [1,2,5,7]
legendLabels = ['Ground Truth','NN Search','Filter','Optimize']
ylabels = ['Center of Twist','Writhe','Center of Twist','Center of Twist']
titles = ['Collar-Head','Collar-Body','Left Sleeve-Arm','Right Sleeve-Arm']
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times,testData[mKey][:,plotInds[0]],c='k',linewidth=2)
ax.plot(times,results['nnsearch']['pred'][:,plotInds[0]],c='r',linewidth=2)
ax.plot(times,results['filter']['pred'][:,plotInds[0]],c='g',linewidth=2)
ax.plot(times,results['optimize']['pred'][:,plotInds[0]],c='b',linewidth=2)
plt.title(titles[0], fontsize=25)
plt.ylabel(ylabels[0], fontsize=20)
plt.xlabel('Time (secs)', fontsize=20)
plt.tick_params(axis='both', labelsize=20)
plt.legend(legendLabels,fontsize=20,loc=4)
plt.savefig('../Results/Exp3/Inference%s%d%d.pdf' % (exp,nS+1,nP+1), format='pdf')
In [ ]:
experiments = ['Pos', 'Shr']
functions = ['filter','nnsearch','optimize']
for nS in range(nShr):
for nP in range(nPos):
print 'Cycle:%d,%d' % (nS+1,nP+1)
testInd = {'Pos':[nS*nPos+nP], 'Shr':[((nS+1)%nShr)*nPos+nP,((nS+1)%nShr)*nPos+(nP+1)%nPos,
((nS+1)%nShr)*nPos+(nP+2)%nPos]}
# load mrd model and obtain plot indices
mrdModel = pickle.load(open('../Models/Model%d%d.p' % (nS+1,nP+1), 'rb'))
scales = 1./np.asarray(mrdModel.Y0.kern.lengthscale)
indices = scales.argsort()[-2:][::-1]
for exp in ['Shr']:
for ind in testInd[exp]:
testData = {}
for key,dim in zip(keys,dims):
testData[key] = Data[names[ind]][key]
results = {}
results['filter'] = filterFunc(mrdModel, testData)
results['nnsearch'] = nnsearchFunc(mrdModel, testData)
results['optimize'] = optimizeFunc(mrdModel, testData)
pickle.dump(results,open('../Results/Exp3/Res%s%d%d%02d.p' % (exp,nS+1,nP+1,ind+1), 'wb'))
# latent space plot
mrdX = []
for func in functions:
mrdX.append(results[func]['latent'])
maxPoints = [300,300,300]
options = {'title':'',
'indices':indices,
'plotVariance':True,
'model':mrdModel.Y0,
'maxPoints':maxPoints,
'legend':['hybrid','nnsearch','optimize'],
'color':[Tango.colorsHex['mediumRed'],Tango.colorsHex['darkGreen'],Tango.colorsHex['mediumBlue']]}
plotLatent(mrdX, options)
plt.savefig('../Results/Exp3/Latent%s%d%d%02d.pdf' % (exp,nS+1,nP+1,ind), format='pdf')
# inference plot
mKey = 'TopCoord'
times = testData['Time'][:,1]
nSamplesTest = testData[mKey].shape[0]
legendInds = [0]
plotInds = [1,2,5,7]
legendLabels = ['Ground Truth','NN Search','Filter','Optimize']
ylabels = ['Center of Twist','Writhe','Center of Twist','Center of Twist']
titles = ['Collar-Head','Collar-Body','Left Sleeve-Arm','Right Sleeve-Arm']
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(times,testData[mKey][:,plotInds[0]],c='k',linewidth=2)
ax.plot(times,results['nnsearch']['pred'][:,plotInds[0]],c='r',linewidth=2)
ax.plot(times,results['filter']['pred'][:,plotInds[0]],c='g',linewidth=2)
ax.plot(times,results['optimize']['pred'][:,plotInds[0]],c='b',linewidth=2)
plt.title(titles[0], fontsize=25)
plt.ylabel(ylabels[0], fontsize=20)
plt.xlabel('Time (secs)', fontsize=20)
plt.tick_params(axis='both', labelsize=20)
plt.legend(legendLabels,fontsize=20,loc=4)
plt.savefig('../Results/Exp3/Inference%s%d%d%02d.pdf' % (exp,nS+1,nP+1,ind+1), format='pdf')