In memory one off

src/Makefile

@@ -135,6 +135,9 @@ install: libssu.so ssu faithpd
 	mkdir -p ${PREFIX}/include/unifrac
 	cp task_parameters.hpp ${PREFIX}/include/unifrac/
 	cp api.hpp ${PREFIX}/include/unifrac/
+	cp biom.hpp ${PREFIX}/include/unifrac/
+	cp tree.hpp ${PREFIX}/include/unifrac/
+	cp biom_interface.hpp ${PREFIX}/include/unifrac/
 
 rapi_test: main
 	mkdir -p ~/.R

src/api.cpp

@@ -48,22 +48,25 @@
                                               return err;                                                      \
                                           }
 
-#define PARSE_SYNC_TREE_TABLE(tree_filename, table_filename) std::ifstream ifs(tree_filename);                                        \
-                                                             std::string content = std::string(std::istreambuf_iterator<char>(ifs),   \
-                                                                                               std::istreambuf_iterator<char>());     \
-                                                             su::BPTree tree = su::BPTree(content);                                   \
-                                                             su::biom table = su::biom(biom_filename);                                \
-                                                             if(table.n_samples <= 0 | table.n_obs <= 0) {                            \
-                                                                 return table_empty;                                                  \
-                                                             }                                                                        \
-                                                             std::string bad_id = su::test_table_ids_are_subset_of_tree(table, tree); \
-                                                             if(bad_id != "") {                                                       \
-                                                                 return table_and_tree_do_not_overlap;                                \
-                                                             }                                                                        \
-                                                             std::unordered_set<std::string> to_keep(table.obs_ids.begin(),           \
-                                                                                                     table.obs_ids.end());            \
-                                                             su::BPTree tree_sheared = tree.shear(to_keep).collapse();
-
+#define SYNC_TREE_TABLE(tree, table) std::unordered_set<std::string> to_keep(table.obs_ids.begin(),           \
+                                                                             table.obs_ids.end());            \
+                                     su::BPTree tree_sheared = tree.shear(to_keep).collapse();
+
+#define PARSE_TREE_TABLE(tree_filename, table_filename) std::ifstream ifs(tree_filename);                                        \
+                                                        std::string content = std::string(std::istreambuf_iterator<char>(ifs),   \
+                                                                                          std::istreambuf_iterator<char>());     \
+                                                        su::BPTree tree = su::BPTree(content);                                   \
+                                                        su::biom table = su::biom(biom_filename);                                \
+                                                        if(table.n_samples <= 0 | table.n_obs <= 0) {                            \
+                                                            return table_empty;                                                  \
+                                                        }                                                                        \
+                                                        std::string bad_id = su::test_table_ids_are_subset_of_tree(table, tree); \
+                                                        if(bad_id != "") {                                                       \
+                                                            return table_and_tree_do_not_overlap;                                \
+                                                        }   
+
+#define PARSE_SYNC_TREE_TABLE(tree_filename, table_filename) PARSE_TREE_TABLE(tree_filename, table_filename) \
+                                                             SYNC_TREE_TABLE(tree, table)
 
 using namespace su;
 using namespace std;
@@ -115,7 +118,7 @@ void initialize_mat(mat_t* &result, biom &table, bool is_upper_triangle) {
     }
 }
 
-void initialize_results_vec(r_vec* &result, biom& table){
+void initialize_results_vec(r_vec* &result, biom &table){
     // Stores results for Faith PD
     result = (r_vec*)malloc(sizeof(results_vec));
     result->n_samples = table.n_samples;
@@ -401,16 +404,13 @@ compute_status faith_pd_one_off(const char* biom_filename, const char* tree_file
     return okay;
 }
 
-compute_status one_off(const char* biom_filename, const char* tree_filename,
-                       const char* unifrac_method, bool variance_adjust, double alpha,
-                       bool bypass_tips, unsigned int nthreads, mat_t** result) {
-
-    CHECK_FILE(biom_filename, table_missing)
-    CHECK_FILE(tree_filename, tree_missing)
+compute_status one_off_inmem(su::biom *table, su::BPTree *tree,
+                             const char* unifrac_method, bool variance_adjust, double alpha,
+                             bool bypass_tips, unsigned int nthreads, mat_t** result) {
     SET_METHOD(unifrac_method, unknown_method)
-    PARSE_SYNC_TREE_TABLE(tree_filename, table_filename)
+    SYNC_TREE_TABLE((*tree), (*table))
 
-    const unsigned int stripe_stop = (table.n_samples + 1) / 2;
+    const unsigned int stripe_stop = (table->n_samples + 1) / 2;
     std::vector<double*> dm_stripes(stripe_stop);
     std::vector<double*> dm_stripes_total(stripe_stop);
 
@@ -422,19 +422,28 @@ compute_status one_off(const char* biom_filename, const char* tree_filename,
     std::vector<su::task_parameters> tasks(nthreads);
     std::vector<std::thread> threads(nthreads);
 
-    set_tasks(tasks, alpha, table.n_samples, 0, stripe_stop, bypass_tips, nthreads);
-    su::process_stripes(table, tree_sheared, method, variance_adjust, dm_stripes, dm_stripes_total, threads, tasks);
+    set_tasks(tasks, alpha, table->n_samples, 0, stripe_stop, bypass_tips, nthreads);
+    su::process_stripes(*table, tree_sheared, method, variance_adjust, dm_stripes, dm_stripes_total, threads, tasks);
 
-    initialize_mat(*result, table, true);  // true -> is_upper_triangle
+    initialize_mat(*result, *table, true);  // true -> is_upper_triangle
     for(unsigned int tid = 0; tid < threads.size(); tid++) {
-        su::stripes_to_condensed_form(dm_stripes,table.n_samples,(*result)->condensed_form,tasks[tid].start,tasks[tid].stop);
+        su::stripes_to_condensed_form(dm_stripes,table->n_samples,(*result)->condensed_form,tasks[tid].start,tasks[tid].stop);
     }
 
-    destroy_stripes(dm_stripes, dm_stripes_total, table.n_samples, 0, 0);
+    destroy_stripes(dm_stripes, dm_stripes_total, table->n_samples, 0, 0);
 
     return okay;
 }
 
+compute_status one_off(const char* biom_filename, const char* tree_filename,
+                       const char* unifrac_method, bool variance_adjust, double alpha,
+                       bool bypass_tips, unsigned int nthreads, mat_t** result) {
+    CHECK_FILE(biom_filename, table_missing)
+    CHECK_FILE(tree_filename, tree_missing)
+    PARSE_TREE_TABLE(tree_filename, table_filename)
+    return one_off_inmem(&table, &tree, unifrac_method, variance_adjust, alpha, bypass_tips, nthreads, result);
+}
+
 // TMat mat_full_fp32_t
 template<class TReal, class TMat>
 compute_status one_off_matrix_T(const char* biom_filename, const char* tree_filename,

src/api.hpp

@@ -2,6 +2,8 @@
 
 #ifdef __cplusplus
 #include <vector>
+#include "tree.hpp"
+#include "biom.hpp"
 #define EXTERN extern "C"
 
 
@@ -154,6 +156,30 @@ EXTERN ComputeStatus one_off(const char* biom_filename, const char* tree_filenam
                              const char* unifrac_method, bool variance_adjust, double alpha,
                              bool bypass_tips, unsigned int threads, mat_t** result);
 
+
+/* Compute UniFrac - condensed form, in-memory
+ *
+ * table <biom> a constructed BIOM object
+ * tree <BPTree> a constructed BPTree object
+ * unifrac_method <const char*> the requested unifrac method.
+ * variance_adjust <bool> whether to apply variance adjustment.
+ * alpha <double> GUniFrac alpha, only relevant if method == generalized.
+ * bypass_tips <bool> disregard tips, reduces compute by about 50%
+ * threads <uint> the number of threads to use.
+ * result <mat_t**> the resulting distance matrix in condensed form, this is initialized within the method so using **
+ *
+ * one_off_inmem returns the following error codes:
+ *
+ * okay           : no problems encountered
+ * unknown_method : the requested method is unknown.
+ * table_empty    : the table does not have any entries
+ */
+#ifdef __cplusplus
+ComputeStatus one_off_inmem(su::biom *table, su::BPTree *tree,
+                            const char* unifrac_method, bool variance_adjust, double alpha,
+                            bool bypass_tips, unsigned int threads, mat_t** result);
+#endif
+
 /* Compute UniFrac - matrix form
  *
  * biom_filename <const char*> the filename to the biom table.

src/biom.cpp

@@ -27,6 +27,8 @@ const std::string SAMPLE_DATA = std::string("/sample/matrix/data");
 const std::string SAMPLE_IDS = std::string("/sample/ids");
 
 biom::biom(std::string filename) {
+    has_hdf5_backing = true;
+
     file = H5File(filename.c_str(), H5F_ACC_RDONLY);
 
     /* establish the datasets */
@@ -58,6 +60,39 @@ biom::biom(std::string filename) {
     create_id_index(obs_ids, obs_id_index);
     create_id_index(sample_ids, sample_id_index);
 
+    malloc_resident(n_obs);
+
+    uint32_t *current_indices = NULL;
+    double *current_data = NULL;
+    for(unsigned int i = 0; i < obs_ids.size(); i++)  {
+        std::string id_ = obs_ids[i];
+        unsigned int n = get_obs_data_direct(id_, current_indices, current_data);
+        obs_counts_resident[i] = n;
+        obs_indices_resident[i] = current_indices;
+        obs_data_resident[i] = current_data;
+    }
+    sample_counts = get_sample_counts();
+}
+
+biom::~biom() {
+    if(obs_indices_resident != NULL && obs_data_resident != NULL) {
+        for(unsigned int i = 0; i < n_obs; i++) {
+            if(obs_indices_resident[i] != NULL)
+                free(obs_indices_resident[i]);
+            if(obs_data_resident[i] != NULL)
+                free(obs_data_resident[i]);
+        }
+    }
+
+    if(obs_indices_resident != NULL)
+        free(obs_indices_resident);
+    if(obs_data_resident != NULL)
+        free(obs_data_resident);
+    if(obs_counts_resident != NULL)
+        free(obs_counts_resident);
+}
+
+void biom::malloc_resident(uint32_t n_obs) { 
     /* load obs sparse data */
     obs_indices_resident = (uint32_t**)malloc(sizeof(uint32_t**) * n_obs);
     if(obs_indices_resident == NULL) {
@@ -77,30 +112,55 @@ biom::biom(std::string filename) {
                 sizeof(unsigned int) * n_obs, __FILE__, __LINE__);
         exit(EXIT_FAILURE);
     }
+}
+
+biom::biom() { 
+    n_obs = 0;
+    malloc_resident(0);
+}
+
+biom::biom(const std::vector<std::string> &obs_ids_in,
+           const std::vector<std::string> &samp_ids_in,
+           const std::vector<uint32_t> &index,
+           const std::vector<uint32_t> &indptr,
+           const std::vector<double> &data) {
+    nnz = data.size();
+    n_samples = samp_ids_in.size();
+    n_obs = obs_ids_in.size();
+
+    sample_ids = std::vector<std::string>();
+    obs_ids = std::vector<std::string>();
+    sample_ids.assign(samp_ids_in.begin(), samp_ids_in.end());
+    obs_ids.assign(obs_ids_in.begin(), obs_ids_in.end());
+
+    /* define a mapping between an ID and its corresponding offset */
+    obs_id_index = std::unordered_map<std::string, uint32_t>();
+    sample_id_index = std::unordered_map<std::string, uint32_t>();
+
+    create_id_index(obs_ids, obs_id_index);
+    create_id_index(sample_ids, sample_id_index);
+
+    malloc_resident(n_obs);
 
     uint32_t *current_indices = NULL;
     double *current_data = NULL;
     for(unsigned int i = 0; i < obs_ids.size(); i++)  {
         std::string id_ = obs_ids[i];
-        unsigned int n = get_obs_data_direct(id_, current_indices, current_data);
+        unsigned int n = get_obs_data_from_sparse(id_, index, indptr, data, current_indices, current_data);
         obs_counts_resident[i] = n;
         obs_indices_resident[i] = current_indices;
         obs_data_resident[i] = current_data;
     }
     sample_counts = get_sample_counts();
 }
 
-biom::~biom() {
-    for(unsigned int i = 0; i < n_obs; i++) {
-        free(obs_indices_resident[i]);
-        free(obs_data_resident[i]);
+void biom::set_nnz() {
+    if(!has_hdf5_backing) {
+        fprintf(stderr, "Lacks HDF5 backing; [%s]:%d\n", 
+                __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
     }
-    free(obs_indices_resident);
-    free(obs_data_resident);
-    free(obs_counts_resident);
-}
 
-void biom::set_nnz() {
     // should these be cached?
     DataType dtype = obs_data.getDataType();
     DataSpace dataspace = obs_data.getSpace();
@@ -111,6 +171,12 @@ void biom::set_nnz() {
 }
 
 void biom::load_ids(const char *path, std::vector<std::string> &ids) {
+    if(!has_hdf5_backing) {
+        fprintf(stderr, "Lacks HDF5 backing; [%s]:%d\n", 
+                __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
     DataSet ds_ids = file.openDataSet(path);
     DataType dtype = ds_ids.getDataType();
     DataSpace dataspace = ds_ids.getSpace();
@@ -138,6 +204,12 @@ void biom::load_ids(const char *path, std::vector<std::string> &ids) {
 }
 
 void biom::load_indptr(const char *path, std::vector<uint32_t> &indptr) {
+    if(!has_hdf5_backing) {
+        fprintf(stderr, "Lacks HDF5 backing; [%s]:%d\n", 
+                __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
     DataSet ds = file.openDataSet(path);
     DataType dtype = ds.getDataType();
     DataSpace dataspace = ds.getSpace();
@@ -159,7 +231,7 @@ void biom::load_indptr(const char *path, std::vector<uint32_t> &indptr) {
     free(dataout);
 }
 
-void biom::create_id_index(std::vector<std::string> &ids, 
+void biom::create_id_index(const std::vector<std::string> &ids, 
                            std::unordered_map<std::string, uint32_t> &map) {
     uint32_t count = 0;
     map.reserve(ids.size());
@@ -168,7 +240,46 @@ void biom::create_id_index(std::vector<std::string> &ids,
     }
 }
 
+unsigned int biom::get_obs_data_from_sparse(const std::string &id_, 
+                                            const std::vector<uint32_t> &index, 
+                                            const std::vector<uint32_t> &indptr, 
+                                            const std::vector<double> &data, 
+                                            uint32_t *& current_indices_out, 
+                                            double *& current_data_out) {
+    uint32_t idx = obs_id_index.at(id_);
+    uint32_t start = indptr[idx];
+    uint32_t end = indptr[idx + 1];
+    unsigned int count = end - start;
+
+    current_indices_out = (uint32_t*)malloc(sizeof(uint32_t) * count);
+    if(current_indices_out == NULL) {
+        fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", 
+                sizeof(uint32_t) * count, __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
+    current_data_out = (double*)malloc(sizeof(double) * count);
+    if(current_data_out == NULL) {
+        fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", 
+                sizeof(double) * count, __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
+    for(int i = 0, offset = start ; offset < end; i++, offset++) {
+        current_indices_out[i] = index[offset];
+        current_data_out[i] = data[offset];
+    }
+
+    return count;
+}
+
 unsigned int biom::get_obs_data_direct(const std::string &id, uint32_t *& current_indices_out, double *& current_data_out) {
+    if(!has_hdf5_backing) {
+        fprintf(stderr, "Lacks HDF5 backing; [%s]:%d\n", 
+                __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
     uint32_t idx = obs_id_index.at(id);
     uint32_t start = obs_indptr[idx];
     uint32_t end = obs_indptr[idx + 1];
@@ -270,6 +381,12 @@ void biom::get_obs_data_range(const std::string &id, unsigned int start, unsigne
 }
 
 unsigned int biom::get_sample_data_direct(const std::string &id, uint32_t *& current_indices_out, double *& current_data_out) {
+    if(!has_hdf5_backing) {
+        fprintf(stderr, "Lacks HDF5 backing; [%s]:%d\n", 
+                __FILE__, __LINE__);
+        exit(EXIT_FAILURE);
+    }
+
     uint32_t idx = sample_id_index.at(id);
     uint32_t start = sample_indptr[idx];
     uint32_t end = sample_indptr[idx + 1];