In [1]:
%matplotlib inline
from matplotlib import pylab as pl
import cPickle as pickle
import pandas as pd
import numpy as np
import os
import random

In [80]:
import scipy.io
W=10
good_jump = []
bad_jump = []
for target in ['Dog_1', 'Dog_2', 'Dog_3', 'Dog_4', 'Dog_5']: #, 'Patient_1', 'Patient_2']:
    print target
    for data_type in ['preictal', 'interictal']:
        last_data = None
        for segment in range(1,1000000):
            fname = '../seizure-data/%s/%s_%s_segment_%04d.mat'%(target,target,data_type,segment)
            try:
                data = scipy.io.loadmat(fname)
            except:
                break
            k = '%s_segment_%d'%(data_type,segment)
            sequence = data[k]['sequence'][0,0][0,0]
            d = data[k]['data'][0,0]

            if last_data is not None:
                if last_sequence+1 == sequence:
                    good_jump.append((last_data,d[:,:W].astype(float))) # copy is critical becayse d[:] is just a window
                else:
                    bad_jump.append((last_data,d[:,:W].astype(float)))
            last_sequence = sequence
            last_data = d[:,-W:].astype(float)


Dog_1
Dog_2
Dog_3
Dog_4
Dog_5

In [66]:
len(good_jump),len(bad_jump),np.vstack((good_jump[0][0][:,-1],good_jump[0][1][:,0])).shape


Out[66]:
(3281, 648, (2, 16))

In [67]:
scipy.spatial.distance.pdist(np.vstack((good_jump[0][0][:,-1],good_jump[0][1][:,0]))).shape


Out[67]:
(1,)

In [157]:
def dotscore(j1, j2):
    return -2.*np.dot(j1,j2)/(np.dot(j1,j1) + np.dot(j2,j2))
def mydist(v1, v2):
    j1 = v1[:,-1]
    j2 = v2[:,0]
#     score = dotscore(j1,j2)
    q1 = v1[:,-1] - v1[:,-2]
    q2 = v2[:,1] - v2[:,0]
    score = -np.dot(j1,j2)/np.sqrt(np.dot(j1,j1)*np.dot(j2,j2))
    return score #np.sqrt(score)

In [158]:
def proc(jumps, alpha=1):
    import scipy.spatial.distance
    jp = np.array([])
    for jump in jumps:
        j = mydist(jump[0], jump[1])
#         j = scipy.spatial.distance.pdist(np.vstack((jump[0][:,-1],jump[1][:,0])),metric='cosine')
        jp = np.hstack((jp,j))
    print jp.shape
    print 'mean',jp.mean()
#     jp = np.abs(jp - jp.mean())
    lm = np.log(2.)/np.median(jp) # or 1/jp.mean()
    print 'lambda',lm
    y,x,_=pl.hist(jp, bins=50, alpha=alpha,normed=True);
    #pl.plot(x,(x[1]-x[0])*len(jp)*lm*np.exp(-lm*x))

In [159]:
proc(good_jump)
proc(bad_jump, alpha=0.5)


(3281,)
mean -0.867352991073
lambda -0.738084615073
(648,)
mean -0.0200832723617
lambda -22.531309286

In [135]:


In [ ]: