Manually run this notebooks on all targets
In [1110]:
target = 'Patient_2'
In [1111]:
%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
In [1112]:
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 [1113]:
FEATURES = 'gen-8_medianwindow-bands2-usf-w60-b0.2-b4-b8-b12-b30-b70-0.1-0.5-0.9'
In [1114]:
LOCAL_HOSTNAME=!hostname
if LOCAL_HOSTNAME[0].startswith('ip-'):
!aws s3 cp s3://udikaggle/sezure/{FEATURES}.tgz data-cache/
!tar xfz data-cache/{FEATURES}.tgz
each target receive a different model
positive examples. The positive examles were upsampled (using gen_ictal=-8)
In [1115]:
pdata = read_data(target, 'preictal', FEATURES)
Np, NF = pdata.X.shape
Np, NF
Out[1115]:
negative examples
In [1116]:
ndata = read_data(target, 'interictal', FEATURES)
Nn, NF = ndata.X.shape
Nn, NF
Out[1116]:
data is broken into segments that should be taken together when splitting to training and validation
In [1117]:
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
In [1118]:
psegments = getsegments(pdata)
Nps = len(psegments)
nsegments = getsegments(ndata)
Nns = len(nsegments)
statistics
In [1119]:
npratio = float(Nn)/Np
print target,1/(1+npratio),Np,Nn
npsratio = float(Nns)/Nps
print target,1/(1+npsratio),Nps,Nns
Ntrainps = 1
Ntrainns = int(Ntrainps*npsratio)
In [1120]:
X = np.concatenate((pdata.X, ndata.X))
y = np.zeros(X.shape[0])
y[:pdata.X.shape[0]] = 1
nsegments = [[s+Np for s in ns] for ns in nsegments]
In [1121]:
latencies = np.concatenate((pdata.latencies,ndata.latencies))
In [1122]:
N, NF = X.shape
N, NF
Out[1122]:
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 [1123]:
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=1000, n_jobs=-1, oob_score=True, max_depth=2)
rf.fit(X, y)
Out[1123]:
In [1124]:
rf.oob_score_
Out[1124]:
In [1125]:
pnweights = rf.feature_importances_
pn_importance_order = pnweights.argsort()[::-1]
pl.plot(rf.feature_importances_[pn_importance_order])
Out[1125]:
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 [1126]:
!rm /tmp/X.mmap
mmX = np.memmap('/tmp/X.mmap', shape=X.shape, dtype=np.float32, mode='w+')
mmX[:,:] = X
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 [1127]:
#!sleep 30
In [1128]:
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[1128]:
copy some information to all engines
In [1129]:
clients['X_shape'] = X.shape
clients['y'] = y
clients['psegments'] = psegments
clients['nsegments'] = nsegments
load the shared memory on all engines
In [1130]:
%%px
import numpy as np
N, NF = X_shape
X = np.memmap('/tmp/X.mmap', shape=X_shape, dtype=np.float32, mode='r')
In [1131]:
%%px
import random, itertools
def random_train_validation_split(psegments=psegments, nsegments=nsegments, N=N):
"""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 = np.array(list(itertools.chain(*ns))) # .ravel does not work - elements of nsegments are not of equal length
sample_validate = np.concatenate((s, n))
return np.setdiff1d(xrange(N),sample_validate,True), sample_validate
We will optimize AUC
In [1132]:
from sklearn.metrics import roc_auc_score
myscore = roc_auc_score
clients['myscore'] = myscore
In [1133]:
%%px
def hyperopt_work(args):
from lockfile import LockFile
sample_train, sample_validate = random_train_validation_split()
X_trn = X[sample_train,:]
y_trn = y[sample_train]
X_val = X[sample_validate,:]
y_val = y[sample_validate]
def t_est(params):
try:
args = {}
for h, p in params.iteritems():
if h.startswith('I'):
args[h[1:]] = int(p)
else:
args[h] = p
from sklearn.ensemble import GradientBoostingClassifier
est = GradientBoostingClassifier()
est.set_params(**dict((k,v) for k,v in args.iteritems() if k not in ['nf']))
nf = args['nf']
est.fit(X_trn[:,:nf], y_trn)
y_train = est.predict_proba(X_trn[:,:nf])
y_validate = est.predict_proba(X_val[:,:nf])
#from sklearn.metrics import average_precision_score
train_score = myscore(y_trn, y_train[:,1])
validate_score = myscore(y_val, y_validate[:,1])
lock = LockFile('.lock')
with lock:
with open('/tmp/hyperopt.txt','a') as fp:
print >>fp, validate_score, train_score, args
from hyperopt import STATUS_OK
return {'loss':1.-validate_score, 'status':STATUS_OK, 'train_score':train_score, 'validate_score':validate_score}
except:
lock = LockFile('.lock')
with lock:
with open('/tmp/hyperopt.txt','a') as fp:
print >>fp, 'failed', args
from hyperopt import STATUS_FAIL
return {'status':STATUS_FAIL, 'loss':1.} # 'loss' is mandatory
max_evals = args.get('max_evals', 10)
from hyperopt import fmin, tpe, Trials
space = args.get('space')
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('/tmp/hyperopt.spkl','ab') as fp:
pickle.dump(trials, fp, -1)
return best
Define statistical space in which we will do our hper-parameter search. An "I" at the start indicates that the value is an int
In [1134]:
%%px
from hyperopt import hp
from math import log
space = {
'n_estimators': 100,
'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(1.)),
'Inf': hp.quniform( 'nf', 10, NF, 1),
'Imax_features': hp.quniform( 'max_features', 10, 30, 1),
'Imax_depth': hp.quniform( 'max_depth', 2, 18, 1),
'Imin_samples_leaf': hp.quniform( 'min_samples_leaf', 2, 30, 1),
'subsample': hp.uniform( 'subsample', 0.2, 0.9),
}
In [1135]:
!rm /tmp/hyperopt.*
run hyperopt searches in parallel on all cores. Each hyperopt search will do 100 evaluations of the hyper parameters.
In [1136]:
%%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 [1137]:
!sort -n -r /tmp/hyperopt.txt > /tmp/hyperopt.sort.txt
!head -n 5 /tmp/hyperopt.sort.txt
In [1138]:
fp = open('/tmp/hyperopt.sort.txt')
hyeropt_results = []
for l in fp:
if l.startswith('failed'):
continue
l = l.split()
validate_score = l[0]
train_score = l[1]
args = eval(''.join(l[2:]))
hyeropt_results.append((validate_score, train_score, args))
fp.close()
len(hyeropt_results)
Out[1138]:
In [1139]:
tdata = read_data(target, 'test', FEATURES)
Nt, NFF = tdata.X.shape
Nt, NFF
Out[1139]:
In [1140]:
!rm /tmp/Xt.mmap
Xt_shape = (tdata.X.shape[0], NF)
mmXt = np.memmap('/tmp/Xt.mmap', shape=Xt_shape, dtype=np.float32, mode='w+')
mmXt[:,:] = tdata.X
del mmXt # flush to disk
In [1141]:
clients['Xt_shape'] = Xt_shape
clients['target'] = target
In [1142]:
%%px
Xt = np.memmap('/tmp/Xt.mmap', shape=Xt_shape, dtype=np.float32, mode='r')
In [1143]:
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]))
# while True:
# use out-of-bag samples to estimate the generalization error
sample_train, sample_validate = random_train_validation_split()
X_trn = X[sample_train,:NF]
y_trn = y[sample_train]
# if np.any(y_trn == 1) and np.any(y_trn == 0):
# break
X_val = X[sample_validate,:NF]
y_val = y[sample_validate]
X_test = Xt[:,:NF]
from sklearn.ensemble import GradientBoostingClassifier
est = GradientBoostingClassifier(n_estimators=1000)
est.set_params(**dict((k,v) for k, v in args.iteritems() if k not in ['nf']))
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('/tmp/predict.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('/tmp/%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
# (their p_ratio will be different) so this error measure is not completely accurate
return myscore(y_val, y_val_est)
In [1144]:
!rm /tmp/predict.spkl
!rm /tmp/{target}_test.spkl
In [1145]:
args_list = []
for res in hyeropt_results:
params = res[2]
args = {}
for k,v in params.iteritems():
if k.startswith('I'):
args[k[1:]] = int(v)
else:
args[k] = v
print args
args_list.append(args)
if len(args_list) >= 16:
break
len(args_list)
Out[1145]:
In [1146]:
results = lv.map(predict_work, args_list*Ncores)
In [1147]:
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 [1148]:
fp = open('/tmp/predict.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
y_est[y_count == 0] = y.mean()
In [1149]:
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]
myscore(y_no_overlap, y_est_no_overlap)
Out[1149]:
In [1150]:
with open('/tmp/%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 140926-GBC-combine
In [1150]: