added working H3, refactored files into tree, adp and utilities

.ipynb_checkpoints/var-checkpoint.py

+import numpy as np
+
+d = np.fromfile("../norm_data/std_LR_091_0000", dtype=np.float32)
+print(d.shape)
+d = d.reshape((d.shape[0]//5,5))
+print(np.cov(d.T))

Makefile

@@ -4,7 +4,7 @@ LDFLAGS=-lm
 
 all: main
 
-obj=src/main/main.c src/tree/tree.c src/common/common.c src/tree/kdtreeV2.c src/tree/heap.c
+obj=src/main/main.c src/tree/tree.c src/common/common.c src/tree/kdtreeV2.c src/tree/heap.c src/adp/adp.c
 
 main: ${obj} 
 	${CC} ${CFLAGS} ${LDFLAGS} ${obj} -o $@

src/adp/adp.h

+#pragma once
+#include "../tree/tree.h"
+
+#define PREALLOC_BORDERS 100
+#define MAX_SERIAL_MERGING 100
+#define PRINT_H2_COMM_SCHEME
+
+typedef struct border_t 
+{
+    float_t density;
+    float_t error;
+    idx_t idx;
+} border_t;
+
+typedef struct sparse_border_t 
+{
+    idx_t i;
+    idx_t j;
+    idx_t idx;
+    float_t density;
+    float_t error;
+} sparse_border_t;
+
+typedef struct adj_list_t 
+{
+    idx_t count;
+    idx_t size;
+    struct sparse_border_t* data;
+} adj_list_t;
+
+
+typedef struct clusters_t 
+{
+    int use_sparse_borders;
+    struct adj_list_t *sparse_borders;
+    struct lu_dynamic_array_t centers;
+    struct border_t **borders;
+    struct border_t *__borders_data;
+    idx_t n;
+} clusters_t;
+
+typedef struct merge_t 
+{
+    idx_t source;
+    idx_t target;
+    float_t density;
+} merge_t;
+
+
+
+void compute_density_kstarnn_rma(global_context_t* ctx, const float_t d, int verbose);
+void compute_density_kstarnn_rma_v2(global_context_t* ctx, const float_t d, int verbose);
+float_t compute_ID_two_NN_ML(global_context_t* ctx, datapoint_info_t* dp_info, idx_t n, int verbose);
+void clusters_allocate(clusters_t * c, int s);
+
+clusters_t Heuristic1(global_context_t *ctx, int verbose);
+void Heuristic2(global_context_t* ctx, clusters_t* cluster);
+void Heuristic3(global_context_t* ctx, clusters_t* cluster, float_t Z, int halo);

src/common/common.c

@@ -202,3 +202,54 @@ void lu_dynamic_array_init(lu_dynamic_array_t * a)
     a -> size = 0;
 }
 
+void ordered_buffer_to_file(global_context_t* ctx, void* buffer, size_t el_size, uint64_t n, const char* fname)
+{
+    //MPI_Barrier(ctx -> mpi_communicator);
+    MPI_DB_PRINT("[MASTER] writing to file %s\n", fname);
+    void* tmp_data; 
+    int* ppp; 
+    int* displs;
+
+    MPI_Barrier(ctx -> mpi_communicator);
+    
+    uint64_t tot_n = 0;
+    MPI_Reduce(&n, &tot_n, 1, MPI_UINT64_T , MPI_SUM, 0, ctx -> mpi_communicator);
+
+    if(I_AM_MASTER) 
+    {
+        tmp_data = (void*)MY_MALLOC(el_size * tot_n );
+        ppp      = (int*)MY_MALLOC(ctx -> world_size * sizeof(int));
+        displs   = (int*)MY_MALLOC(ctx -> world_size * sizeof(int));
+
+    }
+    
+    int nn = (int)n;
+    MPI_Gather(&nn, 1, MPI_INT, ppp, 1, MPI_INT, 0, ctx -> mpi_communicator);
+
+    if(I_AM_MASTER)
+    {
+        displs[0] = 0;
+        for(int i = 0; i < ctx -> world_size; ++i) ppp[i]    = el_size  * ppp[i];
+        for(int i = 1; i < ctx -> world_size; ++i) displs[i] = displs[i - 1] + ppp[i - 1];
+            
+    }
+
+    MPI_Gatherv(buffer, (int)(el_size * n), 
+            MPI_BYTE, tmp_data, ppp, displs, MPI_BYTE, 0, ctx -> mpi_communicator);
+
+    if(I_AM_MASTER)
+    {
+        FILE* file = fopen(fname,"w");
+        if(!file)
+        {
+            printf("Cannot open file %s ! Aborting \n", fname);
+        }
+        fwrite(tmp_data, 1, el_size * tot_n, file);
+        fclose(file);
+        free(tmp_data);
+        free(ppp);
+        free(displs);
+
+    }
+    MPI_Barrier(ctx -> mpi_communicator);
+}

src/common/common.h

@@ -29,6 +29,15 @@ typedef struct datapoint_info_t {
 	#define float_t double
 #endif
 
+#ifdef USE_FLOAT32
+#define MPI_MY_FLOAT MPI_FLOAT
+#else
+#define MPI_MY_FLOAT MPI_DOUBLE
+#endif
+
+
+#define I_AM_MASTER ctx->mpi_rank == 0
+
 
 #define MY_TRUE  1
 #define MY_FALSE 0
@@ -173,4 +182,4 @@ void lu_dynamic_array_reserve(lu_dynamic_array_t * a, idx_t n);
 void lu_dynamic_array_init(lu_dynamic_array_t * a);
 void print_error_code(int err);
 
-
+void ordered_data_to_file(global_context_t* ctx);

src/main/main.c

@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include "../common/common.h"
 #include "../tree/tree.h"
+#include "../adp/adp.h"
 
 //
 
@@ -89,3 +90,206 @@ int main(int argc, char** argv) {
 
 
 
+void simulate_master_read_and_scatter(int dims, size_t n, global_context_t *ctx) 
+{
+    float_t *data;
+    TIME_DEF
+    double elapsed_time;
+
+    if (ctx->mpi_rank == 0) 
+    {
+        //data = read_data_file(ctx, "../norm_data/50_blobs_more_var.npy", MY_TRUE);
+        //ctx->dims = 2;
+        //data = read_data_file(ctx, "../norm_data/50_blobs.npy", MY_TRUE);
+        // std_g0163178_Me14_091_0000
+    
+        // 190M points
+        // std_g2980844_091_0000
+        // data = read_data_file(ctx,"../norm_data/std_g2980844_091_0000",MY_TRUE);
+        
+        /* 1M points ca.*/
+        data = read_data_file(ctx,"../norm_data/std_LR_091_0001",MY_TRUE);
+
+        /* BOX */
+        // data = read_data_file(ctx,"../norm_data/std_Box_256_30_092_0000",MY_TRUE);
+
+        /* 8M points */
+        
+        // data = read_data_file(ctx,"../norm_data/std_g0144846_Me14_091_0001",MY_TRUE);
+
+        //88M 
+        //data = read_data_file(ctx,"../norm_data/std_g5503149_091_0000",MY_TRUE);
+
+        //
+        //34 M
+        //data = read_data_file(ctx,"../norm_data/std_g1212639_091_0001",MY_TRUE);
+        ctx->dims = 5;
+
+        // ctx -> n_points = 5 * 100000;
+        ctx->n_points = ctx->n_points / ctx->dims;
+        //ctx->n_points = (ctx->n_points * 5) / 10;
+        // ctx -> n_points = ctx -> world_size * 1000;
+
+        //ctx -> n_points = 10000000 * ctx -> world_size;
+        //generate_random_matrix(&data, ctx -> dims, ctx -> n_points, ctx);
+        //mpi_printf(ctx, "Read %lu points in %u dims\n", ctx->n_points, ctx->dims);
+    }
+    //MPI_DB_PRINT("[MASTER] Reading file and scattering\n");
+    
+    /* communicate the total number of points*/
+    MPI_Bcast(&(ctx->dims), 1, MPI_UINT32_T, 0, ctx->mpi_communicator);
+    MPI_Bcast(&(ctx->n_points), 1, MPI_UINT64_T, 0, ctx->mpi_communicator);
+
+    /* compute the number of elements to recieve for each processor */
+    int *send_counts = (int *)MY_MALLOC(ctx->world_size * sizeof(int));
+    int *displacements = (int *)MY_MALLOC(ctx->world_size * sizeof(int));
+
+    displacements[0] = 0;
+    send_counts[0] = ctx->n_points / ctx->world_size;
+    send_counts[0] += (ctx->n_points % ctx->world_size) > 0 ? 1 : 0;
+    send_counts[0] = send_counts[0] * ctx->dims;
+
+    for (int p = 1; p < ctx->world_size; ++p) 
+    {
+        send_counts[p] = (ctx->n_points / ctx->world_size);
+        send_counts[p] += (ctx->n_points % ctx->world_size) > p ? 1 : 0;
+        send_counts[p] = send_counts[p] * ctx->dims;
+        displacements[p] = displacements[p - 1] + send_counts[p - 1];
+    }
+
+
+    ctx->local_n_points = send_counts[ctx->mpi_rank] / ctx->dims;
+
+    float_t *pvt_data = (float_t *)MY_MALLOC(send_counts[ctx->mpi_rank] * sizeof(float_t));
+
+    MPI_Scatterv(data, send_counts, displacements, MPI_MY_FLOAT, pvt_data, send_counts[ctx->mpi_rank], MPI_MY_FLOAT, 0, ctx->mpi_communicator);
+
+    ctx->local_data = pvt_data;
+
+    int k_local = 20;
+    int k_global = 20;
+
+    uint64_t *global_bin_counts_int = (uint64_t *)MY_MALLOC(k_global * sizeof(uint64_t));
+
+    pointset_t original_ps;
+    original_ps.data = ctx->local_data;
+    original_ps.dims = ctx->dims;
+    original_ps.n_points = ctx->local_n_points;
+    original_ps.lb_box = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t));
+    original_ps.ub_box = (float_t*)MY_MALLOC(ctx -> dims * sizeof(float_t));
+
+    float_t tol = 0.002;
+
+    top_kdtree_t tree;
+    TIME_START;
+    top_tree_init(ctx, &tree);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("Initializing global kdtree", elapsed_time);
+
+    TIME_START;
+    build_top_kdtree(ctx, &original_ps, &tree, k_global, tol);
+    exchange_points(ctx, &tree);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("Top kdtree build and domain decomposition", elapsed_time);
+    //test_the_idea(ctx);
+
+    TIME_START;
+    kdtree_v2 local_tree;
+    kdtree_v2_init( &local_tree, ctx -> local_data, ctx -> local_n_points, (unsigned int)ctx -> dims);
+    int k = 300;
+    //int k = 30;
+
+    datapoint_info_t* dp_info = (datapoint_info_t*)MY_MALLOC(ctx -> local_n_points * sizeof(datapoint_info_t));            
+    /* initialize, to cope with valgrind */
+    for(uint64_t i = 0; i < ctx -> local_n_points; ++i)
+    {
+        dp_info[i].ngbh.data = NULL;
+        dp_info[i].ngbh.N = 0;
+        dp_info[i].ngbh.count = 0;
+        dp_info[i].g = 0.f;
+        dp_info[i].log_rho = 0.f;
+        dp_info[i].log_rho_c = 0.f;
+        dp_info[i].log_rho_err = 0.f;
+        dp_info[i].array_idx = -1;
+        dp_info[i].kstar = -1;
+        dp_info[i].is_center = -1;
+        dp_info[i].cluster_idx = -1;
+    }
+
+    build_local_tree(ctx, &local_tree);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("Local trees init and build", elapsed_time);
+
+    TIME_START;
+    MPI_DB_PRINT("----- Performing ngbh search -----\n");
+    MPI_Barrier(ctx -> mpi_communicator);
+
+    mpi_ngbh_search(ctx, dp_info, &tree, &local_tree, ctx -> local_data, k);
+
+    MPI_Barrier(ctx -> mpi_communicator);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("Total time for all knn search", elapsed_time)
+
+
+
+    TIME_START;
+    //float_t id = id_estimate(ctx, dp_info, ctx -> local_n_points, 0.9, MY_FALSE);
+    float_t id = compute_ID_two_NN_ML(ctx, dp_info, ctx -> local_n_points, MY_FALSE);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("ID estimate", elapsed_time)
+
+    MPI_DB_PRINT("ID %lf \n",id);
+
+    TIME_START;
+
+    float_t  z = 2;
+
+    ctx -> local_datapoints = dp_info;
+    //compute_density_kstarnn_rma(ctx, id, MY_FALSE);
+    compute_density_kstarnn_rma_v2(ctx, id, MY_FALSE);
+    compute_correction(ctx, z);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("Density estimate", elapsed_time)
+
+    TIME_START;
+    clusters_t clusters = Heuristic1(ctx, MY_FALSE);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("H1", elapsed_time)
+
+    TIME_START;
+    clusters_allocate(&clusters, 1);
+    Heuristic2(ctx, &clusters);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("H2", elapsed_time)
+
+
+    TIME_START;
+    int halo = 0;
+    Heuristic3(ctx, &clusters, z, halo);
+    elapsed_time = TIME_STOP;
+    LOG_WRITE("H3", elapsed_time)
+
+
+    /* find density */ 
+    #if defined (WRITE_NGBH)
+        ordered_data_to_file(ctx);
+    #endif
+
+    /*
+    free(foreign_dp_info);
+    */
+    
+    
+    top_tree_free(ctx, &tree);
+    kdtree_v2_free(&local_tree);
+
+    free(send_counts);
+    free(displacements);
+    //free(dp_info);
+
+    if (ctx->mpi_rank == 0) free(data);
+
+    original_ps.data = NULL;
+    free_pointset(&original_ps);
+    free(global_bin_counts_int);
+}

src/tree/heap.h

@@ -8,7 +8,7 @@
 #define T double
 #define DATA_DIMS 0 
 
-#define HERE printf("%d reached line %d\n", ctx -> mpi_rank, __LINE__);
+#define HERE printf("%d reached line %d\n", ctx -> mpi_rank, __LINE__); MPI_Barrier(ctx -> mpi_communicator);
 
 #ifdef USE_FLOAT32
     #define FLOAT_TYPE float

src/tree/tree.h

@@ -2,6 +2,7 @@
 #include "../common/common.h"
 #include "heap.h"
 #include "kdtreeV2.h"
+
 #include <stdlib.h>
 #include <stdio.h>
 #include <memory.h>
@@ -84,41 +85,24 @@ typedef struct top_kdtree_t
 	size_t _capacity;
 	struct top_kdtree_node_t* _nodes;
 	struct top_kdtree_node_t* root;
-
 } top_kdtree_t;
 
