# from __future__ import exam_success
from __future__ import absolute_import 
from __future__ import print_function  

# Standard imports
%matplotlib inline
import os
import sklearn
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import random
import pandas as pd
import scipy.stats as stats

# Sk cheats
from sklearn.cross_validation import cross_val_score
from sklearn import grid_search
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import SVR
#from sklearn.preprocessing import Imputer   # get rid of nan
from sklearn.decomposition import NMF        # to add features based on the latent representation
from sklearn.decomposition import ProjectedGradientNMF

# Faster gradient boosting
import xgboost as xgb

# For neural networks models
from keras.models import Sequential
from keras.layers.core import Dense, Dropout, Activation 
from keras.optimizers import SGD, RMSprop

"""
Explanation: FDMS TME3
Kaggle How Much Did It Rain? II
Florian Toque & Paul Willot  
Notes
We tried different regressor model, like GBR, SVM, MLP, Random Forest and KNN as recommanded by the winning team of the Kaggle on taxi trajectories. So far GBR seems to be the best, slightly better than the RF.
The new features we exctracted only made a small impact on predictions but still improved them consistently.
We tried to use a LSTM to take advantage of the sequential structure of the data but it didn't work too well, probably because there is not enought data (13M lines divided per the average length of sequences (15), less the 30% of fully empty data)
End of explanation
"""


%%time
#filename = "data/train.csv"
#filename = "data/reduced_train_100000.csv"
filename = "data/reduced_train_5000000.csv"
raw = pd.read_csv(filename)
raw = raw.set_index('Id')

raw.columns

raw['Expected'].describe()

"""
Explanation: 13.765.202 lines in train.csv  
8.022.757 lines in test.csv  

Few words about the dataset
Predictions is made in the USA corn growing states (mainly Iowa, Illinois, Indiana) during the season with the highest rainfall (as illustrated by Iowa for the april to august months)
The Kaggle page indicate that the dataset have been shuffled, so working on a subset seems acceptable
The test set is not a extracted from the same data as the training set however, which make the evaluation trickier
Load the dataset
End of explanation
"""


# Considering that the gauge may concentrate the rainfall, we set the cap to 1000
# Comment this line to analyse the complete dataset 
l = len(raw)
raw = raw[raw['Expected'] < 300]  #1000
print("Dropped %d (%0.2f%%)"%(l-len(raw),(l-len(raw))/float(l)*100))

raw.head(5)

raw.describe()

"""
Explanation: Per wikipedia, a value of more than 421 mm/h is considered "Extreme/large hail"
If we encounter the value 327.40 meter per hour, we should probably start building Noah's ark
Therefor, it seems reasonable to drop values too large, considered as outliers
End of explanation
"""


# We select all features except for the minutes past,
# because we ignore the time repartition of the sequence for now

features_columns = list([u'Ref', u'Ref_5x5_10th',
       u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
       u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
       u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
       u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
       u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
       u'Kdp_5x5_50th', u'Kdp_5x5_90th'])

