From 51db464d45ca9799eb2ee0840b64d7626dc77a59 Mon Sep 17 00:00:00 2001
From: Apostolof <apotwohd@gmail.com>
Date: Sun, 21 Jan 2018 02:17:07 +0200
Subject: [PATCH] Change meanshift func to recursive Change constants to
 variables Init arguments parsing

---
 serial.c             |  33 +++++-----
 serialDeclarations.c | 145 ++++++++++++++++++++++++-------------------
 serialDeclarations.h |  12 ++--
 3 files changed, 106 insertions(+), 84 deletions(-)

diff --git a/serial.c b/serial.c
index faafbc2..67b263a 100755
--- a/serial.c
+++ b/serial.c
@@ -4,17 +4,18 @@
 
 #include "serialDeclarations.h"
 
+int NUMBER_OF_POINTS = 600;
+int DIMENSIONS = 2;
+char* POINTS_FILENAME = "data/X.bin";
+char* LABELS_FILENAME = "data/L.bin";
+
 struct timeval startwtime, endwtime;
 double seq_time;
 
 int main(int argc, char **argv){
+	int h = 1;
 
-//    if (argc<2){
-//        printf("%s\n", "Specify the k");
-//        return 1;
-//    }
-//    = atoi(argv[1]);  // the k-parameter
-
+	//get_args(argc, argv, &h); commented out while in development
 
     FILE *f;
 //    f = fopen(X, "rb");
@@ -22,19 +23,19 @@ int main(int argc, char **argv){
 //    long int pos = ftell(f);
 //    fclose(f);
 //    int elements = pos / sizeof(double);  // number of total elements (points*dimension)
-//    int points = elements/COLUMNS;
+//    int points = elements/DIMENSIONS;
 //    //printf("points : %d \n", points);
-    f = fopen(X, "rb");
-    double ** vectors;
-    vectors = alloc_2d_double(ROWS, COLUMNS);
-    for (int i=0; i<ROWS; i++){
-        int out = fread(vectors[i], sizeof(double), COLUMNS, f);
+    f = fopen(POINTS_FILENAME, "rb");
+    double **vectors;
+    vectors = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
+    for (int i=0; i<NUMBER_OF_POINTS; i++){
+        int out = fread(vectors[i], sizeof(double), DIMENSIONS, f);
     }
 
     save_matrix(vectors, 0);
 
     // initializing file that will contain the labels (train)
-    f = fopen(L, "rb");
+    f = fopen(LABELS_FILENAME, "rb");
     // NOTE : Labels were classified as <class 'numpy.uint8'>
     // variables of type uint8 are stored as 1-byte (8-bit) unsigned integers
     fseek(f, 0L, SEEK_END);
@@ -48,7 +49,6 @@ int main(int argc, char **argv){
     fclose(f);
 
     // MEAN SHIFT OPTIONS
-    int h = 1;
     parameters params;
     params.epsilon = 0.0001;
     params.verbose = false;
@@ -56,10 +56,11 @@ int main(int argc, char **argv){
     parameters *opt;
     opt = &params;
 
+    double **shiftedPoints;
     // tic
     gettimeofday (&startwtime, NULL);
 
-    meanshift(vectors, h, opt);
+    int iterations = meanshift(vectors, &shiftedPoints, h, opt, 1);
 
     // toc
     gettimeofday (&endwtime, NULL);
@@ -67,5 +68,5 @@ int main(int argc, char **argv){
     printf("%s wall clock time = %f\n","Mean Shift", seq_time);
 
     //TODO write output points to file -> plot later
-
+    //save_matrix(vectors, iterations);
 }
\ No newline at end of file
diff --git a/serialDeclarations.c b/serialDeclarations.c
index 4947958..7cfe344 100644
--- a/serialDeclarations.c
+++ b/serialDeclarations.c
@@ -6,76 +6,94 @@
 
 #include "serialDeclarations.h"
 
-void meanshift(double **originalPoints, int h, parameters *opt){
+void get_args(int argc, char **argv, int *h){
+    if (argc != 6) {
+        printf("Usage: %s h N D Pd Pl\nwhere:\n", argv[0]);
+        printf("\th is the variance\n");
+        printf("\tN is the the number of points\n");
+        printf("\tD is the number of dimensions of each point\n");
+        printf("\tPd is the path of the dataset file\n");
+        printf("\tPl is the path of the labels file\n");
+        exit(1);
+    }
+
+    *h = atoi(argv[1]);
+    NUMBER_OF_POINTS = atoi(argv[2]);
+    DIMENSIONS = atoi(argv[3]);
+    POINTS_FILENAME = argv[4];
+    LABELS_FILENAME = argv[5];
+}
+
+int meanshift(double **originalPoints, double ***shiftedPoints, int h
+    , parameters *opt, int iteration){
 
-    double **y;
-    y = alloc_2d_double(ROWS, COLUMNS);
-    y = duplicate(originalPoints, y, ROWS, COLUMNS);
+    // allocates space and copies original points on first iteration
+    if (iteration == 1){
+        (*shiftedPoints) = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
+        (*shiftedPoints) = duplicate(originalPoints, (*shiftedPoints)
+            , NUMBER_OF_POINTS, DIMENSIONS);
+    }
 
     // mean shift vector
     double **meanShiftVector;
-    meanShiftVector = alloc_2d_double(ROWS, COLUMNS);
+    meanShiftVector = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
     // initialize elements of meanShiftVector to inf
-    for (int i=0;i<ROWS;i++){
-        for (int j=0;j<COLUMNS;j++){
+    for (int i=0;i<NUMBER_OF_POINTS;i++){
+        for (int j=0;j<DIMENSIONS;j++){
             meanShiftVector[i][j] = DBL_MAX;
         }
     }
 
-    // initialize iteration counter
-    int iter = 0;
-
-    // printf("%f \n", opt->epsilon);
-
-    double ** kernelMatrix = alloc_2d_double(ROWS, ROWS);
-    double *denominator = malloc(ROWS * sizeof(double));
-
-    /** iterate until convergence **/
-    // printf("norm : %f \n", norm(m, ROWS, COLUMNS));
-    while (norm(meanShiftVector, ROWS, COLUMNS) > opt->epsilon) {
-        iter = iter +1;
-        // find pairwise distance matrix (inside radius)
-        // [I, D] = rangesearch(x,y,h);
-        for (int i=0; i<ROWS; i++){
-            double sum =0;
-            for (int j=0; j<ROWS; j++){
-                double dist = calculateDistance(y[i],originalPoints[j]);
-
-                if (i == j){
-                    kernelMatrix[i][j] = 1;
-                } else if (dist < h*h){
-                    kernelMatrix[i][j] = dist * dist;
-                    // compute kernel matrix
-                    double pow = ((-1)*(kernelMatrix[i][j]))/(2*(h*h));
-                    kernelMatrix[i][j] = exp(pow);
-                } else {
-                    kernelMatrix[i][j] = 0;
-                }
-                sum = sum + kernelMatrix[i][j];
+    double **kernelMatrix = alloc_2d_double(NUMBER_OF_POINTS, NUMBER_OF_POINTS);
+    double *denominator = malloc(NUMBER_OF_POINTS * sizeof(double));
+
+    // find pairwise distance matrix (inside radius)
+    // [I, D] = rangesearch(x,y,h);
+    for (int i=0; i<NUMBER_OF_POINTS; i++){
+        double sum =0;
+        for (int j=0; j<NUMBER_OF_POINTS; j++){
+            double dist = calculateDistance((*shiftedPoints)[i],originalPoints[j]);
+
+            if (i == j){
+                kernelMatrix[i][j] = 1;
+            } else if (dist < h*h){
+                kernelMatrix[i][j] = dist * dist;
+                // compute kernel matrix
+                double pow = ((-1)*(kernelMatrix[i][j]))/(2*(h*h));
+                kernelMatrix[i][j] = exp(pow);
+            } else {
+                kernelMatrix[i][j] = 0;
             }
-            denominator[i] = sum;
+            sum = sum + kernelMatrix[i][j];
         }
+        denominator[i] = sum;
+    }
 
-        // create new y vector
-        double** y_new = alloc_2d_double(ROWS, COLUMNS);
-
-        multiply(kernelMatrix, originalPoints, y_new);
-        // divide element-wise
-        for (int i=0; i<ROWS; i++){
-            for (int j=0; j<COLUMNS; j++){
-                y_new[i][j] = y_new[i][j] / denominator[i];
-                // calculate mean-shift vector
-                meanShiftVector[i][j] = y_new[i][j] - y[i][j];
-                // update y
-                y[i][j] = y_new[i][j];
-            }
+    // create new y vector
+    double **y_new = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
+
+    multiply(kernelMatrix, originalPoints, y_new);
+    // divide element-wise
+    for (int i=0; i<NUMBER_OF_POINTS; i++){
+        for (int j=0; j<DIMENSIONS; j++){
+            y_new[i][j] = y_new[i][j] / denominator[i];
+            // calculate mean-shift vector
+            meanShiftVector[i][j] = y_new[i][j] - (*shiftedPoints)[i][j];
         }
+    }
+    shiftedPoints = &y_new;
+
+    save_matrix((*shiftedPoints), iteration);
 
-        save_matrix(y, iter);
+    double current_norm = norm(meanShiftVector, NUMBER_OF_POINTS, DIMENSIONS);
+    printf("Iteration n. %d, error %f \n", iteration, current_norm);
 
-        printf("Iteration n. %d, error %f \n", iter, norm(meanShiftVector, ROWS, COLUMNS));
-        // TODO maybe keep y for live display later?
+    /** iterate until convergence **/
+    if (current_norm > opt->epsilon) {
+        return meanshift(originalPoints, shiftedPoints, h, opt, ++iteration);
     }
+
+    return iteration;
 }
 
 // TODO check why there's is a difference in the norm calculate in matlab
@@ -92,12 +110,13 @@ double norm(double **matrix, int rows, int cols){
 }
 
 void multiply(double **matrix1, double **matrix2, double **output){
-    // W dims are ROWS ROWS and x dims are ROWS COLUMNS
+    // W dims are NUMBER_OF_POINTS NUMBER_OF_POINTS
+    // and x dims are NUMBER_OF_POINTS DIMENSIONS
 
-    for (int i=0; i<ROWS; i++){
-        for (int j=0; j<COLUMNS; j++){
+    for (int i=0; i<NUMBER_OF_POINTS; i++){
+        for (int j=0; j<DIMENSIONS; j++){
             output[i][j] = 0;
-            for (int k=0; k<ROWS; k++){
+            for (int k=0; k<NUMBER_OF_POINTS; k++){
                 output[i][j] += matrix1[i][k] * matrix2[k][j];
             }
         }
@@ -106,7 +125,7 @@ void multiply(double **matrix1, double **matrix2, double **output){
 
 double calculateDistance(double *y, double *x){
     double sum = 0, dif;
-    for (int i=0; i<COLUMNS; i++){
+    for (int i=0; i<DIMENSIONS; i++){
         dif = y[i]-x[i];
         sum += dif * dif;
     }
@@ -142,15 +161,15 @@ void print_matrix(double **array, int rows, int cols){
     }
 }
 
-void save_matrix(double **matrix,int iteration){
+void save_matrix(double **matrix, int iteration){
     char filename[18];
     snprintf(filename, sizeof(filename), "%s%d", "output/output_", iteration);
     FILE *iterOutput;
     iterOutput = fopen(filename, "w");
-    for (int rows=0; rows<ROWS; ++rows){
-        for (int cols=0; cols<COLUMNS; ++cols){
+    for (int rows=0; rows<NUMBER_OF_POINTS; ++rows){
+        for (int cols=0; cols<DIMENSIONS; ++cols){
             fprintf(iterOutput, "%f", matrix[rows][cols]);
-            if (cols != COLUMNS - 1){
+            if (cols != DIMENSIONS - 1){
                 fprintf(iterOutput, ",");
             }
         }
diff --git a/serialDeclarations.h b/serialDeclarations.h
index a43679e..b70ad47 100644
--- a/serialDeclarations.h
+++ b/serialDeclarations.h
@@ -3,10 +3,10 @@
 
 #include <stdbool.h>
 
-#define X "data/X.bin"
-#define L "data/L.bin"
-#define COLUMNS     2
-#define ROWS        600
+extern int NUMBER_OF_POINTS;
+extern int DIMENSIONS;
+extern char* POINTS_FILENAME;
+extern char* LABELS_FILENAME;
 
 typedef struct parameters {
     double epsilon;
@@ -14,7 +14,9 @@ typedef struct parameters {
     bool display;
 } parameters;
 
-void meanshift(double **x, int h, struct parameters *opt);
+void get_args(int argc, char **argv, int *h);
+int meanshift(double **originalPoints, double ***shiftedPoints, int h
+    , parameters *opt, int iteration);
 double norm(double ** m, int rows, int cols);
 void multiply(double ** matrix1, double ** matrix2, double ** output);
 double calculateDistance(double *, double *);