#IMPORTING NECESSARY LIBRARIES

import matplotlib.pyplot as plt #import matplotlib

from matplotlib.pyplot import figure #import figure from matplotlib.pyplot

import numpy as np #import numpy

import tensorflow as tf #import tensorflow

import tensorflow_datasets as tfds #import tensorflow_datasets

from tensorflow import keras #import keras

import pandas as pd #import pandas

from scipy.fftpack import irfft, rfft #import rfft and irfft functions

print("importing datasets")

#IMPORTING THE DATASET pure.xlsx

df = pd.read_excel('pure.xlsx',header=None) #import pure dataset

pure_acc = df.to_numpy()

#IMPORTING THE DATASET noisy.xlsx

df = pd.read_excel('noisy.xlsx',header=None)#import noisy dataset

noisy_acc = df.to_numpy() #convert to numpy array

#IMPORTING THE DATASET test_pure.xlsx

df = pd.read_excel('test_pure.xlsx',header=None)#import test_pure dataset

test_pure_acc = df.to_numpy()#convert to numpy array

#IMPORTING THE DATASET test_noisy.xlsx

df = pd.read_excel('test_noisy.xlsx',header=None)#import test_noisy dataset

test_noisy_acc = df.to_numpy()#convert to numpy array

print("RANDOMIZING ORDER OF DATASET TO AVOID BIAS")

#RANDOMIZING ORDER OF DATASET TO AVOID BIAS

n = np.linspace(0,noisy_acc.shape[0]-1,noisy_acc.shape[0],dtype=int) #create array from 0 to 700

np.random.shuffle(n) #shauffle the ARRAY

temp_noisy = noisy_acc #save noisy dataset in temporary array

temp_pure = pure_acc #save pure dataset in temporary array

for i in n:

  noisy_acc[i] = temp_noisy[n[i]]#randomizing noisy datasets

  pure_acc[i] = temp_pure[n[i]]#randomizing pure dataset

print("DATA PREPROCESSING, CONVERTING TRAINING DATA TO FREQUENCY DOMAIN")

#DATA PREPROCESSING, CONVERTING TRAINING DATA TO FREQUENCY DOMAIN

noisy_acc_freq = np.zeros((noisy_acc.shape[0],noisy_acc.shape[1]),np.float64) #create array of zeros of same shape as noisy_acc

for i in range(noisy_acc.shape[0]): #looping in noisy_acc ARRAY

  noisy_acc_freq[i] = rfft(noisy_acc[i]) #converting all signals to freq domin

pure_acc_freq = np.zeros((pure_acc.shape[0],pure_acc.shape[1]),np.float64) #create array of zeros of same shape as pure_acc

for i in range(pure_acc.shape[0]): #looping in pure_acc ARRAY

  pure_acc_freq[i] = rfft(pure_acc[i])#converting all signals to freq domin

noisy_acc_freq = noisy_acc_freq.reshape(noisy_acc_freq.shape[0],noisy_acc_freq.shape[1],1)#RESHAPING THE NOISY ACCELEROMETER ARRAY

print("MODELING THE CNN FOR DENOISING")

#MODELING THE CNN FOR DENOISING

# linear activations are used as unlike regression it is a pattern to pattern matching

