{"cells":[{"metadata":{},"cell_type":"markdown","source":"A light weight time to event regression using differintegrated features."},{"metadata":{},"cell_type":"markdown","source":"This is an approximation of the workflow that I followed to develop the final submission for LANL Earthquake prediction challenge. Sadly the best submission that I have was not the one that I selected for the final submission, however, the performance was really good. Hope it helps to develop new and more accurate models."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport differint.differint as df\n\nfrom sklearn import preprocessing as pr\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.model_selection import train_test_split\n\nfrom sklearn.ensemble import RandomForestRegressor\n\nimport csv\nfrom os import listdir\nfrom os.path import isfile, join","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Several functions used along with the kernel"},{"metadata":{"trusted":true},"cell_type":"code","source":"#General plot style\ndef PlotStyle(Axes,Title,x_label,y_label):\n    \n    Axes.spines['top'].set_visible(False)\n    Axes.spines['right'].set_visible(False)\n    Axes.spines['bottom'].set_visible(True)\n    Axes.spines['left'].set_visible(True)\n    Axes.xaxis.set_tick_params(labelsize=12)\n    Axes.yaxis.set_tick_params(labelsize=12)\n    Axes.set_ylabel(y_label,fontsize=14)\n    Axes.set_xlabel(x_label,fontsize=14)\n    Axes.set_title(Title)\n\ndef MinimalLoader(filename, delimiter=',', dtype=float):\n  \n  \"\"\"\n  modified from SO answer by Joe Kington\n  \"\"\"\n  def IterFunc():\n    with open(filename, 'r') as infile:\n      for line in infile:\n        line = line.rstrip().split(delimiter)\n        for item in line:\n          yield dtype(item)\n    MinimalLoader.rowlength = len(line)\n\n  data = np.fromiter(IterFunc(), dtype=dtype)\n  data = data.reshape((-1, MinimalLoader.rowlength))\n\n  return data\n\n#Mean and standard deviation of a time series \ndef GetScalerParameters(TimeSeries):\n  return np.mean(TimeSeries),np.std(TimeSeries)\n\n#Generates a Zero mean and unit variance signal \ndef MakeScaledSeries(Signal,MeanValue,StdValue):\n  StandardSignal=[(val-MeanValue)/StdValue for val in Signal]\n  return StandardSignal\n\n#Makes a matrix of time series samples \ndef MakeSamplesMatrix(TimeSeries,TimeToFailure,FragmentSize,delay):\n    \n  \"\"\"\n  TimeSeries -> Data to be sampled\n  TimeToFailure-> Time to failure \n  FragmentSize-> Size of the time series sample\n  delay-> Number of steps to wait to get a new sample, manages the amount of overlay between samples\n  \"\"\"\n  \n  cData=TimeSeries\n  cTim=TimeToFailure\n  cFrag=FragmentSize\n  container=[]\n  time=[]\n  nData=len(cData)\n  counter=0\n  \n  for k in range(nData-cFrag):\n    \n    if counter==delay:\n      \n      cSample=list(cData[k:k+cFrag])\n      container.append(cSample)\n      time.append(cTim[k+cFrag])\n      counter=0\n      \n    else:\n      counter=counter+1\n\n  return np.array(container),np.array(time)\n\n#Data features\ndef MakeFeaturesRF(DataMatrix):\n  \n  cont=[]\n  featuresList=[np.mean,np.std,np.ptp,np.min,np.sum]\n  \n  for sample in DataMatrix:\n    sampCont=[]\n    for feature in featuresList:\n      sampCont.append(feature(sample))\n      sampCont.append(np.abs(feature(sample)))\n    \n    sampCont.append(np.mean(np.diff(sample)))\n    sampCont.append(np.mean(np.abs(np.diff(sample))))\n    cont.append(sampCont)\n    \n  return np.array(cont)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Data=MinimalLoader(r'../input/lanl15/train15.csv',delimiter=',')\n\nAcData=Data[:,0]\nTimeData=Data[:,1]\n\nGlobalMean,GlobalStd=GetScalerParameters(AcData)\nScaledData=MakeScaledSeries(AcData,GlobalMean,GlobalStd)\n\ndel Data,AcData","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Consider the following random forest regression "},{"metadata":{"trusted":true},"cell_type":"code","source":"SamplesData,TimeTE=MakeSamplesMatrix(ScaledData,TimeData,10000,10000)\nFeaturesData=MakeFeaturesRF(SamplesData)\nRFScaler=pr.MinMaxScaler()\nRFScaler.fit(FeaturesData)\nFeaturesData=RFScaler.transform(FeaturesData)\n\nRFR=RandomForestRegressor(n_estimators=100)\nRFR.fit(FeaturesData,TimeTE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the following feature importances"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(1)\nplt.plot(RFR.feature_importances_)\nax=plt.gca()\nPlotStyle(ax,'','Features','Importance')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see that the first and the third most important features have what it looks like a linear correlation with the time to event, and the second one returns a central tendency of the samples. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(2,figsize=(15,5))\nplt.subplot(131)\nplt.plot(FeaturesData[:,-1])\nax=plt.gca()\nPlotStyle(ax,'','Samples','Abs Sum Of Changes')\nplt.subplot(132)\nplt.plot(FeaturesData[:,2])\nax=plt.gca()\nPlotStyle(ax,'','Samples','Standard Deviation')\nplt.subplot(133)\nplt.plot(FeaturesData[:,-2])\nax=plt.gca()\nPlotStyle(ax,'','Samples','Mean')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Both features, standard deviation and the absolute sum of changes, measure the dissimilarity between consecutive points or to the sample mean. That difference could be seen as an approximation of the first derivative of the sample. However, a full integer differentiation could lead to an excessive loss of information.\n\nWe can see that by plotting the autocorrelation as a measure of information/memory, a full integer differentiation removes almost all the correlation inside the data. However, we can perform a differintegration between [-0.1,0.25] and retain most of the information."},{"metadata":{"trusted":true},"cell_type":"code","source":"Corr=[]\norder=[]\ndf1=ScaledData[0:1500000]\n\nfor d in np.linspace(-1,1,20): \n  df2=df.GL(d,df1,num_points=len(df1)) \n  df2=MakeScaledSeries(df2,df2.mean(),df2.std())\n  corr=np.corrcoef(df1,df2)[0,1] \n  Corr.append(corr)\n  order.append(d)\n\nplt.figure(3)\nplt.plot(order,Corr)\nax=plt.gca()\nPlotStyle(ax,'','Differintegration order','Correlation')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A visual inspection of the absolute sum of changes of the differintegrated features shows a slight change in the slope. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(4,figsize=(15,5))\n\nfor order in [-0.1,0,0.1,0.2]:\n    Der0=[np.sum(np.abs(df.GL(order,df1[k:k+10000],num_points=10000))) for k in range(300000,len(df1),10000)]\n    Der1=[np.mean(df.GL(order,df1[k:k+10000],num_points=10000)) for k in range(300000,len(df1),10000)]\n    \n    plt.subplot(121)\n    plt.plot(Der0,label='Order ='+str(order))\n    plt.legend(loc=3)\n    ax=plt.gca()\n    PlotStyle(ax,'','Time','Abs Sum Of Changes')\n    \n    plt.subplot(122)\n    plt.plot(Der1,label='Order ='+str(order))\n    plt.legend(loc=3)\n    ax=plt.gca()\n    PlotStyle(ax,'','Time','Mean')\n    \ndel df1,Der0,FeaturesData,SamplesData","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The final feature engineering performed to the selected features was to take the logarithm of the absolute sum of changes, trying to get a more steep slope between earthquakes. Box-Cox and Yeo-Jhonson also were considered, but at least I was not able to get any performance from those transformations. "},{"metadata":{"trusted":true},"cell_type":"code","source":"SeriesFragment=10000\nDelay=int(0.1*SeriesFragment)\nDerOrders=np.linspace(-0.1, 0.25, 6, endpoint=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the features already selected the next problem to tackle is the highly unbalanced dataset, to equalize the samples, first, the original signal is resampled using an overlapping scheme, allowing a 90% overlap between samples. Then all the samples are shuffled and the time to failure is segmented in 1000 bins and only 15 samples are taken for each bin. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Location function \ndef GetSampleLoc(SampleTime,boundaries):\n  \n  \"\"\"\n  \n  Returns the bin index of a time to the next eartquake sample \n  \n  SampleTime: Time To the next eartquake sample\n  boundaries: list of the boundaries of the bined time to the next earquake distribution \n  \n  \"\"\"\n  \n  for k in range(len(boundaries)-1):\n      \n    if SampleTime>=boundaries[k] and SampleTime<=boundaries[k+1]:\n        \n      cLoc=k\n      break\n      \n  return cLoc\n\n#Equalizes the samples over the range of time to the next earthquake\ndef MakeEqualizedSamples(DataSamples,TimeSamples):\n  \n  \"\"\"\n  \n  DataSamples:  Matrix of size (SampleSize,NumberOfSamples), contains the time \n                series samples\n  Time Samples: Array of size (NumberOfSamples), contains the time to the next \n                earthquake\n  \n  \"\"\"\n  \n  cData=DataSamples\n  cTime=TimeSamples\n  nData=len(cTime)\n  nBins=1000\n  \n  cMin,cMax=np.min(cTime),np.max(cTime)\n  bins=np.linspace(cMin,cMax,num=nBins+1)\n  \n  SamplesCount=[0 for k in range(nBins)]\n  \n  Xcont=[]\n  Ycont=[]\n  \n  index=[k for k in range(len(cTime))]\n  np.random.shuffle(index)\n  \n  for k in range(nData):\n    \n    cXval=cData[index[k]]\n    cYval=cTime[index[k]]\n    \n    cLoc=GetSampleLoc(cYval,bins)\n    \n    if SamplesCount[cLoc]<=15:\n      \n      Xcont.append(list(cXval))\n      Ycont.append(cYval)\n      SamplesCount[cLoc]=SamplesCount[cLoc]+1\n      \n  return np.array(Xcont),np.array(Ycont)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Samples,Times=MakeSamplesMatrix(ScaledData,TimeData,SeriesFragment,Delay)\nSamplesE,TimesE=MakeEqualizedSamples(Samples,Times)\n\ndel Samples,Times","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(5,figsize=(15,5))\nplt.subplot(121)\nn, bins, patches=plt.hist(TimeTE,bins=1000)\nax=plt.gca()\nPlotStyle(ax,'Normal Sampling','','')\nplt.subplot(122)\nn, bins, patches=plt.hist(TimesE,bins=1000)\nax=plt.gca()\nPlotStyle(ax,'Random Sampling','','')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With the equalized samples the logarithm of the mean absolute sum of changes and the mean of the differintegrated features is calculated over the selected interval. And such interval is divided into six values, resulting in twelve features per sample. "},{"metadata":{"trusted":true},"cell_type":"code","source":"#Calculate the features for each sample \ndef CalculateFeatures(Sample,Orders):\n  \n  \"\"\"\n  Sample: Time series fragment\n  Orders: Array of non integer differentiation orders \n  \"\"\"\n\n  container=[]\n  nSample=len(Sample)\n  \n  for order in Orders:\n      \n    derSample=df.GL(order,Sample,num_points=nSample)\n    absSample=np.abs(derSample)\n\n    container.append(np.log(1+np.mean(absSample)))\n    container.append(np.mean(derSample))\n\n  return container\n\n#A brief description \ndef MakeDataMatrix(Samples,Orders):\n  \n  \"\"\"\n  Samples: Matrix of time series samples \n  Orders: Array of non integer differentiation orders\n  \"\"\"\n  \n  container=[]\n  \n  for samp in Samples:\n    \n    container.append(CalculateFeatures(samp,Orders))\n    \n  return np.array(container)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All the features are scaled and the data is divided into two sets, one for grid search cross-validation (90%) and the rest to visualize the performance of the best solution "},{"metadata":{"trusted":true},"cell_type":"code","source":"Xtrain0=MakeDataMatrix(SamplesE,DerOrders)\nToMinMax=pr.MinMaxScaler()\nToMinMax.fit(Xtrain0)\nMMData=ToMinMax.transform(Xtrain0)\n\nXtrain,Xtest,Ytrain,Ytest=train_test_split(MMData,TimesE, train_size = 0.9,test_size=0.1,shuffle=True)\n\ndel Xtrain0,MMData","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For the hyperparameter optimization a grid search is performed with 2 folds. "},{"metadata":{"trusted":true},"cell_type":"code","source":"params={'n_estimators':[10,100,150,200],\n        'max_depth':[2,4,8,16,32,None],\n        'min_samples_split':[0.1,0.5,1.0],\n        'min_samples_leaf':[1,2,4],\n        'bootstrap':[True,False]}\n\n\nRFR=RandomForestRegressor() \nFinalModel=GridSearchCV(RFR,params,cv=2,verbose=1,n_jobs=2)\nFinalModel.fit(Xtrain,Ytrain)\npreds4 = FinalModel.predict(Xtest)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Performance of the best model"},{"metadata":{"trusted":true},"cell_type":"code","source":"MAE='Mean Absolute Error = ' +str(sum(np.abs(preds4-Ytest))/len(Ytest))\n\nplt.figure(3)\nplt.plot(preds4,Ytest,'bo',alpha=0.15)\nplt.plot([0,17],[0,17],'r')\nplt.xlim([0,17])\nplt.ylim([0,17])\nax=plt.gca()\nPlotStyle(ax,MAE,'Predicted','Real')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Loading the test data "},{"metadata":{"trusted":true},"cell_type":"code","source":"TestSamples=np.genfromtxt(r'../input/tested/test.csv',delimiter=',')\nTestIds=np.genfromtxt(r'../input/tested/ord.csv',delimiter=',')\n\nSamplesFeatures=MakeDataMatrix(TestSamples,DerOrders)\nScaledTest=ToMinMax.transform(SamplesFeatures)\nfinal=FinalModel.predict(ScaledTest)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Saving the predictions"},{"metadata":{"trusted":true},"cell_type":"code","source":"PredictionDir=r'../predictions.csv'\nfirstRow=['seg_id','time_to_failure']\n\nwith open(PredictionDir,'w',newline='') as output:\n        \n    writer=csv.writer(output)\n    nData=len(final)\n    writer.writerow(firstRow)\n            \n    for k in range(nData):\n      cRow=[TestIds[k],final[k]]\n      writer.writerow(cRow)\n        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":""}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}