Exercise 3 for the course "Parallel and distributed systems" of THMMY in AUTH university.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

305 lines
9.4 KiB

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include "serial_declarations.h"
#include <stdbool.h>
extern int NUMBER_OF_POINTS;
extern int DIMENSIONS;
extern char* POINTS_FILENAME;
extern char* LABELS_FILENAME;
typedef struct parameters {
double epsilon;
bool verbose;
bool display;
} parameters;
//Function get_args parses command line arguments.
void get_args(int argc, char **argv, int *h);
//Function meanshift recursively shifts original points according to th
//mean-shift algorithm saving the result to shiftedPoints. Struct opt has user
//options, h is the desirable deviation, iteration is this call's iteration
//number.
int meanshift(double **original_points, double ***shifted_points, int h
, parameters *opt, int iteration);
//Function norm returns the second norm of matrix of dimensions rowsXcols.
double norm(double **matrix, int rows, int cols);
//Function multiply calculates the product of matrices 1 and 2 into output.
void multiply(double **matrix1, double **matrix2, double **output);
//Function calculateDistance returns the distance between x and y vectors.
double calculateDistance(double *y, double *x);
//Function alloc_2d_double allocates rows*cols bytes of continuous memory.
double **alloc_2d_double(int rows, int cols);
//Function duplicate copies the values of source array to dest array.
void duplicate(double **source, int rows, int cols, double ***dest);
//Function print_matrix prints array of dimensions rowsXcols to the console.
void print_matrix(double **array, int rows, int cols);
//Function save_matrix prints matrix in a csv file with path/filename
//"output/output_iteration". If a file already exists new lines are concatenated.
void save_matrix(double **matrix
, int iteration);
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;
//get_args(argc, argv, &h); commented out while in development
FILE *f;
// f = fopen(X, "rb");
// fseek(f, 0L, SEEK_END);
// long int pos = ftell(f);
// fclose(f);
// int elements = pos / sizeof(double); // number of total elements (points*dimension)
// int points = elements/DIMENSIONS;
// //printf("points : %d \n", points);
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(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);
long int pos = ftell(f);
rewind(f);
//printf("position : %ld \n", pos);
int label_elements = pos/ sizeof(char);
char *labels = (char*)malloc(label_elements* sizeof(char));
fseek(f, 0L, SEEK_SET);
int out = fread(labels, sizeof(char), label_elements, f);
fclose(f);
// MEAN SHIFT OPTIONS
parameters params;
params.epsilon = 0.0001;
params.verbose = false;
params.display = false;
parameters *opt;
opt = &params;
double **shifted_points;
// tic
gettimeofday (&startwtime, NULL);
int iterations = meanshift(vectors, &shifted_points, h, opt, 1);
// toc
gettimeofday (&endwtime, NULL);
seq_time = (double)((endwtime.tv_usec - startwtime.tv_usec)/1.0e6 + endwtime.tv_sec - startwtime.tv_sec);
printf("%s wall clock time = %f\n","Mean Shift", seq_time);
//TODO write output points to file -> plot later
//save_matrix(shifted_points, iterations);
}
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 **original_points, double ***shifted_points, int h
, parameters *opt, int iteration){
// allocates space and copies original points on first iteration
if (iteration == 1){
(*shifted_points) = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
duplicate(original_points, NUMBER_OF_POINTS, DIMENSIONS, shifted_points);
}
// mean shift vector
double **mean_shift_vector;
mean_shift_vector = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
// initialize elements of mean_shift_vector to inf
for (int i=0;i<NUMBER_OF_POINTS;i++){
for (int j=0;j<DIMENSIONS;j++){
mean_shift_vector[i][j] = DBL_MAX;
}
}
double **kernel_matrix = 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((*shifted_points)[i]
, original_points[j]);
if (i == j){
kernel_matrix[i][j] = 1;
} else if (dist < h*h){
kernel_matrix[i][j] = dist * dist;
// compute kernel matrix
double pow = ((-1)*(kernel_matrix[i][j]))/(2*(h*h));
kernel_matrix[i][j] = exp(pow);
} else {
kernel_matrix[i][j] = 0;
}
sum = sum + kernel_matrix[i][j];
}
denominator[i] = sum;
}
// create new y vector
double **new_shift = alloc_2d_double(NUMBER_OF_POINTS, DIMENSIONS);
// build nominator
multiply(kernel_matrix, original_points, new_shift);
// divide element-wise
for (int i=0; i<NUMBER_OF_POINTS; i++){
for (int j=0; j<DIMENSIONS; j++){
new_shift[i][j] = new_shift[i][j] / denominator[i];
// calculate mean-shift vector at the same time
mean_shift_vector[i][j] = new_shift[i][j] - (*shifted_points)[i][j];
}
}
// frees previously shifted points, they're now garbage
free((*shifted_points)[0]);
// updates shifted points pointer to the new array address
shifted_points = &new_shift;
save_matrix((*shifted_points), iteration);
double current_norm = norm(mean_shift_vector, NUMBER_OF_POINTS, DIMENSIONS);
printf("Iteration n. %d, error %f \n", iteration, current_norm);
// clean up this iteration's allocates
free(mean_shift_vector[0]);
free(mean_shift_vector);
free(kernel_matrix[0]);
free(kernel_matrix);
free(denominator);
/** iterate until convergence **/
if (current_norm > opt->epsilon) {
return meanshift(original_points, shifted_points, h, opt, ++iteration);
}
return iteration;
}
// TODO check why there's is a difference in the norm calculate in matlab
double norm(double **matrix, int rows, int cols){
double sum=0, temp_mul=0;
for (int i=0; i<rows; i++) {
for (int j=0; j<cols; j++) {
temp_mul = matrix[i][j] * matrix[i][j];
sum = sum + temp_mul;
}
}
double norm = sqrt(sum);
return norm;
}
void multiply(double **matrix1, double **matrix2, double **output){
// W dims are NUMBER_OF_POINTS NUMBER_OF_POINTS
// and x dims are NUMBER_OF_POINTS DIMENSIONS
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<NUMBER_OF_POINTS; k++){
output[i][j] += matrix1[i][k] * matrix2[k][j];
}
}
}
}
double calculateDistance(double *y, double *x){
double sum = 0, dif;
for (int i=0; i<DIMENSIONS; i++){
dif = y[i]-x[i];
sum += dif * dif;
}
double distance = sqrt(sum);
return distance;
}
double **alloc_2d_double(int rows, int cols) {
double *data = (double *) malloc(rows*cols*sizeof(double));
double **array = (double **) malloc(rows*sizeof(double*));
for (int i=0; i<rows; i++)
array[i] = &(data[cols*i]);
return array;
}
void duplicate(double **source, int rows, int cols, double ***dest){
for (int i=0; i<rows; i++){
for (int j=0; j<cols; j++){
(*dest)[i][j] = source[i][j];
}
}
}
void print_matrix(double **array, int rows, int cols){
for (int i=0; i<cols; i++){
for (int j=0; j<rows; j++){
printf("%f ", array[j][i]);
}
printf("\n");
}
}
void save_matrix(double **matrix, int iteration){
char filename[18];
snprintf(filename, sizeof(filename), "%s%d", "output/output_", iteration);
FILE *file;
file = fopen(filename, "w");
for (int rows=0; rows<NUMBER_OF_POINTS; ++rows){
for (int cols=0; cols<DIMENSIONS; ++cols){
fprintf(file, "%f", matrix[rows][cols]);
if (cols != DIMENSIONS - 1){
fprintf(file, ",");
}
}
fprintf(file, "\n");
}
}