-typedef struct border_t {
-  float_t density;
-  float_t error;
-  idx_t idx;
-} border_t;
-
-typedef struct sparse_border_t {
-  idx_t i;
-  idx_t j;
-  idx_t idx;
-  float_t density;
-  float_t error;
-} sparse_border_t;
-
-typedef struct adj_list_t {
-  idx_t count;
-  idx_t size;
-  struct sparse_border_t* data;
-} adj_list_t;
 
 
-typedef struct clusters_t {
-  int use_sparse_borders;
-  struct adj_list_t *sparse_borders;
-  struct lu_dynamic_array_t centers;
-  struct border_t **borders;
-  struct border_t *__borders_data;
-  idx_t n;
-} clusters_t;
+void     simulate_master_read_and_scatter(int, size_t, global_context_t* ); 
+float_t* read_data_file(global_context_t *ctx, const char *fname, const int file_in_float32);
+void     top_tree_init(global_context_t *ctx, top_kdtree_t *tree); 
+void     build_top_kdtree(global_context_t *ctx, pointset_t *og_pointset, top_kdtree_t *tree, int n_bins, float_t tolerance);
+void     exchange_points(global_context_t* ctx, top_kdtree_t* tree);
+void     build_local_tree(global_context_t* ctx, kdtree_v2* local_tree);
+void     mpi_ngbh_search(global_context_t* ctx, datapoint_info_t* dp_info, top_kdtree_t* top_tree, kdtree_v2* local_tree, float_t* data, int k);
+float_t  compute_ID_two_NN_ML(global_context_t* ctx, datapoint_info_t* dp_info, idx_t n, int verbose);
+void     compute_density_kstarnn_rma_v2(global_context_t* ctx, const float_t d, int verbose);
+void     compute_correction(global_context_t* ctx, float_t Z);
+int      foreign_owner(global_context_t*, idx_t);
+void     top_tree_free(global_context_t *ctx, top_kdtree_t *tree);
 
 
-void simulate_master_read_and_scatter(int, size_t, global_context_t* );