In [ ]:

    

from __future__ import division
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt

from pybasicbayes.distributions import Regression
from pybasicbayes.util.text import progprint_xrange
from autoregressive.distributions import AutoRegression

import pyximport
pyximport.install()
import os
os.chdir('../pylds')
from models import LDS, DefaultLDS
from distributions import Regression_diag, AutoRegression_input
from obs_scheme import ObservationScheme
from user_util import gen_pars, rand_rotation_matrix, init_LDS_model, collect_LDS_stats
#npr.seed(0)

save_file = '/home/mackelab/Desktop/Projects/Stitching/results/LDS_save'
py3res = np.load(save_file + '.npz')

data = py3res['data']
T,p = data.shape
stateseq = py3res['stateseq']
n   = stateseq.shape[1]
obs_scheme = py3res['obs_scheme'].reshape(1,)[0]
obs_scheme = ObservationScheme(p=p, T=T, 
                               sub_pops=obs_scheme['sub_pops'], 
                               obs_pops=obs_scheme['obs_pops'], 
                               obs_time=obs_scheme['obs_time'])

m = len(py3res['likes'])-1
eps = py3res['eps']

pars_true = py3res['pars_true'].reshape(1,)[0] # some numpy (cross-version 3.x-2.x?) 
tmp = {}
for i in range(len(pars_true.keys())):
    tmp[str(pars_true.keys()[i])] =  pars_true.values()[i]
pars_true = tmp
pars_init = py3res['pars_init'].reshape(1,)[0] # savez wierdness
tmp = {}
for i in range(len(pars_init.keys())):
    tmp[str(pars_init.keys()[i])] =  pars_init.values()[i]
pars_init = tmp

print('(T, p, n, m) = ', (T, p, n, m))


    

In [ ]:

    

def update(model):
    model.EM_step()
    return model.log_likelihood()                    



###################
#    EM cycles    #
###################


# get E-step results for init pars
model = init_LDS_model(pars_init, data, obs_scheme) # set to initialisation
model.E_step()
stats_init,_ = collect_LDS_stats(model)

# get EM-step results after one iteration
model.M_step()
model.E_step()
stats_first,pars_first = collect_LDS_stats(model)

# get EM-step results after m iterations                    
model = init_LDS_model(pars_init, data, obs_scheme) # reset to initialisation                    
print 'fitting'
likes = [update(model) for _ in progprint_xrange(m)]
stats_hat,pars_hat = collect_LDS_stats(model)
                                                        
# get EM-step results from true parameters
model = init_LDS_model(pars_true, data, obs_scheme) # reset to true pars
model.E_step()
stats_true,_ = collect_LDS_stats(model)
model.M_step()                                        

###################
#  store results  #
###################

from scipy.io import savemat # store results for comparison with Matlab code   
from scipy.linalg import solve_discrete_lyapunov as dtlyap # solve discrete-time Lyapunov equation
save_file_m = {'x': model.states_list[0].stateseq, 
               'y': model.states_list[0].data,
               'u' : [], 
               'll' : likes, 
               'T' : model.states_list[0].T, 
               'Trial': len(model.states_list), 
               'elapsedTime' : 0,
               'ifUseB':False, 
               'ifUseA':True, 
               'epsilon':0,
               'ifRDiagonal':False,
               'covConvEps':0,        
               'truePars':pars_true,
               'initPars':pars_init,
               'firstPars':pars_first,
               'estPars': pars_hat,
               'stats_0': stats_init,
               'stats_1': stats_first,
               'stats_h': stats_hat,
               'stats_true': stats_true,
               'Pi':dtlyap(pars_true['A'], pars_true['Q']),
               'Pi_h':dtlyap(pars_hat['A'], pars_hat['Q']),
               'Pi_t':pars_true['A'].dot(dtlyap(pars_true['A'], pars_true['Q'])),
               'Pi_t_h':pars_hat['A'].dot(dtlyap(pars_hat['A'], pars_hat['Q'])),
               'obsScheme' : []}

savemat(save_file+'_mattjj',save_file_m) # does the actual saving


    

In [ ]:

    




    

In [ ]: