Experiment 3: TRo Journal


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

Inference Functions



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

Plotting Functions



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

Main Analysis



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')

Shirt Evaluation


The p-values for T-shirt test inference are not sufficient. This experiment is to perform test inference for more trials and perform significance test on larger dataset.


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')