In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
import pyrem as pr
import numpy as np
import pylab as pl
import pandas as pd
from scipy import signal


/home/quentin/.virtualenvs/pyrem/lib/python2.7/site-packages/pandas/io/excel.py:626: UserWarning: Installed openpyxl is not supported at this time. Use >=1.6.1 and <2.0.0.
  .format(openpyxl_compat.start_ver, openpyxl_compat.stop_ver))

In [3]:
seq = np.genfromtxt("/home/quentin/Desktop/rem_chunking/test_eeg_no_header_small.txt")
to_remove =10 * int(seq.size / 100.0)
seq = seq[to_remove: -to_remove]
# seq = np.sin(np.linspace(0, 200, 100000)) *10
# seq = seq + np.sin(np.linspace(0, 20, 100000)) * 2
# seq = seq + np.random.normal(0,10,100000)
sig = pr.Signal(seq, 200.0)

In [4]:
sig = sig.resample(256)
print sig.duration
sig[0:300000]


1:06:39.996094
Out[4]:

In [5]:
n = 25
c1 = sig[sig.sampling_freq * 20 * n:sig.sampling_freq * 20 * (n+1)]

f, pxx = signal.welch(c1, sig.sampling_freq, window="blackmanharris")
pl.plot(f, pxx)

n = 20
c2 = sig[sig.sampling_freq * 20 * n:sig.sampling_freq * 20 * (n+1)]
#boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann,
f, pxx = signal.welch(c2, sig.sampling_freq, window="bartlett",nperseg=256)
pl.plot(f, pxx)
pl.xscale("log")

print f[np.argmax(pxx)]


7.0

In [6]:
ff = pr.features.FeatureFactory()

df = pd.concat([ff(1, c1), ff(2, c2)])
df


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-8deff7184fb5> in <module>()
      1 ff = pr.features.FeatureFactory()
      2 
----> 3 df = pd.concat([ff(1, c1), ff(2, c2)])
      4 df

/home/quentin/comput/pyrem-git/src/pyrem/features/feature_factory.pyc in __call__(self, t, signal)
     17         ]
     18     def __call__(self, t, signal):
---> 19         dfs = [ group(t,signal) for group in self.feature_groups]
     20         return pd.concat(dfs, axis=1)

/home/quentin/comput/pyrem-git/src/pyrem/features/feature_base.pyc in __call__(self, t, signal)
      9     prefix = None
     10     def __call__(self, t, signal):
---> 11         feature_dict = self._make_feature_vec(signal)
     12 
     13         data_frame = pd.DataFrame(feature_dict, index=[t])

/home/quentin/comput/pyrem-git/src/pyrem/features/entropy.pyc in _make_feature_vec(self, signal)
     40 
     41         out = dict()
---> 42         out["spectral"] = self.spectral_entropy(signal,[0,50]) # fixme magic number here
     43         return out

/home/quentin/comput/pyrem-git/src/pyrem/features/entropy.pyc in spectral_entropy(self, signal, Band)
     30         for i in xrange(0, len(Power_Ratio) - 1):
     31             Spectral_Entropy += Power_Ratio[i] * log(Power_Ratio[i])
---> 32         Spectral_Entropy /= log(len(Power_Ratio))    # to save time, minus one is omitted
     33         return -1 * Spectral_Entropy
     34 

NameError: global name 'log' is not defined

In [7]:
ff = pr.features.FeatureFactory()
dfs = [ff(t,s) for t, s in sig.embed_seq(sig.sampling_freq * 30, sig.sampling_freq * 20 )]
features = pd.concat(dfs)
features.shape


Out[7]:
(2364, 11)

In [50]:
plot(features["welch_mean"], features["welch_sd"],".")
pl.figure()
plot(features["welch_skew"], features["welch_sd"],".")
pl.figure()
plot(features["welch_skew"], features["welch_median"],".")
pl.figure()
plot(features["welch_median"], features["welch_mean"],".")
pl.figure()


Out[50]:
<matplotlib.figure.Figure at 0x7f9fa40b9f90>
<matplotlib.figure.Figure at 0x7f9fa40b9f90>

In [33]:
pl.hist(np.array(features["welch_median"]), bins=100)


