diff --git a/public/c++/VectorPairHMM/LoadTimeInitializer.cc b/public/c++/VectorPairHMM/LoadTimeInitializer.cc index 4188ea3fb..a1843141e 100644 --- a/public/c++/VectorPairHMM/LoadTimeInitializer.cc +++ b/public/c++/VectorPairHMM/LoadTimeInitializer.cc @@ -86,3 +86,25 @@ void LoadTimeInitializer::debug_close() m_filename_to_fptr.clear(); } +void LoadTimeInitializer::dump_sandbox(unsigned haplotypeLength, unsigned readLength, char* haplotypeBasesArray, testcase& tc) +{ + ofstream& dumpFptr = m_sandbox_fptr; + for(unsigned k=0;k +#include "template.h" class LoadTimeInitializer { public: @@ -10,6 +11,10 @@ class LoadTimeInitializer void debug_dump(std::string filename, std::string s, bool to_append, bool add_newline=true); void debug_close(); + void dump_sandbox(unsigned haplotypeLength, unsigned readLength, char* haplotypeBasesArray, testcase& tc); + void open_sandbox() { m_sandbox_fptr.open("sandbox.txt", std::ios::app); } + void close_sandbox() { m_sandbox_fptr.close(); } + jfieldID m_readBasesFID; jfieldID m_readQualsFID; jfieldID m_insertionGOPFID; @@ -37,6 +42,7 @@ class LoadTimeInitializer uint64_t m_bytes_copied; private: std::map m_filename_to_fptr; + std::ofstream m_sandbox_fptr; }; extern LoadTimeInitializer g_load_time_initializer; diff --git a/public/c++/VectorPairHMM/jni_common.h b/public/c++/VectorPairHMM/jni_common.h index d1ce86e9c..1e23f66ee 100644 --- a/public/c++/VectorPairHMM/jni_common.h +++ b/public/c++/VectorPairHMM/jni_common.h @@ -7,6 +7,7 @@ //#define DEBUG 1 //#define DEBUG0_1 1 //#define DEBUG3 1 +/*#define DUMP_TO_SANDBOX 1*/ #define DIRECT_ACCESS_TO_JAVA_HEAP_MEMORY 1 diff --git a/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc b/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc index e837e5765..ba31da420 100644 --- a/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc +++ b/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc @@ -5,6 +5,16 @@ #include "utils.h" #include "LoadTimeInitializer.h" +char* all_ptrs[] = { + "TCAAACCGAAATAAAGGCCAGTATATCCATATCCTTCCCATAAATGTTGATGGAAGAATTATTTGGAAGCCATATAGAATGAAATGACTCTATACACAAATTAAAACACAAAAACGTACTCAAAATAGTCCAGAGACTACAACTTCAAATGCAAAACTATAAATAATCTAAAAGAAAACCTAAGAGACATTC", + "GTCCAGAGACTACAACTTCAAATGCAAAACTATAAATAATCTAACAGAAAACCTAAGAGACATTC", + ">D?@BAEEEEDEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE@EEEEEEEEEEDEEEEEE?", + "IIIIIIIIIIIIIIIIIIIIIIIIHHHHIIIIIIIIIIIIIIIIIIHHHHIIIIIIIIIIIIIIN", + "IIIIIIIIIIIIIIIIIIIIIIIIHHHHIIIIIIIIIIIIIIIIIIHHHHIIIIIIIIIIIIIIN", + "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" +}; +char all_arrays[6][16384]; + using namespace std; JNIEXPORT jlong JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniGetMachineType @@ -46,6 +56,17 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless initialize_function_pointers((uint64_t)mask); cout.flush(); } +#if 0 + for(unsigned i=0;i<6;++i) + { + unsigned length = strlen(all_ptrs[i]); + for(unsigned j=0;j<16384;++j) + all_arrays[i][j] = all_ptrs[i][j%length]; + } + for(unsigned i=2;i<6;++i) + for(unsigned j=0;j<16384;++j) + all_arrays[i][j] = all_arrays[i][j]-33; +#endif } //Since the list of haplotypes against which the reads are evaluated in PairHMM is the same for a region, @@ -127,6 +148,9 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless readBasesArrayVector.resize(numReads); #ifdef DO_PROFILING start_time = get_time(); +#endif +#ifdef DUMP_TO_SANDBOX + g_load_time_initializer.open_sandbox(); #endif for(unsigned i=0;i tc_vector; tc_vector.clear(); - testcase tc; + tc_vector.resize(BATCH_SIZE+4); + vector results_vec; + results_vec.clear(); + results_vec.resize(tc_vector.size()); + vector baseline_results; + baseline_results.clear(); + baseline_results.resize(tc_vector.size()); + + bool all_ok = true; uint64_t total_time = 0; + uint64_t baseline_time = 0; + unsigned total_count = 0; + unsigned num_testcases = 0; + //unsigned curr_batch_size = rand()%BATCH_SIZE + 4; //min batch size + unsigned curr_batch_size = BATCH_SIZE; while(1) { - int break_value = use_old_read_testcase ? read_testcase(&tc, fptr) : read_mod_testcase(ifptr,&tc,true); + int break_value = use_old_read_testcase ? read_testcase(&(tc_vector[num_testcases]), fptr) : + read_mod_testcase(ifptr,&(tc_vector[num_testcases]),true); if(break_value >= 0) - tc_vector.push_back(tc); - if(tc_vector.size() == BATCH_SIZE || (break_value < 0 && tc_vector.size() > 0)) + ++num_testcases; + if(num_testcases == curr_batch_size || (break_value < 0 && num_testcases > 0)) { - vector results_vec; - results_vec.clear(); - results_vec.resize(tc_vector.size()); get_time(); #pragma omp parallel for schedule(dynamic,chunk_size) num_threads(12) - for(unsigned i=0;i(&tc); - baseline_result = log10(baseline_result) - log10(ldexp(1.0, 1020.0)); + testcase& tc = tc_vector[i]; + float result_avxf = compute_full_prob(&tc); + double result = 0; + if (result_avxf < MIN_ACCEPTED) { + double result_avxd = compute_full_prob(&tc); + result = log10(result_avxd) - log10(ldexp(1.0, 1020.0)); + } + else + result = (double)(log10f(result_avxf) - log10f(ldexpf(1.f, 120.f))); + baseline_results[i] = result; + } + baseline_time += get_time(); + for(unsigned i=0;i 1e-5 && rel_error > 1e-5) - cout << std::scientific << baseline_result << " "<hap, tc->rs, q, i, d, c) != 6) return -1; + tc->haplen = strlen(tc->hap); tc->rslen = strlen(tc->rs); + assert(strlen(q) == tc->rslen); + assert(strlen(i) == tc->rslen); + assert(strlen(d) == tc->rslen); + assert(strlen(c) == tc->rslen); //assert(tc->rslen < MROWS); tc->ihap = (int *) malloc(tc->haplen*sizeof(int)); tc->irs = (int *) malloc(tc->rslen*sizeof(int)); @@ -111,7 +116,8 @@ int read_testcase(testcase *tc, FILE* ifp) _i = normalize(i[x]); _d = normalize(d[x]); _c = normalize(c[x]); - tc->q[x] = (_q < 6) ? 6 : _q; + tc->q[x] = (_q < 6) ? 6 : _q; + //tc->q[x] = _q; tc->i[x] = _i; tc->d[x] = _d; tc->c[x] = _c; @@ -120,13 +126,15 @@ int read_testcase(testcase *tc, FILE* ifp) for (x = 0; x < tc->haplen; x++) tc->ihap[x] = tc->hap[x]; - + free(q); free(i); free(d); free(c); free(line); + + return 0; }