def getXy(raw):
    selected_columns = list([ u'minutes_past',u'radardist_km', u'Ref', u'Ref_5x5_10th',
       u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
       u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
       u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
       u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
       u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
       u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
    
    data = raw[selected_columns]
    
    docX, docY = [], []
    for i in data.index.unique():
        if isinstance(data.loc[i],pd.core.series.Series):
            m = [data.loc[i].as_matrix()]
            docX.append(m)
            docY.append(float(raw.loc[i]["Expected"]))
        else:
            m = data.loc[i].as_matrix()
            docX.append(m)
            docY.append(float(raw.loc[i][:1]["Expected"]))
    X , y = np.array(docX) , np.array(docY)
    return X,y

"""
Explanation: We regroup the data by ID
End of explanation
"""


#noAnyNan = raw.loc[raw[features_columns].dropna(how='any').index.unique()]
noAnyNan = raw.dropna()

noFullNan = raw.loc[raw[features_columns].dropna(how='all').index.unique()]

fullNan = raw.drop(raw[features_columns].dropna(how='all').index)

print(len(raw))
print(len(noAnyNan))
print(len(noFullNan))
print(len(fullNan))

"""
Explanation: On fully filled dataset
End of explanation
"""


%%time
#X,y=getXy(noAnyNan)
X,y=getXy(noFullNan)

#%%time
#XX = [np.array(t).mean(0) for t in X]
#XX = np.array([np.append(np.nanmean(np.array(t),0),(np.array(t)[1:] - np.array(t)[:-1]).sum(0) ) for t in X])

%%time
global_means = np.nanmean(noFullNan,0)

XX=[]
for t in X:
    nm = np.nanmean(t,0)
    for idx,j in enumerate(nm):
        if np.isnan(j):
            nm[idx]=global_means[idx]
    XX.append(nm)
XX=np.array(XX)

# rescale to clip min at 0 (for non negative matrix factorization)
XX_rescaled=XX[:,:]-np.min(XX,0)

%%time
nmf = NMF(max_iter=1000)
W = nmf.fit_transform(XX_rescaled)
#H = nn.components_

# used to fill fully empty datas
global_means = np.nanmean(noFullNan,0)

# reduce the sequence structure of the data and produce
# new hopefully informatives features
def addFeatures(X,mf=0):
    # used to fill fully empty datas
    #global_means = np.nanmean(X,0)
    
    XX=[]
    nbFeatures=float(len(X[0][0]))
    for idxt,t in enumerate(X):
        
        # compute means, ignoring nan when possible, marking it when fully filled with nan
        nm = np.nanmean(t,0)
        tt=[]
        for idx,j in enumerate(nm):
            if np.isnan(j):
                nm[idx]=global_means[idx]
                tt.append(1)
            else:
                tt.append(0)
        tmp = np.append(nm,np.append(tt,tt.count(0)/nbFeatures))
        
        # faster if working on fully filled data:
        #tmp = np.append(np.nanmean(np.array(t),0),(np.array(t)[1:] - np.array(t)[:-1]).sum(0) )
        
        # add the percentiles
        tmp = np.append(tmp,np.nanpercentile(t,10,axis=0))
        tmp = np.append(tmp,np.nanpercentile(t,50,axis=0))
        tmp = np.append(tmp,np.nanpercentile(t,90,axis=0))
        
        for idx,i in enumerate(tmp):
            if np.isnan(i):
                tmp[idx]=0

        # adding the dbz as a feature
        test = t
        try:
            taa=test[:,0]
        except TypeError:
            taa=[test[0][0]]
        valid_time = np.zeros_like(taa)
        valid_time[0] = taa[0]
        for n in xrange(1,len(taa)):
            valid_time[n] = taa[n] - taa[n-1]
        valid_time[-1] = valid_time[-1] + 60 - np.sum(valid_time)
        valid_time = valid_time / 60.0


        sum=0
        try:
            column_ref=test[:,2]
        except TypeError:
            column_ref=[test[0][2]]
        for dbz, hours in zip(column_ref, valid_time):
            # See: https://en.wikipedia.org/wiki/DBZ_(meteorology)
            if np.isfinite(dbz):
                mmperhr = pow(pow(10, dbz/10)/200, 0.625)
                sum = sum + mmperhr * hours
        
        if not(mf is 0):
            tmp = np.append(tmp,mf[idxt])

        XX.append(np.append(np.array(sum),tmp))
        #XX.append(np.array([sum]))
        #XX.append(tmp)
    return np.array(XX)

%%time
XX=addFeatures(X,mf=W)
#XX=addFeatures(X)

os.system('say "終わりだ"') #its over!

def splitTrainTest(X, y, split=0.2):
    tmp1, tmp2 = [], []
    ps = int(len(X) * (1-split))
    index_shuf = range(len(X))
    random.shuffle(index_shuf)
    for i in index_shuf:
        tmp1.append(X[i])
        tmp2.append(y[i])
    return tmp1[:ps], tmp2[:ps], tmp1[ps:], tmp2[ps:]

X_train,y_train, X_test, y_test = splitTrainTest(XX,y)

"""
Explanation: Predicitons
As a first try, we make predictions on the complete data, and return the 50th percentile and uncomplete and fully empty data
End of explanation
"""


def manualScorer(estimator, X, y):
    err = (estimator.predict(X_test)-y_test)**2
    return -err.sum()/len(err)

"""
Explanation: 
End of explanation
"""


svr = SVR(kernel='rbf', C=800.0)

%%time
srv = svr.fit(X_train,y_train)

print(svr.score(X_train,y_train))
print(svr.score(X_test,y_test))

err = (svr.predict(X_train)-y_train)**2
err.sum()/len(err)

err = (svr.predict(X_test)-y_test)**2
err.sum()/len(err)

%%time
svr_score = cross_val_score(svr, XX, y, cv=5)
print("Score: %s\nMean: %.03f"%(svr_score,svr_score.mean()))

"""
Explanation: max prof 24
nb trees 84
min sample per leaf 17
min sample to split 51
End of explanation
"""


knn = KNeighborsRegressor(n_neighbors=6,weights='distance',algorithm='ball_tree')

#parameters = {'weights':('distance','uniform'),'algorithm':('auto', 'ball_tree', 'kd_tree', 'brute')}
parameters = {'n_neighbors':range(1,10,1)}
grid_knn = grid_search.GridSearchCV(knn, parameters,scoring=manualScorer)

%%time
grid_knn.fit(X_train,y_train)

print(grid_knn.grid_scores_)
print("Best: ",grid_knn.best_params_)

knn = grid_knn.best_estimator_

%%time
knn= knn.fit(XX,y)

knn= knn.fit(X_train,y_train)

print(knn.score(X_train,y_train))
print(knn.score(X_test,y_test))

err = (knn.predict(X_train)-y_train)**2
err.sum()/len(err)

err = (knn.predict(X_test)-y_test)**2
err.sum()/len(err)

"""
Explanation: 
End of explanation
"""


etreg = ExtraTreesRegressor(n_estimators=200, max_depth=None, min_samples_split=1, random_state=0)

parameters = {'n_estimators':range(100,200,20)}
grid_rf = grid_search.GridSearchCV(etreg, parameters,n_jobs=2,scoring=manualScorer)

%%time
grid_rf.fit(X_train,y_train)

print(grid_rf.grid_scores_)
print("Best: ",grid_rf.best_params_)

grid_rf.best_params_

#etreg = grid_rf.best_estimator_

%%time
etreg = etreg.fit(XX,y)

%%time
etreg = etreg.fit(X_train,y_train)

print(etreg.score(X_train,y_train))
print(etreg.score(X_test,y_test))

err = (etreg.predict(X_train)-y_train)**2
err.sum()/len(err)

err = (etreg.predict(X_test)-y_test)**2
err.sum()/len(err)

"""
Explanation: 
End of explanation
"""


rfr = RandomForestRegressor(n_estimators=200, criterion='mse', max_depth=None, min_samples_split=2,
                            min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_features='auto',
                            max_leaf_nodes=None, bootstrap=True, oob_score=False, n_jobs=-1,
                            random_state=None, verbose=0, warm_start=False)

%%time
rfr = rfr.fit(XX,y)

%%time
rfr = rfr.fit(X_train,y_train)

print(rfr.score(X_train,y_train))
print(rfr.score(X_test,y_test))

"""
Explanation: 
End of explanation
"""


# the dbz feature does not influence xgbr so much
xgbr = xgb.XGBRegressor(max_depth=6, learning_rate=0.1, n_estimators=700, silent=True,
                        objective='reg:linear', nthread=-1, gamma=0, min_child_weight=1,
                        max_delta_step=0, subsample=1, colsample_bytree=1, colsample_bylevel=1,
                        reg_alpha=0, reg_lambda=1, scale_pos_weight=1, base_score=0.5,
                        seed=0, missing=None)

%%time
xgbr = xgbr.fit(XX,y)

%%time
xgbr = xgbr.fit(X_train,y_train)

# without the nmf features
# print(xgbr.score(X_train,y_train))
## 0.993948231144
# print(xgbr.score(X_test,y_test))
## 0.613931733332

# with nmf features
print(xgbr.score(X_train,y_train))
print(xgbr.score(X_test,y_test))

"""
Explanation: 
End of explanation
"""


gbr = GradientBoostingRegressor(loss='ls', learning_rate=0.1, n_estimators=900,
                                subsample=1.0, min_samples_split=2, min_samples_leaf=1, max_depth=4, init=None,
                                random_state=None, max_features=None, alpha=0.5,
                                verbose=0, max_leaf_nodes=None, warm_start=False)

%%time
gbr = gbr.fit(XX,y)
#os.system('say "終わりだ"') #its over!

os.system('say "終わりだ"') #its over!

%%time
gbr = gbr.fit(X_train,y_train)
#os.system('say "終わりだ"') #its over!

#parameters = {'max_depth':range(2,5,1),'alpha':[0.5,0.6,0.7,0.8,0.9]}
#parameters = {'subsample':[0.2,0.4,0.5,0.6,0.8,1]}
#parameters = {'subsample':[0.2,0.5,0.6,0.8,1],'n_estimators':[800,1000,1200]}
#parameters = {'max_depth':range(2,4,1)}
parameters = {'n_estimators':[400,800,1100]}
#parameters = {'loss':['ls', 'lad', 'huber', 'quantile'],'alpha':[0.3,0.5,0.8,0.9]}
#parameters = {'learning_rate':[0.1,0.5,0.9]}



grid_gbr = grid_search.GridSearchCV(gbr, parameters,n_jobs=2,scoring=manualScorer)

%%time
grid_gbr = grid_gbr.fit(X_train,y_train)

print(grid_gbr.grid_scores_)
print("Best: ",grid_gbr.best_params_)

print(gbr.score(X_train,y_train))
print(gbr.score(X_test,y_test))

err = (gbr.predict(X_train)-y_train)**2
print(err.sum()/len(err))
err = (gbr.predict(X_test)-y_test)**2
print(err.sum()/len(err))

err = (gbr.predict(X_train)-y_train)**2
print(err.sum()/len(err))
err = (gbr.predict(X_test)-y_test)**2
print(err.sum()/len(err))

"""
Explanation: 
End of explanation
"""


t = []
for i in XX:
    t.append(np.count_nonzero(~np.isnan(i)) / float(i.size))
pd.DataFrame(np.array(t)).describe()

import pickle

with open('data/knn.pkl', 'wb') as fid:
    pickle.dump(knn, fid)  
with open('data/etreg.pkl', 'wb') as fid:
    pickle.dump(etreg, fid) 
with open('data/rfr.pkl', 'wb') as fid:
    pickle.dump(rfr, fid)
with open('data/xgbr.pkl', 'wb') as fid:
    pickle.dump(xgbr, fid)
with open('data/gbr.pkl', 'wb') as fid:
    pickle.dump(gbr, fid) 

# load it again
#with open('my_dumped_classifier.pkl', 'rb') as fid:
#    gnb_loaded = pickle.load(fid)

#knn,etreg,rfr,xgbr,gbr

"""
Explanation: 
End of explanation
"""


svr.predict(X_test)

s = modelList[0]

t.mean(1)

modelList = [svr,knn,etreg,rfr,xgbr,gbr]
reducedModelList = [knn,etreg,xgbr,gbr]

score_train = [[str(f).split("(")[0],f.score(X_train,y_train)] for f in modelList]
score_test = [[str(f).split("(")[0],f.score(X_test,y_test)] for f in modelList]
for idx,i in enumerate(score_train):
    print(i[0])
    print("   train: %.03f"%i[1])
    print("   test:  %.03f"%score_test[idx][1])


#reducedModelList = [knn,etreg,xgbr,gbr]
globalPred = np.array([f.predict(XX) for f in reducedModelList]).T
#globalPred.mean(1)

globalPred[0]

y[0]

err = (globalPred.mean(1)-y)**2
print(err.sum()/len(err))

err = (globalPred.mean(1)-y)**2
print(err.sum()/len(err))

for f in modelList:
    print(str(f).split("(")[0])
    err = (f.predict(XX)-y)**2
    print(err.sum()/len(err))

err = (XX[:,0]-y)**2
print(err.sum()/len(err))

for f in modelList:
    print(str(f).split("(")[0])
    print(f.score(XX,y))

XX[:10,0] # feature 0 is marshall-palmer

svrMeta = SVR()

%%time
svrMeta = svrMeta.fit(globalPred,y)

err = (svrMeta.predict(globalPred)-y)**2
print(err.sum()/len(err))

"""
Explanation: 
End of explanation
"""




in_dim = len(XX[0])
out_dim = 1  

model = Sequential()
# Dense(64) is a fully-connected layer with 64 hidden units.
# in the first layer, you must specify the expected input data shape:
# here, 20-dimensional vectors.
model.add(Dense(128, input_shape=(in_dim,)))
model.add(Activation('tanh'))
model.add(Dropout(0.5))
model.add(Dense(1, init='uniform'))
model.add(Activation('linear'))

#sgd = SGD(lr=0.1, decay=1e-6, momentum=0.9, nesterov=True)
#model.compile(loss='mean_squared_error', optimizer=sgd)

rms = RMSprop()
model.compile(loss='mean_squared_error', optimizer=rms)

#model.fit(X_train, y_train, nb_epoch=20, batch_size=16)
#score = model.evaluate(X_test, y_test, batch_size=16)

prep = []
for i in y_train:
    prep.append(min(i,20))

prep=np.array(prep)
mi,ma = prep.min(),prep.max()
fy = (prep-mi) / (ma-mi)
#my = fy.max()
#fy = fy/fy.max()

model.fit(np.array(X_train), fy, batch_size=10, nb_epoch=10, validation_split=0.1)  

pred = model.predict(np.array(X_test))*ma+mi

err = (pred-y_test)**2
err.sum()/len(err)

r = random.randrange(len(X_train))
print("(Train) Prediction %0.4f, True: %0.4f"%(model.predict(np.array([X_train[r]]))[0][0]*ma+mi,y_train[r]))

r = random.randrange(len(X_test))
print("(Test)  Prediction %0.4f, True: %0.4f"%(model.predict(np.array([X_test[r]]))[0][0]*ma+mi,y_test[r]))

"""
Explanation: Here for legacy
End of explanation
"""


%%time
#filename = "data/reduced_test_5000.csv"
filename = "data/test.csv"
test = pd.read_csv(filename)
test = test.set_index('Id')

features_columns = list([u'Ref', u'Ref_5x5_10th',
       u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
       u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
       u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
       u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
       u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
       u'Kdp_5x5_50th', u'Kdp_5x5_90th'])

def getX(raw):
    selected_columns = list([ u'minutes_past',u'radardist_km', u'Ref', u'Ref_5x5_10th',
       u'Ref_5x5_50th', u'Ref_5x5_90th', u'RefComposite',
       u'RefComposite_5x5_10th', u'RefComposite_5x5_50th',
       u'RefComposite_5x5_90th', u'RhoHV', u'RhoHV_5x5_10th',
       u'RhoHV_5x5_50th', u'RhoHV_5x5_90th', u'Zdr', u'Zdr_5x5_10th',
       u'Zdr_5x5_50th', u'Zdr_5x5_90th', u'Kdp', u'Kdp_5x5_10th',
       u'Kdp_5x5_50th', u'Kdp_5x5_90th'])
    
    data = raw[selected_columns]
    
    docX= []
    for i in data.index.unique():
        if isinstance(data.loc[i],pd.core.series.Series):
            m = [data.loc[i].as_matrix()]
            docX.append(m)
        else:
            m = data.loc[i].as_matrix()
            docX.append(m)
    X = np.array(docX)
    return X

#%%time
#X=getX(test)

#tmp = []
#for i in X:
#    tmp.append(len(i))
#tmp = np.array(tmp)
#sns.countplot(tmp,order=range(tmp.min(),tmp.max()+1))
#plt.title("Number of ID per number of observations\n(On test dataset)")
#plt.plot()

#testFull = test.dropna()
testNoFullNan = test.loc[test[features_columns].dropna(how='all').index.unique()]


%%time
test_X=getX(testNoFullNan)  # 1min
#XX = [np.array(t).mean(0) for t in X]  # 10s


test_XX=[]
for t in test_X:
    nm = np.nanmean(t,0)
    for idx,j in enumerate(nm):
        if np.isnan(j):
            nm[idx]=global_means[idx]
    test_XX.append(nm)
test_XX=np.array(test_XX)

# rescale to clip min at 0 (for non negative matrix factorization)
test_XX_rescaled=test_XX[:,:]-np.min(test_XX,0)


%%time
test_W = nmf.transform(test_XX_rescaled)


%%time
test_XX=addFeatures(test_X,mf=test_W)


#pd.DataFrame(xgbr.predict(XX)).describe()


%%time
reducedModelList = [knn,etreg,xgbr,gbr]
test_globalPred = np.array([f.predict(test_XX) for f in reducedModelList]).T
predTest = test_globalPred.mean(1)


predFull = zip(testNoFullNan.index.unique(),predTest)


testNan = test.drop(test[features_columns].dropna(how='all').index)


tmp = np.empty(len(testNan))
tmp.fill(0.445000)   # 50th percentile of full Nan dataset
predNan = zip(testNan.index.unique(),tmp)


pred = predFull + predNan


testLeft = test.drop(testNan.index.unique()).drop(testFull.index.unique())


tmp = np.empty(len(testLeft))
tmp.fill(1.27)   # 50th percentile of full Nan dataset
predLeft = zip(testLeft.index.unique(),tmp)


len(testFull.index.unique())


len(testNan.index.unique())

len(testLeft.index.unique())

pred = predFull + predNan + predLeft


pred.sort(key=lambda x: x[0], reverse=False)


#reducedModelList = [knn,etreg,xgbr,gbr]
#globalPred = np.array([f.predict(XX) for f in reducedModelList]).T
#globalPred.mean(1)

submission = pd.DataFrame(pred)
submission.columns = ["Id","Expected"]
submission.head()


submission.loc[submission['Expected']<0,'Expected'] = 0.445


submission.to_csv("submit5.csv",index=False)


print("ok")

filename = "data/sample_solution.csv"
sol = pd.read_csv(filename)


sol

ss = np.array(sol)

%%time
for a,b in predFull:
    ss[a-1][1]=b

ss

sub = pd.DataFrame(pred)
sub.columns = ["Id","Expected"]
sub.Id = sub.Id.astype(int)
sub.head()

sub.to_csv("submit3.csv",index=False)

"""
Explanation: Predict on testset
End of explanation
"""


