Manually run this notebooks on all targets
In [38]:
target = 'Patient_2'
In [39]:
%matplotlib inline
from matplotlib import pylab as pl
import cPickle as pickle
import pandas as pd
import numpy as np
import os
import re
import math
import sys
import random
from collections import defaultdict
In [40]:
import sys
sys.path.append('..')
from common.data import CachedDataLoader
cached_data_loader = CachedDataLoader('../data-cache')
def read_data(target, data_type, features):
return cached_data_loader.load('data_%s_%s_%s'%(data_type,target,features),None)
In [41]:
FEATURES = 'gen-8_allbands2-usf-w60-b0.2-b4-b8-b12-b30-b70'
FEATURES1 = 'gen-8_maxdiff-60'
nbands = 0
nwindows = 0
for p in FEATURES.split('-'):
if p[0] == 'b':
nbands += 1
elif p[0] == 'w':
nwindows = int(p[1:])
nbands -= 1
nbands, nwindows
Out[41]:
each target receive a different model
positive examples. The positive examles were upsampled (using gen_ictal=-8)
In [42]:
from common.data import CachedDataLoader
cached_data_loader = CachedDataLoader('../data-cache')
def read_data(target, data_type, features=FEATURES):
fname = 'data_%s_%s_%s'%(data_type,target,features)
print fname
return cached_data_loader.load(fname,None)
def process(X, X1, percentile=[0.05,0.95], nunits=2, mask_level=7000):
N, Nf = X.shape
print '# samples',N,'# power points', Nf
nchannels = Nf / (nbands*nwindows)
print '# channels', nchannels
fix = defaultdict(int)
newX = []
for i in range(N):
nw = nwindows//nunits
windows = X[i,:].reshape((nunits,nw,-1))
mask = X1[i,:].reshape((nunits,nw,-1)) # max value for each channel
for j in range(nunits):
for k in range(nchannels):
m = mask[j,:,k] > mask_level # find large windows
if np.any(m):
# print 'FIX', sum(m)
fix[sum(m)] += 1
if not np.all(m): # make sure we had at least one good window so we can re use its values
# replace the bands of a large windows with the mean of the bands in all other windows
windows[j,m,k*nbands:(k+1)*nbands] = np.mean(windows[j,~m,k*nbands:(k+1)*nbands], axis=0)
sorted_windows = np.sort(windows, axis=1)
features = np.concatenate([sorted_windows[:,int(p*nw),:] for p in percentile], axis=-1)
newX.append(features.ravel())
newX = np.array(newX)
print sorted(fix.items())
return newX
def getsegments(pdata):
segments = []
start = 0
last_l = 0
for i,l in enumerate(pdata.latencies):
if l<last_l:
segments.append(range(start,i))
start = i
last_l = l
segments.append(range(start,i+1))
return segments
def getdata():
pdata = read_data(target, 'preictal') # positive examples
Np, _ = pdata.X.shape
print 'Positive examples',Np,
psegments = getsegments(pdata)
Nps = len(psegments)
print 'sequences:',Nps
ndata = read_data(target, 'interictal') # negative examples
Nn, _ = ndata.X.shape
print 'Negative',Nn,
nsegments = getsegments(ndata)
Nns = len(nsegments)
print 'sequences:',Nns
X = np.concatenate((pdata.X, ndata.X))
print 'p-ratio',float(Np)/(Np+Nn), 'sequences p-ratio:',float(Nps)/(Nps+Nns)
nsegments = [[s+Np for s in ns] for ns in nsegments]
latencies = np.concatenate((pdata.latencies,ndata.latencies))
pdata1 = read_data(target, 'preictal', FEATURES1) # positive examples
ndata1 = read_data(target, 'interictal', FEATURES1) # negative examples
X1 = np.concatenate((pdata1.X, ndata1.X))
print 'Training:'
X = process(X, X1)
y = np.zeros(X.shape[0])
y[:Np] = 1
print 'Test:'
tdata = read_data(target, 'test') # test examples
tdata1 = read_data(target, 'test', FEATURES1) # test examples
Xt = process(tdata.X, tdata1.X)
return X,y,Xt,psegments,nsegments,latencies
data is broken into segments that should be taken together when splitting to training and validation
In [43]:
X,y,Xt,psegments,nsegments,latencies = getdata()
In [44]:
N, NF = X.shape
N, NF
Out[44]:
I am using RF because it needs very little (hyper) parameter tuning. On purpose I am using a small depth, because I am not interested in the best prediction (which is already high) but with the feature importance after taking into account pair interactions.
In [45]:
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=1000, n_jobs=-1, oob_score=True, max_depth=2)
rf.fit(X, y)
Out[45]:
In [46]:
rf.oob_score_
Out[46]:
In [47]:
pnweights = rf.feature_importances_
pn_importance_order = pnweights.argsort()[::-1]
pl.plot(rf.feature_importances_[pn_importance_order])
Out[47]:
we will use Gradient Boosting Classifier which usually gives better results than L1 or RF. In addition, like RF, it does not require normalization or PCA. However, unlike RF or L1 it has many hyper parameters that can effect its performance. In addition we need to decide how many features we want to use which is another hyper-parameter. Instead of manually guessing, we can use the hyperopt to do the guesing for us
we will perform several hyperopt search in parallel each running on a different bootstrap sample of the data
The data itself is identical so there is no need to duplicate it for each process and we will use shared memory (shmem)
In [48]:
!rm /tmp/X.mmap
mmX = np.memmap('/tmp/X.mmap', shape=X.shape, dtype=np.float32, mode='w+')
mmX[:,:] = X[:,pn_importance_order]
del mmX # flush to disk
We will use ipython parallel processing infrastructure. Visit Clusters tab in the Home page of ipython and start 8 engines (or as many cores you have on your machine) from the default profile
OR you can run the command line:
ipcluster start --n=8
wait a little bit (otherwise you will get an error on next cell)
In [49]:
#!sleep 30
In [50]:
from IPython.parallel import Client
client = Client()
lv = client.load_balanced_view()
#lv.set_flags(block = False, retries = 0)
clients=client[:]
Ncores = len(clients)
Ncores
Out[50]:
copy some information to all engines
In [51]:
clients['X_shape'] = X.shape
clients['y'] = y
clients['psegments'] = psegments
clients['nsegments'] = nsegments
load the shared memory on all engines
In [52]:
%%px
import numpy as np
N, NF = X_shape
X = np.memmap('/tmp/X.mmap', shape=X_shape, dtype=np.float32, mode='r')
In [53]:
%%px
import random, itertools
def random_train_validation_split(psegments=psegments, nsegments=nsegments, N=N, pratio=1):
"""Randomly pick one positive segment for validation and a matching number of negative segments"""
Nps = len(psegments)
assert Nps > 1
Nns = len(nsegments)
assert Nns > 1
npsratio = float(Nns)/Nps
Ntrainps = 1
Ntrainns = min(max(1,int(Ntrainps*npsratio+0.5)), Nns-1) # make sure we have something to train
s = random.choice(psegments)
ns = random.sample(nsegments,Ntrainns) # sequence based
n = list(itertools.chain(*ns)) # .ravel does not work - elements of nsegments are not of equal length
sample_validate = s + n
random.shuffle(sample_validate)
all_p = list(itertools.chain(*psegments))
all_n = list(itertools.chain(*nsegments))
testp = list(set(all_p) - set(s))
if pratio != 1:
testp *= pratio
# ntestp = len(testp)
# boot_ntestp = int(ntestp*pratio)
# w = np.ones(ntestp)/float(ntestp)
# testp = [testp[i] for i, n in enumerate(np.random.multinomial(boot_ntestp, w))
# for k in xrange(n)]
testn = list(set(all_n) - set(n))
sample_test = testp + testn
random.shuffle(sample_test)
return sample_test, sample_validate
We will optimize AUC
In [54]:
%%px
from gradient_boost_loss import NAUC
import os
def hyperopt_work(args):
from lockfile import LockFile
space = args.get('space')
pratio = int(space['pratio'])
sample_train, sample_validate = random_train_validation_split(pratio=pratio)
X_trn = X[sample_train,:]
y_trn = y[sample_train]
assert y_trn.mean() > 0.01 and y_trn.mean() < 0.99
X_val = X[sample_validate,:]
y_val = y[sample_validate]
def t_est(args):
try:
from sklearn.ensemble import GradientBoostingClassifier
est = GradientBoostingClassifier()
params = dict((k,int(v) if isinstance(v, float) and abs(v - int(v)) < 1e-3 else v)
for k, v in args.iteritems() if k not in ['nf', 'pratio'])
est.set_params(**params)
# est.loss__ = NAUC(2)
nf = int(args['nf'])
est.fit(X_trn[:,:nf], y_trn)
y_train = est.predict_proba(X_trn[:,:nf])[:,1]
y_validate = est.predict_proba(X_val[:,:nf])[:,1]
from sklearn.metrics import roc_auc_score
train_score = roc_auc_score(y_trn, y_train)
validate_score = roc_auc_score(y_val, y_validate)
lock = LockFile('.lock')
with lock:
with open('../data-cache/hyperopt.txt','a') as fp:
print >>fp, validate_score, train_score, os.getpid(), args
from hyperopt import STATUS_OK
return {'loss':1.-validate_score, 'status':STATUS_OK, 'pid':os.getpid(),
'train_score':train_score, 'validate_score':validate_score}
except Exception as e:
lock = LockFile('.lock')
with lock:
with open('../data-cache/hyperopt.txt','a') as fp:
print >>fp, 'failed', e, args
from hyperopt import STATUS_FAIL
return {'status':STATUS_FAIL, 'loss':1.} # 'loss' is mandatory
max_evals = args.get('max_evals', 100)
from hyperopt import fmin, tpe, Trials
# trials = Trials()
best = fmin( t_est, space, algo=tpe.suggest, max_evals=max_evals) #, trials=trials)
# import cPickle as pickle
# lock = LockFile('.lock')
# with lock:
# with open('../data-cache/hyperopt.spkl','ab') as fp:
# pickle.dump(trials, fp, -1)
return best
Define statistical space in which we will do our hyper-parameter search
In [55]:
%%px
from hyperopt import hp
from math import log
space = {
'n_estimators': 1000,
'pratio': 1,
'learning_rate': hp.loguniform('learning_rate', np.log(0.001), np.log(1.)),
'nf': hp.quniform( 'nf', 10, NF, 1),
'max_features': hp.quniform( 'max_features', 10, 30, 1),
'max_depth': hp.quniform( 'max_depth', 2, 18, 1),
'min_samples_leaf': hp.quniform( 'min_samples_leaf', 2, 30, 1),
# 'subsample': hp.uniform( 'subsample', 0.2, 0.9),
}
In [56]:
!rm ../data-cache/hyperopt.*
run hyperopt searches in parallel on all cores. Each hyperopt search will do 100 evaluations of the hyper parameters.
In [57]:
%%px
hyperopt_work({'space':space, 'max_evals':100})
wait for the jobs to end. This will take some time. Also your computer can get really hot, so use the time to arange some cooling to it.
In [58]:
!sort -n -r ../data-cache/hyperopt.txt > ../data-cache/hyperopt.sort.txt
!head -n 5 ../data-cache/hyperopt.sort.txt
In [59]:
fp = open('../data-cache/hyperopt.sort.txt')
hyeropt_results = []
for l in fp:
if l.startswith('failed'):
continue
l = l.split()
validate_score = float(l[0])
train_score = float(l[1])
pid = int(l[2])
args = eval(''.join(l[3:]))
hyeropt_results.append((validate_score, train_score, pid, args))
fp.close()
len(hyeropt_results)
Out[59]:
In [60]:
Nt, NFF = Xt.shape
assert NF == NFF
Nt, NFF
Out[60]:
In [61]:
!rm /tmp/Xt.mmap
Xt_shape = (Nt, NF)
mmXt = np.memmap('/tmp/Xt.mmap', shape=Xt_shape, dtype=np.float32, mode='w+')
mmXt[:,:] = Xt[:,pn_importance_order]
del mmXt # flush to disk
In [62]:
clients['Xt_shape'] = Xt_shape
clients['target'] = target
In [63]:
%%px
Xt = np.memmap('/tmp/Xt.mmap', shape=Xt_shape, dtype=np.float32, mode='r')
In [64]:
def predict_work(args):
from lockfile import LockFile
from sklearn.ensemble import GradientBoostingClassifier
import cPickle as pickle
N = X_shape[0]
NF = int(args.get('nf', X_shape[1]))
pratio = int(args.get('pratio',1))
# use out-of-bag samples to estimate the generalization error
sample_train, sample_validate = random_train_validation_split(pratio=pratio)
X_trn = X[sample_train,:NF]
y_trn = y[sample_train]
X_val = X[sample_validate,:NF]
y_val = y[sample_validate]
X_test = Xt[:,:NF]
from sklearn.ensemble import GradientBoostingClassifier
est = GradientBoostingClassifier()
params = dict((k,int(v) if isinstance(v, float) and abs(v - int(v)) < 1e-3 else v)
for k, v in args.iteritems() if k not in ['nf', 'pratio'])
est.set_params(**params)
# est.loss__ = NAUC(2)
try:
est.fit(X_trn, y_trn)
y_val_est = est.predict_proba(X_val)[:,1]
y_test_est = est.predict_proba(X_test)[:,1]
lock = LockFile('.lock')
with lock:
with open('../data-cache/validate.spkl','ab') as fp:
# in a later stage we will use the OOB samples to make predictions on all samples
# so keep a record of the index of each OOB sample
pickle.dump((sample_validate,y_val_est), fp, -1)
with open('../data-cache/%s_test.spkl'%target,'ab') as fp:
# in a later stage we will use the OOB samples to make predictions on all samples
# so keep a record of the index of each OOB sample
pickle.dump(y_test_est, fp, -1)
except Exception as e:
return e
from sklearn.metrics import roc_auc_score
# (their p_ratio will be different) so this error measure is not completely accurate
return roc_auc_score(y_val, y_val_est)
In [65]:
!rm ../data-cache/validate.spkl
!rm ../data-cache/{target}_test.spkl
In [66]:
from collections import defaultdict
pid2args_list = defaultdict(list)
for res in hyeropt_results:
validation_score = res[0]
test_score = res[1]
pid = res[2]
args = res[3]
pid2args_list[pid].append((validation_score, args))
In [67]:
args_list = []
for v in pid2args_list.values():
print v[1]
for vv in v[:4]:
args_list.append(vv[1])
In [68]:
args_list
Out[68]:
In [69]:
results = lv.map(predict_work, args_list*Ncores)
In [70]:
import IPython
itr = results.__iter__()
while True:
try:
r = itr.next()
except StopIteration:
print 'stopped'
break
except IPython.parallel.error.RemoteError as e:
print e
continue
except Exception as e:
print e.__class__
continue
print r
In [71]:
fp = open('../data-cache/validate.spkl','rb')
count = 0
y_est = np.zeros(N)
y_count = np.zeros(N)
vals = []
while True:
try:
sample_validate,y_val_est = pickle.load(fp)
except:
break
count += 1
y_est[sample_validate] += y_val_est
y_count[sample_validate] += 1
idx = y_val_est.argsort()[::-1]
n = len(y_val_est)
val_recall_support = np.zeros(n)
p_sum = 0.
for i,j in enumerate(idx):
p_sum += float(y[sample_validate[j]])
val_recall_support[i] = p_sum
val_x = np.linspace(0.,100.,n)
vals.append((val_x, val_recall_support))
y_est /= y_count
In [72]:
np.sum(y_count == 0)
Out[72]:
In [73]:
from sklearn.metrics import roc_auc_score
y_no_overlap = [r for r,l in zip(y, latencies) if abs(l-int(l)) < 0.01]
y_est_no_overlap = [r for r,l in zip(y_est, latencies) if abs(l-int(l)) < 0.01]
roc_auc_score(y_no_overlap, y_est_no_overlap)
Out[73]:
In [74]:
with open('../data-cache/%s_predict.spkl'%target,'wb') as fp:
pickle.dump((y_no_overlap, y_est_no_overlap), fp, -1)
after running the entire notebook, again and again, on all targets - continue to 141107-GBC-combine
In [ ]: