refactored tool chain setup, able to debug and run from vscode nowHEADmaster

10 files changed, 1134 insertions, 0 deletions

diff --git a/OASIS/c/src/ekd_blas.c b/OASIS/c/src/ekd_blas.c
new file mode 100644
index 0000000..8555fbf
--- /dev/null
+++ b/OASIS/c/src/ekd_blas.c
@@ -0,0 +1,191 @@
+#include "ekd_blas.h"

+

+void print_vector(float* vector, uint32_t size) {

+    for (uint32_t i=0; i<size; i++) printf("%.2f ", vector[i]);

+}

+

+void print_matrix(float** matrix, uint32_t rows, uint32_t cols) {

+    for (uint32_t i=0; i<rows; i++) {

+        printf("\n");

+        for (uint32_t j=0; j<cols; j++) printf("%.2f ", matrix[i][j]);

+    }

+}

+

+void copy_vector(float* input_vector, float* resultant_vector, uint32_t size) {

+    for (uint32_t i=0; i<size; i++) resultant_vector[i] = input_vector[i];

+}

+

+void copy_matrix(float** input_matrix, float** resultant_matrix, uint32_t rows, uint32_t cols) {

+    for (uint32_t i=0; i<rows; i++){

+        for (uint32_t j=0; j<cols; j++) resultant_matrix[i][j] = input_matrix[i][j];

+    }

+}

+

+float trace_matrix(float** input_matrix, uint32_t size) {

+    float trace = 0.0f;

+    for (uint32_t i=0; i<size; i++) {

+        for (uint32_t j=0; j< size; j++) {

+            if (i == j) trace += input_matrix[i][j];

+        }

+   }

+

+   return trace;

+}

+

+void shift_matrix(float** input_matrix, uint32_t size, float trace, shift_matrix_operation_e operation) {

+    for (uint32_t i=0; i<size; i++) {

+        for (uint32_t j=0; j<size; j++) {

+            if (i == j) {

+                if (operation == SHIFT_ADD) {

+                    input_matrix[i][j] += trace;

+                } else if (operation == SHIFT_SUBTRACT) {

+                    input_matrix[i][j] -= trace;

+                }

+            }

+        }

+    }

+}

+

+void eye(float** input_matrix, uint32_t size) {

+    for (uint32_t i=0; i<size; i++) {

+        for (uint32_t j=0; j<size; j++) input_matrix[i][j] = (i == j) ? 1 : 0;

+    }

+}

+

+void outer_product(float* input_vector, uint32_t size, float** resultant_matrix) {

+    // for (uint32_t i=0; i<size; i++) {

+    //     for (uint32_t j=0; j<size; j++) resultant_matrix[i][j] = input_vector[i] * input_vector[j]; 

+    // }

+    for (uint32_t i = 0; i < size; i++) {

+        for (uint32_t j = i; j < size; j++) {

+            float value = input_vector[i] * input_vector[j];

+            resultant_matrix[i][j] = value;

+            resultant_matrix[j][i] = value;

+        }

+    }

+}

+

+float inner_product(float* input_vector_1, float* input_vector_2, uint32_t size) {

+    float resultant_scalar = 0.0;

+    for (uint32_t i=0; i<size; i++) resultant_scalar += input_vector_1[i] * input_vector_2[i];

+    return resultant_scalar;

+}

+

+void transpose(float** input_matrix, uint32_t rows, uint32_t cols, float** transposed_matrix) {

+    for (uint32_t i=0; i<rows; i++){

+        for (uint32_t j=0; j<cols; j++) transposed_matrix[j][i] = input_matrix[i][j];

+    }

+}

+

+void matrix_operation(float** input_matrix_A, uint32_t rows_A, uint32_t cols_A, float** input_matrix_B, uint32_t rows_B, uint32_t cols_B, float** resultant_matrix, matrix_matrix_operation_e operation) {

+        if (operation == MM_ADD) {

+            for (uint32_t i=0; i<rows_A; i++) {

+                for (uint32_t j=0; j<cols_A; j++) resultant_matrix[i][j] = input_matrix_A[i][j] + input_matrix_B[i][j];

+            }

+        } else if (operation == MM_SUBTRACT) {

+            for (uint32_t i=0; i<rows_A; i++) {

+                for (uint32_t j=0; j<cols_A; j++) resultant_matrix[i][j] = input_matrix_A[i][j] - input_matrix_B[i][j];

+            }

+        } else if (operation == MM_MULTIPLY) {

+            for (uint32_t i=0; i<rows_A; i++) {

+                for (uint32_t j=0; j<cols_B; j++) {

+                    resultant_matrix[i][j] = 0;

+                    for (uint32_t k=0; k<cols_A; k++) resultant_matrix[i][j] += input_matrix_A[i][k] * input_matrix_B[k][j];

+                }

+            }

+        } else if (operation == MM_MULTIPLY_OMP) {

+            uint32_t a, b, c;

+            #pragma omp parallel for private(a, b, c) shared(input_matrix_A, input_matrix_B, resultant_matrix)

+            for (a=0; a<rows_A; a++) {

+                for (b=0; b<cols_B; b++) {

+                resultant_matrix[a][b] = 0;

+                    for (c=0; c<cols_A; c++) resultant_matrix[a][b] += input_matrix_A[a][c] * input_matrix_B[c][b];

+                }

+             }

+        } else if (operation == MM_MULTIPLY_GEMM) {

+            float* temp_matrix_A = malloc(rows_A * cols_A * sizeof(float));

+            float* temp_matrix_B = malloc(rows_B * cols_B * sizeof(float));

+            float* temp_resultant_matrix = malloc(rows_A * cols_B * sizeof(float));

+

+            for (uint32_t i = 0; i < rows_A; i++) {

+                for (uint32_t j = 0; j < cols_A; j++) {

+                    temp_matrix_A[i * cols_A + j] = input_matrix_A[i][j];

+                }

+            }

+

+            for (uint32_t i = 0; i < rows_B; i++) {

+                for (uint32_t j = 0; j < cols_B; j++) {

+                    temp_matrix_B[i * cols_B + j] = input_matrix_B[i][j];

+                }

+            }

+

+            cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, rows_A, cols_B, cols_A, 1.0, temp_matrix_A, cols_A, temp_matrix_B, cols_B, 0.0, temp_resultant_matrix, cols_B);

+

+            for (uint32_t i = 0; i < rows_A; i++) {

+                for (uint32_t j = 0; j < cols_B; j++) {

+                    resultant_matrix[i][j] = temp_resultant_matrix[i * cols_B + j];

+                }

+            }

+

+            free(temp_matrix_A);

+            free(temp_matrix_B);

+            free(temp_resultant_matrix);

+        }

+}

+

+uint8_t sign(float num) {

+    return (num > 0) - (num < 0);

+}

+

+void vector_normalize(float* vector, uint32_t size) {

+    float norm = 0.0;

+    for (uint32_t i=0; i<size; i++) norm += vector[i] * vector[i];

+    norm = sqrt(norm);

+    for (uint32_t i=0; i<size; i++) vector[i] = (norm == 0) ? 0 : vector[i]/norm;

+}

+

+float vector_norm(float* vector, uint32_t size) {

+    float norm = 0.0;

+    for (uint32_t i=0; i<size; i++) norm += vector[i] * vector[i];

+

+    return sqrt(norm);

+}

+

+float matrix_norm(float** input_matrix_A1, float** input_matrix_A, uint32_t rows, uint32_t cols) {

+    float norm = 0.0;

+    for (uint32_t i=0; i<rows; i++) {

+        for (uint32_t j=0; j<cols; j++) norm += fabs(input_matrix_A1[i][j] - input_matrix_A[i][j]) * fabs(input_matrix_A1[i][j] - input_matrix_A[i][j]);

+    }

+

+    return sqrt(norm);

+}

+

+float cosine_similarity(float* curr_vector, float* new_vector, uint32_t size) {

+    float dot_product = 0.0f;

+    for (uint32_t i=0; i<size; i++) dot_product += curr_vector[i] * new_vector[i];

+

+    return (vector_norm(curr_vector, size) == 0 || vector_norm(new_vector, size) == 0) ? 1e-6 : dot_product / (vector_norm(curr_vector, size) * vector_norm(new_vector, size));

+}

+

+void matrix_vector_multiply(float** input_matrix, uint32_t rows, uint32_t cols, float* input_vector, float* resultant_vector, matrix_vector_operation_e operation) {

+    if (operation == MV_MULTIPLY) {

+        for (uint32_t i=0; i<rows; i++) {

+            resultant_vector[i] = 0;

+            for (uint32_t j=0; j<cols; j++) {

+                resultant_vector[i] += input_matrix[i][j] * input_vector[j];

+            }

+        }

+    } else if (operation == MV_MULTIPLY_GEMV) {

+        float* temp_matrix = malloc(rows * cols * sizeof(float));

+

+        for (uint32_t i = 0; i < rows; i++) {

+            for (uint32_t j = 0; j < cols; j++) {

+                temp_matrix[i * cols + j] = input_matrix[i][j];

+            }

+        }

+

+        cblas_sgemv(CblasRowMajor, CblasNoTrans, rows, cols, 1.0f, temp_matrix, cols, input_vector, 1, 0.0f, resultant_vector, 1);

+

+        free(temp_matrix);

+    }

+}
\ No newline at end of file
diff --git a/OASIS/c/src/ekd_blas.h b/OASIS/c/src/ekd_blas.h
new file mode 100644
index 0000000..4c9919d
--- /dev/null
+++ b/OASIS/c/src/ekd_blas.h
@@ -0,0 +1,49 @@
+#ifndef EKD_BLAS_H

+#define EKD_BLAS_H

+

+#include <stdint.h>

+#include <stdlib.h>

+#include <stdio.h>

+#include <string.h>

+#include <math.h>

+#include <omp.h>

+#include <cblas.h>

+

+typedef enum {

+    MM_ADD,

+    MM_SUBTRACT,

+    MM_MULTIPLY,

+    MM_MULTIPLY_OMP,

+    MM_MULTIPLY_GEMM

+} matrix_matrix_operation_e;

+

+typedef enum {

+    MV_MULTIPLY,

+    MV_MULTIPLY_GEMV

+} matrix_vector_operation_e;

+

+typedef enum {

+    SHIFT_ADD,

+    SHIFT_SUBTRACT

+} shift_matrix_operation_e;

+

+// Prototypes

+void print_matrix(float** matrix, uint32_t rows, uint32_t cols);

+void print_vector(float* vector, uint32_t size);

+void copy_matrix(float** input_matrix, float** resultant_matrix, uint32_t rows, uint32_t cols);

+void copy_vector(float* input_vector, float* resultant_vector, uint32_t size);

+void shift_matrix(float** input_matrix, uint32_t size, float trace, shift_matrix_operation_e operation);

+float trace_matrix(float** input_matrix, uint32_t size);

+void eye(float** input_matrix, uint32_t size);

+void outer_product(float* input_vector, uint32_t size, float** resultant_matrix); // outer_productes itself

+float inner_product(float* input_vector_1, float* input_vector_2, uint32_t size); // Dots itself

+void matrix_operation(float** input_matrix_A, uint32_t rows_A, uint32_t cols_A, float** input_matrix_B, uint32_t rows_B, uint32_t cols_B, float** resultant_matrix, matrix_matrix_operation_e operation);

+void transpose(float** input_matrix, uint32_t rows, uint32_t cols, float** transposed_matrix);

+uint8_t sign(float num);

+void vector_normalize(float* vector, uint32_t size);

+float vector_norm(float* vector, uint32_t size);

+float matrix_norm(float** input_matrix_A1, float** input_matrix_A, uint32_t rows, uint32_t cols);

+float cosine_similarity(float* curr_vector, float* new_vector, uint32_t size);

+void matrix_vector_multiply(float** input_matrix, uint32_t rows, uint32_t cols, float* input_vector, float* resultant_vector, matrix_vector_operation_e operation);

+

+#endif /* EKD_BLAS_H */

diff --git a/OASIS/c/src/ekd_cv.c b/OASIS/c/src/ekd_cv.c
new file mode 100644
index 0000000..1f9afad
--- /dev/null
+++ b/OASIS/c/src/ekd_cv.c
@@ -0,0 +1,90 @@
+#include "ekd_cv.h"

+

+void read_bmp(char* filename, uint8_t** data) {

+    bmp_header_t header;

+    bmp_info_t info;

+

+    FILE* file = fopen(filename, "rb");

+    if (file == NULL) {

+        printf("Unable to open file.\n");

+        return;

+    }

+

+    fread(&header.signature, sizeof(uint16_t), 1, file);

+    fread(&header.file_size, sizeof(uint32_t), 1, file);

+    fread(&header.reserved, sizeof(uint32_t), 1, file);

+    fread(&header.data_offset, sizeof(uint32_t), 1, file);

+

+    fread(&info.size, sizeof(uint32_t), 1, file);

+    fread(&info.width, sizeof(uint32_t), 1, file);

+    fread(&info.height, sizeof(uint32_t), 1, file);

+    fread(&info.planes, sizeof(uint16_t), 1, file);

+    fread(&info.bit_depth, sizeof(uint16_t), 1, file);

+    fread(&info.compression, sizeof(uint32_t), 1, file);

+    fread(&info.image_size, sizeof(uint32_t), 1, file);

+    fread(&info.x_pixels_per_meter, sizeof(uint32_t), 1, file);

+    fread(&info.y_pixels_per_meter, sizeof(uint32_t), 1, file);

+    fread(&info.colours_used, sizeof(uint32_t), 1, file);

+    fread(&info.important_colours, sizeof(uint32_t), 1, file);

+

+    uint8_t* raw_data = (uint8_t*) malloc(info.width * info.height * sizeof(uint8_t));

+

+    fseek(file, header.data_offset, SEEK_SET);

+    fread(raw_data, sizeof(uint8_t), info.width * info.height, file);

+

+    fclose(file);

+

+    for (uint8_t i=0; i<info.height; i++) {

+        for (uint8_t j=0; j < info.width; j++) data[i][j] = raw_data[(i * info.width) + j];

+    }

+

+    free(raw_data);

+}

+

+void read_dir(char* dneg_path, char* dpos_path, char** file_name_arr) {

+    DIR* dneg_dir = opendir(dneg_path);

+    DIR* dpos_dir = opendir(dpos_path);

+    struct dirent *entry;

+    uint8_t item_id = 0; 

+

+    while ((entry = readdir(dneg_dir)) != NULL) {

+        if (strcmp(entry->d_name, ".") == 0 || strcmp(entry->d_name, "..") == 0) {

+            continue;

+        }

+        strcat(file_name_arr[item_id], entry->d_name);

+        item_id++;

+    }

+

+    while ((entry = readdir(dpos_dir)) != NULL) {

+        if (strcmp(entry->d_name, ".") == 0 || strcmp(entry->d_name, "..") == 0) {

+            continue;

+        }

+        strcat(file_name_arr[item_id], entry->d_name);

+        item_id++;

+    }

+ 

+    closedir(dneg_dir);

+    closedir(dpos_dir);

+}

+

+void flatten_image(uint8_t** image, float* image_flattened) {

+    for (uint8_t i=0; i<IMAGE_HEIGHT; i++) {

+        for (uint8_t j=0; j<IMAGE_WIDTH; j++) image_flattened[(i*IMAGE_WIDTH)+j] = (float)image[i][j];

+    }

+}

+

+void shuffle_data_labels(float** data, int32_t* labels, uint8_t num_images) {

+    srand(time(NULL));  // initialize random seed

+    for (int i = num_images - 1; i > 0; i--) {

+        // Generate a random index between 0 and i (inclusive)

+        int j = rand() % (i + 1);

+

+        // Swap the rows of data_flatten and labels at indices i and j

+        float* temp_data = data[i];

+        data[i] = data[j];

+        data[j] = temp_data;

+        int32_t temp_label = labels[i];

+        labels[i] = labels[j];

+        labels[j] = temp_label;

+    }

+}
\ No newline at end of file
diff --git a/OASIS/c/src/ekd_cv.h b/OASIS/c/src/ekd_cv.h
new file mode 100644
index 0000000..146144c
--- /dev/null
+++ b/OASIS/c/src/ekd_cv.h
@@ -0,0 +1,45 @@
+#ifndef EKD_CV_H

+#define EKD_CV_H

+

+#include <stdlib.h>

+#include <stdio.h>

+#include <stdint.h>

+#include <dirent.h>

+#include <string.h>

+#include <time.h>

+

+#define NUM_DNEG_IMAGES 134

+#define NUM_DPOS_IMAGES 99

+#define NUM_IMAGES 233

+#define IMAGE_WIDTH 144

+#define IMAGE_HEIGHT 176

+#define FEATURE_VECTOR_SIZE 25344

+

+typedef struct {

+    uint16_t signature;

+    uint32_t file_size; 

+    uint32_t reserved; 

+    uint32_t data_offset; 

+} bmp_header_t;

+

+typedef struct {

+    uint32_t size;

+    uint32_t width;

+    uint32_t height;

+    uint16_t planes;

+    uint16_t bit_depth;

+    uint32_t compression;

+    uint32_t image_size;

+    uint32_t x_pixels_per_meter;

+    uint32_t y_pixels_per_meter;

+    uint32_t colours_used;

+    uint32_t important_colours;

+} bmp_info_t;

+

+// Prototypes

+void read_bmp(char* filename, uint8_t** data_flatten);

+void read_dir(char* dneg_path, char* dpos_path, char** file_name_arr);

+void flatten_image(uint8_t** image, float* image_flattened);

+void shuffle_data_labels(float** data, int32_t* labels, uint8_t num_images);

+

+#endif /* EKD_CV_H */
\ No newline at end of file
diff --git a/OASIS/c/src/ekd_pca.c b/OASIS/c/src/ekd_pca.c
new file mode 100644
index 0000000..f76d20e
--- /dev/null
+++ b/OASIS/c/src/ekd_pca.c
@@ -0,0 +1,286 @@
+#include "ekd_pca.h"

+

+void centre_scale(float** input_matrix, uint32_t rows, uint32_t cols, centre_scale_operation_e operation) {

+    float* means = (float*) malloc(cols * sizeof(float));

+    float* std_devs = (float*) malloc(cols * sizeof(float));

+

+    // Iterating column by column

+    for (uint32_t j = 0; j<cols; j++) {

+        float col_sum = 0.0;

+        for (uint32_t i = 0; i<rows; i++){

+            col_sum += input_matrix[i][j];

+        }

+        means[j] = col_sum / (float)rows;

+

+        float variance = 0.0;

+        for (uint32_t i = 0; i<rows; i++) variance += (input_matrix[i][j] - means[j]) * (input_matrix[i][j] - means[j]);

+        std_devs[j] = (variance < 1e-6) ? 1e-6 : sqrtf(variance / (float)rows);

+    }

+

+    // Centering -> subtract means, Scaling -> Divide by std_devs

+    for (uint32_t j=0; j<cols; j++) {

+        for (uint32_t i=0; i<rows; i++) {

+            if (operation == MEAN_UNIT_VARIANCE) {

+                input_matrix[i][j] = (input_matrix[i][j] - means[j]) / std_devs[j];

+            } else if (operation == MEAN) {

+                input_matrix[i][j] = (input_matrix[i][j] - means[j]);

+            }

+        }

+    }

+

+    free(means);

+    free(std_devs);

+}

+

+void covariance(float** input_matrix, uint32_t rows, uint32_t cols, float** covariance_matrix) {

+    float** transposed_matrix = (float**) malloc(cols * sizeof(float*));

+    for (uint32_t i=0; i<cols; i++) transposed_matrix[i] = (float*) malloc(rows * sizeof(float));

+

+    transpose(input_matrix, rows, cols, transposed_matrix);

+

+    matrix_operation(transposed_matrix, cols, rows, input_matrix, rows, cols, covariance_matrix, MM_MULTIPLY_GEMM);

+

+    // Divide by n since we have full population

+    for (uint32_t i=0; i<cols; i++) {

+        for (uint32_t j=0; j<cols; j++) covariance_matrix[i][j] /= rows-1;

+    }

+    

+    for (uint32_t i=0; i<cols; i++) free(transposed_matrix[i]);

+    free(transposed_matrix);

+}

+

+void house_holder_vector(float* vector, uint32_t size) {

+

+    // Computes the norm of the input vector

+    float norm_input_vector = 0.0;

+    for (uint32_t i=0; i<size; i++) norm_input_vector += vector[i] * vector[i];

+    norm_input_vector = sqrt(norm_input_vector);

+

+    // Since the rest of output vector already is the input vector, change the first index

+    vector[0] = vector[0] + sign(vector[0]) * norm_input_vector;

+}

+

+void house_holder_transformation(float* house_holder_vector, uint32_t size, float** house_holder_matrix) {

+    float** identity_matrix = (float**) malloc(size * sizeof(float*));

+    for (uint32_t i=0; i<size; i++) identity_matrix[i] = (float*) malloc(size * sizeof(float));

+

+    float** numerator = (float**) malloc(size * sizeof(float*));

+    for (uint32_t i=0; i<size; i++) numerator[i] = (float*) malloc(size * sizeof(float));

+

+    eye(identity_matrix, size);

+

+    float denominator = inner_product(house_holder_vector, house_holder_vector, size);

+

+    outer_product(house_holder_vector, size, numerator);

+

+    for (uint32_t i=0; i<size; i++) {

+        for (uint32_t j=0; j<size; j++) numerator[i][j] = (denominator == 0) ? 0 : (2 / denominator) * numerator[i][j];

+    }

+

+    matrix_operation(identity_matrix, size, size, numerator, size, size, house_holder_matrix, MM_SUBTRACT);

+

+    for (uint32_t i=0; i<size; i++){

+        free(numerator[i]);

+        free(identity_matrix[i]);

+    } 

+    free(numerator);

+    free(identity_matrix);

+}

+

+void qr_decomposition(float** input_matrix, uint32_t size, float** resultant_matrix) {

+    float** matrix_Q = (float**) malloc(size * sizeof(float*)); // Has size x size dimensions

+    float** matrix_Q_prime = (float**) malloc(size * sizeof(float*)); // Has size x size dimensions

+    float** matrix_R = (float**) malloc(size * sizeof(float*)); // Has size x size dimensions

+    float** matrix_R_prime = (float**) malloc(size * sizeof(float*)); // Has size x size dimensions

+    for(uint32_t i=0; i<size; i++) {

+        matrix_Q[i] = (float*) malloc(size * sizeof(float));

+        matrix_Q_prime[i] = (float*) malloc(size * sizeof(float));

+        matrix_R[i] = (float*) malloc(size * sizeof(float));

+        matrix_R_prime[i] = (float*) malloc(size * sizeof(float));

+    }

+

+    for (uint32_t i=0; i<size; i++) memcpy(matrix_R[i], input_matrix[i], size * sizeof(float));

+

+    eye(matrix_Q, size);

+

+    for (uint32_t i=0; i<size; i++) {

+        float* current_house_holder_vector = (float*) malloc((size - i) * sizeof(float));

+

+        float** current_house_holder_matrix = (float**) malloc((size - i) * sizeof(float*)); // Has (size-i)x(size-i) dimensions

+        for (uint32_t j=0; j<(size - i); j++) current_house_holder_matrix[j] = (float*) malloc((size - i) * sizeof(float));

+

+        float** matrix_H = (float**) malloc(size * sizeof(float*)); // Has size x size dimensions

+        for (uint32_t j=0; j<size; j++) matrix_H[j] = (float*) malloc(size * sizeof(float));

+

+        eye(matrix_H, size);

+

+        for(uint32_t j=0; j<(size - i); j++) current_house_holder_vector[j] = matrix_R[j+i][i];

+

+        house_holder_vector(current_house_holder_vector, (size - i));

+

+        house_holder_transformation(current_house_holder_vector, (size - i), current_house_holder_matrix);

+

+        // Starting from the ith ith sub matrix

+        for (uint32_t j=i; j<size; j++) {

+            for (uint32_t k=i; k<size; k++) matrix_H[j][k] = current_house_holder_matrix[j-i][k-i];

+        }

+

+        matrix_operation(matrix_H, size, size, matrix_R, size, size, matrix_R_prime, MM_MULTIPLY_GEMM);

+        matrix_operation(matrix_Q, size, size, matrix_H, size, size, matrix_Q_prime, MM_MULTIPLY_GEMM);

+

+        // print_matrix(matrix_H, size, size);

+        

+        for (uint32_t j=0; j<size; j++) memcpy(matrix_Q[j], matrix_Q_prime[j], size * sizeof(float));

+        for (uint32_t j=0; j<size; j++) memcpy(matrix_R[j], matrix_R_prime[j], size * sizeof(float));

+

+        for (uint32_t j=0; j<size; j++) free(matrix_H[j]);

+        free(matrix_H);

+        for (uint32_t j=0; j<(size - i); j++) free(current_house_holder_matrix[j]);

+        free(current_house_holder_matrix);

+        free(current_house_holder_vector);

+    }

+

+    matrix_operation(matrix_R, size, size, matrix_Q, size, size, resultant_matrix, MM_MULTIPLY_GEMM);

+    matrix_operation(matrix_Q, size, size, matrix_R, size, size, resultant_matrix, MM_MULTIPLY_GEMM);

+

+    print_matrix(matrix_R, size, size);

+

+    for (uint32_t i=0; i<size; i++) {

+        free(matrix_R[i]);

+        free(matrix_R_prime[i]);

+        free(matrix_Q_prime[i]);

+        free(matrix_Q[i]);

+    }

+    free(matrix_R_prime);

+    free(matrix_R);

+    free(matrix_Q_prime);

+    free(matrix_Q);

+}

+

+void qr_algorithm(float** input_matrix, uint32_t rows, uint32_t cols, float** resultant_matrix) {

+    float** covariance_matrix = (float**) malloc(cols * sizeof(float*)); // Has feature x feature dimensions

+    for (uint32_t i=0; i<cols; i++) covariance_matrix[i] = (float*) malloc(cols * sizeof(float));

+

+    float** covariance_matrix_prime = (float**) malloc(cols * sizeof(float*)); // Has feature x feature dimensions

+    for (uint32_t i=0; i<cols; i++) covariance_matrix_prime[i] = (float*) malloc(cols * sizeof(float));

+

+    centre_scale(input_matrix, rows, cols, MEAN);

+    covariance(input_matrix, rows, cols, covariance_matrix);

+

+    qr_decomposition(covariance_matrix, cols, covariance_matrix_prime);

+    while (matrix_norm(covariance_matrix_prime, covariance_matrix, cols, cols) > TOLERANCE) {

+        for (uint32_t i=0; i<cols; i++) {

+            for (uint32_t j=0; j<cols; j++) {

+                if (i == j) printf("%.2f ", covariance_matrix[i][j]);

+            } 

+        }

+

+        memcpy(covariance_matrix, covariance_matrix_prime, cols);

+        qr_decomposition(covariance_matrix, cols, covariance_matrix_prime);

+

+    for (uint32_t i=0; i<cols; i++) free(covariance_matrix[i]);

+    free(covariance_matrix);

+    }

+}

+

+void power_iteration(float** input_matrix, uint32_t size, float* eigen_vector) {

+    // Initialize first guess vector

+    float* new_vector = (float*) malloc(size * sizeof(float));

+    float* curr_vector = (float*) malloc(size * sizeof(float));

+    for (uint32_t i=0; i<size; i++) curr_vector[i] = (float) rand() / RAND_MAX;

+

+    // MAX_INTERATIONS in case the singular vector doesn't converge

+    for (uint32_t i=0; i<MAX_ITERATIONS; i++) {

+        matrix_vector_multiply(input_matrix, size, size, curr_vector, new_vector, MV_MULTIPLY);

+        vector_normalize(curr_vector, size);

+        vector_normalize(new_vector, size);

+        if (fabs(fabs(cosine_similarity(new_vector, curr_vector, size))-1.0) < TOLERANCE) break;

+        memcpy(curr_vector, new_vector, size * sizeof(float));

+    }

+

+    memcpy(eigen_vector, new_vector, size * sizeof(float));

+    free(curr_vector);

+    free(new_vector);

+}

+

+float eigen_value_compute(float* eigen_vector, float** input_matrix, uint32_t size) {

+    if (vector_norm(eigen_vector, size) == 0) return 0;

+

+    float lambda = 0.0f;

+    for (uint32_t i=0; i<size; i++) {

+        lambda += input_matrix[i][i] * eigen_vector[i] * eigen_vector[i];

+        for (uint32_t j=0; j<size; j++) {

+            if (i != j) {

+                lambda += input_matrix[i][j] * eigen_vector[i] * eigen_vector[j];

+            }

+        }

+    }

+    lambda /= vector_norm(eigen_vector, size) * vector_norm(eigen_vector, size);

+    return lambda;

+}

+

+void subtract_eigen(float* eigen_vector, float** input_matrix, float** resultant_matrix, float eigen_value, uint32_t size) {

+    float** outer_product_matrix = (float**) malloc(size * sizeof(float*));

+    for (uint32_t i=0; i<size; i++) outer_product_matrix[i] = (float*) malloc(size * sizeof(float));

+

+    outer_product(eigen_vector, size, outer_product_matrix);

+

+    for (uint32_t i=0; i<size; i++) {

+        for (uint32_t j=0; j<size; j++) outer_product_matrix[i][j] *= eigen_value;

+    }

+

+    matrix_operation(input_matrix, size, size, outer_product_matrix, size, size, resultant_matrix, MM_SUBTRACT);

+

+    for (uint32_t i=0 ; i<size; i++) free(outer_product_matrix[i]);

+    free(outer_product_matrix);

+}

+

+void pca_covariance_method(float** input_matrix, uint32_t rows, uint32_t cols, float** resultant_matrix, uint32_t components) {

+    float** input_matrix_original = (float**) malloc(rows * sizeof(float*)); // Has sample x feature dimensions

+    for (uint32_t i=0; i<rows; i++) input_matrix_original[i] = (float*) malloc(cols * sizeof(float));

+

+    float* eigen_vector = (float*) malloc(cols * sizeof(float)); // Has feature dimensions

+

+    float** covariance_matrix = (float**) malloc(cols * sizeof(float*)); // Has feature x feature dimensions

+    float** covariance_matrix_prime = (float**) malloc(cols * sizeof(float*));

+    float** projection_matrix = (float**) malloc(cols * sizeof(float*)); // Has features x components dimensions

+

+    for (uint32_t i=0; i<cols; i++) {

+        projection_matrix[i] = (float*) malloc(components * sizeof(float));

+        covariance_matrix[i] = (float*) malloc(cols * sizeof(float));

+        covariance_matrix_prime[i] = (float*) malloc(cols * sizeof(float));

+    }

+

+    copy_matrix(input_matrix, input_matrix_original, rows, cols);

+

+    centre_scale(input_matrix, rows, cols, MEAN);

+    covariance(input_matrix, rows, cols, covariance_matrix);

+

+    for (uint32_t i=0; i<components; i++) {

+        power_iteration(covariance_matrix, cols, eigen_vector);

+        float eigen_value = eigen_value_compute(eigen_vector, covariance_matrix, cols);

+        subtract_eigen(eigen_vector, covariance_matrix, covariance_matrix_prime, eigen_value, cols);

+        for (uint32_t j=0; j<cols; j++) projection_matrix[j][i] = eigen_vector[j];

+       

+        copy_matrix(covariance_matrix_prime, covariance_matrix, cols, cols);

+    }

+

+    matrix_operation(input_matrix_original, rows, cols, projection_matrix, cols, components, resultant_matrix, MM_MULTIPLY_GEMM); // Has sample x components dimensions

+

+    free(eigen_vector);

+

+    for (uint32_t i=0; i<cols; i++) {

+        free(projection_matrix[i]);

+        free(covariance_matrix_prime[i]);

+        free(covariance_matrix[i]);

+    }

+

+    for (uint32_t i=0; i<rows; i++) {

+        free(input_matrix_original[i]);

+    }

+

+    free(covariance_matrix_prime);

+    free(covariance_matrix); 

+    free(input_matrix_original);

+    free(projection_matrix);

+}

diff --git a/OASIS/c/src/ekd_pca.h b/OASIS/c/src/ekd_pca.h
new file mode 100644
index 0000000..842e485
--- /dev/null
+++ b/OASIS/c/src/ekd_pca.h
@@ -0,0 +1,38 @@
+#ifndef EKD_PCA_H

+#define EKD_PCA_H

+

+#include <stdint.h>

+#include <stdlib.h>

+#include <stdio.h>

+#include <string.h>

+#include <math.h>

+

+#include "ekd_blas.h"

+

+#define NUM_DNEG_IMAGES 134

+#define NUM_DPOS_IMAGES 99

+#define NUM_IMAGES 233

+#define IMAGE_WIDTH 144

+#define IMAGE_HEIGHT 176

+#define FEATURE_VECTOR_SIZE 25344

+#define MAX_ITERATIONS 100

+#define TOLERANCE 1e-3

+

+typedef enum {

+    MEAN_UNIT_VARIANCE,

+    MEAN

+} centre_scale_operation_e;

+

+// Prototypes

+void centre_scale(float** input_matrix, uint32_t rows, uint32_t cols, centre_scale_operation_e operation);

+void covariance(float** input_matrix, uint32_t rows, uint32_t cols, float** covariance_matrix);

+void house_holder_vector(float* input_vector, uint32_t size);

+void house_holder_transformation(float* house_holder_vector, uint32_t size, float** house_holder_matrix);

+void qr_decomposition(float** input_matrix, uint32_t size, float** resultant_matrix);

+void qr_algorithm(float** input_matrix, uint32_t rows, uint32_t cols, float** resultant_matrix);

+void power_iteration(float** input_matrix, uint32_t size, float* eigen_vector);

+float eigen_value_compute(float* eigen_vector, float** input_matrix, uint32_t size);

+void subtract_eigen(float* eigen_vector, float** input_matrix, float** resultant_matrix, float eigen_value, uint32_t size);

+void pca_covariance_method(float** input_matrix, uint32_t rows, uint32_t cols, float** resultant_matrix, uint32_t components);

+

+#endif /* EKD_PCA_H */
\ No newline at end of file
diff --git a/OASIS/c/src/ekd_svm.c b/OASIS/c/src/ekd_svm.c
new file mode 100644
index 0000000..13c155d
--- /dev/null
+++ b/OASIS/c/src/ekd_svm.c
@@ -0,0 +1,259 @@
+#include "ekd_svm.h"
+
+void svm_preprocess(float** input_matrix, uint32_t rows, uint32_t cols) {
+    centre_scale(input_matrix, rows, cols, MEAN_UNIT_VARIANCE);
+}
+
+float rbf_kernel(float* sample_1, float* sample_2, uint32_t features, float gamma) {
+    float distance = 0.0f;
+    for (uint32_t i=0; i<features; i++) distance += pow(sample_1[i] - sample_2[i], 2);
+    return exp(-gamma * distance);
+}
+
+void kernel_matrix_transform(float** input_matrix, float** resultant_matrix, uint32_t samples, uint32_t features, float gamma) {
+    for (uint32_t i=0; i<samples; i++) {
+        for (uint32_t j=0; j<samples; j++) {
+            float kernel = rbf_kernel(input_matrix[i], input_matrix[j], features, gamma);
+            resultant_matrix[i][j] = kernel;
+        }
+    }
+}
+
+void kernel_row_transform(float** input_matrix, float* query_row, float* resultant_row, uint32_t samples, uint32_t features, float gamma) {
+    for (uint32_t i=0; i<samples; i++) {
+        float kernel = rbf_kernel(input_matrix[i], query_row, features, gamma);
+        resultant_row[i] = kernel;
+    }
+}
+
+float get_error(uint32_t index, float* alphas, float** kernel_matrix, int32_t* labels, uint32_t samples) {
+    float error = 0.0f;
+
+    for (uint32_t i=0; i<samples; i++) error += alphas[i] * labels[i] * kernel_matrix[i][index];
+
+    error -= labels[index];
+    return error;
+}
+
+uint32_t find_second_alpha(uint32_t first_index, float first_error_val, float* alphas, float** kernel_matrix, int32_t* labels, uint32_t samples) {
+    int32_t second_index = -1;
+    float max_error_diff = 0;
+
+    for (uint32_t i=0; i<samples; i++) {
+        if (i == first_index) continue;
+
+        float second_error_val = get_error(i, alphas, kernel_matrix, labels, samples);
+        float error_diff = fabs(first_error_val - second_error_val);
+        if (error_diff > max_error_diff) {
+            max_error_diff = error_diff;
+            second_index = i;
+        }
+    }
+    return second_index;
+}
+
+void calculate_bounds(float first_label, float second_label, float first_alpha, float second_alpha, float penalty, float* lower_bound, float* upper_bound) {
+    if (first_label != second_label) {
+        *lower_bound = fmax(0, second_alpha - first_alpha);
+        *upper_bound = fmin(penalty, penalty + second_alpha - first_alpha);
+    } else {
+        *lower_bound = fmax(0, second_alpha + first_alpha - penalty);
+        *upper_bound = fmin(penalty, second_alpha + first_alpha);
+    }
+}
+
+float clip_second_alpha_new(float second_alpha_new, float lower_bound, float upper_bound) {
+    if (second_alpha_new < lower_bound) {
+        return lower_bound;
+    } else if (second_alpha_new > upper_bound) {
+        return upper_bound;
+    } else {
+        return second_alpha_new;
+    }
+}
+
+void update_alphas(float first_alpha_new, float second_alpha_new, float* alphas, uint32_t first_index, uint32_t second_index) {
+    alphas[first_index] = first_alpha_new;
+    alphas[second_index] = second_alpha_new;
+}
+
+float calculate_bias(float first_error_val, float second_error_val, float first_alpha_new, float first_alpha, float second_alpha_new, float second_alpha, float kernel_val_ii, float kernel_val_ij, float kernel_val_jj, float first_label, float second_label, float penalty) {
+    float first_bias = first_error_val + first_label * (first_alpha_new - first_alpha) * kernel_val_ii + second_label * (second_alpha_new - second_alpha) * kernel_val_ij;
+    float second_bias = first_error_val + first_label * (first_alpha_new - first_alpha) * kernel_val_ij + second_label * (second_alpha_new - second_alpha) * kernel_val_jj;
+    float bias;
+    if (first_alpha_new > 0 && first_alpha_new < penalty) {
+        bias = first_bias;
+    } else if (second_alpha_new > 0 && second_alpha_new < penalty) {
+        bias = second_bias;
+    } else {
+        bias = (first_bias + second_bias) / 2;
+    }
+    return bias;
+}
+
+void update_errors(float* errors, uint32_t samples, uint32_t first_index, uint32_t second_index, float first_label, float second_label, float first_alpha_new, float first_alpha, float second_alpha_new, float second_alpha, float* bias, float** kernel_matrix) {
+    float bias_temp = *bias;
+    for (uint32_t i=0; i<samples; i++) {
+        if (i == first_index || i == second_index) continue;
+        float error = errors[i];
+        error += first_label * (first_alpha_new - first_alpha) * kernel_matrix[first_index][i];
+        error -= bias_temp - *bias;
+        errors[i] = error;
+    }
+}
+
+uint8_t optimize_alphas(uint32_t first_index, uint32_t second_index, float* alphas, float* bias, float* errors, float** kernel_matrix, int32_t* labels, uint32_t samples, float penalty) {
+    if (first_index == second_index) return 0;
+
+    float first_label = (float) labels[first_index];
+    float second_label = (float) labels[second_index];
+
+    float first_alpha = alphas[first_index];
+    float second_alpha = alphas[second_index];
+
+    float first_error_val = errors[first_index];
+    float second_error_val = errors[second_index];
+
+    float kernel_val_ii = kernel_matrix[first_index][first_index];
+    float kernel_val_jj = kernel_matrix[second_index][second_index];
+    float kernel_val_ij = kernel_matrix[first_index][second_index];
+
+    // Second derivative of objective function
+    float learning_rate = 2 * kernel_val_ij - kernel_val_ii - kernel_val_jj;
+    if (learning_rate >= 0) return 0;
+
+    float second_alpha_new = second_alpha - (second_label * (first_error_val - second_error_val) / learning_rate);
+
+    float lower_bound, upper_bound;
+    calculate_bounds(first_label, second_label, first_alpha, second_alpha, penalty, &lower_bound, &upper_bound);
+    second_alpha_new = clip_second_alpha_new(second_alpha_new, lower_bound, upper_bound);
+
+    if (fabs(second_alpha_new - second_alpha) < SVM_TOLERANCE * (second_alpha_new + second_alpha + SVM_TOLERANCE)) return 0;
+    
+    float first_alpha_new = first_alpha + first_label * second_label * (second_alpha - second_alpha_new);
+
+    update_alphas(first_alpha_new, second_alpha_new, alphas, first_index, second_index);
+
+    float bias_new = calculate_bias(first_error_val, second_error_val, first_alpha_new, first_alpha, second_alpha_new, second_alpha, kernel_val_ii, kernel_val_ij, kernel_val_jj, first_label, second_label, penalty);
+    *bias = bias_new;
+
+    update_errors(errors, samples, first_index, second_index, first_label, second_label, first_alpha_new, first_alpha, second_alpha_new, second_alpha, bias, kernel_matrix);
+
+    return 1;
+}
+
+uint8_t verify_kkt(uint32_t first_index, float tolerance, float* alphas, float* bias, float* errors, float** kernel_matrix, int32_t* labels, uint32_t samples, float penalty) {
+    float first_label = (float) labels[first_index];
+    float first_alpha =  alphas[first_index];
+    float first_error_val = errors[first_index];
+    float first_residual_error = first_error_val * first_label;
+    
+    if ((first_residual_error < -tolerance && first_alpha < penalty) || (first_residual_error > tolerance && first_alpha > 0)) {
+        // Alpha is able be optimized, now find second alpha using heuristic method
+        uint32_t second_index = find_second_alpha(first_index, first_error_val, alphas, kernel_matrix, labels, samples);
+        // If the alphas have been optimized, return
+        if (optimize_alphas(first_index, second_index, alphas, bias, errors, kernel_matrix, labels, samples, penalty)) return 1;
+
+        // Could not find second alpha or couldn't optimize alphas values
+        uint32_t max_index = -1;
+        float max_error = -FLT_MAX;
+        for (uint32_t i=0; i<samples; i++) {
+            if (alphas[i] > 0 && alphas[i] < penalty) {
+                float error = errors[i];
+                if (fabs(error - first_error_val) > max_error) {
+                    max_index = i;
+                    max_error = fabs(error - first_error_val);
+                }
+            }
+        }
+
+        if (max_index != -1 && optimize_alphas(first_index, max_index, alphas, bias, errors, kernel_matrix, labels, samples, penalty) == 1) return 1;
+        
+        // Still couldn't optimize alpha values:w
+        return 0;
+    }
+
+    // Alpha values do not need to be optimized
+    return 0;
+}
+
+void train_svm(float** data_matrix, uint32_t samples, uint32_t features, float* alphas, float* bias, int32_t* labels, float penalty) {
+    // Centre and scale reduced_dimension_matrix
+    svm_preprocess(data_matrix, samples, features);
+
+    float** kernel_matrix = (float**) malloc(samples * sizeof(float*));
+    for (uint32_t i=0; i<samples; i++) kernel_matrix[i] = (float*) malloc(samples * sizeof(float));
+
+    // Compute kernel matrix using RBF
+    kernel_matrix_transform(data_matrix, kernel_matrix, samples, features, GAMMA);
+
+    // Initialize bias to 0, alphas/errors to random numbers between -1 and 1
+    *bias = 0;
+
+    float* errors = (float*) malloc(samples * sizeof(float));
+    for (uint32_t i=0; i<samples; i++) {
+        errors[i] = (float) rand() / (float) RAND_MAX * 2 - 1;
+        alphas[i] = (float) rand() / (float) RAND_MAX * 2 - 1;
+    }
+
+    uint32_t num_changed_alphas = 0;
+    uint32_t search_passes = 0;
+    uint8_t searched_all_flag = 1;
+
+    while ((search_passes < MAX_SVM_PASSES && num_changed_alphas > 0) || searched_all_flag) {
+        num_changed_alphas = 0;
+
+        if (searched_all_flag) {
+            // Loops over all samples
+            for (uint32_t sample_index=0; sample_index<samples; sample_index++) 
+                num_changed_alphas += verify_kkt(sample_index, SVM_TOLERANCE, alphas, bias, errors, kernel_matrix, labels, samples, penalty);
+        } else {
+            // Loops over non support vector samples (not 0 and not C -> penalty parameter)
+            for (uint32_t sample_index=0; sample_index<samples; sample_index++) {
+                if (alphas[sample_index] > 0 && alphas[sample_index] < penalty)
+                    num_changed_alphas += verify_kkt(sample_index, SVM_TOLERANCE, alphas, bias, errors, kernel_matrix, labels, samples, penalty);
+            }
+        }
+
+        if (searched_all_flag) {
+            searched_all_flag = 0;
+        } else if (num_changed_alphas == 0) {
+            searched_all_flag = 1;
+        }
+        search_passes++;
+    }
+
+    for (uint32_t i=0; i<samples; i++) free (kernel_matrix[i]);
+    free(kernel_matrix);
+    free(errors);
+}
+
+int32_t svm_predict(float** test_data, int32_t* labels_test, uint32_t num_test, uint32_t num_components, float* alphas, float bias, float gamma) {
+    int num_correct = 0;
+    for (int i=0; i<num_test; i++) {
+        float transformed_row[num_test];
+        kernel_row_transform(test_data, test_data[i], transformed_row, num_test, num_components, gamma);
+        float prediction = bias;
+        for (int j = 0; j<num_test; j++) {
+            prediction += alphas[j] * transformed_row[j];
+        }
+        if (prediction >= 0 && labels_test[i] == 1) {
+            num_correct++;
+        } else if (prediction < 0 && labels_test[i] == -1) {
+            num_correct++;
+        }
+    }
+    return num_correct;
+}
+
+float test_svm(float** data_test, int32_t* labels_test, uint32_t num_test, uint32_t num_components, float* alphas, float bias) {
+    uint32_t num_correct = 0;
+    float* prediction = (float*) malloc(num_test * sizeof(float));
+    for (uint8_t i=0; i<num_test; i++) {
+        prediction[i] = svm_predict(data_test, labels_test, num_test, num_components, alphas, bias, GAMMA);
+        if (prediction[i] * labels_test[i] > 0) num_correct++;
+    }
+    float accuracy = ((float) num_correct / num_test) * 100;
+    free(prediction);
+    return accuracy;
+}
\ No newline at end of file
diff --git a/OASIS/c/src/ekd_svm.h b/OASIS/c/src/ekd_svm.h
new file mode 100644
index 0000000..259eb5c
--- /dev/null
+++ b/OASIS/c/src/ekd_svm.h
@@ -0,0 +1,40 @@
+#ifndef EKD_SVM_H
+#define EKD_SVM_H
+
+#include <stdint.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <float.h>
+
+#include "ekd_pca.h"
+#include "ekd_blas.h"
+
+#define NUM_DNEG_IMAGES 134
+#define NUM_DPOS_IMAGES 99
+#define NUM_IMAGES 233
+#define IMAGE_WIDTH 144
+#define IMAGE_HEIGHT 176
+#define FEATURE_VECTOR_SIZE 25344
+#define MAX_SVM_PASSES 100
+#define GAMMA 0.9
+#define SVM_TOLERANCE 1e-3
+
+void svm_preprocess(float** input_matrix, uint32_t rows, uint32_t cols);
+float rbf_kernel(float* sample_1, float* sample_2, uint32_t features, float gamma);
+void kernel_matrix_transform(float** input_matrix, float** resultant_matrix, uint32_t samples, uint32_t features, float gamma);
+void kernel_row_transform(float** input_matrix, float* query_row, float* resultant_row, uint32_t samples, uint32_t features, float gamma);
+float get_error(uint32_t index, float* alphas, float** kernel_matrix, int32_t* labels, uint32_t num_samples);
+uint32_t find_second_alpha(uint32_t first_index, float first_error_val, float* alphas, float** kernel_matrix, int32_t* labels, uint32_t num_samples);
+void calculate_bounds(float first_label, float second_label, float first_alpha, float second_alpha, float penalty, float* lower_bound, float* upper_bound);
+float clip_second_alpha_new(float second_alpha_new, float lower_bound, float upper_bound);
+void update_alphas(float first_alpha_new, float second_alpha_new, float* alphas, uint32_t first_index, uint32_t second_index);
+float calculate_bias(float first_error_val, float second_error_val, float first_alpha_new, float first_alpha, float second_alpha_new, float second_alpha, float kernel_val_ii, float kernel_val_ij, float kernel_val_jj, float first_label, float second_label, float penalty);
+void update_errors(float* errors, uint32_t num_samples, uint32_t first_index, uint32_t second_index, float first_label, float second_label, float first_alpha_new, float first_alpha, float second_alpha_new, float second_alpha, float* bias, float** kernel_matrix);
+uint8_t optimize_alphas(uint32_t first_index, uint32_t second_index, float* alphas, float* bias, float* errors, float** kernel_matrix, int32_t* labels, uint32_t num_samples, float penalty);
+uint8_t verify_kkt(uint32_t first_index, float tolerance, float* alphas, float* bias, float* errors, float** kernel_matrix, int32_t* labels, uint32_t num_samples, float penalty);
+void train_svm(float** data_matrix, uint32_t samples, uint32_t features, float* alphas, float* bias, int32_t* labels, float penalty);
+float test_svm(float** data_test, int32_t* labels_test, uint32_t num_test, uint32_t num_components, float* alphas, float bias);
+
+#endif /* EKD_SVM_H */
\ No newline at end of file
diff --git a/OASIS/c/src/main.c b/OASIS/c/src/main.c
new file mode 100644
index 0000000..ef8918d
--- /dev/null
+++ b/OASIS/c/src/main.c
@@ -0,0 +1,108 @@
+#include "main.h"

+

+// uint64_t seconds() {

+//     struct timespec start;

+//     clock_gettime(CLOCK_MONOTONIC_RAW, &start);

+//     return start.tv_sec;

+// }

+

+int main() {

+    // uint64_t start = seconds();

+    // Allocate memory

+    char** file_names = (char**) malloc((NUM_IMAGES) * sizeof(char*));

+    for (uint8_t i=0; i<NUM_IMAGES; i++) {

+        file_names[i] = (char*) malloc(FILE_NAME_SIZE * sizeof(char));

+        // Populate file name array with file path

+        file_names[i] = (i < NUM_DNEG_IMAGES) ? strcpy(file_names[i], DNEG_FOLDER_PATH) : strcpy(file_names[i], DPOS_FOLDER_PATH);

+    }   

+    uint8_t*** data = (uint8_t***) malloc((NUM_IMAGES) * sizeof(uint8_t**));

+    for (uint8_t i=0; i<NUM_IMAGES; i++) {

+        data[i] = (uint8_t**) malloc(IMAGE_HEIGHT * sizeof(uint8_t*)); 

+        for (uint8_t j=0; j<IMAGE_HEIGHT; j++) data[i][j] = (uint8_t*)malloc(IMAGE_WIDTH * sizeof(uint8_t));

+    }

+

+    float** data_flatten = (float**) malloc((NUM_IMAGES) * sizeof(float*));

+    for (uint8_t i=0; i<NUM_IMAGES; i++) data_flatten[i] = (float*) malloc(FEATURE_VECTOR_SIZE * sizeof(float));

+

+    int32_t* labels = (int32_t*) malloc((NUM_IMAGES) * sizeof(int32_t));

+

+    uint8_t num_train = (uint8_t) (0.8 * NUM_IMAGES);

+    uint8_t num_test = NUM_IMAGES - num_train;

+

+    float** data_train = (float**) malloc (num_train * sizeof(float*));

+    float** data_test = (float**) malloc (num_test * sizeof(float*));

+

+    int32_t* labels_train = (int32_t*) malloc(num_train * sizeof(int32_t));

+    int32_t* labels_test = (int32_t*) malloc(num_test * sizeof(int32_t));

+

+    float** data_reduced_dimensions = (float**) malloc((NUM_IMAGES) * sizeof(float*));

+    for (uint32_t i=0; i<NUM_IMAGES; i++) data_reduced_dimensions[i] = (float*) malloc(NUM_COMPONENTS * sizeof(float));

+

+    float* alphas = (float*) malloc((NUM_IMAGES) * sizeof(float));

+

+    // Concatenate the file name array with the image name

+    read_dir(DNEG_FOLDER_PATH, DPOS_FOLDER_PATH, file_names);

+    printf("read\n");

+

+    // Read in an array of 2D images, then flatten

+    for (uint8_t i=0; i<NUM_IMAGES; i++) {

+        labels[i] = (i < NUM_DNEG_IMAGES) ? -1 : 1;

+        read_bmp(file_names[i], data[i]);

+        flatten_image(data[i], data_flatten[i]);

+    }

+

+    printf("flattened\n");

+

+    // PCA Transform

+    pca_covariance_method(data_flatten, NUM_IMAGES, FEATURE_VECTOR_SIZE, data_reduced_dimensions, NUM_COMPONENTS);

+    printf("pca\n");

+

+    // Shuffle data and labels

+    shuffle_data_labels(data_reduced_dimensions, labels, NUM_IMAGES);

+    printf("shuffled\n");

+

+    // Populate the training and testing data

+    for (uint8_t i=0; i<num_train; i++) {

+        data_train[i] = data_reduced_dimensions[i];

+        labels_train[i] = labels[i];

+    }

+

+    for (uint8_t i=num_train; i<NUM_IMAGES; i++) {

+        data_test[i - num_train] = data_reduced_dimensions[i];

+        labels_test[i - num_train] = labels[i];

+    }

+

+    printf("split\n");

+

+    // Train SVM

+    float bias;

+    train_svm(data_train, num_train, NUM_COMPONENTS, alphas, &bias, labels, PENALTY_PARAM);

+    printf("trained\n");

+

+    // Test SVM

+    float accuracy = test_svm(data_test, labels_test, num_test, NUM_COMPONENTS, alphas, bias);

+    printf("Accuracy: %.2f%%\n", accuracy);

+ 

+    // Free memory

+    for (uint8_t i=0; i<NUM_IMAGES; i++) {

+        free(file_names[i]);

+        free(data_flatten[i]);

+        free(data_reduced_dimensions[i]);

+        for (uint8_t j=0; j<IMAGE_HEIGHT; j++) free(data[i][j]);

+        free(data[i]);

+    }

+    free(file_names);

+    free(data);

+    free(data_flatten);

+    free(data_reduced_dimensions);

+    free(labels);

+    free(labels_train);

+    free(labels_test);

+    free(data_train);

+    free(data_test);

+    free(alphas);

+

+    // uint64_t end = seconds();

+    // printf("Time Elapsed: %lds\n", end-start);

+    return 0;

+}
\ No newline at end of file
diff --git a/OASIS/c/src/main.h b/OASIS/c/src/main.h
new file mode 100644
index 0000000..1aaee54
--- /dev/null
+++ b/OASIS/c/src/main.h
@@ -0,0 +1,28 @@
+#ifndef MAIN_H

+#define MAIN_H

+

+#include <stdint.h>

+#include <string.h>

+#include <stddef.h>

+#include <stdlib.h>

+#include <stdio.h>

+#include <time.h>

+

+#include "ekd_cv.h"

+#include "ekd_pca.h"

+#include "ekd_blas.h"

+#include "ekd_svm.h"

+

+#define NUM_DNEG_IMAGES 134

+#define NUM_DPOS_IMAGES 99

+#define NUM_IMAGES 233

+#define DNEG_FOLDER_PATH "/home/eric/URA/OASIS/dataset/bmp/dneg/"

+#define DPOS_FOLDER_PATH "/home/eric/URA/OASIS/dataset/bmp/dpos/"

+#define IMAGE_WIDTH 144

+#define IMAGE_HEIGHT 176

+#define FEATURE_VECTOR_SIZE 25344

+#define FILE_NAME_SIZE 150

+#define NUM_COMPONENTS 2

+#define PENALTY_PARAM 1

+

+#endif /* MAIN_H */
\ No newline at end of file