Feature Selection to Optimise Human Activity Recognition Model

debug = False

#defining the root directory of the pressure mat data
pmRoot = "data/"

#timing information - start timer
start = timer()

#for subjects 1 to 30
for i in range(1, 31):
    
    #for exercises 1 to 7
    for j in range(1, 8):
        
        #read in a pressure mat data frame
        
        #formatting the numbers to match file system
        if(i<10): folder = "0" + str(i) + "/"
        else: folder = str(i) + "/"
        file = "0" + str(j) + "_pm_1.csv"
        
        #read in the .csv files
        if (debug == True): print(pmRoot + folder + file)
        pm_df_in = pd.read_csv(pmRoot + folder + file)
        
        #calculate the total pressure for each time increment
        pm_temp_totals = pm_df_in.drop(pm_df_in.columns[0], axis = 1).sum(axis = 1)
        
        #reformatting the time string into an actual datetime object
        for k in range(len(pm_df_in)):
            timeString = pm_df_in.iloc[k,0]
            if len(timeString) < 20:
                timeString = timeString + "."
            while len(timeString) < 26:
                timeString = timeString + "0"
            pm_df_in.iat[k,0] = datetime.strptime(timeString, "%Y-%m-%d %H:%M:%S.%f")
        
        #create a measure of absolute time for each data point
        pm_temp_a_Time = pm_df_in.iloc[:,0]
        
        #create a measure of relative time for each data point
        #have to do this in a for loop for some reason?
        pm_temp_r_Time = [None] * len(pm_temp_a_Time)
        for time in range(len(pm_temp_a_Time)):
            pm_temp_r_Time[time] = pm_temp_a_Time[time] - pm_temp_a_Time[0]
        
        #create a dataframe to store our bad boys in
        pm_temp_df = pd.DataFrame(columns = colNames)
        
        #add in our data
        pm_temp_df["a_Time"] = pm_temp_a_Time
        pm_temp_df["r_Time"] = pm_temp_r_Time
        pm_temp_df["Total_p"] = pm_temp_totals
        
        pm_temp_df[coords] = pm_df_in.iloc[:,1:len(pm_df_in.columns)]
        
        #label the data with the subject and exercise
        pm_temp_df["Subject"] = pd.Series([i for x in range(len(pm_temp_df.index))], index=pm_temp_df.index)
        pm_temp_df["Exercise"] = pd.Series([j for x in range(len(pm_temp_df.index))], index=pm_temp_df.index)
        
        #append to our dataframe
        pm_df = pm_df.append(pm_temp_df)
        
#         break
#     break

#timing information - end timer
end = timer()
dt = end - start
print("Importing time elapsed: {:.2f}".format(dt))

pm_df

#1520 is just a nice one, can alter the row number to anything
rowNum = 1520

#create an array using the pressure point readings in the specified row 
heatmapData = pm_df.iloc[rowNum, 6:len(pm_df.columns)].to_numpy() #start from column index 6 until the end
heatmapData = heatmapData.reshape(32, 16).astype(float) #reshape the datapoints into a 32x16 grid

plt.figure(figsize=(16, 8)) #defining the size and dimensions of the graph
sns.heatmap(heatmapData.T, vmin=0, vmax=1550) #display the heatmap with defined scale
#sns.heatmap(heatmapData.T) #display the heatmap with relative scale

# plotting total pressure against time

whichsubject = 1
whichexercise = 1
numexercise = 7

fig, ax = plt.subplots(numexercise)

for i in range(1,numexercise+1):
    exercise = pm_df.loc[(pm_df['Subject'] == whichsubject) & (pm_df['Exercise'] == i)]
    ax[i-1].plot(exercise.loc[:,'r_Time'], exercise.loc[:,'Total_p']) #plot relative time against total pressure
    #ax[i-1].get_yaxis().set_visible(False) #hide y axis
    ax[i-1].get_xaxis().set_visible(False) #hide x axis
    ax[i-1].title.set_text("Exercise {}".format(i)) #set titles per subplot according to exercise number
    plt.subplots_adjust(top=2, hspace=0.75) #adjust spacing between subplots

# define a function that calculates the approximate entropy
# from https://en.wikipedia.org/wiki/Approximate_entropy

def ApEn(U, m, r) -> float:
    """Approximate_entropy."""

    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [
            len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0)
            for x_i in x
        ]
        return (N - m + 1.0) ** (-1) * sum(np.log(C))

    N = len(U)

    return abs(_phi(m + 1) - _phi(m))

# calculate the ApEn for each feature

column_apen = []

for i in tqdm(range(len(sub1ex2.columns))):
    column_apen.append(ApEn(sub1ex2.iloc[:,i].values, 2, 3))
    
#plot the distribution of ApEn values
sns.displot(column_apen)

#calculating mutual information between each individual sensor result and the total pressure

from sklearn.feature_selection import mutual_info_classif

mutInf = sklearn.feature_selection.mutual_info_regression(sub1ex2.values, sub1ex2_df['Total_p'].values)

#plotting mutual information against number of features
sns.displot(mutInf)

#create a list of centiles that will be used to reduce the volume of features in the df

#100% of dataset reduced by 30% each time:
reduce = 0.3
current = 100
centiles = []

while current > 1:
    centiles.append(current)
    current = current * (1 - reduce)

#print(centiles)

from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.feature_selection import SelectPercentile, mutual_info_classif

X = pm_df_fs #all the pressure sensor features
Y = pm_df['Exercise'] #the target i.e. exercise number

print(X.shape)

start = timer()

X_list = []

for i in tqdm(range(len(centiles))):
    X_list.append(SelectPercentile(mutual_info_classif, centiles[i]).fit_transform(X, Y))
    print("Number of features retained: {}".format(X_list[i].shape[1]))

end = timer()
dt = end - start
print("Feature selecting time elapsed: {:.2f}s".format(dt))

from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report

neigh = KNeighborsClassifier(n_neighbors=3)

fit_time = []
predict_time = []
reports = []

for i in range(len(centiles)): #for each volume of data
    start = timer()
    x_train, x_test, y_train, y_test = train_test_split(X_list[i], Y) #for each dataset in x_list (30% stepwise feature reduction), split the dataset into training and testing sets.
    neigh.fit(x_train, y_train) #run the training data through the KNN
    end = timer()
    fit_time.append(end - start) #record the training time. this will help us compare model performance
    print("Fitting time elapsed for {} columns: {:.2f}s".format(x_train.shape[1], fit_time[i]))
    
    start = timer()
    pred = neigh.predict(x_test) #run test data through our trained KNN to predict exercise performed
    end = timer()
    predict_time.append(end - start) #record the prediction time. this will help us compare model performance
    print("Predicting time elapsed: {:.2f}s".format(predict_time[i]))

    c_mat = confusion_matrix(y_test, pred) #compare the test predictions against the actual target (hidden to trained model). produce a confusion matrix.
    print("Confusion matrix:\n{}".format(c_mat))
    
    reports.append(classification_report(y_test, pred, output_dict = True)) #writes a report of results
    print(reports[i])
    
    print("\n\n")   

#pull macro average precision results from the above reports into a list
macro_avg = []

for i in range(len(reports)):
    macro_avg.append(reports[i]['macro avg']['precision'])

t = list(range(len(centiles)))

numCol = []
for i in range(len(X_list)):
    numCol.append(X_list[i].shape[1])
    
fig, ax1 = plt.subplots()

color = 'tab:red'
ax1.set_xlabel('Number of Columns Used')
ax1.set_ylabel('Macro Average Precision', color=color)
ax1.plot(t, macro_avg, color=color)
ax1.tick_params(axis='y', labelcolor=color)

plt.xticks(t, numCol)

ax2 = ax1.twinx() #instantiate a second axes that shares the same x-axis

color = 'tab:blue'
ax2.set_ylabel('Prediction Time (s)', color=color)  #already handled the x-label with ax1
ax2.plot(t, predict_time, color=color)
ax2.tick_params(axis='y', labelcolor=color)

fig.tight_layout()  #otherwise the right y-label is slightly clipped
plt.show()