{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\nimport os\nimport random\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\n\n# (try to) use a GPU for computation?\nuse_cuda=True\nif use_cuda and torch.cuda.is_available():\n  mydevice=torch.device('cuda')\nelse:\n  mydevice=torch.device('cpu')\n\n# training .npy files base dir\ntrain_dir='../input/g2net-gravitational-wave-detection/train'\n# location to store preprocessed file\ncache_dir='/kaggle/temp/preprocessed/'\ntrain_file='../input/g2net-gravitational-wave-detection/training_labels.csv'\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-16T21:13:07.598727Z","iopub.execute_input":"2021-07-16T21:13:07.599358Z","iopub.status.idle":"2021-07-16T21:13:09.118513Z","shell.execute_reply.started":"2021-07-16T21:13:07.599244Z","shell.execute_reply":"2021-07-16T21:13:09.117359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use the Karhunen Loeve transform, bordered autocorrelation method [Maccone, 2012](https://www.researchgate.net/publication/301171057_A_simple_introduction_to_the_KLT_and_BAM-KLT) to convert the time series data to be fed into a classifier.","metadata":{}},{"cell_type":"code","source":"# KLT BAM method\ndef klt_bam(infile, K, M):\n  # infile: input .npy file of time series\n  # K: how many largest eigenvalues ?\n  # M: max window size (lags) for correlation, must be > K (e.g. M=256+K)\n  # output: K eigenvalues of largest magnitude \n  # input size: ndim x N\n  # output size: M x K x ndim (removing zeros at beginning)\n  a=np.load(infile)\n  # get dimensions\n  ndim,N=a.shape\n\n  # normalize \n  for chan in range(ndim):\n    a[chan]=a[chan]/np.linalg.norm(a[chan])\n\n  assert(M>=K)\n\n  # storage \n  eigs=np.zeros((ndim,M,K))\n\n  # iterate over each channel\n  for chan in range(ndim):\n     x=a[chan]\n     # find covariance\n     # loop over window size\n     for m in range(1,M+1):\n       R=np.zeros((m,m))\n       # first row\n       r=np.zeros((m))\n       for ci in range(m):\n          r[ci]=np.dot(x[ci:N],x[0:N-ci])/(N-ci)\n\n       # fill matrix\n       for ci in range(m):\n         R[ci,ci:m]=r[0:m-ci]\n         R[ci+1:m,ci]=r[1:m-ci].transpose()\n\n       # find eigenvalues (in ascending order)\n       d=np.flip(np.linalg.eigvalsh(R))\n       \n       # select K largest\n       if K>m:\n         eigs[chan,m-1,K-m:K]=d[:m]\n       else:\n         eigs[chan,m-1,:]=d[:K]\n\n  # calculate difference along the lag axis\n  eigs=np.diff(eigs,axis=1)\n  # take log(abs())\n  eigs=np.log(np.abs(eigs)+1e-9)\n  # select lags K to M\n  eigs=eigs[:,-(M-K):,:]\n  # make data zero mean unit variance\n  eigs=(eigs-np.mean(eigs))/(np.sqrt(np.var(eigs)))\n\n  return eigs\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-16T21:13:09.120077Z","iopub.execute_input":"2021-07-16T21:13:09.120443Z","iopub.status.idle":"2021-07-16T21:13:09.137421Z","shell.execute_reply.started":"2021-07-16T21:13:09.120406Z","shell.execute_reply":"2021-07-16T21:13:09.135952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use a 1D CNN classifier that reads in the rows .and the columns of the KLT-BAM output","metadata":{}},{"cell_type":"code","source":"# 1D CNN model (both rows/columns are fed separately)\nclass Class1D(nn.Module):\n   def __init__(self,M=256,K=12,channels=3):\n     super(Class1D,self).__init__()\n     self.K=K\n     self.M=M\n     self.channels=channels\n     self.xconv0=nn.Conv1d(channels,8,3,stride=2,padding=1)\n     self.xconv1=nn.Conv1d(8,12,3,stride=2,padding=1)\n     self.xconv2=nn.Conv1d(12,24,3,stride=2,padding=1)\n     self.xconv3=nn.Conv1d(24,48,3,stride=2,padding=1)\n     self.xconv4=nn.Conv1d(48,96,3,stride=2,padding=1)\n     self.yconv0=nn.Conv1d(channels,8,3,stride=2,padding=1)\n     self.yconv1=nn.Conv1d(8,12,3,stride=2,padding=1)\n     self.yconv2=nn.Conv1d(12,24,3,stride=2,padding=1)\n     self.yconv3=nn.Conv1d(24,48,3,stride=2,padding=1)\n     self.yconv4=nn.Conv1d(48,96,3,stride=2,padding=1)\n     self.yconv5=nn.Conv1d(96,120,3,stride=2,padding=1)\n     self.xpool0=nn.AvgPool1d(2,stride=2)\n     self.xpool1=nn.AvgPool1d(2,stride=2)\n     self.ypool0=nn.AvgPool1d(2,stride=2)\n     self.ypool1=nn.AvgPool1d(2,stride=2)\n\n     self.fc1=nn.Linear(96*24*2,2400)\n     self.fc2=nn.Linear(2400,120)\n     self.fc3=nn.Linear(120,1)\n\n\n   def forward(self,E):\n     # E: batch,channel,M,K\n     x=torch.flatten(E,start_dim=2,end_dim=3)\n     x=self.xpool0(F.elu(self.xconv0(x)))\n     x=self.xpool1(F.elu(self.xconv1(x)))\n     x=(F.elu(self.xconv2(x)))\n     x=(F.elu(self.xconv3(x)))\n     x=(F.elu(self.xconv4(x)))\n     y=torch.flatten(torch.transpose(E,2,3),start_dim=2,end_dim=3)\n     y=self.ypool0(F.elu(self.yconv0(y)))\n     y=self.ypool1(F.elu(self.yconv1(y)))\n     y=(F.elu(self.yconv2(y)))\n     y=(F.elu(self.yconv3(y)))\n     y=(F.elu(self.yconv4(y)))\n     x=torch.flatten(x,start_dim=1)\n     y=torch.flatten(y,start_dim=1)\n     z=torch.cat((x,y),1)\n     z=F.elu(self.fc1(z))\n     z=F.elu(self.fc2(z))\n     z=torch.sigmoid(self.fc3(z))\n     return z","metadata":{"execution":{"iopub.status.busy":"2021-07-16T21:13:09.139340Z","iopub.execute_input":"2021-07-16T21:13:09.139847Z","iopub.status.idle":"2021-07-16T21:13:09.163135Z","shell.execute_reply.started":"2021-07-16T21:13:09.139796Z","shell.execute_reply":"2021-07-16T21:13:09.161823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The training loop follows below, it randomly selects a minibatch of time series data, calculates the KLT for that minibatch and feeds it into the CNN classifier. It also saves the pre-processed data in a cache, to avoid calculating the KLT again if that time series is used again.","metadata":{}},{"cell_type":"code","source":"# how many epochs\nNepoch=800\nbatchsize=40\n# how many largest eigenvalues\nK=12\n# max window size for correlation (excluding K)\nM=256\n\nnet=Class1D(M,K).to(mydevice)\n\nload_model=False\nsave_model=True\n# update from a saved model \nif load_model:\n  checkpoint=torch.load('./net.model',map_location=mydevice)\n  net.load_state_dict(checkpoint['model_state_dict'])\n  net.train() # initialize for training (BN,dropout)\n\n# get available file size\ntrain_df=pd.read_csv(train_file)\nnfiles,ncols=train_df.shape\nassert(ncols==2)\ndel train_df\n","metadata":{"execution":{"iopub.status.busy":"2021-07-16T21:13:09.164938Z","iopub.execute_input":"2021-07-16T21:13:09.165243Z","iopub.status.idle":"2021-07-16T21:13:09.855455Z","shell.execute_reply.started":"2021-07-16T21:13:09.165215Z","shell.execute_reply":"2021-07-16T21:13:09.854467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch.optim as optim\ncriterion=torch.nn.BCELoss()\noptimizer=optim.Adam(net.parameters(), lr=0.0001)\n","metadata":{"execution":{"iopub.status.busy":"2021-07-16T21:13:09.859589Z","iopub.execute_input":"2021-07-16T21:13:09.861124Z","iopub.status.idle":"2021-07-16T21:13:09.868400Z","shell.execute_reply.started":"2021-07-16T21:13:09.861060Z","shell.execute_reply":"2021-07-16T21:13:09.866989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for epoch in range(Nepoch):\n  ######################### data generation\n  # randomly select a subset of rows equal to batch size\n  skip=sorted(random.sample(range(nfiles),nfiles-batchsize))\n  train_df=pd.read_csv(train_file,skiprows=skip)\n\n  train_batch=train_df.to_numpy()\n  filenames=train_batch[:,0]\n  filelabels=train_batch[:,1]\n\n  X=torch.zeros((batchsize,3,M,K)).to(mydevice)\n  labels=torch.zeros((batchsize,1)).to(mydevice)\n  nbatch=0\n  for fname,flab in zip(filenames,filelabels):\n    labels[nbatch,0]=float(flab)\n    fullname=train_dir+os.sep+fname[0]+os.sep+fname[1]+os.sep+fname[2]+os.sep+fname+'.npy'\n    # preprocessed file\n    cachedfile=cache_dir+os.sep+fname[0]+os.sep+fname[1]+os.sep+fname[2]+os.sep+fname+'.npy'\n    fulldir=cache_dir+os.sep+fname[0]+os.sep+fname[1]+os.sep+fname[2]\n    # check if preprocessed file exists\n    if not os.path.exists(fulldir):\n      os.makedirs(fulldir)\n    if not os.path.isfile(cachedfile):\n      # preprocess file\n      E=klt_bam(fullname,K,M+K)\n      # create file\n      np.save(cachedfile,E)\n    else:\n      # now we have a file that is preprocessed\n      # size MxKx3\n      E=np.load(cachedfile)\n\n    # save as image\n    #tv.utils.save_image(torch.from_numpy(E),'Img_'+fname+'.png',normalize=True)\n    X[nbatch]=torch.from_numpy(E).to(mydevice)\n    nbatch+=1\n\n\n  ############################# training phase\n  def closure():\n    if torch.is_grad_enabled():\n      optimizer.zero_grad()\n    y=net(X)\n    #print(y-labels)\n    loss=criterion(y,labels)\n    if loss.requires_grad:\n       loss.backward()\n       print('%d loss %f'%(epoch,loss.data.item()))\n    return loss\n\n\n  optimizer.step(closure)\n\n\n# save model (and other extra items)\nif save_model:\n   torch.save({\n            'model_state_dict':net.state_dict(),\n           },'./net.model')\n","metadata":{"execution":{"iopub.status.busy":"2021-07-16T21:13:09.869829Z","iopub.execute_input":"2021-07-16T21:13:09.870139Z","iopub.status.idle":"2021-07-16T21:21:24.413538Z","shell.execute_reply.started":"2021-07-16T21:13:09.870109Z","shell.execute_reply":"2021-07-16T21:21:24.411189Z"},"trusted":true},"execution_count":null,"outputs":[]}]}