Out[33]:
(array([  31.,  104.,  178.,  191.,  217.,  151.,  137.,  104.,   84.,
          63.,   53.,   58.,   59.,   47.,   36.,   30.,   31.,   32.,
          28.,   27.,   38.,   20.,   23.,   24.,   19.,   18.,   15.,
          14.,   11.,    9.,   13.,   11.,   15.,   18.,   13.,   16.,
          15.,   20.,   14.,   16.,   20.,   24.,   11.,   15.,   17.,
          19.,   18.,   15.,   10.,   12.,    7.,    8.,   13.,    5.,
           7.,    8.,    3.,    3.,    3.,    4.,    3.,    4.,    2.,
           5.,    2.,    0.,    3.,    4.,    3.,    3.,    3.,    4.,
           1.,    4.,    2.,    7.,    3.,    6.,    8.,    5.,    6.,
           3.,    6.,    2.,    5.,    7.,   10.,    0.,    3.,    3.,
           3.,    0.,    4.,    4.,    6.,    0.,    1.,    3.,    0.,    1.]),
 array([  1.31459704,   1.57313057,   1.83166411,   2.09019765,
          2.34873118,   2.60726472,   2.86579825,   3.12433179,
          3.38286533,   3.64139886,   3.8999324 ,   4.15846593,
          4.41699947,   4.67553301,   4.93406654,   5.19260008,
          5.45113361,   5.70966715,   5.96820068,   6.22673422,
          6.48526776,   6.74380129,   7.00233483,   7.26086836,
          7.5194019 ,   7.77793544,   8.03646897,   8.29500251,
          8.55353604,   8.81206958,   9.07060312,   9.32913665,
          9.58767019,   9.84620372,  10.10473726,  10.3632708 ,
         10.62180433,  10.88033787,  11.1388714 ,  11.39740494,
         11.65593847,  11.91447201,  12.17300555,  12.43153908,
         12.69007262,  12.94860615,  13.20713969,  13.46567323,
         13.72420676,  13.9827403 ,  14.24127383,  14.49980737,
         14.75834091,  15.01687444,  15.27540798,  15.53394151,
         15.79247505,  16.05100858,  16.30954212,  16.56807566,
         16.82660919,  17.08514273,  17.34367626,  17.6022098 ,
         17.86074334,  18.11927687,  18.37781041,  18.63634394,
         18.89487748,  19.15341102,  19.41194455,  19.67047809,
         19.92901162,  20.18754516,  20.44607869,  20.70461223,
         20.96314577,  21.2216793 ,  21.48021284,  21.73874637,
         21.99727991,  22.25581345,  22.51434698,  22.77288052,
         23.03141405,  23.28994759,  23.54848113,  23.80701466,
         24.0655482 ,  24.32408173,  24.58261527,  24.84114881,
         25.09968234,  25.35821588,  25.61674941,  25.87528295,
         26.13381648,  26.39235002,  26.65088356,  26.90941709,  27.16795063]),
 <a list of 100 Patch objects>)

In [9]:
from sklearn import cluster

In [45]:
f = np.array(features[["welch_median","welch_sd","welch_kurtosis","welch_skew"]] )

km = cluster.KMeans(2)
fit = km.fit(f)

In [42]:
fit


Out[42]:
KMeans(copy_x=True, init='k-means++', max_iter=300, n_clusters=2, n_init=10,
    n_jobs=1, precompute_distances=True, random_state=None, tol=0.0001,
    verbose=0)

In [49]:
p = km.predict(f)
color_map = ["k","b","r","g"]
pl.figure()
for i, (c, x, y) in enumerate(zip(p, features["welch_median"], features["welch_sd"])):
    pl.plot(x, y,'.',color=color_map[c])
pl.figure()
for i, (c, x, y) in enumerate(zip(p, features["welch_median"], features["welch_kurtosis"])):
    pl.plot(x, y,'.',color=color_map[c])
pl.figure()
for i, (c, x, y) in enumerate(zip(p, features["welch_skew"], features["welch_kurtosis"])):
    pl.plot(x, y,'.',color=color_map[c])



In [12]:


In [ ]: