import numpy as np
from msmbuilder.Trajectory import Trajectory
from msmbuilder.Project import Project
from euclid.clustering import concatenate_trajectories
from time import time
from euclid import dihedral
from Geometry import DihedralTools

proj = Project.LoadFromHDF('/Users/rmcgibbo/Desktop/msmb_tutorial/ProjectInfo.h5')
trjs = [proj.LoadTraj(i) for i in range(proj['NumTrajs'])]
trj = concatenate_trajectories(trjs)

phi_indices = DihedralTools.GetAtomIndicesForPhi(trj)
print 'phi indices'
print phi_indices

xyzlist = trj['XYZList']
start = time()
phi = np.zeros((len(xyzlist), len(phi_indices)-1))
for i in range(len(xyzlist)):
    phi_xyz = xyzlist[i,phi_indices]
    phi[i] = DihedralTools.ExtractDihedralsFromArray(phi_xyz)
print 'Python coe time', time() - start

phi_indices = phi_indices[0:-1,:]
start = time()
newphi = dihedral.compute_dihedrals(trj, phi_indices, degrees=True)
print 'C code time', time() - start

diff = newphi - phi

import matplotlib.pyplot as pp
pp.title('Difference between codes')
pp.hist(np.ravel(diff))
pp.show()