From 79f8ab57399827eaf84eaaeff0562bf690fae82f Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Tue, 15 Aug 2017 16:24:40 -0700 Subject: [PATCH 01/23] Test for NULL on malloc... --- sucpp/biom.cpp | 46 ++++++++++++++++++++++++++++++++++++++++++++++ sucpp/su.cpp | 12 ++++++++++++ sucpp/unifrac.cpp | 37 +++++++++++++++++++++++++++++++++++-- 3 files changed, 93 insertions(+), 2 deletions(-) diff --git a/sucpp/biom.cpp b/sucpp/biom.cpp index 8c21c539..a058fd2a 100644 --- a/sucpp/biom.cpp +++ b/sucpp/biom.cpp @@ -51,8 +51,24 @@ biom::biom(std::string filename) { /* load obs sparse data */ obs_indices_resident = (uint32_t**)malloc(sizeof(uint32_t**) * n_obs); + if(obs_indices_resident == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(uint32_t**) * n_obs, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } obs_data_resident = (double**)malloc(sizeof(double**) * n_obs); + if(obs_data_resident == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double**) * n_obs, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } obs_counts_resident = (unsigned int*)malloc(sizeof(unsigned int) * n_obs); + if(obs_counts_resident == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(unsigned int) * n_obs, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } + uint32_t *current_indices = NULL; double *current_data = NULL; for(unsigned int i = 0; i < obs_ids.size(); i++) { @@ -95,6 +111,11 @@ void biom::load_ids(const char *path, std::vector &ids) { /* the IDs are a dataset of variable length strings */ char **dataout = (char**)malloc(sizeof(char*) * dims[0]); + if(dataout == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(char*) * dims[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } ds_ids.read((void*)dataout, dtype); ids.reserve(dims[0]); @@ -116,6 +137,11 @@ void biom::load_indptr(const char *path, std::vector &indptr) { dataspace.getSimpleExtentDims(dims, NULL); uint32_t *dataout = (uint32_t*)malloc(sizeof(uint32_t) * dims[0]); + if(dataout == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(uint32_t) * dims[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } ds.read((void*)dataout, dtype); indptr.reserve(dims[0]); @@ -154,7 +180,17 @@ unsigned int biom::get_obs_data_direct(std::string id, uint32_t *& current_indic data_dataspace.selectHyperslab(H5S_SELECT_SET, count, offset); current_indices_out = (uint32_t*)malloc(sizeof(uint32_t) * count[0]); + if(current_indices_out == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(uint32_t) * count[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } current_data_out = (double*)malloc(sizeof(double) * count[0]); + if(current_data_out == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * count[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } obs_indices.read((void*)current_indices_out, indices_dtype, indices_memspace, indices_dataspace); obs_data.read((void*)current_data_out, data_dtype, data_memspace, data_dataspace); @@ -198,7 +234,17 @@ unsigned int biom::get_sample_data_direct(std::string id, uint32_t *& current_in data_dataspace.selectHyperslab(H5S_SELECT_SET, count, offset); current_indices_out = (uint32_t*)malloc(sizeof(uint32_t) * count[0]); + if(current_indices_out == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(uint32_t) * count[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } current_data_out = (double*)malloc(sizeof(double) * count[0]); + if(current_data_out == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * count[0], __FILE__, __LINE__); + exit(EXIT_FAILURE); + } sample_indices.read((void*)current_indices_out, indices_dtype, indices_memspace, indices_dataspace); sample_data.read((void*)current_data_out, data_dtype, data_memspace, data_dataspace); diff --git a/sucpp/su.cpp b/sucpp/su.cpp index 06f90b63..d45ec16c 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -108,8 +108,20 @@ int main(int argc, char **argv){ unsigned int chunksize = dm_stripes.size() / nthreads; unsigned int start = 0; unsigned int end = dm_stripes.size(); + unsigned int *starts = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); + if(starts == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(unsigned int) * nthreads, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } unsigned int *ends = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); + if(ends == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(unsigned int) * nthreads, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } + std::vector threads(nthreads); for(unsigned int tid = 0; tid < threads.size(); tid++) { starts[tid] = start; diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 6eee14c2..c6168116 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -54,6 +54,11 @@ double* PropStack::pop(uint32_t node) { double *vec; if(prop_stack.empty()) { posix_memalign((void **)&vec, 32, sizeof(double) * defaultsize); + if(vec == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * defaultsize, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } } else { vec = prop_stack.top(); @@ -68,8 +73,18 @@ double** su::deconvolute_stripes(std::vector &stripes, uint32_t n) { // would be better to just do striped_to_condensed_form double **dm; dm = (double**)malloc(sizeof(double*) * n); + if(dm == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double*) * n, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } for(unsigned int i = 0; i < n; i++) { dm[i] = (double*)malloc(sizeof(double) * n); + if(dm[i] == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * n, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } dm[i][i] = 0; } @@ -317,15 +332,30 @@ void su::unifrac(biom &table, double *embedded_proportions; double length; posix_memalign((void **)&embedded_proportions, 32, sizeof(double) * table.n_samples * 2); + if(embedded_proportions == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * table.n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } // thread local memory for(int i = start; i < end; i++){ posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * table.n_samples); + if(dm_stripes[i] == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * table.n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } for(int j = 0; j < table.n_samples; j++) dm_stripes[i][j] = 0.; if(unifrac_method == unweighted || unifrac_method == weighted_normalized) { posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * table.n_samples); + if(dm_stripes_total[i] == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * table.n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } for(int j = 0; j < table.n_samples; j++) dm_stripes_total[i][j] = 0.; } @@ -336,8 +366,6 @@ void su::unifrac(biom &table, node = tree.postorderselect(k); length = tree.lengths[node]; - // optional: embed this block into a prefetch thread. very easy. double buffer - // already established for embedded_proportions node_proportions = propstack.pop(node); set_proportions(node_proportions, tree, node, table, propstack); embed_proportions(embedded_proportions, node_proportions, table.n_samples); @@ -442,6 +470,11 @@ std::vector su::make_strides(unsigned int n_samples) { for(unsigned int i = 0; i < n_rotations; i++) { double* tmp; posix_memalign((void **)&tmp, 32, sizeof(double) * n_samples); + if(tmp == NULL) { + fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + sizeof(double) * n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } for(int j = 0; j < n_samples; j++) tmp[j] = 0.0; dm_stripes[i] = tmp; From 96db81ccfe2bad55d7c2023bd06c416cc26780a2 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Tue, 15 Aug 2017 16:52:59 -0700 Subject: [PATCH 02/23] Test for file existence --- sucpp/su.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/sucpp/su.cpp b/sucpp/su.cpp index d45ec16c..56cdcc85 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -24,6 +24,12 @@ void err(std::string msg) { usage(); } +// https://stackoverflow.com/a/19841704/19741 +bool is_file_exists(const char *fileName) { + std::ifstream infile(fileName); + return infile.good(); +} + int main(int argc, char **argv){ InputParser input(argc, argv); if(input.cmdOptionExists("-h")){ @@ -61,20 +67,28 @@ int main(int argc, char **argv){ err("table filename missing"); return EXIT_FAILURE; } + if(!is_file_exists(table_filename.c_str())) { + return EXIT_FAILURE; + } if(tree_filename.empty()) { err("tree filename missing"); return EXIT_FAILURE; } + if(!is_file_exists(tree_filename.c_str())) { + return EXIT_FAILURE; + } if(output_filename.empty()) { err("output filename missing"); return EXIT_FAILURE; } + if(method_string.empty()) { err("method missing"); return EXIT_FAILURE; } + if(nthreads_arg.empty()) { nthreads = 1; } else { From b82cb315a00df9da5f94e8ae9f890bd4a678ee58 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Tue, 15 Aug 2017 17:19:59 -0700 Subject: [PATCH 03/23] Cleaning up warning messages --- sucpp/biom.cpp | 20 ++++++++++---------- sucpp/su.cpp | 4 ++-- sucpp/tree.cpp | 2 +- sucpp/unifrac.cpp | 42 +++++++++++++++++++++++------------------- 4 files changed, 36 insertions(+), 32 deletions(-) diff --git a/sucpp/biom.cpp b/sucpp/biom.cpp index a058fd2a..f62c4aff 100644 --- a/sucpp/biom.cpp +++ b/sucpp/biom.cpp @@ -52,19 +52,19 @@ biom::biom(std::string filename) { /* load obs sparse data */ obs_indices_resident = (uint32_t**)malloc(sizeof(uint32_t**) * n_obs); if(obs_indices_resident == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(uint32_t**) * n_obs, __FILE__, __LINE__); exit(EXIT_FAILURE); } obs_data_resident = (double**)malloc(sizeof(double**) * n_obs); if(obs_data_resident == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double**) * n_obs, __FILE__, __LINE__); exit(EXIT_FAILURE); } obs_counts_resident = (unsigned int*)malloc(sizeof(unsigned int) * n_obs); if(obs_counts_resident == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(unsigned int) * n_obs, __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -82,7 +82,7 @@ biom::biom(std::string filename) { } biom::~biom() { - for(int i = 0; i < n_obs; i++) { + for(unsigned int i = 0; i < n_obs; i++) { free(obs_indices_resident[i]); free(obs_data_resident[i]); } @@ -112,7 +112,7 @@ void biom::load_ids(const char *path, std::vector &ids) { /* the IDs are a dataset of variable length strings */ char **dataout = (char**)malloc(sizeof(char*) * dims[0]); if(dataout == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(char*) * dims[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -138,7 +138,7 @@ void biom::load_indptr(const char *path, std::vector &indptr) { uint32_t *dataout = (uint32_t*)malloc(sizeof(uint32_t) * dims[0]); if(dataout == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(uint32_t) * dims[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -181,13 +181,13 @@ unsigned int biom::get_obs_data_direct(std::string id, uint32_t *& current_indic current_indices_out = (uint32_t*)malloc(sizeof(uint32_t) * count[0]); if(current_indices_out == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(uint32_t) * count[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } current_data_out = (double*)malloc(sizeof(double) * count[0]); if(current_data_out == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * count[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -235,13 +235,13 @@ unsigned int biom::get_sample_data_direct(std::string id, uint32_t *& current_in current_indices_out = (uint32_t*)malloc(sizeof(uint32_t) * count[0]); if(current_indices_out == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(uint32_t) * count[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } current_data_out = (double*)malloc(sizeof(double) * count[0]); if(current_data_out == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * count[0], __FILE__, __LINE__); exit(EXIT_FAILURE); } diff --git a/sucpp/su.cpp b/sucpp/su.cpp index 56cdcc85..75fdf625 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -125,13 +125,13 @@ int main(int argc, char **argv){ unsigned int *starts = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); if(starts == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(unsigned int) * nthreads, __FILE__, __LINE__); exit(EXIT_FAILURE); } unsigned int *ends = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); if(ends == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(unsigned int) * nthreads, __FILE__, __LINE__); exit(EXIT_FAILURE); } diff --git a/sucpp/tree.cpp b/sucpp/tree.cpp index f91b0540..a2ed967e 100644 --- a/sucpp/tree.cpp +++ b/sucpp/tree.cpp @@ -244,7 +244,7 @@ int32_t BPTree::enclose(uint32_t i) { } int32_t BPTree::bwd(uint32_t i, int d) { - int32_t target_excess = excess[i] + d; + uint32_t target_excess = excess[i] + d; for(int current_idx = i - 1; current_idx >= 0; current_idx--) { if(excess[current_idx] == target_excess) return current_idx; diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index c6168116..41eee88e 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -55,7 +55,7 @@ double* PropStack::pop(uint32_t node) { if(prop_stack.empty()) { posix_memalign((void **)&vec, 32, sizeof(double) * defaultsize); if(vec == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * defaultsize, __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -74,14 +74,14 @@ double** su::deconvolute_stripes(std::vector &stripes, uint32_t n) { double **dm; dm = (double**)malloc(sizeof(double*) * n); if(dm == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double*) * n, __FILE__, __LINE__); exit(EXIT_FAILURE); } for(unsigned int i = 0; i < n; i++) { dm[i] = (double*)malloc(sizeof(double) * n); if(dm[i] == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * n, __FILE__, __LINE__); exit(EXIT_FAILURE); } @@ -113,9 +113,8 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s unsigned int start, unsigned int stop) { double *dm_stripe; - int j; //__m256d ymm0, ymm1; - for(int stripe=start; stripe < stop; stripe++) { + for(unsigned int stripe=start; stripe < stop; stripe++) { dm_stripe = dm_stripes[stripe]; /* intrinsics yield about a 2x reduction in runtime on llvm. they @@ -157,16 +156,18 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s // double v = embedded_proportions[k + stripe + 1]; // dm_stripe[k] += fabs(u - v) * length; //} - for(j = 0; j < n_samples / 4; j++) { + for(unsigned int j = 0; j < n_samples / 4; j++) { int k = j * 4; double u1 = embedded_proportions[k]; double u2 = embedded_proportions[k + 1]; double u3 = embedded_proportions[k + 2]; double u4 = embedded_proportions[k + 3]; + double v1 = embedded_proportions[k + stripe + 1]; double v2 = embedded_proportions[k + stripe + 2]; double v3 = embedded_proportions[k + stripe + 3]; double v4 = embedded_proportions[k + stripe + 4]; + dm_stripe[k] += fabs(u1 - v1) * length; dm_stripe[k + 1] += fabs(u2 - v2) * length; dm_stripe[k + 2] += fabs(u3 - v3) * length; @@ -174,9 +175,10 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s } if((n_samples % 4) != 0) { - for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { + for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { double u = embedded_proportions[k]; double v = embedded_proportions[k + stripe + 1]; + dm_stripe[k] += fabs(u - v) * length; } } @@ -222,7 +224,7 @@ void _normalized_weighted_unifrac_task(std::vector &__restrict__ dm_str } if((n_samples % 4) != 0) { - for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { + for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { double u = embedded_proportions[k]; double v = embedded_proportions[k + stripe + 1]; @@ -243,16 +245,17 @@ void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, double *dm_stripe; double *dm_stripe_total; - for(int stripe = start; stripe < stop; stripe++) { + for(unsigned int stripe = start; stripe < stop; stripe++) { dm_stripe = dm_stripes[stripe]; dm_stripe_total = dm_stripes_total[stripe]; - for(int j = 0; j < n_samples / 4; j++) { + for(unsigned int j = 0; j < n_samples / 4; j++) { int k = j * 4; int32_t u1 = embedded_proportions[k] > 0; int32_t u2 = embedded_proportions[k + 1] > 0; int32_t u3 = embedded_proportions[k + 2] > 0; int32_t u4 = embedded_proportions[k + 3] > 0; + int32_t v1 = embedded_proportions[k + stripe + 1] > 0; int32_t v2 = embedded_proportions[k + stripe + 2] > 0; int32_t v3 = embedded_proportions[k + stripe + 3] > 0; @@ -262,6 +265,7 @@ void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, dm_stripe[k + 1] += (u2 ^ v2) * length; dm_stripe[k + 2] += (u3 ^ v3) * length; dm_stripe[k + 3] += (u4 ^ v4) * length; + dm_stripe_total[k] += (u1 | v1) * length; dm_stripe_total[k + 1] += (u2 | v2) * length; dm_stripe_total[k + 2] += (u3 | v3) * length; @@ -269,7 +273,7 @@ void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, } if((n_samples % 4) != 0) { - for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { + for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { int32_t u = embedded_proportions[k] > 0; int32_t v = embedded_proportions[k + stripe + 1] > 0; @@ -333,30 +337,30 @@ void su::unifrac(biom &table, double length; posix_memalign((void **)&embedded_proportions, 32, sizeof(double) * table.n_samples * 2); if(embedded_proportions == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * table.n_samples, __FILE__, __LINE__); exit(EXIT_FAILURE); } // thread local memory - for(int i = start; i < end; i++){ + for(unsigned int i = start; i < end; i++){ posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * table.n_samples); if(dm_stripes[i] == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * table.n_samples, __FILE__, __LINE__); exit(EXIT_FAILURE); } - for(int j = 0; j < table.n_samples; j++) + for(unsigned int j = 0; j < table.n_samples; j++) dm_stripes[i][j] = 0.; if(unifrac_method == unweighted || unifrac_method == weighted_normalized) { posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * table.n_samples); if(dm_stripes_total[i] == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * table.n_samples, __FILE__, __LINE__); exit(EXIT_FAILURE); } - for(int j = 0; j < table.n_samples; j++) + for(unsigned int j = 0; j < table.n_samples; j++) dm_stripes_total[i][j] = 0.; } } @@ -471,11 +475,11 @@ std::vector su::make_strides(unsigned int n_samples) { double* tmp; posix_memalign((void **)&tmp, 32, sizeof(double) * n_samples); if(tmp == NULL) { - fprintf(stderr, "Failed to allocate %d bytes; [%s]:%d\n", + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", sizeof(double) * n_samples, __FILE__, __LINE__); exit(EXIT_FAILURE); } - for(int j = 0; j < n_samples; j++) + for(unsigned int j = 0; j < n_samples; j++) tmp[j] = 0.0; dm_stripes[i] = tmp; } From fe30abf1c1ddf9fb31fbb39bccba4f2bcf3f03db Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 10:30:34 -0700 Subject: [PATCH 04/23] ENH: generalized unifrac --- sucpp/su.cpp | 36 +++++++++++++---- sucpp/test_su.cpp | 49 +++++++++++++++++++++-- sucpp/unifrac.cpp | 100 ++++++++++++++++++++++++++++------------------ sucpp/unifrac.hpp | 30 ++++++++------ 4 files changed, 152 insertions(+), 63 deletions(-) diff --git a/sucpp/su.cpp b/sucpp/su.cpp index 75fdf625..0d0d494f 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -9,13 +9,14 @@ #include void usage() { - std::cout << "usage: ssu -i -o -m [METHOD] -t [-n threads]" << std::endl; + std::cout << "usage: ssu -i -o -m [METHOD] -t [-n threads] [-a alpha]" << std::endl; std::cout << std::endl; std::cout << " -i\tThe input BIOM table." << std::endl; std::cout << " -t\tThe input phylogeny in newick." << std::endl; - std::cout << " -m\tThe method, [unweighted | weighted_normalized | weighted_unnormalized]." << std::endl; + std::cout << " -m\tThe method, [unweighted | weighted_normalized | weighted_unnormalized | generalized]." << std::endl; std::cout << " -o\tThe output distance matrix." << std::endl; std::cout << " -n\t[OPTIONAL] The number of threads, default is 1." << std::endl; + std::cout << " -a\t[OPTIONAL] Generalized UniFrac alpha, default is 1." << std::endl; std::cout << std::endl; } @@ -60,6 +61,7 @@ int main(int argc, char **argv){ const std::string &output_filename = input.getCmdOption("-o"); const std::string &method_string = input.getCmdOption("-m"); const std::string &nthreads_arg = input.getCmdOption("-n"); + const std::string &gunifrac_arg = input.getCmdOption("-a"); su::Method method; @@ -95,12 +97,21 @@ int main(int argc, char **argv){ nthreads = atoi(nthreads_arg.c_str()); } + double g_unifrac_alpha; + if(gunifrac_arg.empty()) { + g_unifrac_alpha = 1.0; + } else { + g_unifrac_alpha = atof(gunifrac_arg.c_str()); + } + if(method_string == "unweighted") method = su::unweighted; else if(method_string == "weighted_normalized") method = su::weighted_normalized; else if(method_string == "weighted_unnormalized") method = su::weighted_unnormalized; + else if(method_string == "generalized") + method = su::generalized; else { err("Unknown method"); return EXIT_FAILURE; @@ -114,8 +125,8 @@ int main(int argc, char **argv){ std::unordered_set to_keep(table.obs_ids.begin(), table.obs_ids.end()); su::BPTree tree_sheared = tree.shear(to_keep).collapse(); - std::vector dm_stripes; // = su::make_strides(table.n_samples); - std::vector dm_stripes_total; // = su::make_strides(table.n_samples); + std::vector dm_stripes; + std::vector dm_stripes_total; dm_stripes.resize((table.n_samples + 1) / 2); dm_stripes_total.resize((table.n_samples + 1) / 2); @@ -135,15 +146,27 @@ int main(int argc, char **argv){ sizeof(unsigned int) * nthreads, __FILE__, __LINE__); exit(EXIT_FAILURE); } + su::task_parameters *tasks = (su::task_parameters*)malloc(sizeof(su::task_parameters) * nthreads); + if(tasks == NULL) { + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", + sizeof(tasks) * nthreads, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } std::vector threads(nthreads); for(unsigned int tid = 0; tid < threads.size(); tid++) { starts[tid] = start; ends[tid] = start + chunksize; + tasks[tid].tid = tid; + tasks[tid].start = start; + tasks[tid].stop = end; + tasks[tid].n_samples = table.n_samples; + tasks[tid].g_unifrac_alpha = g_unifrac_alpha; start = start + chunksize; } // the last thread gets any trailing bits ends[threads.size() - 1] = end; + tasks[threads.size() - 1].stop = end; for(unsigned int tid = 0; tid < threads.size(); tid++) { threads[tid] = std::thread(su::unifrac, @@ -152,9 +175,7 @@ int main(int argc, char **argv){ method, std::ref(dm_stripes), std::ref(dm_stripes_total), - starts[tid], - ends[tid], - tid); + &tasks[tid]); } for(unsigned int tid = 0; tid < threads.size(); tid++) { @@ -163,6 +184,7 @@ int main(int argc, char **argv){ free(starts); free(ends); + free(tasks); double **dm = su::deconvolute_stripes(dm_stripes, table.n_samples); diff --git a/sucpp/test_su.cpp b/sucpp/test_su.cpp index 4f91b4ea..64178fe7 100644 --- a/sucpp/test_su.cpp +++ b/sucpp/test_su.cpp @@ -651,7 +651,10 @@ void test_unnormalized_weighted_unifrac() { std::vector strides = su::make_strides(6); std::vector strides_total = su::make_strides(6); - su::unifrac(table, tree, su::weighted_unnormalized, strides, strides_total, 0, 3, 0); + su::task_parameters task_p; + task_p.start = 0; task_p.stop = 3; task_p.tid = 0; task_p.n_samples = 6; + + su::unifrac(table, tree, su::weighted_unnormalized, strides, strides_total, &task_p); for(unsigned int i = 0; i < 3; i++) { for(unsigned int j = 0; j < 6; j++) { ASSERT(fabs(strides[i][j] - exp[i][j]) < 0.000001); @@ -661,6 +664,37 @@ void test_unnormalized_weighted_unifrac() { SUITE_END(); } +void test_generalized_unifrac() { + SUITE_START("test generalized unifrac"); + + std::vector threads(1); + su::BPTree tree = su::BPTree("(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1);"); + su::biom table = su::biom("test.biom"); + + std::vector w_exp; + double w_stride1[] = {0.38095238, 0.33333333, 0.73333333, 0.33333333, 0.5, 0.26785714}; + double w_stride2[] = {0.58095238, 0.66666667, 0.86666667, 0.25, 0.28571429, 0.45833333}; + double w_stride3[] = {0.47619048, 0.66666667, 0.46666667, 0.47619048, 0.66666667, 0.46666667}; + w_exp.push_back(w_stride1); + w_exp.push_back(w_stride2); + w_exp.push_back(w_stride3); + std::vector w_strides = su::make_strides(6); + std::vector w_strides_total = su::make_strides(6); + su::task_parameters w_task_p; + w_task_p.start = 0; w_task_p.stop = 3; w_task_p.tid = 0; w_task_p.n_samples = 6; + w_task_p.g_unifrac_alpha = 1.0; + + su::unifrac(table, tree, su::generalized, w_strides, w_strides_total, &w_task_p); + + for(unsigned int i = 0; i < 3; i++) { + for(unsigned int j = 0; j < 6; j++) { + ASSERT(fabs(w_strides[i][j] - w_exp[i][j]) < 0.000001); + } + free(w_strides[i]); + } + SUITE_END(); +} + void test_make_strides() { SUITE_START("test make stripes"); std::vector exp; @@ -695,7 +729,10 @@ void test_unweighted_unifrac() { std::vector strides = su::make_strides(6); std::vector strides_total = su::make_strides(6); - su::unifrac(table, tree, su::unweighted, strides, strides_total, 0, 3, 0); + su::task_parameters task_p; + task_p.start = 0; task_p.stop = 3; task_p.tid = 0; task_p.n_samples = 6; + + su::unifrac(table, tree, su::unweighted, strides, strides_total, &task_p); for(unsigned int i = 0; i < 3; i++) { for(unsigned int j = 0; j < 6; j++) { @@ -722,8 +759,11 @@ void test_normalized_weighted_unifrac() { exp.push_back(stride3); std::vector strides = su::make_strides(6); std::vector strides_total = su::make_strides(6); - - su::unifrac(table, tree, su::weighted_normalized, strides, strides_total, 0, 3, 0); + + su::task_parameters task_p; + task_p.start = 0; task_p.stop = 3; task_p.tid = 0; task_p.n_samples = 6; + + su::unifrac(table, tree, su::weighted_normalized, strides, strides_total, &task_p); for(unsigned int i = 0; i < 3; i++) { for(unsigned int j = 0; j < 6; j++) { ASSERT(fabs(strides[i][j] - exp[i][j]) < 0.000001); @@ -856,6 +896,7 @@ int main(int argc, char** argv) { test_unweighted_unifrac(); test_unnormalized_weighted_unifrac(); test_normalized_weighted_unifrac(); + test_generalized_unifrac(); test_unifrac_sample_counts(); printf("\n"); diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 41eee88e..fd8c1f54 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -109,12 +109,10 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, double length, - uint32_t n_samples, - unsigned int start, - unsigned int stop) { + const su::task_parameters* task_p) { double *dm_stripe; //__m256d ymm0, ymm1; - for(unsigned int stripe=start; stripe < stop; stripe++) { + for(unsigned int stripe=task_p->start; stripe < task_p->stop; stripe++) { dm_stripe = dm_stripes[stripe]; /* intrinsics yield about a 2x reduction in runtime on llvm. they @@ -156,7 +154,7 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s // double v = embedded_proportions[k + stripe + 1]; // dm_stripe[k] += fabs(u - v) * length; //} - for(unsigned int j = 0; j < n_samples / 4; j++) { + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { int k = j * 4; double u1 = embedded_proportions[k]; double u2 = embedded_proportions[k + 1]; @@ -174,8 +172,8 @@ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_s dm_stripe[k + 3] += fabs(u4 - v4) * length; } - if((n_samples % 4) != 0) { - for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { double u = embedded_proportions[k]; double v = embedded_proportions[k + stripe + 1]; @@ -189,18 +187,16 @@ void _normalized_weighted_unifrac_task(std::vector &__restrict__ dm_str std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, double length, - uint32_t n_samples, - unsigned int start, - unsigned int stop) { + const su::task_parameters* task_p) { double *dm_stripe; double *dm_stripe_total; // point of thread - for(unsigned int stripe = start; stripe < stop; stripe++) { + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { dm_stripe = dm_stripes[stripe]; dm_stripe_total = dm_stripes_total[stripe]; - for(unsigned int j = 0; j < n_samples / 4; j++) { + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { int k = j * 4; double u1 = embedded_proportions[k]; double u2 = embedded_proportions[k + 1]; @@ -217,39 +213,62 @@ void _normalized_weighted_unifrac_task(std::vector &__restrict__ dm_str dm_stripe[k + 2] += fabs(u3 - v3) * length; dm_stripe[k + 3] += fabs(u4 - v4) * length; - dm_stripe_total[k] += fabs(u1 + v1) * length; - dm_stripe_total[k + 1] += fabs(u2 + v2) * length; - dm_stripe_total[k + 2] += fabs(u3 + v3) * length; - dm_stripe_total[k + 3] += fabs(u4 + v4) * length; + dm_stripe_total[k] += (u1 + v1) * length; + dm_stripe_total[k + 1] += (u2 + v2) * length; + dm_stripe_total[k + 2] += (u3 + v3) * length; + dm_stripe_total[k + 3] += (u4 + v4) * length; } - if((n_samples % 4) != 0) { - for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { double u = embedded_proportions[k]; double v = embedded_proportions[k + stripe + 1]; dm_stripe[k] += fabs(u - v) * length; - dm_stripe_total[k] += fabs(u + v) * length; + dm_stripe_total[k] += (u + v) * length; } } } } +void _generalized_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + // point of thread + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + double u1 = embedded_proportions[j]; + double v1 = embedded_proportions[j + stripe + 1]; + double sub1 = fabs(u1 - v1); + double sum1 = u1 + v1; + double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; + dm_stripe[j] += (sum1 != 0.0) ? (sum_pow1 * (sub1 / sum1)) : 0.0; + dm_stripe_total[j] += sum_pow1; + } + } +} + void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, double length, - uint32_t n_samples, - unsigned int start, - unsigned int stop) { + const su::task_parameters* task_p) { double *dm_stripe; double *dm_stripe_total; - for(unsigned int stripe = start; stripe < stop; stripe++) { + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { dm_stripe = dm_stripes[stripe]; dm_stripe_total = dm_stripes_total[stripe]; - for(unsigned int j = 0; j < n_samples / 4; j++) { + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { int k = j * 4; int32_t u1 = embedded_proportions[k] > 0; int32_t u2 = embedded_proportions[k + 1] > 0; @@ -272,8 +291,8 @@ void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, dm_stripe_total[k + 3] += (u4 | v4) * length; } - if((n_samples % 4) != 0) { - for(unsigned int k = n_samples - (n_samples % 4); k < n_samples; k++) { + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { int32_t u = embedded_proportions[k] > 0; int32_t v = embedded_proportions[k + stripe + 1] > 0; @@ -306,17 +325,18 @@ void su::unifrac(biom &table, Method unifrac_method, std::vector &dm_stripes, std::vector &dm_stripes_total, - unsigned int start, - unsigned int end, - unsigned int tid) { + const su::task_parameters* task_p) { + + if(table.n_samples != task_p->n_samples) { + fprintf(stderr, "Task and table n_samples not equal\n"); + exit(EXIT_FAILURE); + } void (*func)(std::vector&, // dm_stripes std::vector&, // dm_stripes_total double*, // embedded_proportions double, // length - uint32_t, // number of samples - unsigned int, // stripe start - unsigned int); // stripe stop + const su::task_parameters*); switch(unifrac_method) { case unweighted: @@ -328,6 +348,9 @@ void su::unifrac(biom &table, case weighted_unnormalized: func = &_unnormalized_weighted_unifrac_task; break; + case generalized: + func = &_generalized_unifrac_task; + break; } PropStack propstack(table.n_samples); @@ -343,7 +366,7 @@ void su::unifrac(biom &table, } // thread local memory - for(unsigned int i = start; i < end; i++){ + for(unsigned int i = task_p->start; i < task_p->stop; i++){ posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * table.n_samples); if(dm_stripes[i] == NULL) { fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", @@ -372,7 +395,7 @@ void su::unifrac(biom &table, node_proportions = propstack.pop(node); set_proportions(node_proportions, tree, node, table, propstack); - embed_proportions(embedded_proportions, node_proportions, table.n_samples); + embed_proportions(embedded_proportions, node_proportions, task_p->n_samples); /* * The values in the example vectors correspond to index positions of an @@ -418,17 +441,16 @@ void su::unifrac(biom &table, * We end up performing N / 2 redundant calculations on the last stripe * (see C) but that is small over large N. */ - func(dm_stripes, dm_stripes_total, embedded_proportions, length, - table.n_samples, start, end); + func(dm_stripes, dm_stripes_total, embedded_proportions, length, task_p); // should make this compile-time support //if((tid == 0) && ((k % 1000) == 0)) // progressbar((float)k / (float)(tree.nparens / 2)); } - if(unifrac_method == weighted_normalized || unifrac_method == unweighted) { - for(unsigned int i = start; i < end; i++) { - for(unsigned int j = 0; j < table.n_samples; j++) { + if(unifrac_method == weighted_normalized || unifrac_method == unweighted || unifrac_method == generalized) { + for(unsigned int i = task_p->start; i < task_p->stop; i++) { + for(unsigned int j = 0; j < task_p->n_samples; j++) { dm_stripes[i][j] = dm_stripes[i][j] / dm_stripes_total[i][j]; } } diff --git a/sucpp/unifrac.hpp b/sucpp/unifrac.hpp index ab4335d8..f4a8159c 100644 --- a/sucpp/unifrac.hpp +++ b/sucpp/unifrac.hpp @@ -4,7 +4,15 @@ #include namespace su { - enum Method {unweighted, weighted_normalized, weighted_unnormalized}; + enum Method {unweighted, weighted_normalized, weighted_unnormalized, generalized}; + + struct task_parameters { + uint32_t n_samples; // number of samples + unsigned int start; // starting stripe + unsigned int stop; // stopping stripe + unsigned int tid; // thread ID + double g_unifrac_alpha; // generalized unifrac alpha + }; class PropStack { private: @@ -24,17 +32,13 @@ namespace su { Method unifrac_method, std::vector &dm_stripes, std::vector &dm_stripes_total, - unsigned int start, - unsigned int end, - unsigned int tid); - void unweighted_unifrac(biom &table, - BPTree &tree, - Method unifrac_method, - std::vector &dm_stripes, - std::vector &dm_stripes_total, - unsigned int start, - unsigned int end, - unsigned int tid); + const task_parameters* task_p); + //void unweighted_unifrac(biom &table, + // BPTree &tree, + // Method unifrac_method, + // std::vector &dm_stripes, + // std::vector &dm_stripes_total, + // const task_parameters* task_p); double** deconvolute_stripes(std::vector &stripes, uint32_t n); void set_proportions(double* props, BPTree &tree, uint32_t node, @@ -49,5 +53,5 @@ namespace su { out[i + n] = val; } } -} +} From a99808cd47698b1cd015351dd34672fe2ed2d3bf Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 10:55:57 -0700 Subject: [PATCH 05/23] Expose via q2 --- q2_state_unifrac/_methods.py | 13 +++++++++++ q2_state_unifrac/plugin_setup.py | 40 ++++++++++++++++++++++++++++---- 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/q2_state_unifrac/_methods.py b/q2_state_unifrac/_methods.py index f9dab867..dbe015b9 100644 --- a/q2_state_unifrac/_methods.py +++ b/q2_state_unifrac/_methods.py @@ -68,3 +68,16 @@ def weighted_unnormalized(table: BIOMV210Format, _run(str(table), str(phylogeny), output_fp, str(threads), 'weighted_unnormalized') return skbio.DistanceMatrix.read(output_fp) + + +def generalized(table: BIOMV210Format, + phylogeny: NewickFormat, + threads: int=1, + alpha: float=1.0)-> skbio.DistanceMatrix: + _sanity() + + with tempfile.TemporaryDirectory() as tmp: + output_fp = os.path.join(tmp, 'foo.dm') + _run(str(table), str(phylogeny), output_fp, str(threads), + 'generalized') + return skbio.DistanceMatrix.read(output_fp) diff --git a/q2_state_unifrac/plugin_setup.py b/q2_state_unifrac/plugin_setup.py index a60c1c33..511134d3 100644 --- a/q2_state_unifrac/plugin_setup.py +++ b/q2_state_unifrac/plugin_setup.py @@ -6,7 +6,7 @@ # The full license is in the file LICENSE, distributed with this software. # ---------------------------------------------------------------------------- -from qiime2.plugin import (Plugin, Properties, Int) +from qiime2.plugin import (Plugin, Properties, Int, Float) from q2_types.feature_table import FeatureTable, Frequency from q2_types.distance_matrix import DistanceMatrix from q2_types.tree import Phylogeny, Rooted @@ -37,7 +37,9 @@ }, outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], name='Unweighted UniFrac', - description=('This method computes Unweighted UniFrac') + description=('This method computes unweighted UniFrac as described in ' + 'Lozupone and Knight 2005 Appl Environ Microbiol; ' + 'DOI: 10.1128/AEM.71.12.8228-8235.2005') ) @@ -54,7 +56,9 @@ }, outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], name='Weighted unnormalizd UniFrac', - description=('This method computes weighted unnormalized UniFrac') + description=('This method computes weighted unnormalized UniFrac as ' + 'described in Lozupone et al. 2007 Appl Environ Microbiol; ' + 'DOI: 10.1128/AEM.01996-06') ) @@ -71,5 +75,33 @@ }, outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], name='Wweighted normalized UniFrac', - description=('This method computes weighted normalized UniFrac') + description=('This method computes weighted normalized UniFrac as ' + 'described in Lozupone et al. 2007 Appl Environ Microbiol; ' + 'DOI: 10.1128/AEM.01996-06') +) + + +plugin.methods.register_function( + function=q2_state_unifrac.generalized, + inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), + 'phylogeny': Phylogeny[Rooted]}, + parameters={'threads': Int, + 'alpha': Float}, + parameter_descriptions={'threads': 'The number of threads to use.', + 'alpha': ('The value of alpha controls importance ' + 'of sample proportions. 1.0 is ' + 'weighted normalized UniFrac. 0.0 is ' + 'close to unweighted UniFrac, but only ' + 'if the sample proportions are ' + 'dichotomized.')}, + input_descriptions={ + 'table': 'A rarefied FeatureTable', + 'phylogeny': ('A rooted phylogeny which relates the observations in ' + 'the table.') + }, + outputs=[('distance_matrix', DistanceMatrix % Properties('phylogenetic'))], + name='Wweighted normalized UniFrac', + description=('This method computes Generalized UniFrac as described in ' + 'Chen et al. 2012 Bioinformatics; ' + 'DOI: 10.1093/bioinformatics/bts342') ) From 23af923941ca49e6f9fdc811d273e75f791d9013 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 11:49:11 -0700 Subject: [PATCH 06/23] TST: d^0 and d^0.5 --- sucpp/test_su.cpp | 52 +++++++++++++++++++++++++++++++++++++++++++++-- sucpp/unifrac.cpp | 10 +++++---- 2 files changed, 56 insertions(+), 6 deletions(-) diff --git a/sucpp/test_su.cpp b/sucpp/test_su.cpp index 64178fe7..ee25440f 100644 --- a/sucpp/test_su.cpp +++ b/sucpp/test_su.cpp @@ -671,6 +671,7 @@ void test_generalized_unifrac() { su::BPTree tree = su::BPTree("(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1);"); su::biom table = su::biom("test.biom"); + // weighted normalized unifrac as computed above std::vector w_exp; double w_stride1[] = {0.38095238, 0.33333333, 0.73333333, 0.33333333, 0.5, 0.26785714}; double w_stride2[] = {0.58095238, 0.66666667, 0.86666667, 0.25, 0.28571429, 0.45833333}; @@ -683,14 +684,61 @@ void test_generalized_unifrac() { su::task_parameters w_task_p; w_task_p.start = 0; w_task_p.stop = 3; w_task_p.tid = 0; w_task_p.n_samples = 6; w_task_p.g_unifrac_alpha = 1.0; - su::unifrac(table, tree, su::generalized, w_strides, w_strides_total, &w_task_p); - + + // as computed by GUniFrac + // Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 + //Sample1 0.0000000 0.4408392 0.6886965 0.7060606 0.5833333 0.3278410 + //Sample2 0.4408392 0.0000000 0.5102041 0.7500000 0.8000000 0.5208125 + //Sample3 0.6886965 0.5102041 0.0000000 0.8649351 0.9428571 0.5952381 + //Sample4 0.7060606 0.7500000 0.8649351 0.0000000 0.5000000 0.4857143 + //Sample5 0.5833333 0.8000000 0.9428571 0.5000000 0.0000000 0.7485714 + //Sample6 0.3278410 0.5208125 0.5952381 0.4857143 0.7485714 0.0000000 + std::vector d0_exp; + double d0_stride1[] = {0.4408392, 0.5102041, 0.8649351, 0.5000000, 0.7485714, 0.3278410}; + double d0_stride2[] = {0.6886965, 0.7500000, 0.9428571, 0.4857143, 0.5833333, 0.5208125}; + double d0_stride3[] = {0.7060606, 0.8000000, 0.5952381, 0.7060606, 0.8000000, 0.5952381}; + d0_exp.push_back(d0_stride1); + d0_exp.push_back(d0_stride2); + d0_exp.push_back(d0_stride3); + std::vector d0_strides = su::make_strides(6); + std::vector d0_strides_total = su::make_strides(6); + su::task_parameters d0_task_p; + d0_task_p.start = 0; d0_task_p.stop = 3; d0_task_p.tid = 0; d0_task_p.n_samples = 6; + d0_task_p.g_unifrac_alpha = 0.0; + su::unifrac(table, tree, su::generalized, d0_strides, d0_strides_total, &d0_task_p); + + // as computed by GUniFrac + // Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 + //Sample1 0.0000000 0.4040518 0.6285560 0.5869439 0.4082483 0.2995673 + //Sample2 0.4040518 0.0000000 0.4160597 0.7071068 0.7302479 0.4860856 + //Sample3 0.6285560 0.4160597 0.0000000 0.8005220 0.9073159 0.5218198 + //Sample4 0.5869439 0.7071068 0.8005220 0.0000000 0.4117216 0.3485667 + //Sample5 0.4082483 0.7302479 0.9073159 0.4117216 0.0000000 0.6188282 + //Sample6 0.2995673 0.4860856 0.5218198 0.3485667 0.6188282 0.0000000 + std::vector d05_exp; + double d05_stride1[] = {0.4040518, 0.4160597, 0.8005220, 0.4117216, 0.6188282, 0.2995673}; + double d05_stride2[] = {0.6285560, 0.7071068, 0.9073159, 0.3485667, 0.4082483, 0.4860856}; + double d05_stride3[] = {0.5869439, 0.7302479, 0.5218198, 0.5869439, 0.7302479, 0.5218198}; + d05_exp.push_back(d05_stride1); + d05_exp.push_back(d05_stride2); + d05_exp.push_back(d05_stride3); + std::vector d05_strides = su::make_strides(6); + std::vector d05_strides_total = su::make_strides(6); + su::task_parameters d05_task_p; + d05_task_p.start = 0; d05_task_p.stop = 3; d05_task_p.tid = 0; d05_task_p.n_samples = 6; + d05_task_p.g_unifrac_alpha = 0.5; + su::unifrac(table, tree, su::generalized, d05_strides, d05_strides_total, &d05_task_p); + for(unsigned int i = 0; i < 3; i++) { for(unsigned int j = 0; j < 6; j++) { ASSERT(fabs(w_strides[i][j] - w_exp[i][j]) < 0.000001); + ASSERT(fabs(d0_strides[i][j] - d0_exp[i][j]) < 0.000001); + ASSERT(fabs(d05_strides[i][j] - d05_exp[i][j]) < 0.000001); } free(w_strides[i]); + free(d0_strides[i]); + free(d05_strides[i]); } SUITE_END(); } diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index fd8c1f54..74619a3b 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -247,11 +247,13 @@ void _generalized_unifrac_task(std::vector &__restrict__ dm_stripes, for(unsigned int j = 0; j < task_p->n_samples; j++) { double u1 = embedded_proportions[j]; double v1 = embedded_proportions[j + stripe + 1]; - double sub1 = fabs(u1 - v1); double sum1 = u1 + v1; - double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; - dm_stripe[j] += (sum1 != 0.0) ? (sum_pow1 * (sub1 / sum1)) : 0.0; - dm_stripe_total[j] += sum_pow1; + if(sum1 != 0.0) { + double sub1 = fabs(u1 - v1); + double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; + dm_stripe[j] += sum_pow1 * (sub1 / sum1); + dm_stripe_total[j] += sum_pow1; + } } } } From b82b7e58b946bd5fdf683e8d73abf6d27ef9dda6 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 13:26:14 -0700 Subject: [PATCH 07/23] Missed import --- q2_state_unifrac/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/q2_state_unifrac/__init__.py b/q2_state_unifrac/__init__.py index bab503ba..8772f077 100644 --- a/q2_state_unifrac/__init__.py +++ b/q2_state_unifrac/__init__.py @@ -10,8 +10,10 @@ from q2_state_unifrac._methods import (unweighted, weighted_normalized, - weighted_unnormalized) + weighted_unnormalized, + generalized) __version__ = pkg_resources.get_distribution('q2-state-unifrac').version -__all__ = ['unweighted', 'weighted_normalized', 'weighted_unnormalized'] +__all__ = ['unweighted', 'weighted_normalized', 'weighted_unnormalized', + 'generalized'] From a98bf69bf80997fe2e6a48d2ab434b7bbe3c48d5 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 14:00:20 -0700 Subject: [PATCH 08/23] Refactoring --- sucpp/unifrac.cpp | 68 ++++++++++++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 28 deletions(-) diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 74619a3b..034e3190 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -322,6 +322,44 @@ void progressbar(float progress) { std::cout.flush(); } +void initialize_embedded(double*& prop, const su::task_parameters* task_p) { + posix_memalign((void **)&prop, 32, sizeof(double) * task_p->n_samples * 2); + if(prop == NULL) { + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", + sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } +} + + +void initialize_stripes(std::vector &dm_stripes, + std::vector &dm_stripes_total, + Method unifrac_method, + const su::task_parameters* task_p) { + for(unsigned int i = task_p->start; i < task_p->stop; i++){ + posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * task_p->n_samples); + if(dm_stripes[i] == NULL) { + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", + sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } + for(unsigned int j = 0; j < task_p->n_samples; j++) + dm_stripes[i][j] = 0.; + + if(unifrac_method == unweighted || unifrac_method == weighted_normalized) { + posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * task_p->n_samples); + if(dm_stripes_total[i] == NULL) { + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", + sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } + for(unsigned int j = 0; j < task_p->n_samples; j++) + dm_stripes_total[i][j] = 0.; + } + } +} + + void su::unifrac(biom &table, BPTree &tree, Method unifrac_method, @@ -360,35 +398,9 @@ void su::unifrac(biom &table, double *node_proportions; double *embedded_proportions; double length; - posix_memalign((void **)&embedded_proportions, 32, sizeof(double) * table.n_samples * 2); - if(embedded_proportions == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * table.n_samples, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } - // thread local memory - for(unsigned int i = task_p->start; i < task_p->stop; i++){ - posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * table.n_samples); - if(dm_stripes[i] == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * table.n_samples, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } - for(unsigned int j = 0; j < table.n_samples; j++) - dm_stripes[i][j] = 0.; - - if(unifrac_method == unweighted || unifrac_method == weighted_normalized) { - posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * table.n_samples); - if(dm_stripes_total[i] == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * table.n_samples, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } - for(unsigned int j = 0; j < table.n_samples; j++) - dm_stripes_total[i][j] = 0.; - } - } + initialize_embedded(embedded_proportions, task_p); + initialize_stripes(std::ref(dm_stripes), std::ref(dm_stripes_total), unifrac_method, task_p); // - 1 to avoid root for(unsigned int k = 0; k < (tree.nparens / 2) - 1; k++) { From 2c083f2ea7f5ac1282cf9f83d3cb06e059c0b1fe Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 14:52:22 -0700 Subject: [PATCH 09/23] regression: task param boundaries were set wrong --- .travis.yml | 2 ++ sucpp/su.cpp | 31 ++++--------------------------- 2 files changed, 6 insertions(+), 27 deletions(-) diff --git a/.travis.yml b/.travis.yml index 83919720..90d7ad7f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -43,5 +43,7 @@ script: - for i in $(find ./ -maxdepth 1 -name 'core*' -print); do gdb $(pwd)/sucpp/ssu core* -ex "thread apply all bt" -ex "set pagination 0" -batch; done; - ls -lrt - qiime state-unifrac unweighted --i-table table.qza --i-phylogeny tree.qza --o-distance-matrix ssu_result.qza --verbose + - qiime state-unifrac unweighted --i-table table.qza --i-phylogeny tree.qza --o-distance-matrix ssu_result_thread.qza --verbose --p-threads 2 - qiime diversity beta-phylogenetic --p-metric unweighted_unifrac --i-table table.qza --i-phylogeny tree.qza --o-distance-matrix sk_result.qza - python compare_dms.py ssu_result.qza sk_result.qza + - python compare_dms.py ssu_result_thread.qza sk_result.qza diff --git a/sucpp/su.cpp b/sucpp/su.cpp index 0d0d494f..f59edd8d 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -134,38 +134,19 @@ int main(int argc, char **argv){ unsigned int start = 0; unsigned int end = dm_stripes.size(); - unsigned int *starts = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); - if(starts == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(unsigned int) * nthreads, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } - unsigned int *ends = (unsigned int*)malloc(sizeof(unsigned int) * nthreads); - if(ends == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(unsigned int) * nthreads, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } - su::task_parameters *tasks = (su::task_parameters*)malloc(sizeof(su::task_parameters) * nthreads); - if(tasks == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(tasks) * nthreads, __FILE__, __LINE__); - exit(EXIT_FAILURE); - } + std::vector tasks; + tasks.resize(nthreads); std::vector threads(nthreads); for(unsigned int tid = 0; tid < threads.size(); tid++) { - starts[tid] = start; - ends[tid] = start + chunksize; tasks[tid].tid = tid; - tasks[tid].start = start; - tasks[tid].stop = end; + tasks[tid].start = start; // stripe start + tasks[tid].stop = start + chunksize; // stripe end tasks[tid].n_samples = table.n_samples; tasks[tid].g_unifrac_alpha = g_unifrac_alpha; start = start + chunksize; } // the last thread gets any trailing bits - ends[threads.size() - 1] = end; tasks[threads.size() - 1].stop = end; for(unsigned int tid = 0; tid < threads.size(); tid++) { @@ -182,10 +163,6 @@ int main(int argc, char **argv){ threads[tid].join(); } - free(starts); - free(ends); - free(tasks); - double **dm = su::deconvolute_stripes(dm_stripes, table.n_samples); std::ofstream output; From e3804bd1a891a5d382a539599786c5ab6d3ed0e2 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 15:15:03 -0700 Subject: [PATCH 10/23] Refactor tasks to separate file --- sucpp/Makefile | 8 +- sucpp/unifrac.cpp | 210 +---------------------------------------- sucpp/unifrac.hpp | 16 +--- sucpp/unifrac_task.cpp | 203 +++++++++++++++++++++++++++++++++++++++ sucpp/unifrac_task.hpp | 33 +++++++ 5 files changed, 246 insertions(+), 224 deletions(-) create mode 100644 sucpp/unifrac_task.cpp create mode 100644 sucpp/unifrac_task.hpp diff --git a/sucpp/Makefile b/sucpp/Makefile index a5be99c5..2cd5ad8d 100644 --- a/sucpp/Makefile +++ b/sucpp/Makefile @@ -29,11 +29,11 @@ else endif -test: tree.o test_su.cpp biom.o unifrac.o - $(CXX) $(CPPFLAGS) -Wno-unused-parameter test_su.cpp -o test_su tree.o biom.o unifrac.o +test: tree.o test_su.cpp biom.o unifrac.o unifrac_task.o + $(CXX) $(CPPFLAGS) -Wno-unused-parameter test_su.cpp -o test_su tree.o biom.o unifrac.o unifrac_task.o -main: tree.o biom.o unifrac.o cmd.o - $(CXX) $(CPPFLAGS) su.cpp -o ssu tree.o biom.o unifrac.o cmd.o +main: tree.o biom.o unifrac.o cmd.o unifrac_task.o + $(CXX) $(CPPFLAGS) su.cpp -o ssu tree.o biom.o unifrac.o cmd.o unifrac_task.o ocltest: opencl.o clang++ -framework OpenCL opencl.cpp -c diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 034e3190..762a3974 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -3,7 +3,6 @@ #include "unifrac.hpp" #include #include -#include #include #include @@ -105,207 +104,6 @@ double** su::deconvolute_stripes(std::vector &stripes, uint32_t n) { return dm; } -void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, - std::vector &__restrict__ dm_stripes_total, - double* __restrict__ embedded_proportions, - double length, - const su::task_parameters* task_p) { - double *dm_stripe; - //__m256d ymm0, ymm1; - for(unsigned int stripe=task_p->start; stripe < task_p->stop; stripe++) { - dm_stripe = dm_stripes[stripe]; - - /* intrinsics yield about a 2x reduction in runtime on llvm. they - * were not effective on linux gcc 4.9.1 or 4.9.2. it is unclear - * if they would be effective on other versions of gcc. - * - * one reason they help is that these for loops are not easily - * autovectorizable. using the intrinsics effectively gets around - * this. ...although, it also appears that loop unrolling works. - * - * it may make sense to revisit the inclusion of intriniscs, however - * support must be tested at compile time, so it's rather annoying - * at the moment. basically, we can't assume the presence of avx2. - */ - //for(j = 0; j < n_samples / 4; j++) { - // int offset = j * 4; - // ymm0 = _mm256_sub_pd(_mm256_load_pd(embedded_proportions + offset), - // _mm256_load_pd(embedded_proportions + (offset + stripe + 1))); // u - v - // ymm1 = _mm256_sub_pd(_mm256_set1_pd(0.0), - // ymm0); // 0.0 - (u - v) - // ymm0 = _mm256_max_pd(ymm0, ymm1); // abs(u - v) - // ymm0 = _mm256_fmadd_pd(ymm0, - // _mm256_set1_pd(length), - // _mm256_load_pd(dm_stripe + offset)); // fused mutiply add, dm_stripe[offset] += (abs(u - v) * length) - // _mm256_store_pd(dm_stripe + offset, ymm0); //store - //} - - //if((n_samples % 4) != 0) { - // for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { - // double u = embedded_proportions[k]; - // ////double v = alt_proportions[j]; - // double v = embedded_proportions[k + stripe + 1]; - // dm_stripe[k] += fabs(u - v) * length; - // } - //} - - //for(int k = 0; k < n_samples; k++) { - // double u = embedded_proportions[k]; - // double v = embedded_proportions[k + stripe + 1]; - // dm_stripe[k] += fabs(u - v) * length; - //} - for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { - int k = j * 4; - double u1 = embedded_proportions[k]; - double u2 = embedded_proportions[k + 1]; - double u3 = embedded_proportions[k + 2]; - double u4 = embedded_proportions[k + 3]; - - double v1 = embedded_proportions[k + stripe + 1]; - double v2 = embedded_proportions[k + stripe + 2]; - double v3 = embedded_proportions[k + stripe + 3]; - double v4 = embedded_proportions[k + stripe + 4]; - - dm_stripe[k] += fabs(u1 - v1) * length; - dm_stripe[k + 1] += fabs(u2 - v2) * length; - dm_stripe[k + 2] += fabs(u3 - v3) * length; - dm_stripe[k + 3] += fabs(u4 - v4) * length; - } - - if((task_p->n_samples % 4) != 0) { - for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { - double u = embedded_proportions[k]; - double v = embedded_proportions[k + stripe + 1]; - - dm_stripe[k] += fabs(u - v) * length; - } - } - } -} - -void _normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, - std::vector &__restrict__ dm_stripes_total, - double* __restrict__ embedded_proportions, - double length, - const su::task_parameters* task_p) { - double *dm_stripe; - double *dm_stripe_total; - - // point of thread - for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { - dm_stripe = dm_stripes[stripe]; - dm_stripe_total = dm_stripes_total[stripe]; - - for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { - int k = j * 4; - double u1 = embedded_proportions[k]; - double u2 = embedded_proportions[k + 1]; - double u3 = embedded_proportions[k + 2]; - double u4 = embedded_proportions[k + 3]; - - double v1 = embedded_proportions[k + stripe + 1]; - double v2 = embedded_proportions[k + stripe + 2]; - double v3 = embedded_proportions[k + stripe + 3]; - double v4 = embedded_proportions[k + stripe + 4]; - - dm_stripe[k] += fabs(u1 - v1) * length; - dm_stripe[k + 1] += fabs(u2 - v2) * length; - dm_stripe[k + 2] += fabs(u3 - v3) * length; - dm_stripe[k + 3] += fabs(u4 - v4) * length; - - dm_stripe_total[k] += (u1 + v1) * length; - dm_stripe_total[k + 1] += (u2 + v2) * length; - dm_stripe_total[k + 2] += (u3 + v3) * length; - dm_stripe_total[k + 3] += (u4 + v4) * length; - } - - if((task_p->n_samples % 4) != 0) { - for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { - double u = embedded_proportions[k]; - double v = embedded_proportions[k + stripe + 1]; - - dm_stripe[k] += fabs(u - v) * length; - dm_stripe_total[k] += (u + v) * length; - } - } - } -} - -void _generalized_unifrac_task(std::vector &__restrict__ dm_stripes, - std::vector &__restrict__ dm_stripes_total, - double* __restrict__ embedded_proportions, - double length, - const su::task_parameters* task_p) { - double *dm_stripe; - double *dm_stripe_total; - - // point of thread - for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { - dm_stripe = dm_stripes[stripe]; - dm_stripe_total = dm_stripes_total[stripe]; - - for(unsigned int j = 0; j < task_p->n_samples; j++) { - double u1 = embedded_proportions[j]; - double v1 = embedded_proportions[j + stripe + 1]; - double sum1 = u1 + v1; - if(sum1 != 0.0) { - double sub1 = fabs(u1 - v1); - double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; - dm_stripe[j] += sum_pow1 * (sub1 / sum1); - dm_stripe_total[j] += sum_pow1; - } - } - } -} - -void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, - std::vector &__restrict__ dm_stripes_total, - double* __restrict__ embedded_proportions, - double length, - const su::task_parameters* task_p) { - double *dm_stripe; - double *dm_stripe_total; - - for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { - dm_stripe = dm_stripes[stripe]; - dm_stripe_total = dm_stripes_total[stripe]; - - for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { - int k = j * 4; - int32_t u1 = embedded_proportions[k] > 0; - int32_t u2 = embedded_proportions[k + 1] > 0; - int32_t u3 = embedded_proportions[k + 2] > 0; - int32_t u4 = embedded_proportions[k + 3] > 0; - - int32_t v1 = embedded_proportions[k + stripe + 1] > 0; - int32_t v2 = embedded_proportions[k + stripe + 2] > 0; - int32_t v3 = embedded_proportions[k + stripe + 3] > 0; - int32_t v4 = embedded_proportions[k + stripe + 4] > 0; - - dm_stripe[k] += (u1 ^ v1) * length; - dm_stripe[k + 1] += (u2 ^ v2) * length; - dm_stripe[k + 2] += (u3 ^ v3) * length; - dm_stripe[k + 3] += (u4 ^ v4) * length; - - dm_stripe_total[k] += (u1 | v1) * length; - dm_stripe_total[k + 1] += (u2 | v2) * length; - dm_stripe_total[k + 2] += (u3 | v3) * length; - dm_stripe_total[k + 3] += (u4 | v4) * length; - } - - if((task_p->n_samples % 4) != 0) { - for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { - int32_t u = embedded_proportions[k] > 0; - int32_t v = embedded_proportions[k + stripe + 1] > 0; - - dm_stripe[k] += (u ^ v) * length; - dm_stripe_total[k] += (u | v) * length; - } - } - } -} - - void progressbar(float progress) { // from http://stackoverflow.com/a/14539953 // @@ -380,16 +178,16 @@ void su::unifrac(biom &table, switch(unifrac_method) { case unweighted: - func = &_unweighted_unifrac_task; + func = &su::_unweighted_unifrac_task; break; case weighted_normalized: - func = &_normalized_weighted_unifrac_task; + func = &su::_normalized_weighted_unifrac_task; break; case weighted_unnormalized: - func = &_unnormalized_weighted_unifrac_task; + func = &su::_unnormalized_weighted_unifrac_task; break; case generalized: - func = &_generalized_unifrac_task; + func = &su::_generalized_unifrac_task; break; } PropStack propstack(table.n_samples); diff --git a/sucpp/unifrac.hpp b/sucpp/unifrac.hpp index f4a8159c..56c923b0 100644 --- a/sucpp/unifrac.hpp +++ b/sucpp/unifrac.hpp @@ -2,18 +2,11 @@ #include #include #include +#include "unifrac_task.hpp" namespace su { enum Method {unweighted, weighted_normalized, weighted_unnormalized, generalized}; - struct task_parameters { - uint32_t n_samples; // number of samples - unsigned int start; // starting stripe - unsigned int stop; // stopping stripe - unsigned int tid; // thread ID - double g_unifrac_alpha; // generalized unifrac alpha - }; - class PropStack { private: std::stack prop_stack; @@ -33,12 +26,7 @@ namespace su { std::vector &dm_stripes, std::vector &dm_stripes_total, const task_parameters* task_p); - //void unweighted_unifrac(biom &table, - // BPTree &tree, - // Method unifrac_method, - // std::vector &dm_stripes, - // std::vector &dm_stripes_total, - // const task_parameters* task_p); + double** deconvolute_stripes(std::vector &stripes, uint32_t n); void set_proportions(double* props, BPTree &tree, uint32_t node, diff --git a/sucpp/unifrac_task.cpp b/sucpp/unifrac_task.cpp new file mode 100644 index 00000000..1437a0b3 --- /dev/null +++ b/sucpp/unifrac_task.cpp @@ -0,0 +1,203 @@ +#include "unifrac_task.hpp" + + +void su::_unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + //__m256d ymm0, ymm1; + for(unsigned int stripe=task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + + /* intrinsics yield about a 2x reduction in runtime on llvm. they + * were not effective on linux gcc 4.9.1 or 4.9.2. it is unclear + * if they would be effective on other versions of gcc. + * + * one reason they help is that these for loops are not easily + * autovectorizable. using the intrinsics effectively gets around + * this. ...although, it also appears that loop unrolling works. + * + * it may make sense to revisit the inclusion of intriniscs, however + * support must be tested at compile time, so it's rather annoying + * at the moment. basically, we can't assume the presence of avx2. + */ + //for(j = 0; j < n_samples / 4; j++) { + // int offset = j * 4; + // ymm0 = _mm256_sub_pd(_mm256_load_pd(embedded_proportions + offset), + // _mm256_load_pd(embedded_proportions + (offset + stripe + 1))); // u - v + // ymm1 = _mm256_sub_pd(_mm256_set1_pd(0.0), + // ymm0); // 0.0 - (u - v) + // ymm0 = _mm256_max_pd(ymm0, ymm1); // abs(u - v) + // ymm0 = _mm256_fmadd_pd(ymm0, + // _mm256_set1_pd(length), + // _mm256_load_pd(dm_stripe + offset)); // fused mutiply add, dm_stripe[offset] += (abs(u - v) * length) + // _mm256_store_pd(dm_stripe + offset, ymm0); //store + //} + + //if((n_samples % 4) != 0) { + // for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { + // double u = embedded_proportions[k]; + // ////double v = alt_proportions[j]; + // double v = embedded_proportions[k + stripe + 1]; + // dm_stripe[k] += fabs(u - v) * length; + // } + //} + + //for(int k = 0; k < n_samples; k++) { + // double u = embedded_proportions[k]; + // double v = embedded_proportions[k + stripe + 1]; + // dm_stripe[k] += fabs(u - v) * length; + //} + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { + int k = j * 4; + double u1 = embedded_proportions[k]; + double u2 = embedded_proportions[k + 1]; + double u3 = embedded_proportions[k + 2]; + double u4 = embedded_proportions[k + 3]; + + double v1 = embedded_proportions[k + stripe + 1]; + double v2 = embedded_proportions[k + stripe + 2]; + double v3 = embedded_proportions[k + stripe + 3]; + double v4 = embedded_proportions[k + stripe + 4]; + + dm_stripe[k] += fabs(u1 - v1) * length; + dm_stripe[k + 1] += fabs(u2 - v2) * length; + dm_stripe[k + 2] += fabs(u3 - v3) * length; + dm_stripe[k + 3] += fabs(u4 - v4) * length; + } + + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { + double u = embedded_proportions[k]; + double v = embedded_proportions[k + stripe + 1]; + + dm_stripe[k] += fabs(u - v) * length; + } + } + } +} + +void su::_normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + // point of thread + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { + int k = j * 4; + double u1 = embedded_proportions[k]; + double u2 = embedded_proportions[k + 1]; + double u3 = embedded_proportions[k + 2]; + double u4 = embedded_proportions[k + 3]; + + double v1 = embedded_proportions[k + stripe + 1]; + double v2 = embedded_proportions[k + stripe + 2]; + double v3 = embedded_proportions[k + stripe + 3]; + double v4 = embedded_proportions[k + stripe + 4]; + + dm_stripe[k] += fabs(u1 - v1) * length; + dm_stripe[k + 1] += fabs(u2 - v2) * length; + dm_stripe[k + 2] += fabs(u3 - v3) * length; + dm_stripe[k + 3] += fabs(u4 - v4) * length; + + dm_stripe_total[k] += (u1 + v1) * length; + dm_stripe_total[k + 1] += (u2 + v2) * length; + dm_stripe_total[k + 2] += (u3 + v3) * length; + dm_stripe_total[k + 3] += (u4 + v4) * length; + } + + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { + double u = embedded_proportions[k]; + double v = embedded_proportions[k + stripe + 1]; + + dm_stripe[k] += fabs(u - v) * length; + dm_stripe_total[k] += (u + v) * length; + } + } + } +} + +void su::_generalized_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + // point of thread + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + double u1 = embedded_proportions[j]; + double v1 = embedded_proportions[j + stripe + 1]; + double sum1 = u1 + v1; + if(sum1 != 0.0) { + double sub1 = fabs(u1 - v1); + double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; + dm_stripe[j] += sum_pow1 * (sub1 / sum1); + dm_stripe_total[j] += sum_pow1; + } + } + } +} + +void su::_unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { + int k = j * 4; + int32_t u1 = embedded_proportions[k] > 0; + int32_t u2 = embedded_proportions[k + 1] > 0; + int32_t u3 = embedded_proportions[k + 2] > 0; + int32_t u4 = embedded_proportions[k + 3] > 0; + + int32_t v1 = embedded_proportions[k + stripe + 1] > 0; + int32_t v2 = embedded_proportions[k + stripe + 2] > 0; + int32_t v3 = embedded_proportions[k + stripe + 3] > 0; + int32_t v4 = embedded_proportions[k + stripe + 4] > 0; + + dm_stripe[k] += (u1 ^ v1) * length; + dm_stripe[k + 1] += (u2 ^ v2) * length; + dm_stripe[k + 2] += (u3 ^ v3) * length; + dm_stripe[k + 3] += (u4 ^ v4) * length; + + dm_stripe_total[k] += (u1 | v1) * length; + dm_stripe_total[k + 1] += (u2 | v2) * length; + dm_stripe_total[k + 2] += (u3 | v3) * length; + dm_stripe_total[k + 3] += (u4 | v4) * length; + } + + if((task_p->n_samples % 4) != 0) { + for(unsigned int k = task_p->n_samples - (task_p->n_samples % 4); k < task_p->n_samples; k++) { + int32_t u = embedded_proportions[k] > 0; + int32_t v = embedded_proportions[k + stripe + 1] > 0; + + dm_stripe[k] += (u ^ v) * length; + dm_stripe_total[k] += (u | v) * length; + } + } + } +} + diff --git a/sucpp/unifrac_task.hpp b/sucpp/unifrac_task.hpp new file mode 100644 index 00000000..2879235c --- /dev/null +++ b/sucpp/unifrac_task.hpp @@ -0,0 +1,33 @@ +#include +#include + +namespace su { + struct task_parameters { + uint32_t n_samples; // number of samples + unsigned int start; // starting stripe + unsigned int stop; // stopping stripe + unsigned int tid; // thread ID + double g_unifrac_alpha; // generalized unifrac alpha + }; + + void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p); + void _normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p); + void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p); + void _generalized_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double length, + const su::task_parameters* task_p); +} From 034cc4c8838d7951f76ae533f1278972a94f390b Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 16:39:35 -0700 Subject: [PATCH 11/23] VAW methods --- sucpp/su.cpp | 39 +++++++---- sucpp/test_su.cpp | 43 ++++++++++++ sucpp/unifrac.cpp | 103 +++++++++++++++++++++++++++-- sucpp/unifrac.hpp | 10 ++- sucpp/unifrac_task.cpp | 145 +++++++++++++++++++++++++++++++++-------- sucpp/unifrac_task.hpp | 29 +++++++++ 6 files changed, 322 insertions(+), 47 deletions(-) diff --git a/sucpp/su.cpp b/sucpp/su.cpp index f59edd8d..34cd4028 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -9,14 +9,15 @@ #include void usage() { - std::cout << "usage: ssu -i -o -m [METHOD] -t [-n threads] [-a alpha]" << std::endl; + std::cout << "usage: ssu -i -o -m [METHOD] -t [-n threads] [-a alpha] [--vaw]" << std::endl; std::cout << std::endl; - std::cout << " -i\tThe input BIOM table." << std::endl; - std::cout << " -t\tThe input phylogeny in newick." << std::endl; - std::cout << " -m\tThe method, [unweighted | weighted_normalized | weighted_unnormalized | generalized]." << std::endl; - std::cout << " -o\tThe output distance matrix." << std::endl; - std::cout << " -n\t[OPTIONAL] The number of threads, default is 1." << std::endl; - std::cout << " -a\t[OPTIONAL] Generalized UniFrac alpha, default is 1." << std::endl; + std::cout << " -i\t\tThe input BIOM table." << std::endl; + std::cout << " -t\t\tThe input phylogeny in newick." << std::endl; + std::cout << " -m\t\tThe method, [unweighted | weighted_normalized | weighted_unnormalized | generalized]." << std::endl; + std::cout << " -o\t\tThe output distance matrix." << std::endl; + std::cout << " -n\t\t[OPTIONAL] The number of threads, default is 1." << std::endl; + std::cout << " -a\t\t[OPTIONAL] Generalized UniFrac alpha, default is 1." << std::endl; + std::cout << " --vaw\t[OPTIONAL] Variance adjusted, default is to not adjust for variance." << std::endl; std::cout << std::endl; } @@ -97,6 +98,7 @@ int main(int argc, char **argv){ nthreads = atoi(nthreads_arg.c_str()); } + bool vaw = input.cmdOptionExists("--vaw"); double g_unifrac_alpha; if(gunifrac_arg.empty()) { g_unifrac_alpha = 1.0; @@ -150,13 +152,22 @@ int main(int argc, char **argv){ tasks[threads.size() - 1].stop = end; for(unsigned int tid = 0; tid < threads.size(); tid++) { - threads[tid] = std::thread(su::unifrac, - std::ref(table), - std::ref(tree_sheared), - method, - std::ref(dm_stripes), - std::ref(dm_stripes_total), - &tasks[tid]); + if(vaw) + threads[tid] = std::thread(su::unifrac_vaw, + std::ref(table), + std::ref(tree_sheared), + method, + std::ref(dm_stripes), + std::ref(dm_stripes_total), + &tasks[tid]); + else + threads[tid] = std::thread(su::unifrac, + std::ref(table), + std::ref(tree_sheared), + method, + std::ref(dm_stripes), + std::ref(dm_stripes_total), + &tasks[tid]); } for(unsigned int tid = 0; tid < threads.size(); tid++) { diff --git a/sucpp/test_su.cpp b/sucpp/test_su.cpp index ee25440f..56964236 100644 --- a/sucpp/test_su.cpp +++ b/sucpp/test_su.cpp @@ -743,6 +743,48 @@ void test_generalized_unifrac() { SUITE_END(); } +void test_vaw_unifrac_weighted_normalized() { + SUITE_START("test vaw weighted normalized unifrac"); + + std::vector threads(1); + su::BPTree tree = su::BPTree("(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1);"); + su::biom table = su::biom("test.biom"); + + // as computed by GUniFrac, the original implementation of VAW-UniFrac + // could not be found. + // Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 + //Sample1 0.0000000 0.4086040 0.6240185 0.4639481 0.2857143 0.2766318 + //Sample2 0.4086040 0.0000000 0.3798594 0.6884992 0.6807616 0.4735781 + //Sample3 0.6240185 0.3798594 0.0000000 0.7713254 0.8812897 0.5047114 + //Sample4 0.4639481 0.6884992 0.7713254 0.0000000 0.6666667 0.2709298 + //Sample5 0.2857143 0.6807616 0.8812897 0.6666667 0.0000000 0.4735991 + //Sample6 0.2766318 0.4735781 0.5047114 0.2709298 0.4735991 0.0000000 + // weighted normalized unifrac as computed above + + std::vector w_exp; + double w_stride1[] = {0.4086040, 0.3798594, 0.7713254, 0.6666667, 0.4735991, 0.2766318}; + double w_stride2[] = {0.6240185, 0.6884992, 0.8812897, 0.2709298, 0.2857143, 0.4735781}; + double w_stride3[] = {0.4639481, 0.6807616, 0.5047114, 0.4639481, 0.6807616, 0.5047114}; + w_exp.push_back(w_stride1); + w_exp.push_back(w_stride2); + w_exp.push_back(w_stride3); + std::vector w_strides = su::make_strides(6); + std::vector w_strides_total = su::make_strides(6); + su::task_parameters w_task_p; + w_task_p.start = 0; w_task_p.stop = 3; w_task_p.tid = 0; w_task_p.n_samples = 6; + w_task_p.g_unifrac_alpha = 1.0; + su::unifrac_vaw(table, tree, su::weighted_normalized, w_strides, w_strides_total, &w_task_p); + + for(unsigned int i = 0; i < 3; i++) { + for(unsigned int j = 0; j < 6; j++) { + ASSERT(fabs(w_strides[i][j] - w_exp[i][j]) < 0.000001); + } + free(w_strides[i]); + } + SUITE_END(); +} + + void test_make_strides() { SUITE_START("test make stripes"); std::vector exp; @@ -945,6 +987,7 @@ int main(int argc, char** argv) { test_unnormalized_weighted_unifrac(); test_normalized_weighted_unifrac(); test_generalized_unifrac(); + test_vaw_unifrac_weighted_normalized(); test_unifrac_sample_counts(); printf("\n"); diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 762a3974..def68a0d 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -124,11 +124,23 @@ void initialize_embedded(double*& prop, const su::task_parameters* task_p) { posix_memalign((void **)&prop, 32, sizeof(double) * task_p->n_samples * 2); if(prop == NULL) { fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + sizeof(double) * task_p->n_samples * 2, __FILE__, __LINE__); exit(EXIT_FAILURE); } } +void initialize_sample_counts(double*& counts, const su::task_parameters* task_p, biom &table) { + posix_memalign((void **)&counts, 32, sizeof(double) * task_p->n_samples * 2); + if(counts == NULL) { + fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", + sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + exit(EXIT_FAILURE); + } + for(unsigned int i = 0; i < table.n_samples; i++) { + counts[i] = table.sample_counts[i]; + counts[i + table.n_samples] = table.sample_counts[i]; + } +} void initialize_stripes(std::vector &dm_stripes, std::vector &dm_stripes_total, @@ -271,15 +283,98 @@ void su::unifrac(biom &table, free(embedded_proportions); } +void su::unifrac_vaw(biom &table, + BPTree &tree, + Method unifrac_method, + std::vector &dm_stripes, + std::vector &dm_stripes_total, + const su::task_parameters* task_p) { + + if(table.n_samples != task_p->n_samples) { + fprintf(stderr, "Task and table n_samples not equal\n"); + exit(EXIT_FAILURE); + } + + void (*func)(std::vector&, // dm_stripes + std::vector&, // dm_stripes_total + double*, // embedded_proportions + double*, // embedded_counts + double*, // sample total counts + double, // length + const su::task_parameters*); + + switch(unifrac_method) { + case unweighted: + func = &su::_vaw_unweighted_unifrac_task; + break; + case weighted_normalized: + func = &su::_vaw_normalized_weighted_unifrac_task; + break; + case weighted_unnormalized: + func = &su::_vaw_unnormalized_weighted_unifrac_task; + break; + case generalized: + func = &su::_vaw_generalized_unifrac_task; + break; + } + PropStack propstack(table.n_samples); + PropStack countstack(table.n_samples); + + uint32_t node; + double *node_proportions; + double *node_counts; + double *embedded_proportions; + double *embedded_counts; + double *sample_total_counts; + double length; + + initialize_embedded(embedded_proportions, task_p); + initialize_embedded(embedded_counts, task_p); + initialize_sample_counts(sample_total_counts, task_p, table); + initialize_stripes(std::ref(dm_stripes), std::ref(dm_stripes_total), unifrac_method, task_p); + + // - 1 to avoid root + for(unsigned int k = 0; k < (tree.nparens / 2) - 1; k++) { + node = tree.postorderselect(k); + length = tree.lengths[node]; + + node_proportions = propstack.pop(node); + node_counts = countstack.pop(node); + + set_proportions(node_proportions, tree, node, table, propstack); + set_proportions(node_counts, tree, node, table, countstack, false); + + embed_proportions(embedded_proportions, node_proportions, task_p->n_samples); + embed_proportions(embedded_counts, node_counts, task_p->n_samples); + + func(dm_stripes, dm_stripes_total, embedded_proportions, embedded_counts, sample_total_counts, length, task_p); + } + + if(unifrac_method == weighted_normalized || unifrac_method == unweighted || unifrac_method == generalized) { + for(unsigned int i = task_p->start; i < task_p->stop; i++) { + for(unsigned int j = 0; j < task_p->n_samples; j++) { + dm_stripes[i][j] = dm_stripes[i][j] / dm_stripes_total[i][j]; + } + } + } + + free(embedded_proportions); + free(embedded_counts); + free(sample_total_counts); +} void su::set_proportions(double* props, BPTree &tree, uint32_t node, biom &table, - PropStack &ps) { + PropStack &ps, + bool normalize) { if(tree.isleaf(node)) { table.get_obs_data(tree.names[node], props); - for(unsigned int i = 0; i < table.n_samples; i++) - props[i] = props[i] / table.sample_counts[i]; + for(unsigned int i = 0; i < table.n_samples; i++) { + props[i] = props[i]; + if(normalize) + props[i] /= table.sample_counts[i]; + } } else { unsigned int current = tree.leftchild(node); diff --git a/sucpp/unifrac.hpp b/sucpp/unifrac.hpp index 56c923b0..d5403cf3 100644 --- a/sucpp/unifrac.hpp +++ b/sucpp/unifrac.hpp @@ -27,11 +27,19 @@ namespace su { std::vector &dm_stripes_total, const task_parameters* task_p); + void unifrac_vaw(biom &table, + BPTree &tree, + Method unifrac_method, + std::vector &dm_stripes, + std::vector &dm_stripes_total, + const task_parameters* task_p); + double** deconvolute_stripes(std::vector &stripes, uint32_t n); void set_proportions(double* props, BPTree &tree, uint32_t node, biom &table, - PropStack &ps); + PropStack &ps, + bool normalize = true); std::vector make_strides(unsigned int n_samples); inline void embed_proportions(double* out, double* in, uint32_t n) { double val; diff --git a/sucpp/unifrac_task.cpp b/sucpp/unifrac_task.cpp index 1437a0b3..d84a60b4 100644 --- a/sucpp/unifrac_task.cpp +++ b/sucpp/unifrac_task.cpp @@ -7,7 +7,6 @@ void su::_unnormalized_weighted_unifrac_task(std::vector &__restrict__ double length, const su::task_parameters* task_p) { double *dm_stripe; - //__m256d ymm0, ymm1; for(unsigned int stripe=task_p->start; stripe < task_p->stop; stripe++) { dm_stripe = dm_stripes[stripe]; @@ -23,33 +22,6 @@ void su::_unnormalized_weighted_unifrac_task(std::vector &__restrict__ * support must be tested at compile time, so it's rather annoying * at the moment. basically, we can't assume the presence of avx2. */ - //for(j = 0; j < n_samples / 4; j++) { - // int offset = j * 4; - // ymm0 = _mm256_sub_pd(_mm256_load_pd(embedded_proportions + offset), - // _mm256_load_pd(embedded_proportions + (offset + stripe + 1))); // u - v - // ymm1 = _mm256_sub_pd(_mm256_set1_pd(0.0), - // ymm0); // 0.0 - (u - v) - // ymm0 = _mm256_max_pd(ymm0, ymm1); // abs(u - v) - // ymm0 = _mm256_fmadd_pd(ymm0, - // _mm256_set1_pd(length), - // _mm256_load_pd(dm_stripe + offset)); // fused mutiply add, dm_stripe[offset] += (abs(u - v) * length) - // _mm256_store_pd(dm_stripe + offset, ymm0); //store - //} - - //if((n_samples % 4) != 0) { - // for(int k = n_samples - (n_samples % 4); k < n_samples; k++) { - // double u = embedded_proportions[k]; - // ////double v = alt_proportions[j]; - // double v = embedded_proportions[k + stripe + 1]; - // dm_stripe[k] += fabs(u - v) * length; - // } - //} - - //for(int k = 0; k < n_samples; k++) { - // double u = embedded_proportions[k]; - // double v = embedded_proportions[k + stripe + 1]; - // dm_stripe[k] += fabs(u - v) * length; - //} for(unsigned int j = 0; j < task_p->n_samples / 4; j++) { int k = j * 4; double u1 = embedded_proportions[k]; @@ -79,6 +51,30 @@ void su::_unnormalized_weighted_unifrac_task(std::vector &__restrict__ } } +void su::_vaw_unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + for(unsigned int stripe=task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + double u = embedded_proportions[j]; + double v = embedded_proportions[j + stripe + 1]; + + double m = sample_total_counts[j] + sample_total_counts[j + stripe + 1]; + double mi = embedded_counts[j] + embedded_counts[j + stripe + 1]; + double vaw = sqrt(mi * (m - mi)); + + if(vaw > 0) + dm_stripe[j] += (fabs(u - v) * length) / vaw; + } + } +} void su::_normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, @@ -127,6 +123,36 @@ void su::_normalized_weighted_unifrac_task(std::vector &__restrict__ dm } } +void su::_vaw_normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + // point of thread + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + double u = embedded_proportions[j]; + double v = embedded_proportions[j + stripe + 1]; + + double m = sample_total_counts[j] + sample_total_counts[j + stripe + 1]; + double mi = embedded_counts[j] + embedded_counts[j + stripe + 1]; + double vaw = sqrt(mi * (m - mi)); + + if(vaw > 0) { + dm_stripe[j] += (fabs(u - v) * length) / vaw; + dm_stripe_total[j] += ((u + v) * length) / vaw; + } + } + } +} void su::_generalized_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, @@ -154,6 +180,39 @@ void su::_generalized_unifrac_task(std::vector &__restrict__ dm_stripes } } +void su::_vaw_generalized_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + // point of thread + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + double m = sample_total_counts[j] + sample_total_counts[j + stripe + 1]; + double mi = embedded_counts[j] + embedded_counts[j + stripe + 1]; + double vaw = sqrt(mi * (m - mi)); + + double u1 = embedded_proportions[j]; + double v1 = embedded_proportions[j + stripe + 1]; + + if(vaw > 0.0) { + double sum1 = (u1 + v1) / vaw; + double sub1 = fabs(u1 - v1) / vaw; + double sum_pow1 = pow(sum1, task_p->g_unifrac_alpha) * length; + dm_stripe[j] += sum_pow1 * (sub1 / sum1); + dm_stripe_total[j] += sum_pow1; + } + } + } +} void su::_unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, @@ -201,3 +260,33 @@ void su::_unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, } } +void su::_vaw_unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p) { + double *dm_stripe; + double *dm_stripe_total; + + for(unsigned int stripe = task_p->start; stripe < task_p->stop; stripe++) { + dm_stripe = dm_stripes[stripe]; + dm_stripe_total = dm_stripes_total[stripe]; + + for(unsigned int j = 0; j < task_p->n_samples; j++) { + int32_t u = embedded_proportions[j] > 0; + int32_t v = embedded_proportions[j + stripe + 1] > 0; + + double m = sample_total_counts[j] + sample_total_counts[j + stripe + 1]; + double mi = embedded_counts[j] + embedded_counts[j + stripe + 1]; + double vaw = sqrt(mi * (m - mi)); + + if(vaw > 0) { + dm_stripe[j] += ((u ^ v) * length) / vaw; + dm_stripe_total[j] += ((u | v) * length) / vaw; + } + } + } +} + diff --git a/sucpp/unifrac_task.hpp b/sucpp/unifrac_task.hpp index 2879235c..a2d470f2 100644 --- a/sucpp/unifrac_task.hpp +++ b/sucpp/unifrac_task.hpp @@ -30,4 +30,33 @@ namespace su { double* __restrict__ embedded_proportions, double length, const su::task_parameters* task_p); + + void _vaw_unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p); + void _vaw_normalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p); + void _vaw_unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p); + void _vaw_generalized_unifrac_task(std::vector &__restrict__ dm_stripes, + std::vector &__restrict__ dm_stripes_total, + double* __restrict__ embedded_proportions, + double* __restrict__ embedded_counts, + double* __restrict__ sample_total_counts, + double length, + const su::task_parameters* task_p); } From eb904a38933f05a214a7448dc2afc93d8a8d175d Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 16:46:05 -0700 Subject: [PATCH 12/23] Citations --- sucpp/su.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/sucpp/su.cpp b/sucpp/su.cpp index 34cd4028..9a5c8bfc 100644 --- a/sucpp/su.cpp +++ b/sucpp/su.cpp @@ -19,6 +19,17 @@ void usage() { std::cout << " -a\t\t[OPTIONAL] Generalized UniFrac alpha, default is 1." << std::endl; std::cout << " --vaw\t[OPTIONAL] Variance adjusted, default is to not adjust for variance." << std::endl; std::cout << std::endl; + std::cout << "Citations: " << std::endl; + std::cout << " For UniFrac, please see:" << std::endl; + std::cout << " Lozupone and Knight Appl Environ Microbiol 2005; DOI: 10.1128/AEM.71.12.8228-8235.2005" << std::endl; + std::cout << " Lozupone et al. Appl Environ Microbiol 2007; DOI: 10.1128/AEM.01996-06" << std::endl; + std::cout << " Hamady et al. ISME 2010; DOI: 10.1038/ismej.2009.97" << std::endl; + std::cout << " Lozupone et al. ISME 2011; DOI: 10.1038/ismej.2010.133" << std::endl; + std::cout << " For Generalized UniFrac, please see: " << std::endl; + std::cout << " Chen et al. Bioinformatics 2012; DOI: 10.1093/bioinformatics/bts342" << std::endl; + std::cout << " For Variance Adjusted UniFrac, please see: " << std::endl; + std::cout << " Chang et al. BMC Bioinformatics 2011; DOI: 10.1186/1471-2105-12-118" << std::endl; + std::cout << std::endl; } void err(std::string msg) { From a91850c5afa49cdce39eca567f52e3d40ccc02e6 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 16:52:42 -0700 Subject: [PATCH 13/23] Expose VAW --- q2_state_unifrac/_methods.py | 29 ++++++++++++++++++++--------- q2_state_unifrac/plugin_setup.py | 24 ++++++++++++++++++++---- 2 files changed, 40 insertions(+), 13 deletions(-) diff --git a/q2_state_unifrac/_methods.py b/q2_state_unifrac/_methods.py index dbe015b9..6865de1b 100644 --- a/q2_state_unifrac/_methods.py +++ b/q2_state_unifrac/_methods.py @@ -23,7 +23,7 @@ def _sanity(): raise ValueError("ssu could not be located!") -def _run(table_fp, tree_fp, output_fp, threads, method): +def _run(table_fp, tree_fp, output_fp, threads, method, vaw=False, alpha=None): cmd = [resource_filename(*ARGS), '-i', table_fp, '-t', tree_fp, @@ -31,53 +31,64 @@ def _run(table_fp, tree_fp, output_fp, threads, method): '-n', threads, '-m', method] + if vaw: + cmd.append('--vaw') + + if alpha is not None: + cmd.append('-a') + cmd.append(str(alpha)) + subprocess.run(cmd, check=True) def unweighted(table: BIOMV210Format, phylogeny: NewickFormat, - threads: int=1)-> skbio.DistanceMatrix: + threads: int=1, + vaw: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'unweighted') + 'unweighted', vaw) return skbio.DistanceMatrix.read(output_fp) def weighted_normalized(table: BIOMV210Format, phylogeny: NewickFormat, - threads: int=1)-> skbio.DistanceMatrix: + threads: int=1, + vaw: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'weighted_normalized') + 'weighted_normalized', vaw) return skbio.DistanceMatrix.read(output_fp) def weighted_unnormalized(table: BIOMV210Format, phylogeny: NewickFormat, - threads: int=1)-> skbio.DistanceMatrix: + threads: int=1, + vaw: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'weighted_unnormalized') + 'weighted_unnormalized', vaw) return skbio.DistanceMatrix.read(output_fp) def generalized(table: BIOMV210Format, phylogeny: NewickFormat, threads: int=1, - alpha: float=1.0)-> skbio.DistanceMatrix: + alpha: float=1.0, + vaw: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'generalized') + 'generalized', vaw, alpha) return skbio.DistanceMatrix.read(output_fp) diff --git a/q2_state_unifrac/plugin_setup.py b/q2_state_unifrac/plugin_setup.py index 511134d3..69d3a326 100644 --- a/q2_state_unifrac/plugin_setup.py +++ b/q2_state_unifrac/plugin_setup.py @@ -6,7 +6,7 @@ # The full license is in the file LICENSE, distributed with this software. # ---------------------------------------------------------------------------- -from qiime2.plugin import (Plugin, Properties, Int, Float) +from qiime2.plugin import (Plugin, Properties, Int, Float, Boolean) from q2_types.feature_table import FeatureTable, Frequency from q2_types.distance_matrix import DistanceMatrix from q2_types.tree import Phylogeny, Rooted @@ -28,8 +28,12 @@ function=q2_state_unifrac.unweighted, inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, - parameters={'threads': Int}, + parameters={'threads': Int, + 'variance-adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.'}, + 'variance-adjusted': + ('Perform variance adjustment based on ' + 'Chang et al. BMC Bioinformatics 2011'), input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' @@ -47,8 +51,12 @@ function=q2_state_unifrac.weighted_unnormalized, inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, - parameters={'threads': Int}, + parameters={'threads': Int, + 'variance-adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.'}, + 'variance-adjusted': + ('Perform variance adjustment based on ' + 'Chang et al. BMC Bioinformatics 2011'), input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' @@ -66,8 +74,12 @@ function=q2_state_unifrac.weighted_normalized, inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, - parameters={'threads': Int}, + parameters={'threads': Int, + 'variance-adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.'}, + 'variance-adjusted': + ('Perform variance adjustment based on ' + 'Chang et al. BMC Bioinformatics 2011'), input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' @@ -86,8 +98,12 @@ inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, + 'variance-adjusted': Bool, 'alpha': Float}, parameter_descriptions={'threads': 'The number of threads to use.', + 'variance-adjusted': + ('Perform variance adjustment based on ' + 'Chang et al. BMC Bioinformatics 2011'), 'alpha': ('The value of alpha controls importance ' 'of sample proportions. 1.0 is ' 'weighted normalized UniFrac. 0.0 is ' From 97232c6e2c0c5dfa0bdd15f904b213fecef6d644 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 17:05:44 -0700 Subject: [PATCH 14/23] travis c++ needs another include --- sucpp/unifrac_task.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/sucpp/unifrac_task.hpp b/sucpp/unifrac_task.hpp index a2d470f2..846489ad 100644 --- a/sucpp/unifrac_task.hpp +++ b/sucpp/unifrac_task.hpp @@ -1,5 +1,6 @@ #include #include +#include namespace su { struct task_parameters { From 2b1226182761d6b7ab457abd0092e5e01f0018ae Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 17:15:04 -0700 Subject: [PATCH 15/23] syntax... --- q2_state_unifrac/plugin_setup.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/q2_state_unifrac/plugin_setup.py b/q2_state_unifrac/plugin_setup.py index 69d3a326..9c57c36e 100644 --- a/q2_state_unifrac/plugin_setup.py +++ b/q2_state_unifrac/plugin_setup.py @@ -6,7 +6,7 @@ # The full license is in the file LICENSE, distributed with this software. # ---------------------------------------------------------------------------- -from qiime2.plugin import (Plugin, Properties, Int, Float, Boolean) +from qiime2.plugin import (Plugin, Properties, Int, Float, Bool) from q2_types.feature_table import FeatureTable, Frequency from q2_types.distance_matrix import DistanceMatrix from q2_types.tree import Phylogeny, Rooted @@ -30,10 +30,10 @@ 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, 'variance-adjusted': Bool}, - parameter_descriptions={'threads': 'The number of threads to use.'}, + parameter_descriptions={'threads': 'The number of threads to use.', 'variance-adjusted': ('Perform variance adjustment based on ' - 'Chang et al. BMC Bioinformatics 2011'), + 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' @@ -53,10 +53,10 @@ 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, 'variance-adjusted': Bool}, - parameter_descriptions={'threads': 'The number of threads to use.'}, + parameter_descriptions={'threads': 'The number of threads to use.', 'variance-adjusted': ('Perform variance adjustment based on ' - 'Chang et al. BMC Bioinformatics 2011'), + 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' @@ -76,10 +76,10 @@ 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, 'variance-adjusted': Bool}, - parameter_descriptions={'threads': 'The number of threads to use.'}, + parameter_descriptions={'threads': 'The number of threads to use.', 'variance-adjusted': ('Perform variance adjustment based on ' - 'Chang et al. BMC Bioinformatics 2011'), + 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ 'table': 'A rarefied FeatureTable', 'phylogeny': ('A rooted phylogeny which relates the observations in ' From 099340638d91bc71ef0db29cf438c56db7ee0ac8 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 18:04:27 -0700 Subject: [PATCH 16/23] wrong param name --- q2_state_unifrac/_methods.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/q2_state_unifrac/_methods.py b/q2_state_unifrac/_methods.py index 6865de1b..ca9d4ae9 100644 --- a/q2_state_unifrac/_methods.py +++ b/q2_state_unifrac/_methods.py @@ -23,7 +23,8 @@ def _sanity(): raise ValueError("ssu could not be located!") -def _run(table_fp, tree_fp, output_fp, threads, method, vaw=False, alpha=None): +def _run(table_fp, tree_fp, output_fp, threads, method, + variance_adjusted=False, alpha=None): cmd = [resource_filename(*ARGS), '-i', table_fp, '-t', tree_fp, @@ -31,7 +32,7 @@ def _run(table_fp, tree_fp, output_fp, threads, method, vaw=False, alpha=None): '-n', threads, '-m', method] - if vaw: + if variance_adjusted: cmd.append('--vaw') if alpha is not None: @@ -44,39 +45,39 @@ def _run(table_fp, tree_fp, output_fp, threads, method, vaw=False, alpha=None): def unweighted(table: BIOMV210Format, phylogeny: NewickFormat, threads: int=1, - vaw: bool=False)-> skbio.DistanceMatrix: + variance_adjusted: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'unweighted', vaw) + 'unweighted', variance_adjusted) return skbio.DistanceMatrix.read(output_fp) def weighted_normalized(table: BIOMV210Format, phylogeny: NewickFormat, threads: int=1, - vaw: bool=False)-> skbio.DistanceMatrix: + variance_adjusted: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'weighted_normalized', vaw) + 'weighted_normalized', variance_adjusted) return skbio.DistanceMatrix.read(output_fp) def weighted_unnormalized(table: BIOMV210Format, phylogeny: NewickFormat, threads: int=1, - vaw: bool=False)-> skbio.DistanceMatrix: + variance_adjusted: bool=False) -> skbio.DistanceMatrix: # noqa _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'weighted_unnormalized', vaw) + 'weighted_unnormalized', variance_adjusted) return skbio.DistanceMatrix.read(output_fp) @@ -84,11 +85,11 @@ def generalized(table: BIOMV210Format, phylogeny: NewickFormat, threads: int=1, alpha: float=1.0, - vaw: bool=False)-> skbio.DistanceMatrix: + variance_adjusted: bool=False)-> skbio.DistanceMatrix: _sanity() with tempfile.TemporaryDirectory() as tmp: output_fp = os.path.join(tmp, 'foo.dm') _run(str(table), str(phylogeny), output_fp, str(threads), - 'generalized', vaw, alpha) + 'generalized', variance_adjusted, alpha) return skbio.DistanceMatrix.read(output_fp) From be69604f42c31cee6e94ee606e3846b290d15ac0 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Wed, 16 Aug 2017 22:03:57 -0700 Subject: [PATCH 17/23] Missing a merge conflict --- sucpp/unifrac.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 669fd76d..08457ef1 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -104,8 +104,6 @@ double** su::deconvolute_stripes(std::vector &stripes, uint32_t n) { return dm; } -<<<<<<< HEAD -======= void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, From 2d30c6279bb0cb914a6aa1c4fc7c0a62c55aa464 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Thu, 17 Aug 2017 12:00:24 -0700 Subject: [PATCH 18/23] Fix param --- q2_state_unifrac/plugin_setup.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/q2_state_unifrac/plugin_setup.py b/q2_state_unifrac/plugin_setup.py index 9c57c36e..338ddc37 100644 --- a/q2_state_unifrac/plugin_setup.py +++ b/q2_state_unifrac/plugin_setup.py @@ -29,9 +29,9 @@ inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, - 'variance-adjusted': Bool}, + 'variance_adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.', - 'variance-adjusted': + 'variance_adjusted': ('Perform variance adjustment based on ' 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ @@ -52,9 +52,9 @@ inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, - 'variance-adjusted': Bool}, + 'variance_adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.', - 'variance-adjusted': + 'variance_adjusted': ('Perform variance adjustment based on ' 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ @@ -75,9 +75,9 @@ inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, - 'variance-adjusted': Bool}, + 'variance_adjusted': Bool}, parameter_descriptions={'threads': 'The number of threads to use.', - 'variance-adjusted': + 'variance_adjusted': ('Perform variance adjustment based on ' 'Chang et al. BMC Bioinformatics 2011')}, input_descriptions={ @@ -98,10 +98,10 @@ inputs={'table': FeatureTable[Frequency] % Properties('uniform-sampling'), 'phylogeny': Phylogeny[Rooted]}, parameters={'threads': Int, - 'variance-adjusted': Bool, + 'variance_adjusted': Bool, 'alpha': Float}, parameter_descriptions={'threads': 'The number of threads to use.', - 'variance-adjusted': + 'variance_adjusted': ('Perform variance adjustment based on ' 'Chang et al. BMC Bioinformatics 2011'), 'alpha': ('The value of alpha controls importance ' From b8e79979c7c5aa216241a7c5c3247b469ee0ea9b Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Fri, 18 Aug 2017 17:43:52 -0700 Subject: [PATCH 19/23] Merge conflict error --- sucpp/unifrac.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 4ab6d2e0..3bbbf400 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -402,9 +402,6 @@ void su::unifrac(biom &table, case generalized: func = &su::_generalized_unifrac_task; break; - case generalized: - func = &_generalized_unifrac_task; - break; } PropStack propstack(table.n_samples); From aa72ae8db4f75fd1c7b55f54191488d6782a8dcc Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Fri, 18 Aug 2017 18:00:15 -0700 Subject: [PATCH 20/23] Merge hell --- sucpp/unifrac.hpp | 8 -------- 1 file changed, 8 deletions(-) diff --git a/sucpp/unifrac.hpp b/sucpp/unifrac.hpp index f4182bef..d5403cf3 100644 --- a/sucpp/unifrac.hpp +++ b/sucpp/unifrac.hpp @@ -7,14 +7,6 @@ namespace su { enum Method {unweighted, weighted_normalized, weighted_unnormalized, generalized}; - struct task_parameters { - uint32_t n_samples; // number of samples - unsigned int start; // starting stripe - unsigned int stop; // stopping stripe - unsigned int tid; // thread ID - double g_unifrac_alpha; // generalized unifrac alpha - }; - class PropStack { private: std::stack prop_stack; From a9df5fceb6441986faec90d05c6709708414f652 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Fri, 18 Aug 2017 18:13:52 -0700 Subject: [PATCH 21/23] Not sure --- q2_state_unifrac/_methods.py | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/q2_state_unifrac/_methods.py b/q2_state_unifrac/_methods.py index 7ec379e4..ca9d4ae9 100644 --- a/q2_state_unifrac/_methods.py +++ b/q2_state_unifrac/_methods.py @@ -93,16 +93,3 @@ def generalized(table: BIOMV210Format, _run(str(table), str(phylogeny), output_fp, str(threads), 'generalized', variance_adjusted, alpha) return skbio.DistanceMatrix.read(output_fp) - - -def generalized(table: BIOMV210Format, - phylogeny: NewickFormat, - threads: int=1, - alpha: float=1.0)-> skbio.DistanceMatrix: - _sanity() - - with tempfile.TemporaryDirectory() as tmp: - output_fp = os.path.join(tmp, 'foo.dm') - _run(str(table), str(phylogeny), output_fp, str(threads), - 'generalized') - return skbio.DistanceMatrix.read(output_fp) From 619aeb9ae01678b1058ef699ae0b6ab688674d0a Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Mon, 21 Aug 2017 08:26:53 -0700 Subject: [PATCH 22/23] Addressing @eldeveloper's comments --- sucpp/test_su.cpp | 36 ++++++++++---------- sucpp/unifrac.cpp | 75 +++++++++++++++++++++++++++--------------- sucpp/unifrac_task.hpp | 43 ++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 45 deletions(-) diff --git a/sucpp/test_su.cpp b/sucpp/test_su.cpp index 6f187385..446c89ff 100644 --- a/sucpp/test_su.cpp +++ b/sucpp/test_su.cpp @@ -744,24 +744,24 @@ void test_generalized_unifrac() { } void test_vaw_unifrac_weighted_normalized() { - SUITE_START("test vaw weighted normalized unifrac"); - - std::vector threads(1); - su::BPTree tree = su::BPTree("(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1);"); - su::biom table = su::biom("test.biom"); - + SUITE_START("test vaw weighted normalized unifrac"); + + std::vector threads(1); + su::BPTree tree = su::BPTree("(GG_OTU_1:1,(GG_OTU_2:1,GG_OTU_3:1):1,(GG_OTU_5:1,GG_OTU_4:1):1);"); + su::biom table = su::biom("test.biom"); + // as computed by GUniFrac, the original implementation of VAW-UniFrac // could not be found. - // Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 - //Sample1 0.0000000 0.4086040 0.6240185 0.4639481 0.2857143 0.2766318 - //Sample2 0.4086040 0.0000000 0.3798594 0.6884992 0.6807616 0.4735781 - //Sample3 0.6240185 0.3798594 0.0000000 0.7713254 0.8812897 0.5047114 - //Sample4 0.4639481 0.6884992 0.7713254 0.0000000 0.6666667 0.2709298 - //Sample5 0.2857143 0.6807616 0.8812897 0.6666667 0.0000000 0.4735991 - //Sample6 0.2766318 0.4735781 0.5047114 0.2709298 0.4735991 0.0000000 + // Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 + //Sample1 0.0000000 0.4086040 0.6240185 0.4639481 0.2857143 0.2766318 + //Sample2 0.4086040 0.0000000 0.3798594 0.6884992 0.6807616 0.4735781 + //Sample3 0.6240185 0.3798594 0.0000000 0.7713254 0.8812897 0.5047114 + //Sample4 0.4639481 0.6884992 0.7713254 0.0000000 0.6666667 0.2709298 + //Sample5 0.2857143 0.6807616 0.8812897 0.6666667 0.0000000 0.4735991 + //Sample6 0.2766318 0.4735781 0.5047114 0.2709298 0.4735991 0.0000000 // weighted normalized unifrac as computed above - std::vector w_exp; + std::vector w_exp; double w_stride1[] = {0.4086040, 0.3798594, 0.7713254, 0.6666667, 0.4735991, 0.2766318}; double w_stride2[] = {0.6240185, 0.6884992, 0.8812897, 0.2709298, 0.2857143, 0.4735781}; double w_stride3[] = {0.4639481, 0.6807616, 0.5047114, 0.4639481, 0.6807616, 0.5047114}; @@ -775,7 +775,7 @@ void test_vaw_unifrac_weighted_normalized() { w_task_p.g_unifrac_alpha = 1.0; su::unifrac_vaw(table, tree, su::weighted_normalized, w_strides, w_strides_total, &w_task_p); - for(unsigned int i = 0; i < 3; i++) { + for(unsigned int i = 0; i < 3; i++) { for(unsigned int j = 0; j < 6; j++) { ASSERT(fabs(w_strides[i][j] - w_exp[i][j]) < 0.000001); } @@ -868,7 +868,7 @@ void test_bptree_shear_simple() { su::BPTree tree = su::BPTree("((3:2,4:3,(6:5)5:4)2:1,7:6,((10:9,11:10)9:8)8:7)r"); // simple - std::unordered_set to_keep = {"4", "6", "7", "10", "11"}; + std::unordered_set to_keep = {"4", "6", "7", "10", "11"}; uint32_t exp_nparens = 20; std::vector exp_structure = {true, true, true, false, true, true, false, false, false, true, @@ -890,7 +890,7 @@ void test_bptree_shear_deep() { su::BPTree tree = su::BPTree("((3:2,4:3,(6:5)5:4)2:1,7:6,((10:9,11:10)9:8)8:7)r"); // deep - std::unordered_set to_keep = {"10", "11"}; + std::unordered_set to_keep = {"10", "11"}; uint32_t exp_nparens = 10; std::vector exp_structure = {true, true, true, true, false, true, false, false, false, false}; @@ -935,7 +935,7 @@ void test_bptree_collapse_edge() { ASSERT(obs.names == exp.names); ASSERT(vec_almost_equal(obs.lengths, exp.lengths)); - SUITE_END(); + SUITE_END(); } void test_unifrac_sample_counts() { diff --git a/sucpp/unifrac.cpp b/sucpp/unifrac.cpp index 3bbbf400..528647d9 100644 --- a/sucpp/unifrac.cpp +++ b/sucpp/unifrac.cpp @@ -51,11 +51,12 @@ double* PropStack::pop(uint32_t node) { * add it to our record of known vectors so we can track our mallocs */ double *vec; + int err = 0; if(prop_stack.empty()) { - posix_memalign((void **)&vec, 32, sizeof(double) * defaultsize); - if(vec == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * defaultsize, __FILE__, __LINE__); + err = posix_memalign((void **)&vec, 32, sizeof(double) * defaultsize); + if(vec == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes, err %d; [%s]:%d\n", + sizeof(double) * defaultsize, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } } @@ -306,10 +307,10 @@ void _unweighted_unifrac_task(std::vector &__restrict__ dm_stripes, void progressbar(float progress) { - // from http://stackoverflow.com/a/14539953 + // from http://stackoverflow.com/a/14539953 // // could encapsulate into a classs for displaying time elapsed etc - int barWidth = 70; + int barWidth = 70; std::cout << "["; int pos = barWidth * progress; for (int i = 0; i < barWidth; ++i) { @@ -322,19 +323,21 @@ void progressbar(float progress) { } void initialize_embedded(double*& prop, const su::task_parameters* task_p) { - posix_memalign((void **)&prop, 32, sizeof(double) * task_p->n_samples * 2); - if(prop == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * task_p->n_samples * 2, __FILE__, __LINE__); + int err = 0; + err = posix_memalign((void **)&prop, 32, sizeof(double) * task_p->n_samples * 2); + if(prop == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes, err %d; [%s]:%d\n", + sizeof(double) * task_p->n_samples * 2, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } } void initialize_sample_counts(double*& counts, const su::task_parameters* task_p, biom &table) { - posix_memalign((void **)&counts, 32, sizeof(double) * task_p->n_samples * 2); - if(counts == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + int err = 0; + err = posix_memalign((void **)&counts, 32, sizeof(double) * task_p->n_samples * 2); + if(counts == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes, err %d; [%s]:%d\n", + sizeof(double) * task_p->n_samples, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } for(unsigned int i = 0; i < table.n_samples; i++) { @@ -347,21 +350,22 @@ void initialize_stripes(std::vector &dm_stripes, std::vector &dm_stripes_total, Method unifrac_method, const su::task_parameters* task_p) { + int err = 0; for(unsigned int i = task_p->start; i < task_p->stop; i++){ - posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * task_p->n_samples); - if(dm_stripes[i] == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + err = posix_memalign((void **)&dm_stripes[i], 32, sizeof(double) * task_p->n_samples); + if(dm_stripes[i] == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes, err %d; [%s]:%d\n", + sizeof(double) * task_p->n_samples, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } for(unsigned int j = 0; j < task_p->n_samples; j++) dm_stripes[i][j] = 0.; if(unifrac_method == unweighted || unifrac_method == weighted_normalized) { - posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * task_p->n_samples); - if(dm_stripes_total[i] == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * task_p->n_samples, __FILE__, __LINE__); + err = posix_memalign((void **)&dm_stripes_total[i], 32, sizeof(double) * task_p->n_samples); + if(dm_stripes_total[i] == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes err %d; [%s]:%d\n", + sizeof(double) * task_p->n_samples, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } for(unsigned int j = 0; j < task_p->n_samples; j++) @@ -402,7 +406,16 @@ void su::unifrac(biom &table, case generalized: func = &su::_generalized_unifrac_task; break; + default: + func = NULL; + break; } + + if(func == NULL) { + fprintf(stderr, "Unknown unifrac task\n"); + exit(1); + } + PropStack propstack(table.n_samples); uint32_t node; @@ -470,7 +483,7 @@ void su::unifrac(biom &table, // should make this compile-time support //if((tid == 0) && ((k % 1000) == 0)) - // progressbar((float)k / (float)(tree.nparens / 2)); + // progressbar((float)k / (float)(tree.nparens / 2)); } if(unifrac_method == weighted_normalized || unifrac_method == unweighted || unifrac_method == generalized) { @@ -517,6 +530,13 @@ void su::unifrac_vaw(biom &table, case generalized: func = &su::_vaw_generalized_unifrac_task; break; + default: + func = NULL; + break; + } + if(func == NULL) { + fprintf(stderr, "Unknown unifrac task\n"); + exit(1); } PropStack propstack(table.n_samples); PropStack countstack(table.n_samples); @@ -601,12 +621,13 @@ std::vector su::make_strides(unsigned int n_samples) { uint32_t n_rotations = (n_samples + 1) / 2; std::vector dm_stripes(n_rotations); + int err = 0; for(unsigned int i = 0; i < n_rotations; i++) { double* tmp; - posix_memalign((void **)&tmp, 32, sizeof(double) * n_samples); - if(tmp == NULL) { - fprintf(stderr, "Failed to allocate %zd bytes; [%s]:%d\n", - sizeof(double) * n_samples, __FILE__, __LINE__); + err = posix_memalign((void **)&tmp, 32, sizeof(double) * n_samples); + if(tmp == NULL || err != 0) { + fprintf(stderr, "Failed to allocate %zd bytes, err %d; [%s]:%d\n", + sizeof(double) * n_samples, err, __FILE__, __LINE__); exit(EXIT_FAILURE); } for(unsigned int j = 0; j < n_samples; j++) diff --git a/sucpp/unifrac_task.hpp b/sucpp/unifrac_task.hpp index 846489ad..013c9e7a 100644 --- a/sucpp/unifrac_task.hpp +++ b/sucpp/unifrac_task.hpp @@ -3,14 +3,38 @@ #include namespace su { + /* task specific compute parameters + * + * n_samples the number of samples being processed + * start the first stride to process + * stop the last stride to process + * tid the thread identifier + * g_unifrac_alpha an alpha value for generalized unifrac + */ struct task_parameters { uint32_t n_samples; // number of samples unsigned int start; // starting stripe unsigned int stop; // stopping stripe unsigned int tid; // thread ID + + // task specific arguments below double g_unifrac_alpha; // generalized unifrac alpha }; + /* void su::unifrac tasks + * + * all methods utilize the same function signature. that signature is as follows: + * + * dm_stripes vector the stripes of the distance matrix being accumulated + * into for unique branch length + * dm_stripes vector the stripes of the distance matrix being accumulated + * into for total branch length (e.g., to normalize unweighted unifrac) + * embedded_proportions the proportions vector for a sample, or rather + * the counts vector normalized to 1. this vector is embedded as it is + * duplicated: if A, B and C are proportions for features A, B, and C, the + * vector will look like [A B C A B C]. + * length the branch length of the current node to its parent. + * tasp_p task specific parameters. void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, @@ -32,6 +56,25 @@ namespace su { double length, const su::task_parameters* task_p); + /* void su::unifrac_vaw tasks + * + * all methods utilize the same function signature. that signature is as follows: + * + * dm_stripes vector the stripes of the distance matrix being accumulated + * into for unique branch length + * dm_stripes vector the stripes of the distance matrix being accumulated + * into for total branch length (e.g., to normalize unweighted unifrac) + * embedded_proportions the proportions vector for a sample, or rather + * the counts vector normalized to 1. this vector is embedded as it is + * duplicated: if A, B and C are proportions for features A, B, and C, the + * vector will look like [A B C A B C]. + * embedded_counts the counts vector embedded in the same way and order as + * embedded_proportions. the values of this array are unnormalized feature + * counts for the subtree. + * sample_total_counts the total unnormalized feature counts for all samples + * embedded in the same way and order as embedded_proportions. + * length the branch length of the current node to its parent. + * tasp_p task specific parameters. void _vaw_unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, From 6a4f1a644611a4e9046b20abcee5197661b58425 Mon Sep 17 00:00:00 2001 From: Daniel McDonald Date: Mon, 21 Aug 2017 08:47:29 -0700 Subject: [PATCH 23/23] Missed close of comment --- sucpp/unifrac_task.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/sucpp/unifrac_task.hpp b/sucpp/unifrac_task.hpp index 013c9e7a..4a2dec6c 100644 --- a/sucpp/unifrac_task.hpp +++ b/sucpp/unifrac_task.hpp @@ -34,7 +34,8 @@ namespace su { * duplicated: if A, B and C are proportions for features A, B, and C, the * vector will look like [A B C A B C]. * length the branch length of the current node to its parent. - * tasp_p task specific parameters. + * task_p task specific parameters. + */ void _unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions, @@ -74,7 +75,8 @@ namespace su { * sample_total_counts the total unnormalized feature counts for all samples * embedded in the same way and order as embedded_proportions. * length the branch length of the current node to its parent. - * tasp_p task specific parameters. + * task_p task specific parameters. + */ void _vaw_unnormalized_weighted_unifrac_task(std::vector &__restrict__ dm_stripes, std::vector &__restrict__ dm_stripes_total, double* __restrict__ embedded_proportions,