diff --git a/public/c++/VectorPairHMM/.gitignore b/public/c++/VectorPairHMM/.gitignore index ce903a6fa..b034f9461 100644 --- a/public/c++/VectorPairHMM/.gitignore +++ b/public/c++/VectorPairHMM/.gitignore @@ -6,6 +6,7 @@ tests hmm_Mohammad pairhmm-template-main *.swp +*.class checker reformat subdir_checkout.sh diff --git a/public/c++/VectorPairHMM/LoadTimeInitializer.cc b/public/c++/VectorPairHMM/LoadTimeInitializer.cc index d403ff892..fb640ef88 100644 --- a/public/c++/VectorPairHMM/LoadTimeInitializer.cc +++ b/public/c++/VectorPairHMM/LoadTimeInitializer.cc @@ -86,7 +86,7 @@ void LoadTimeInitializer::debug_close() m_filename_to_fptr.clear(); } -void LoadTimeInitializer::dump_sandbox(testcase& tc) +void LoadTimeInitializer::dump_sandbox(testcase& tc, unsigned tc_idx, unsigned numReads, unsigned numHaplotypes) { unsigned haplotypeLength = tc.haplen; unsigned readLength = tc.rslen; @@ -108,5 +108,7 @@ void LoadTimeInitializer::dump_sandbox(testcase& tc) dumpFptr<<" "; for(unsigned k=0;kGetStringUTFChars(fileNameString, 0); + do_compute((char*)fileName); + env->ReleaseStringUTFChars(fileNameString, fileName); +} + diff --git a/public/c++/VectorPairHMM/Sandbox.h b/public/c++/VectorPairHMM/Sandbox.h new file mode 100644 index 000000000..4ac1ea24c --- /dev/null +++ b/public/c++/VectorPairHMM/Sandbox.h @@ -0,0 +1,71 @@ +/* DO NOT EDIT THIS FILE - it is machine generated */ +#include +/* Header for class Sandbox */ + +#ifndef _Included_Sandbox +#define _Included_Sandbox +#ifdef __cplusplus +extern "C" { +#endif +#undef Sandbox_enableAll +#define Sandbox_enableAll -1LL +/* + * Class: Sandbox + * Method: jniGetMachineType + * Signature: ()J + */ +JNIEXPORT jlong JNICALL Java_Sandbox_jniGetMachineType + (JNIEnv *, jobject); + +/* + * Class: Sandbox + * Method: jniInitializeClassFieldsAndMachineMask + * Signature: (Ljava/lang/Class;Ljava/lang/Class;J)V + */ +JNIEXPORT void JNICALL Java_Sandbox_jniInitializeClassFieldsAndMachineMask + (JNIEnv *, jobject, jclass, jclass, jlong); + +/* + * Class: Sandbox + * Method: jniInitializeHaplotypes + * Signature: (I[LSandbox/JNIHaplotypeDataHolderClass;)V + */ +JNIEXPORT void JNICALL Java_Sandbox_jniInitializeHaplotypes + (JNIEnv *, jobject, jint, jobjectArray); + +/* + * Class: Sandbox + * Method: jniFinalizeRegion + * Signature: ()V + */ +JNIEXPORT void JNICALL Java_Sandbox_jniFinalizeRegion + (JNIEnv *, jobject); + +/* + * Class: Sandbox + * Method: jniComputeLikelihoods + * Signature: (II[LSandbox/JNIReadDataHolderClass;[LSandbox/JNIHaplotypeDataHolderClass;[DI)V + */ +JNIEXPORT void JNICALL Java_Sandbox_jniComputeLikelihoods + (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jdoubleArray, jint); + +/* + * Class: Sandbox + * Method: jniClose + * Signature: ()V + */ +JNIEXPORT void JNICALL Java_Sandbox_jniClose + (JNIEnv *, jobject); + +/* + * Class: Sandbox + * Method: doEverythingNative + * Signature: ([B)V + */ +JNIEXPORT void JNICALL Java_Sandbox_doEverythingNative + (JNIEnv *, jobject, jstring); + +#ifdef __cplusplus +} +#endif +#endif diff --git a/public/c++/VectorPairHMM/Sandbox.java b/public/c++/VectorPairHMM/Sandbox.java new file mode 100644 index 000000000..c41a276c2 --- /dev/null +++ b/public/c++/VectorPairHMM/Sandbox.java @@ -0,0 +1,278 @@ +import java.util.List; +import java.util.LinkedList; +import java.util.Map; +import java.util.HashMap; +import java.io.File; +import java.util.Scanner; +import java.io.IOException; +import java.io.FileNotFoundException; +import java.io.InputStreamReader; + +public class Sandbox { + + private long setupTime = 0; + private long computeTime = 0; + //Used to copy references to byteArrays to JNI from reads + protected class JNIReadDataHolderClass { + public byte[] readBases = null; + public byte[] readQuals = null; + public byte[] insertionGOP = null; + public byte[] deletionGOP = null; + public byte[] overallGCP = null; + } + + //Used to copy references to byteArrays to JNI from haplotypes + protected class JNIHaplotypeDataHolderClass { + public byte[] haplotypeBases = null; + } + + /** + * Return 64-bit mask representing machine capabilities + * Bit 0 is LSB, bit 63 MSB + * Bit 0 represents sse4.2 availability + * Bit 1 represents AVX availability + */ + public native long jniGetMachineType(); + public static final long enableAll = 0xFFFFFFFFFFFFFFFFl; + + + /** + * Function to initialize the fields of JNIReadDataHolderClass and JNIHaplotypeDataHolderClass from JVM. + * C++ codegets FieldIDs for these classes once and re-uses these IDs for the remainder of the program. Field IDs do not + * change per JVM session + * @param readDataHolderClass class type of JNIReadDataHolderClass + * @param haplotypeDataHolderClass class type of JNIHaplotypeDataHolderClass + * @param mask mask is a 64 bit integer identical to the one received from jniGetMachineType(). Users can disable usage of some hardware features by zeroing some bits in the mask + * */ + private native void jniInitializeClassFieldsAndMachineMask(Class readDataHolderClass, Class haplotypeDataHolderClass, long mask); + + private static Boolean isVectorLoglessPairHMMLibraryLoaded = false; + //The constructor is called only once inside PairHMMLikelihoodCalculationEngine + public Sandbox() { + synchronized(isVectorLoglessPairHMMLibraryLoaded) { + //Load the library and initialize the FieldIDs + if(!isVectorLoglessPairHMMLibraryLoaded) { + System.loadLibrary("VectorLoglessPairHMM"); + isVectorLoglessPairHMMLibraryLoaded = true; + jniInitializeClassFieldsAndMachineMask(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, enableAll); //need to do this only once + } + } + } + + private native void jniInitializeHaplotypes(final int numHaplotypes, JNIHaplotypeDataHolderClass[] haplotypeDataArray); + + //Used to transfer data to JNI + //Since the haplotypes are the same for all calls to computeLikelihoods within a region, transfer the haplotypes only once to the JNI per region + public void initialize(final List haplotypes) { + int numHaplotypes = haplotypes.size(); + JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes]; + int idx = 0; + for(final JNIHaplotypeDataHolderClass currHaplotype : haplotypes) + { + haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass(); + haplotypeDataArray[idx].haplotypeBases = currHaplotype.haplotypeBases; + ++idx; + } + jniInitializeHaplotypes(numHaplotypes, haplotypeDataArray); + } + /** + * Tell JNI to release arrays - really important if native code is directly accessing Java memory, if not + * accessing Java memory directly, still important to release memory from C++ + */ + private native void jniFinalizeRegion(); + + + public void finalizeRegion() + { + jniFinalizeRegion(); + } + + /** + * Real compute kernel + */ + private native void jniComputeLikelihoods(int numReads, int numHaplotypes, JNIReadDataHolderClass[] readDataArray, + JNIHaplotypeDataHolderClass[] haplotypeDataArray, double[] likelihoodArray, int maxNumThreadsToUse); + + public void computeLikelihoods(final List reads, final List haplotypes) { + //System.out.println("Region : "+reads.size()+" x "+haplotypes.size()); + long startTime = System.nanoTime(); + int readListSize = reads.size(); + int numHaplotypes = haplotypes.size(); + int numTestcases = readListSize*numHaplotypes; + JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize]; + int idx = 0; + for(JNIReadDataHolderClass read : reads) + { + readDataArray[idx] = new JNIReadDataHolderClass(); + readDataArray[idx].readBases = read.readBases; + readDataArray[idx].readQuals = read.readQuals; + readDataArray[idx].insertionGOP = read.insertionGOP; + readDataArray[idx].deletionGOP = read.deletionGOP; + readDataArray[idx].overallGCP = read.overallGCP; + ++idx; + } + + double[] mLikelihoodArray = new double[readListSize*numHaplotypes]; //to store results + setupTime += (System.nanoTime() - startTime); + //for(reads) + // for(haplotypes) + // compute_full_prob() + jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, null, mLikelihoodArray, 12); + + computeTime += (System.nanoTime() - startTime); + } + + /** + * Print final profiling information from native code + */ + public native void jniClose(); + public void close() + { + System.out.println("Time spent in setup for JNI call : "+(setupTime*1e-9)+" compute time : "+(computeTime*1e-9)); + jniClose(); + } + + public void parseSandboxFile(String filename) + { + File file = new File(filename); + Scanner input = null; + try + { + input = new Scanner(file); + } + catch(FileNotFoundException e) + { + System.err.println("File "+filename+" cannot be found/read"); + return; + } + int idx = 0; + int numReads = 0; + int numHaplotypes = 0; + int readIdx = 0, testCaseIdx = 0, haplotypeIdx = 0; + LinkedList haplotypeList = new LinkedList(); + LinkedList readList = new LinkedList(); + + byte[][] byteArray = new byte[6][]; + boolean firstLine = true; + String[] currTokens = new String[8]; + while(input.hasNextLine()) + { + String line = input.nextLine(); + Scanner lineScanner = new Scanner(line); + idx = 0; + while(lineScanner.hasNext()) + currTokens[idx++] = lineScanner.next(); + if(idx == 0) + break; + assert(idx >= 6); + //start of new region + if(idx == 8) + { + if(!firstLine) + { + initialize(haplotypeList); + computeLikelihoods(readList, haplotypeList); + finalizeRegion(); + } + try + { + numReads = Integer.parseInt(currTokens[6]); + } + catch(NumberFormatException e) + { + numReads = 1; + } + try + { + numHaplotypes = Integer.parseInt(currTokens[7]); + } + catch(NumberFormatException e) + { + numHaplotypes = 1; + } + haplotypeIdx = readIdx = testCaseIdx = 0; + readList.clear(); + haplotypeList.clear(); + } + if(haplotypeIdx < numHaplotypes) + { + JNIHaplotypeDataHolderClass X = new JNIHaplotypeDataHolderClass(); + X.haplotypeBases = currTokens[0].getBytes(); + haplotypeList.add(X); + } + if(testCaseIdx%numHaplotypes == 0) + { + JNIReadDataHolderClass X = new JNIReadDataHolderClass(); + X.readBases = currTokens[1].getBytes(); + for(int i=2;i<6;++i) + { + byteArray[i] = currTokens[i].getBytes(); + for(int j=0;j 0 && readList.size() > 0) + { + initialize(haplotypeList); + computeLikelihoods(readList, haplotypeList); + finalizeRegion(); + } + + close(); + input.close(); + } + + private native void doEverythingNative(String filename); + + public static void main(String[] args) + { + if(args.length <= 0) + { + System.err.println("Needs 1 argument - "); + System.exit(-1); + } + //// Get runtime + //java.lang.Runtime rt = java.lang.Runtime.getRuntime(); + //// Start a new process: UNIX command ls + //String cmd = "/home/karthikg/broad/gsa-unstable/public/c++/VectorPairHMM/checker "+args[0]; + //try + //{ + //System.out.println(cmd); + //java.lang.Process p = rt.exec(cmd); + //try + //{ + //p.waitFor(); + //java.io.InputStream is = p.getInputStream(); + //java.io.BufferedReader reader = new java.io.BufferedReader(new InputStreamReader(is)); + //// And print each line + //String s = null; + //while ((s = reader.readLine()) != null) { + //System.out.println(s); + //} + //is.close(); + //} + //catch(InterruptedException e) + //{ + //System.err.println(e); + //} + //} + //catch(IOException e) + //{ + //System.err.println(e); + //} + Sandbox t = new Sandbox(); + t.doEverythingNative(args[0]); + //t.parseSandboxFile(args[0]); + } +} diff --git a/public/c++/VectorPairHMM/Sandbox_JNIHaplotypeDataHolderClass.h b/public/c++/VectorPairHMM/Sandbox_JNIHaplotypeDataHolderClass.h new file mode 100644 index 000000000..7f78f0178 --- /dev/null +++ b/public/c++/VectorPairHMM/Sandbox_JNIHaplotypeDataHolderClass.h @@ -0,0 +1,13 @@ +/* DO NOT EDIT THIS FILE - it is machine generated */ +#include +/* Header for class Sandbox_JNIHaplotypeDataHolderClass */ + +#ifndef _Included_Sandbox_JNIHaplotypeDataHolderClass +#define _Included_Sandbox_JNIHaplotypeDataHolderClass +#ifdef __cplusplus +extern "C" { +#endif +#ifdef __cplusplus +} +#endif +#endif diff --git a/public/c++/VectorPairHMM/Sandbox_JNIReadDataHolderClass.h b/public/c++/VectorPairHMM/Sandbox_JNIReadDataHolderClass.h new file mode 100644 index 000000000..a9312ff3b --- /dev/null +++ b/public/c++/VectorPairHMM/Sandbox_JNIReadDataHolderClass.h @@ -0,0 +1,13 @@ +/* DO NOT EDIT THIS FILE - it is machine generated */ +#include +/* Header for class Sandbox_JNIReadDataHolderClass */ + +#ifndef _Included_Sandbox_JNIReadDataHolderClass +#define _Included_Sandbox_JNIReadDataHolderClass +#ifdef __cplusplus +extern "C" { +#endif +#ifdef __cplusplus +} +#endif +#endif diff --git a/public/c++/VectorPairHMM/baseline.cc b/public/c++/VectorPairHMM/baseline.cc index 2f80acdb0..d92a21ecf 100644 --- a/public/c++/VectorPairHMM/baseline.cc +++ b/public/c++/VectorPairHMM/baseline.cc @@ -1,28 +1,28 @@ #include "headers.h" #include "template.h" - +extern uint64_t exceptions_array[128]; template NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) { - int r, c; - int ROWS = tc->rslen + 1; - int COLS = tc->haplen + 1; + int r, c; + int ROWS = tc->rslen + 1; + int COLS = tc->haplen + 1; - Context ctx; + Context ctx; - NUMBER M[ROWS][COLS]; - NUMBER X[ROWS][COLS]; - NUMBER Y[ROWS][COLS]; - NUMBER p[ROWS][6]; + NUMBER M[ROWS][COLS]; + NUMBER X[ROWS][COLS]; + NUMBER Y[ROWS][COLS]; + NUMBER p[ROWS][6]; - p[0][MM] = ctx._(0.0); - p[0][GapM] = ctx._(0.0); - p[0][MX] = ctx._(0.0); - p[0][XX] = ctx._(0.0); - p[0][MY] = ctx._(0.0); - p[0][YY] = ctx._(0.0); - for (r = 1; r < ROWS; r++) - { + p[0][MM] = ctx._(0.0); + p[0][GapM] = ctx._(0.0); + p[0][MX] = ctx._(0.0); + p[0][XX] = ctx._(0.0); + p[0][MY] = ctx._(0.0); + p[0][YY] = ctx._(0.0); + for (r = 1; r < ROWS; r++) + { int _i = tc->i[r-1] & 127; int _d = tc->d[r-1] & 127; int _c = tc->c[r-1] & 127; @@ -34,50 +34,62 @@ NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) p[r][YY] = ctx.ph2pr[_c]; //p[r][MY] = (r == ROWS - 1) ? ctx._(1.0) : ctx.ph2pr[_d]; //p[r][YY] = (r == ROWS - 1) ? ctx._(1.0) : ctx.ph2pr[_c]; - } + } - for (c = 0; c < COLS; c++) - { + for (c = 0; c < COLS; c++) + { M[0][c] = ctx._(0.0); X[0][c] = ctx._(0.0); Y[0][c] = ctx.INITIAL_CONSTANT / (tc->haplen); - } + } - for (r = 1; r < ROWS; r++) - { + for (r = 1; r < ROWS; r++) + { M[r][0] = ctx._(0.0); X[r][0] = X[r-1][0] * p[r][XX]; Y[r][0] = ctx._(0.0); - } + } - NUMBER result = ctx._(0.0); + NUMBER result = ctx._(0.0); - for (r = 1; r < ROWS; r++) + for (r = 1; r < ROWS; r++) for (c = 1; c < COLS; c++) { - char _rs = tc->rs[r-1]; - char _hap = tc->hap[c-1]; - int _q = tc->q[r-1] & 127; - NUMBER distm = ctx.ph2pr[_q]; - if (_rs == _hap || _rs == 'N' || _hap == 'N') + fexcept_t flagp; + char _rs = tc->rs[r-1]; + char _hap = tc->hap[c-1]; + int _q = tc->q[r-1] & 127; + NUMBER distm = ctx.ph2pr[_q]; + if (_rs == _hap || _rs == 'N' || _hap == 'N') distm = ctx._(1.0) - distm; - else - distm = distm/3; - M[r][c] = distm * (M[r-1][c-1] * p[r][MM] + X[r-1][c-1] * p[r][GapM] + Y[r-1][c-1] * p[r][GapM]); - X[r][c] = M[r-1][c] * p[r][MX] + X[r-1][c] * p[r][XX]; - Y[r][c] = M[r][c-1] * p[r][MY] + Y[r][c-1] * p[r][YY]; + else + distm = distm/3; + + + //feclearexcept(FE_ALL_EXCEPT); + M[r][c] = distm * (M[r-1][c-1] * p[r][MM] + X[r-1][c-1] * p[r][GapM] + Y[r-1][c-1] * p[r][GapM]); + //M[r][c] = (M[r-1][c-1] * p[r][MM] + X[r-1][c-1] * p[r][GapM] + Y[r-1][c-1] * p[r][GapM]); + //STORE_FP_EXCEPTIONS(flagp, exceptions_array); + + //feclearexcept(FE_ALL_EXCEPT); + X[r][c] = M[r-1][c] * p[r][MX] + X[r-1][c] * p[r][XX]; + //STORE_FP_EXCEPTIONS(flagp, exceptions_array); + + //feclearexcept(FE_ALL_EXCEPT); + Y[r][c] = M[r][c-1] * p[r][MY] + Y[r][c-1] * p[r][YY]; + //STORE_FP_EXCEPTIONS(flagp, exceptions_array); } - for (c = 0; c < COLS; c++) - { + for (c = 0; c < COLS; c++) + { result += M[ROWS-1][c] + X[ROWS-1][c]; - } + } - if (before_last_log != NULL) + if (before_last_log != NULL) *before_last_log = result; - return result; - //return ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT; + return result; + //return ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT; } template double compute_full_prob(testcase* tc, double* nextbuf); diff --git a/public/c++/VectorPairHMM/headers.h b/public/c++/VectorPairHMM/headers.h index f502ff5ce..feb687660 100644 --- a/public/c++/VectorPairHMM/headers.h +++ b/public/c++/VectorPairHMM/headers.h @@ -1,3 +1,6 @@ +#ifndef COMMON_HEADERS_H +#define COMMON_HEADERS_H + #include #include #include @@ -22,4 +25,16 @@ #include #include #include +#include +#define STORE_FP_EXCEPTIONS(flagp, exceptions_array) \ + fegetexceptflag(&flagp, FE_OVERFLOW | FE_UNDERFLOW | FE_DIVBYZERO | FE_INVALID | __FE_DENORM); \ + exceptions_array[FE_INVALID] += ((flagp & FE_INVALID)); \ + exceptions_array[__FE_DENORM] += ((flagp & __FE_DENORM) >> 1); \ + exceptions_array[FE_DIVBYZERO] += ((flagp & FE_DIVBYZERO) >> 2); \ + exceptions_array[FE_OVERFLOW] += ((flagp & FE_OVERFLOW) >> 3); \ + exceptions_array[FE_UNDERFLOW] += ((flagp & FE_UNDERFLOW) >> 4); \ + feclearexcept(FE_ALL_EXCEPT); + + +#endif diff --git a/public/c++/VectorPairHMM/jni_common.h b/public/c++/VectorPairHMM/jni_common.h index 1e23f66ee..4c63a0411 100644 --- a/public/c++/VectorPairHMM/jni_common.h +++ b/public/c++/VectorPairHMM/jni_common.h @@ -2,9 +2,9 @@ #define JNI_COMMON_H #include -#define ENABLE_ASSERTIONS 1 +/*#define ENABLE_ASSERTIONS 1*/ #define DO_PROFILING 1 -//#define DEBUG 1 +/*#define DEBUG 1*/ //#define DEBUG0_1 1 //#define DEBUG3 1 /*#define DUMP_TO_SANDBOX 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 e31eba232..c39d2ec3f 100644 --- a/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc +++ b/public/c++/VectorPairHMM/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc @@ -5,16 +5,6 @@ #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 @@ -56,17 +46,6 @@ 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, @@ -211,44 +190,18 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless jbyte* haplotypeBasesArray = haplotypeBasesArrayVector[j].second; tc_array[tc_idx].rslen = (int)readLength; tc_array[tc_idx].haplen = (int)haplotypeLength; -#if 0 tc_array[tc_idx].hap = (char*)haplotypeBasesArray; tc_array[tc_idx].rs = (char*)readBasesArray; tc_array[tc_idx].q = (char*)readQualsArray; tc_array[tc_idx].i = (char*)insertionGOPArray; tc_array[tc_idx].d = (char*)deletionGOPArray; tc_array[tc_idx].c = (char*)overallGCPArray; -#endif - //#define MEMCPY_HACK -#ifdef MEMCPY_HACK - tc_array[tc_idx].hap = new char[haplotypeLength]; - tc_array[tc_idx].rs = new char[readLength]; - tc_array[tc_idx].q = new char[readLength]; - tc_array[tc_idx].i = new char[readLength]; - tc_array[tc_idx].d = new char[readLength]; - tc_array[tc_idx].c = new char[readLength]; - memcpy(tc_array[tc_idx].hap, haplotypeBasesArray, haplotypeLength); - memcpy(tc_array[tc_idx].rs, readBasesArray, readLength); - memcpy(tc_array[tc_idx].q, readQualsArray, readLength); - memcpy(tc_array[tc_idx].i, insertionGOPArray, readLength); - memcpy(tc_array[tc_idx].d, deletionGOPArray, readLength); - memcpy(tc_array[tc_idx].c, overallGCPArray, readLength); -#endif -#if 0 - tc_array[tc_idx].hap = (char*)all_arrays[0]; - tc_array[tc_idx].rs = (char*)all_arrays[1]; - tc_array[tc_idx].q = (char*)all_arrays[2]; - tc_array[tc_idx].i = (char*)all_arrays[3]; - tc_array[tc_idx].d = (char*)all_arrays[4]; - tc_array[tc_idx].c = (char*)all_arrays[5]; -#endif - -#ifdef DUMP_TO_SANDBOX - g_load_time_initializer.dump_sandbox(tc_array[tc_idx]); -#endif #ifdef DO_PROFILING g_load_time_initializer.m_sumProductReadLengthHaplotypeLength += (readLength*haplotypeLength); g_load_time_initializer.m_sumSquareProductReadLengthHaplotypeLength += ((readLength*haplotypeLength)*(readLength*haplotypeLength)); +#endif +#ifdef DUMP_TO_SANDBOX + g_load_time_initializer.dump_sandbox(tc_array[tc_idx], tc_idx, numReads, numHaplotypes); #endif ++tc_idx; } @@ -265,9 +218,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless g_load_time_initializer.m_sumReadLengths += readLength; #endif } -#ifdef DUMP_TO_SANDBOX - g_load_time_initializer.close_sandbox(); -#endif + #ifdef DO_PROFILING g_load_time_initializer.m_data_transfer_time += get_time(); #endif @@ -296,16 +247,6 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless } else result = (double)(log10f(result_avxf) - log10f(ldexpf(1.f, 120.f))); -#if 0 - double result = 0; - testcase& tc = tc_array[tc_idx]; - for(unsigned k=0;k=0;--i)//note the order - reverse of GET { @@ -360,6 +289,9 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLogless #ifdef DEBUG g_load_time_initializer.debug_close(); #endif +#ifdef DUMP_TO_SANDBOX + g_load_time_initializer.close_sandbox(); +#endif } //Release haplotypes at the end of a region diff --git a/public/c++/VectorPairHMM/pairhmm-1-base.cc b/public/c++/VectorPairHMM/pairhmm-1-base.cc index 811a67b8c..c7a85e647 100644 --- a/public/c++/VectorPairHMM/pairhmm-1-base.cc +++ b/public/c++/VectorPairHMM/pairhmm-1-base.cc @@ -2,16 +2,11 @@ //#define DEBUG0_1 1 //#define DEBUG3 1 #include "headers.h" -#include "template.h" #include "utils.h" #include "LoadTimeInitializer.h" - using namespace std; - #define RUN_HYBRID -vector results_vec; -vector tc_vector; int main(int argc, char** argv) { #define BATCH_SIZE 5 @@ -20,6 +15,8 @@ int main(int argc, char** argv) cerr << "Needs path to input file as argument\n"; exit(0); } + do_compute(argv[1]); +#if 0 bool use_old_read_testcase = false; if(argc >= 3 && string(argv[2]) == "1") use_old_read_testcase = true; @@ -39,7 +36,8 @@ int main(int argc, char** argv) ifptr.open(argv[1]); assert(ifptr.is_open()); } - + vector results_vec; + vector tc_vector; tc_vector.clear(); tc_vector.resize(BATCH_SIZE+4); results_vec.clear(); @@ -58,6 +56,7 @@ int main(int argc, char** argv) testcase tc_in; int break_value = 0; tc_vector.clear(); + g_load_time_initializer.open_sandbox(); while(1) { break_value = use_old_read_testcase ? read_testcase(&tc_in, fptr) : @@ -137,7 +136,7 @@ int main(int argc, char** argv) tc_vector.clear(); if(all_ok) cout << "All outputs acceptable\n"; - cout << "Total vector time "<< ((double)total_time)/1e9 << " baseline time "<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)); + //tc->ihap = (int *) malloc(tc->haplen*sizeof(int)); + //tc->irs = (int *) malloc(tc->rslen*sizeof(int)); tc->q = (char *) malloc(sizeof(char) * tc->rslen); tc->i = (char *) malloc(sizeof(char) * tc->rslen); @@ -121,10 +121,10 @@ int read_testcase(testcase *tc, FILE* ifp) tc->i[x] = _i; tc->d[x] = _d; tc->c[x] = _c; - tc->irs[x] = tc->rs[x]; + //tc->irs[x] = tc->rs[x]; } - for (x = 0; x < tc->haplen; x++) - tc->ihap[x] = tc->hap[x]; + //for (x = 0; x < tc->haplen; x++) + //tc->ihap[x] = tc->hap[x]; free(q); @@ -286,4 +286,182 @@ uint64_t diff_time(struct timespec& prev_time) return (uint64_t)((curr_time.tv_sec-prev_time.tv_sec)*1000000000+(curr_time.tv_nsec-prev_time.tv_nsec)); } +//#define USE_PAPI +#ifdef USE_PAPI +#include "papi.h" +#define NUM_PAPI_COUNTERS 4 +#endif +uint64_t exceptions_array[128]; +void do_compute(char* filename) +{ + memset(exceptions_array, 0, 128*sizeof(uint64_t)); + _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); + //assert(feenableexcept(FE_DIVBYZERO | FE_INVALID) >= 0); +#ifdef USE_PAPI + PAPI_num_counters(); + //int events[NUM_PAPI_COUNTERS] = { PAPI_TOT_INS, PAPI_TOT_CYC, PAPI_L1_DCM, PAPI_L1_ICM, PAPI_L3_TCM, PAPI_TLB_DM, PAPI_TLB_IM }; + //char* eventnames[NUM_PAPI_COUNTERS]= { "instructions", "cycles", "l1d_misses", "l1i_misses", "l3_misses", "dtlb_misses", "itlb_misses" }; + //long long values[NUM_PAPI_COUNTERS] = { 0, 0, 0, 0, 0, 0, 0 }; + //long long accum_values[NUM_PAPI_COUNTERS] = { 0, 0, 0, 0, 0, 0, 0 }; + //int events[NUM_PAPI_COUNTERS] = { PAPI_TOT_INS, PAPI_TOT_CYC, PAPI_L1_ICM }; + //char* eventnames[NUM_PAPI_COUNTERS]= { "instructions", "cycles", "l1i_misses"}; + //assert(PAPI_event_name_to_code("PERF_COUNT_HW_STALLED_CYCLES_FRONTEND",&(events[2])) == PAPI_OK); + int events[NUM_PAPI_COUNTERS] = { 0, 0, 0, 0 }; + //assert(PAPI_event_name_to_code("ICACHE:IFETCH_STALL",&(events[2])) == PAPI_OK); + //assert(PAPI_event_name_to_code("MACHINE_CLEARS:e",&(events[3])) == PAPI_OK); + char* eventnames[NUM_PAPI_COUNTERS]= { "instructions", "cycles", "ifetch_stall", "store_misses" }; + assert(PAPI_event_name_to_code("ix86arch::INSTRUCTION_RETIRED",&(events[0])) == PAPI_OK); + assert(PAPI_event_name_to_code("UNHALTED_REFERENCE_CYCLES",&(events[1])) == PAPI_OK); + assert(PAPI_event_name_to_code("ICACHE:IFETCH_STALL", &(events[2])) == PAPI_OK); + assert(PAPI_event_name_to_code("perf::L1-DCACHE-STORE-MISSES", &(events[3])) == PAPI_OK); + long long values[NUM_PAPI_COUNTERS] = { 0, 0, 0, 0 }; + long long accum_values[NUM_PAPI_COUNTERS] = { 0, 0, 0, 0 }; + +#endif +#define BATCH_SIZE 100000 + bool use_old_read_testcase = true; + unsigned chunk_size = 100; + std::ifstream ifptr; + FILE* fptr = 0; + if(use_old_read_testcase) + { + fptr = fopen(filename,"r"); + if(fptr == 0) + cerr << "Could not open file "< tc_vector; + tc_vector.clear(); + vector results_vec; + results_vec.clear(); + vector baseline_results; + baseline_results.clear(); + + 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; + + testcase tc_in; + int break_value = 0; + while(1) + { + break_value = use_old_read_testcase ? read_testcase(&tc_in, fptr) : + read_mod_testcase(ifptr, &tc_in, true); + tc_vector.push_back(tc_in); + if(break_value >= 0) + ++num_testcases; + if(num_testcases == curr_batch_size || (break_value < 0 && num_testcases > 0)) + { + results_vec.resize(tc_vector.size()); + baseline_results.resize(tc_vector.size()); + + get_time(); +#ifdef USE_PAPI + assert(PAPI_start_counters(events, NUM_PAPI_COUNTERS) == PAPI_OK); +#endif +#pragma omp parallel for schedule(dynamic,chunk_size) num_threads(12) + for(unsigned i=0;i(&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 << "Line "<