model1 = keras.Sequential([ #CREATING sequential model

   keras.layers.ZeroPadding1D(padding=3), #zeropadding layer

   keras.layers.Conv1D(16, 7, strides=1, activation='linear'), #convolutional layer filter size=7, no of filters=16

   keras.layers.ZeroPadding1D(padding=8), #zeropadding layer

   keras.layers.Conv1D(32, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=32

   keras.layers.Conv1D(32, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=32

   keras.layers.Conv1D(32, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=32

   keras.layers.Conv1D(16, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=16

   keras.layers.Conv1D(16, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=16

   keras.layers.Conv1D(16, 3, strides=1, activation='linear'),#convolutional layer filter size=3, no of filters=16

   keras.layers.Flatten(), #flatten layer

   keras.layers.Dense(16, activation='linear'),#fully connected layer, size = 16

   keras.layers.Dense(pure_acc_freq.shape[1], activation=None)#fully connected layer size=size of signal

])

optim = tf.keras.optimizers.Adam(3e-4) #using adam optimizer, learning rate=3 x 10^-4

model1.compile(optimizer=optim,  loss = 'mse',metrics=[tf.keras.metrics.RootMeanSquaredError('rmse')]) #defining loss function and optimizer functions

model1.fit(noisy_acc_freq, pure_acc_freq, epochs=100, batch_size=16) #training model1 for 100 epochs on batch size of 16

print("MODELING THE ANN FOR DENOISING")

#MODELING THE ANN FOR DENOISING

# linear activations are used as unlike regression it is a pattern to pattern matching

model2 = keras.Sequential([#CREATING sequential model

   keras.layers.Flatten(),#flatten layer

   keras.layers.Dense(4096, activation='linear'), #fully connected layer, size = 4096

   keras.layers.Dense(8192, activation='linear'),#fully connected layer, size = 8192

   keras.layers.Dense(4096, activation='linear'),#fully connected layer, size = 4096

   keras.layers.Dense(2048, activation='linear'),#fully connected layer, size = 2048

   keras.layers.Dense(pure_acc_freq.shape[1], activation=None)#fully connected layer, size = size of signal

])

optim = tf.keras.optimizers.SGD(1e-3)#using SGD optimizer, learning rate=1 x 10^-3

#the momentum aspect of Adam caused it to spiral out of control

model2.compile(optimizer=optim, loss = 'mse', metrics=[tf.keras.metrics.RootMeanSquaredError('rmse')])#defining loss function and optimizer functions

model2.fit(noisy_acc_freq, pure_acc_freq, epochs=100, batch_size=12)#training model2 for 100 epochs on batch size 12

print("SAVING THE MODEL")

#SAVING THE MODEL

model1.save('cnn1.h5')#SAVING MODEL1

model2.save('ann1.h5')#SAVING MODEL2

new_model = tf.keras.models.load_model('cnn1.h5') #loading CNN under the variable new_model

print("COMPARING THE PURE, NOISY AND DENOISED SIGNALS USING MATPLOTLIB")

#COMPARING THE PURE, NOISY AND DENOISED SIGNALS USING MATPLOTLIB

z = test_noisy_acc[1] #z has signal no. 1038 from the test_noisy dataset.

z= rfft(z) # converting z to freq DOMAIN

z = z.reshape(1,z.shape[0],1) #reshaping z ARRAY

y_denoised = new_model.predict(z) #denoising z using new_model

y_denoised = irfft(y_denoised) #converting z back to time DOMAIN

figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k') #defining size and shape of the plot

x = np.linspace(start=0,stop=7,num=701) #creating array= {0,0.01,0.02....6.98,6.99,7}

y_noisy = test_noisy_acc[1]  #plotting the noisy SIGNAL

plt.plot(x,y_noisy)             # plotting NOISY SIGNAL

figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')#defining size and shape of the plot

plt.plot(x,pure_acc[1])      #plotting PURE SIGNAL

figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')#defining size and shape of the plot

plt.plot(x,y_denoised.reshape(701))     #plotting DENOISED SIGNAL

print("CALCULATING IDR BY CALCULATING RELATIVE DISPLACEMENTS BY DOUBLE INTEGRATION:pure")

#NOW THAT OUR ML ALGORITHM IS WORKING FINE, WE WILL MOVE ON TO EVALUATING IT USING [IO,LS,CP]

#CALCULATING IDR BY CALCULATING RELATIVE DISPLACEMENTS BY DOUBLE INTEGRATION. AND FINALLY CREATING [IO,LS,CP] ARRAY FOR FULL PURE #DATASET

time = 1/360 #time taken between 2 readings. Sampling rate = 360Hz

pure_classification = np.zeros((int(test_pure_acc.shape[0]*0.5),3),np.float64) #creatng array of zeros for classinfiying the pure signals

v = np.zeros(test_pure_acc.shape,np.float64) #velocity

disp = np.zeros(test_pure_acc.shape,np.float64) #displacement

floor_height = 2.75 #standard floor height in India = 2.75m

for i in range(0,test_pure_acc.shape[0]): #iterating over the test_pure acceleration signals

  for j in range(1,test_pure_acc.shape[1]):   #going over each sample

      v[i][j] = v[i][j-1] + (((test_pure_acc[i][j-1]+test_pure_acc[i][j])/2) * (time)) #integrating for velocity

for i in range(0,test_pure_acc.shape[0]):#iterating over the velocity

  for j in range(1,test_pure_acc.shape[1]):  #going over each sample

      disp[i][j] = disp[i][j-1] + (((v[i][j-1]+v[i][j])/2) * (time)) #integrating for displacement

for i in range(0,disp.shape[0],2): #going over adjacent displacement arrays

  idr = np.zeros(disp.shape[1],np.float64) #creating array of zeros for IDRs

  for j in range(disp.shape[1]): #going over adjacent displacement arrays

    idr[j] = ( np.abs(disp[i][j]-disp[i+1][j]) )/(floor_height) #calculating IDRs

  '''

 if idr < 0.007       => Immediate Occupancy

 if idr 0.007 to 0.05 => Life Safety

 if idr >0.05         => Collapse prevention

 '''

  scores=np.array([0,0,0]) #scores

  '''

 io_score=0th index

 ls_score=1st index

 cp_score=2nd index

 '''

  for k in range (idr.shape[0]):#calculating scores for function/array of IDRs

    if idr[k]<0.007:#if IDR<0.007

      scores[0]+=1 #increment scores[0] by 1

    elif idr[k]>0.05:#if IDR>0.05

      scores[2]+=1 #increment scores[2] by 1

    else:#if 0.007<IDR<0.05

      scores[1]+=1#increment scores[2] by 1

  #most severe score is considered for labeling the dataset

  if scores[2]>0:#if scores[2]>0:

    scores = [0,0,1]#scores becomes [0,0,1]

  elif scores[1]>0:#if scores[1]>0 and scores[2]=0:

    scores = [0,1,0]#scores becomes [0,1,0]

  else: #if scores[0]>0 and scores[1]=0 and scores[2]=0

    scores = [1,0,0]#scores becomes [1,0,0]

  pure_classification[int(i/2)]=scores #adding scores to pure_classification array

  scores=np.array([0,0,0]) #resetting scores array

#FOR NOISY CLASSIFICATION

print("CALCULATING IDR BY CALCULATING RELATIVE DISPLACEMENTS BY DOUBLE INTEGRATION:noisy")

time = 1/360 #time taken between 2 readings. Sampling rate = 360Hz

noisy_classification = np.zeros((int(test_noisy_acc.shape[0]*0.5),3),np.float64)#creatng array of zeros for classinfiying the noisy signals

v = np.zeros(test_noisy_acc.shape,np.float64) #velocity

disp = np.zeros(test_noisy_acc.shape,np.float64) #displacement

for i in range(0,test_noisy_acc.shape[0]):#iterating over the test_noisy acceleration signals

  for j in range(1,test_noisy_acc.shape[1]):  #going over each sample

      v[i][j] = v[i][j-1] + (((test_noisy_acc[i][j-1]+test_noisy_acc[i][j])/2) * (time)) #integrating for velocity

for i in range(0,test_noisy_acc.shape[0]):#iterating over the velocity signals

  for j in range(1,test_noisy_acc.shape[1]):#going over each sample

      disp[i][j] = disp[i][j-1] + (((v[i][j-1]+v[i][j])/2) * (time))#integrating for displacement

for i in range(0,disp.shape[0],2): #going over adjacent displacement arrays

  idr = np.zeros(disp.shape[1],np.float64)#creating array of zeros for IDRs

  for j in range(disp.shape[1]): #going over adjacent displacement arrays

    idr[j] = ( np.abs(disp[i][j]-disp[i+1][j]) )/(floor_height)#calculating IDRs

  '''

 if idr < 0.007       => Immediate Occupancy

 if idr 0.007 to 0.05 => Life Safety

 if idr >0.05         => Collapse prevention

 '''

  scores=np.array([0,0,0])#scores

  '''

 io_score=0th index

 ls_score=1st index

 cp_score=2nd index

 '''

  for k in range (idr.shape[0]):#calculating scores for function/array of IDRs

    if idr[k]<0.007:#if IDR<0.007

      scores[0]+=1 #increment scores[0] by 1

    elif idr[k]>0.05:#if IDR>0.05

      scores[2]+=1 #increment scores[2] by 1

    else:#if 0.007<IDR<0.05

      scores[1]+=1#increment scores[2] by 1

#most severe score is considered for labeling the dataset

  if scores[2]>0:#if scores[2]>0:

    scores = [0,0,1]#scores becomes [0,0,1]

  elif scores[1]>0:#if scores[1]>0 and scores[2]=0:

    scores = [0,1,0]#scores becomes [0,1,0]

  else: #if scores[0]>0 and scores[1]=0 and scores[2]=0

    scores = [1,0,0]#scores becomes [1,0,0]

  noisy_classification[int(i/2)]=scores#adding scores to noisy_classification array

  scores=np.array([0,0,0])#resetting scores array

#calculating accuracy without denoising.

b=0 # b=0

for i in range(pure_classification.shape[0]):#iterating over pure_classification

  if pure_classification[i][0] == noisy_classification[i][0] and pure_classification[i][1] == noisy_classification[i][1] and pure_classification[i][2] == noisy_classification[i][2] : #if pure and noisy classifications match

    b+=1 #increment b by 1

print("accuracy without denoising = ",end="") #priniting text

print(np.float64(b)/pure_classification.shape[0])#printing accuracy. divide b by the total number of classifications

print("denoising noisy test set")

#DENOISING TEST_NOISY_ACC

test_denoised_acc = np.zeros((test_noisy_acc.shape[0],test_noisy_acc.shape[1]),np.float64) # create empty array of zeros to store denoised signals of test_noisy

for i in range(test_noisy_acc.shape[0]): #iterate over test_noisy dataset

  z = rfft(test_noisy_acc[i]) #convert to freq DOMAIN

  z = z.reshape(1,701,1) #reshape the array containing the signal

  test_denoised_acc[i] = irfft(new_model.predict(z))#denoise the signal and convert it back to time domain an store it.

print("CALCULATING IDR BY CALCULATING RELATIVE DISPLACEMENTS BY DOUBLE INTEGRATION:denoised")

#GETTING CLASSIFICATION ARRAY FOR DENOISED SIGNALS

time = 1/360 #time taken between 2 readings. Sampling rate = 360Hz

denoised_classification = np.zeros((int(test_denoised_acc.shape[0]*0.5),3),np.float64)#creatng array of zeros for classinfiying the denoised signals

v = np.zeros(test_denoised_acc.shape,np.float64) #velocity

disp = np.zeros(test_denoised_acc.shape,np.float64) #displacement

for i in range(0,test_denoised_acc.shape[0]):#iterating over the test_denoised signals

  for j in range(1,test_denoised_acc.shape[1]):   #going over each sample

      v[i][j] = v[i][j-1] + (((test_denoised_acc[i][j-1]+test_denoised_acc[i][j])/2) * (time))#integrating for velocity

for i in range(0,test_denoised_acc.shape[0]):#iterating over the velocity signals

  for j in range(1,test_denoised_acc.shape[1]):   #going over each sample

      disp[i][j] = disp[i][j-1] + (((v[i][j-1]+v[i][j])/2) * (time))#integrating for displacement

for i in range(0,disp.shape[0],2):#going over adjacent displacement array

  idr = np.zeros(disp.shape[1],np.float64)#creating array of zeros for IDRs

  for j in range(disp.shape[1]):#going over adjacent displacement arrays

    idr[j] = ( np.abs(disp[i][j]-disp[i+1][j]) )/(floor_height)#calculating IDRs

  '''

 if idr < 0.007       => Immediate Occupancy

 if idr 0.007 to 0.05 => Life Safety

 if idr >0.05         => Collapse prevention

 '''

  scores=np.array([0,0,0])#scores

  '''

 io_score=0th index

 ls_score=1st index

 cp_score=2nd index

 '''

  for k in range (idr.shape[0]):#calculating scores for function/array of IDRs

    if idr[k]<0.007:#if IDR<0.007

      scores[0]+=1 #increment scores[0] by 1

    elif idr[k]>0.05:#if IDR>0.05

      scores[2]+=1 #increment scores[2] by 1

    else:#if 0.007<IDR<0.05

      scores[1]+=1#increment scores[2] by 1

#most severe score is considered for labeling the dataset

  if scores[2]>0:#if scores[2]>0:

    scores = [0,0,1]#scores becomes [0,0,1]

  elif scores[1]>0:#if scores[1]>0 and scores[2]=0:

    scores = [0,1,0]#scores becomes [0,1,0]

  else: #if scores[0]>0 and scores[1]=0 and scores[2]=0

    scores = [1,0,0]#scores becomes [1,0,0]

  denoised_classification[int(i/2)]=scores#adding scores to denoised_classification array

  scores=np.array([0,0,0])#resetting scores array

#calculating accuracy with denoising.

b=0 # b=0

for i in range(pure_classification.shape[0]):#iterating over pure_classification

  if pure_classification[i][0] == denoised_classification[i][0] and pure_classification[i][1] == denoised_classification[i][1] and pure_classification[i][2] == denoised_classification[i][2] : #if pure and denoised classifications match

    b+=1 #increment b by 1

print("accuracy with denoising = ",end="") #priniting text

print(np.float64(b)/pure_classification.shape[0])#printing accuracy. divide b by the total number of classifications

print("DIRECT CLASSIFICATION USING ACCELEROMETER DATA")

##DIRECT CLASSIFICATION USING ACCELEROMETER DATA

#RELATIVE ACCELERATION CALCULATION

rel_acc_noisy = np.zeros((int(noisy_acc.shape[0]*0.5),noisy_acc.shape[1]),dtype=np.float64) # Creating array to store relative accelerations between noisy data

rel_acc_pure = np.zeros((int(pure_acc.shape[0]*0.5),pure_acc.shape[1]),dtype=np.float64) # Creating array to store relative accelerations between pure data

a=0 #a is assigned to 0

for i in range(rel_acc_pure.shape[0]):

  rel_acc_noisy[i] = noisy_acc[a+1] - noisy_acc[a] #evaluating the relative accelerations between adjacent noisy_acc

  rel_acc_pure[i]  = pure_acc[a+1] - pure_acc[a] #evaluating the relative accelerations between adjacent pure_acc

  a+=2 #incrementing a by 2 to get the next 2 accelerations

rel_acc_noisy = rel_acc_noisy.reshape(rel_acc_noisy.shape[0],rel_acc_noisy.shape[1],1) #reshaping the array

print("Getting the pure classifications from pure_acc")

### Getting the pure classifications from pure_acc

#For relative acceleration calculation we assume 2 simultaneous accelerometer-time signals

time = 1/360 #time taken between 2 readings. Sampling rate = 360Hz

pure_classification2 = np.zeros((int(pure_acc.shape[0]*0.5),3),np.float64) #creatng array of zeros for classinfiying the pure training signals

v = np.zeros(pure_acc.shape,np.float64) #velocity

disp = np.zeros(pure_acc.shape,np.float64) #displacement

for i in range(0,pure_acc.shape[0]): #iterating over the training pure signals

  for j in range(1,pure_acc.shape[1]):  #going over each sample

      v[i][j] = v[i][j-1] + (((pure_acc[i][j-1]+pure_acc[i][j])/2) * (time)) #integrating for velocity

for i in range(0,pure_acc.shape[0]): #iterating over velocities

  for j in range(1,pure_acc.shape[1]):   #double integration.

      disp[i][j] = disp[i][j-1] + (((v[i][j-1]+v[i][j])/2) * (time))#integrating for displacement

for i in range(0,disp.shape[0],2):#going over adjacent displacement array

  idr = np.zeros(disp.shape[1],np.float64)#creating array of zeros for IDRs

  for j in range(disp.shape[1]):#going over adjacent displacement arrays

    idr[j] = ( np.abs(disp[i][j]-disp[i+1][j]) )/(floor_height)#calculating IDRs

  '''

 if idr < 0.007       => Immediate Occupancy

 if idr 0.007 to 0.05 => Life Safety

 if idr >0.05         => Collapse prevention

 '''

  scores=np.array([0,0,0])#scores

  '''

 io_score=0th index

 ls_score=1st index

 cp_score=2nd index

 '''

  for k in range (idr.shape[0]):#calculating scores for function/array of IDRs

    if idr[k]<0.007:#if IDR<0.007

      scores[0]+=1 #increment scores[0] by 1

    elif idr[k]>0.05:#if IDR>0.05

      scores[2]+=1 #increment scores[2] by 1

    else:#if 0.007<IDR<0.05

      scores[1]+=1#increment scores[2] by 1

  scores = np.floor(scores/(np.amax(scores))) #putting 1 in position that has maximum score and 0 in other positions

  pure_classification2[int(i/2)]=scores #adding scores to denoised_classification array

  scores=np.array([0,0,0])#resetting scores array

print("CNN FOR DIRECT CLASSIFICATION")

###CNN FOR DIRECT CLASSIFICATION

model3 = keras.Sequential([ #sequential model

   keras.layers.Conv1D(32, 32, strides=1, activation='relu'), #convolutional layer, kernal size=32, number of filters = 32

   keras.layers.Conv1D(32, 32, strides=1, activation='relu'),#convolutional layer, kernal size=32, number of filters = 32

   keras.layers.MaxPool1D(pool_size=8),#maxpool layer, pooling filtersize = 8

   keras.layers.Flatten(), #flatten layer

   keras.layers.Dense(512, activation='relu'),#fully connected layer, size = 512

   keras.layers.Dense(512, activation='relu'), #fully connected layer, size = 512

   keras.layers.Dense(3, activation='softmax')]) #fully connected layer, size = 3 to get output as one of the 3 classes

model3.compile(loss=keras.losses.categorical_crossentropy,optimizer=keras.optimizers.SGD(),metrics=['accuracy']) # creating model using SGD optimizer and catagorical cross-entropy loss function.

model3.fit(rel_acc_noisy, pure_classification2, epochs=5, batch_size=16) # training the model. 5 epochs , batch size of 16

print("ANN FOR DIRECT CLASSIFICATION")

#ANN FOR DIRECT CLASSIFICATION

model4 = keras.Sequential([ #sequential model

   keras.layers.Flatten(),#flatten layer

   keras.layers.Dense(512, activation='relu'),#fully connected layer, size = 512

   keras.layers.Dense(512, activation='relu'),#fully connected layer, size = 512

   keras.layers.Dense(512, activation='relu'),#fully connected layer, size = 512

   keras.layers.Dense(3, activation='softmax')])#fully connected layer, size = 3 to get output as one of the 3 classes

model4.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.SGD(),metrics=['accuracy']) # creating model using SGD optimizer and catagorical cross-entropy loss function.

model4.fit(rel_acc_noisy, pure_classification2, epochs=5, batch_size=4)# training the model. 5 epochs , batch size of 4

print("final accuracy testing on unseen data!")

e = np.zeros((int(test_noisy_acc.shape[0]*0.5),test_noisy_acc.shape[1]),dtype=np.float64) #creating empty array  for storing test_noisy dataset's relative accelerations.

#e is the relative acceleration for the test_noisy_acc

a=0 # a   =   0

for i in range(e.shape[0]): # iterating over test_noisy dataset

  e[i] = test_noisy_acc[a+1] - test_noisy_acc[a] # evaluating relative acceletation and storing it in e

  a+=2 # increment a by 2 to get the next 2 accelerations.

e = e.reshape(e.shape[0],e.shape[1],1) # reshape the array e.

#EVALUATING MODELS

model3.evaluate(e,pure_classification) # evaluating the performance of the CNN on unseen data