From e90405cd1fd1f8d824194f91f971aa32da667cc7 Mon Sep 17 00:00:00 2001 From: Karthik Gururaj Date: Thu, 16 Jan 2014 19:53:50 -0800 Subject: [PATCH] 1. Nested loops over reads and haplotypes moved to C++ through JNI 2. OpenMP support added 3. Using direct access to Java primitive arrays 4. Debug messages disabled --- .gitignore | 2 + Makefile | 17 +- convert_char.cc | 186 ++++++++ hmm_mask.cc | 60 +-- ...e_sting_utils_pairhmm_JNILoglessPairHMM.cc | 417 +++++++++++++++--- ...te_sting_utils_pairhmm_JNILoglessPairHMM.h | 15 + pairhmm-1-base.cc | 324 ++++++-------- pairhmm-template-main.cc | 8 +- template.h | 79 +--- 9 files changed, 737 insertions(+), 371 deletions(-) create mode 100644 convert_char.cc diff --git a/.gitignore b/.gitignore index dcc44daa4..cb70b53a7 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ tests hmm_Mohammad pairhmm-template-main *.swp +checker +reformat diff --git a/Makefile b/Makefile index f858a7241..70bf96dd3 100644 --- a/Makefile +++ b/Makefile @@ -16,11 +16,11 @@ CXX=icc LDFLAGS=-lm $(OMPLDFLAGS) #BIN:=pairhmm-1-base #pairhmm-2-omp pairhmm-3-hybrid-float-double pairhmm-4-hybrid-diagonal pairhmm-5-hybrid-diagonal-homogeneus pairhmm-6-onlythreediags pairhmm-7-presse pairhmm-8-sse #pairhmm-dev -BIN:=libJNILoglessPairHMM.so pairhmm-template-main +BIN:=libJNILoglessPairHMM.so pairhmm-template-main checker #SOURCES=pairhmm-1-base.cc input.cc -LIBSOURCES=org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc hmm_mask.cc -SOURCES=$(LIBSOURCES) pairhmm-template-main.cc +LIBSOURCES=org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc convert_char.cc +SOURCES=$(LIBSOURCES) pairhmm-template-main.cc pairhmm-1-base.cc LIBOBJECTS=$(LIBSOURCES:.cc=.o) DEPDIR=.deps DF=$(DEPDIR)/$(*).d @@ -29,11 +29,14 @@ all: $(BIN) -include $(addprefix $(DEPDIR)/,$(SOURCES:.cc=.d)) -pairhmm-template-main: pairhmm-template-main.o hmm_mask.o - $(CXX) -fopenmp -o $@ $^ $(LDFLAGS) +checker: pairhmm-1-base.o convert_char.o + $(CXX) $(OMPCFLAGS) -o $@ $^ $(LDFLAGS) + +pairhmm-template-main: pairhmm-template-main.o convert_char.o + $(CXX) $(OMPCFLAGS) -o $@ $^ $(LDFLAGS) libJNILoglessPairHMM.so: $(LIBOBJECTS) - $(CXX) -shared -o $@ $(LIBOBJECTS) + $(CXX) $(OMPCFLAGS) -shared -o $@ $(LIBOBJECTS) %.o: %.cc @mkdir -p $(DEPDIR) @@ -41,4 +44,4 @@ libJNILoglessPairHMM.so: $(LIBOBJECTS) clean: - rm -f $(BIN) *.o + rm -rf $(BIN) *.o $(DEPDIR) diff --git a/convert_char.cc b/convert_char.cc new file mode 100644 index 000000000..01c5a5137 --- /dev/null +++ b/convert_char.cc @@ -0,0 +1,186 @@ +#include "template.h" +uint8_t ConvertChar::conversionTable[255]; +using namespace std; + +int normalize(char c) +{ + return ((int) (c - 33)); +} + +int read_testcase(testcase *tc, FILE* ifp) +{ + char *q, *i, *d, *c, *line = NULL; + int _q, _i, _d, _c; + int x, size = 0; + ssize_t read; + + read = getline(&line, (size_t *) &size, ifp == 0 ? stdin : ifp); + if (read == -1) + return -1; + + + tc->hap = (char *) malloc(size); + tc->rs = (char *) malloc(size); + q = (char *) malloc(size); + i = (char *) malloc(size); + d = (char *) malloc(size); + c = (char *) malloc(size); + + if (sscanf(line, "%s %s %s %s %s %s\n", tc->hap, tc->rs, q, i, d, c) != 6) + return -1; + + tc->haplen = strlen(tc->hap); + tc->rslen = strlen(tc->rs); + assert(tc->rslen < MROWS); + tc->ihap = (int *) malloc(tc->haplen*sizeof(int)); + tc->irs = (int *) malloc(tc->rslen*sizeof(int)); + + //tc->q = (int *) malloc(sizeof(int) * tc->rslen); + //tc->i = (int *) malloc(sizeof(int) * tc->rslen); + //tc->d = (int *) malloc(sizeof(int) * tc->rslen); + //tc->c = (int *) malloc(sizeof(int) * tc->rslen); + + for (x = 0; x < tc->rslen; x++) + { + _q = normalize(q[x]); + _i = normalize(i[x]); + _d = normalize(d[x]); + _c = normalize(c[x]); + tc->q[x] = (_q < 6) ? 6 : _q; + tc->i[x] = _i; + tc->d[x] = _d; + tc->c[x] = _c; + tc->irs[x] = tc->rs[x]; + } + 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; +} + +unsigned MAX_LINE_LENGTH = 65536; +int convToInt(std::string s) +{ + int i; + std::istringstream strin(s); + strin >> i; + return i; +} + +void tokenize(std::ifstream& fptr, std::vector& tokens) +{ + int i = 0; + std::string tmp; + std::vector myVec; + vector line; + line.clear(); + line.resize(MAX_LINE_LENGTH); + vector tmpline; + tmpline.clear(); + tmpline.resize(MAX_LINE_LENGTH); + myVec.clear(); + + while(!fptr.eof()) + { + i = 0; + bool still_read_line = true; + unsigned line_position = 0; + while(still_read_line) + { + fptr.getline(&(tmpline[0]), MAX_LINE_LENGTH); + if(line_position + MAX_LINE_LENGTH > line.size()) + line.resize(2*line.size()); + for(unsigned i=0;i> std::skipws >> tmp; + if(tmp != "") + { + myVec.push_back(tmp); + ++i; + //std::cout < 0) + break; + } + tokens.clear(); + //std::cout << "Why "< tokens; + tokens.clear(); + tokenize(fptr, tokens); + if(tokens.size() == 0) + return -1; + tc->hap = new char[tokens[0].size()+2]; + tc->haplen = tokens[0].size(); + memcpy(tc->hap, tokens[0].c_str(), tokens[0].size()); + tc->rs = new char[tokens[1].size()+2]; + tc->rslen = tokens[1].size(); + //cout << "Lengths "<haplen <<" "<rslen<<"\n"; + memcpy(tc->rs, tokens[1].c_str(),tokens[1].size()); + assert(tokens.size() == 2 + 4*(tc->rslen)); + assert(tc->rslen < MROWS); + for(unsigned j=0;jrslen;++j) + tc->q[j] = convToInt(tokens[2+0*tc->rslen+j]); + for(unsigned j=0;jrslen;++j) + tc->i[j] = convToInt(tokens[2+1*tc->rslen+j]); + for(unsigned j=0;jrslen;++j) + tc->d[j] = convToInt(tokens[2+2*tc->rslen+j]); + for(unsigned j=0;jrslen;++j) + tc->c[j] = convToInt(tokens[2+3*tc->rslen+j]); + + if(reformat) + { + ofstream ofptr; + ofptr.open("reformat/debug_dump.txt",first_call ? ios::out : ios::app); + assert(ofptr.is_open()); + ofptr << tokens[0] << " "; + ofptr << tokens[1] << " "; + for(unsigned j=0;jrslen;++j) + ofptr << ((char)(tc->q[j]+33)); + ofptr << " "; + for(unsigned j=0;jrslen;++j) + ofptr << ((char)(tc->i[j]+33)); + ofptr << " "; + for(unsigned j=0;jrslen;++j) + ofptr << ((char)(tc->d[j]+33)); + ofptr << " "; + for(unsigned j=0;jrslen;++j) + ofptr << ((char)(tc->c[j]+33)); + ofptr << " 0 false\n"; + + ofptr.close(); + first_call = false; + } + + + return tokens.size(); +} diff --git a/hmm_mask.cc b/hmm_mask.cc index f128a39d9..c335a0160 100644 --- a/hmm_mask.cc +++ b/hmm_mask.cc @@ -30,9 +30,7 @@ #include #include #include -#include -#include -#include +#include "template.h" using namespace std ; @@ -50,43 +48,14 @@ string getBinaryStr (T val, int numBitsToWrite) { return oss.str() ; } - -typedef struct -{ - int rslen, haplen, *q, *i, *d, *c; - char *hap, *rs; -} testcase; - int normalize(char c) { return ((int) (c - 33)); } -class ConvertChar { - - static uint8_t conversionTable[255] ; - - -public: - - static void init() { - conversionTable['A'] = 0 ; - conversionTable['C'] = 1 ; - conversionTable['T'] = 2 ; - conversionTable['G'] = 3 ; - conversionTable['N'] = 4 ; - } - - static inline uint8_t get(uint8_t input) { - return conversionTable[input] ; - } - -} ; uint8_t ConvertChar::conversionTable[255] ; - - int read_testcase(testcase *tc, FILE* ifp) { char *q, *i, *d, *c, *line = NULL; @@ -94,7 +63,7 @@ int read_testcase(testcase *tc, FILE* ifp) int x, size = 0; ssize_t read; - read = getline(&line, (size_t *) &size, ifp); + read = getline(&line, (size_t *) &size, ifp == 0 ? stdin : ifp); if (read == -1) return -1; @@ -111,20 +80,14 @@ int read_testcase(testcase *tc, FILE* ifp) tc->haplen = strlen(tc->hap); tc->rslen = strlen(tc->rs); - tc->q = (int *) malloc(sizeof(int) * tc->rslen); - tc->i = (int *) malloc(sizeof(int) * tc->rslen); - tc->d = (int *) malloc(sizeof(int) * tc->rslen); - tc->c = (int *) malloc(sizeof(int) * tc->rslen); + assert(tc->rslen < MROWS); + tc->ihap = (int *) malloc(tc->haplen*sizeof(int)); + tc->irs = (int *) malloc(tc->rslen*sizeof(int)); - // Convert hap and rs to 3 bits - - /* - for (int ci=0; ci < tc->haplen; ++ci) - tc->hap[ci] = ConvertChar::get(tc->hap[ci]) ; - - for (int ci=0; ci < tc->rslen; ++ci) - tc->rs[ci] = ConvertChar::get(tc->rs[ci]) ; - */ + //tc->q = (int *) malloc(sizeof(int) * tc->rslen); + //tc->i = (int *) malloc(sizeof(int) * tc->rslen); + //tc->d = (int *) malloc(sizeof(int) * tc->rslen); + //tc->c = (int *) malloc(sizeof(int) * tc->rslen); for (x = 0; x < tc->rslen; x++) { @@ -136,8 +99,11 @@ 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]; } - + for (x = 0; x < tc->haplen; x++) + tc->ihap[x] = tc->hap[x]; + free(q); free(i); diff --git a/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc b/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc index 29668a306..b8298161e 100644 --- a/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc +++ b/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc @@ -2,6 +2,7 @@ #include #include #include "org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h" +#include #include #include @@ -12,12 +13,14 @@ #include #include + #include #include #include using namespace std; +//#define DEBUG 1 +//#define DEBUG0_1 1 //#define DEBUG3 1 -#define DEBUG 1 #include "template.h" @@ -42,6 +45,12 @@ class LoadTimeInitializer { ConvertChar::init(); } + jfieldID m_readBasesFID; + jfieldID m_readQualsFID; + jfieldID m_insertionGOPFID; + jfieldID m_deletionGOPFID; + jfieldID m_overallGCPFID; + jfieldID m_haplotypeBasesFID; }; LoadTimeInitializer g_load_time_initializer; @@ -55,73 +64,6 @@ void debug_dump(string filename, string s, bool to_append, bool add_newline) fptr.close(); } -#define INT_STORE_ARRAY_SIZE 2048 -int insertionGOPIntArray[INT_STORE_ARRAY_SIZE]; -int deletionGOPIntArray[INT_STORE_ARRAY_SIZE]; -int overallGCPIntArray[INT_STORE_ARRAY_SIZE]; -int readQualsIntArray[INT_STORE_ARRAY_SIZE]; - -JNIEXPORT jdouble JNICALL -Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniSubComputeReadLikelihoodGivenHaplotypeLog10( - JNIEnv* env, jobject thisObject, - jint readLength, jint haplotypeLength, - jbyteArray readBases, jbyteArray haplotypeBases, jbyteArray readQuals, - jbyteArray insertionGOP, jbyteArray deletionGOP, jbyteArray overallGCP, - jint hapStartIndex - ) -{ - jboolean is_copy = JNI_FALSE; - jbyte* readBasesArray = (env)->GetByteArrayElements(readBases, &is_copy); - jbyte* haplotypeBasesArray = (env)->GetByteArrayElements(haplotypeBases, &is_copy); - jbyte* readQualsArray = (env)->GetByteArrayElements(readQuals, &is_copy); - jbyte* insertionGOPArray = (env)->GetByteArrayElements(insertionGOP, &is_copy); - jbyte* deletionGOPArray = (env)->GetByteArrayElements(deletionGOP, &is_copy); - jbyte* overallGCPArray = (env)->GetByteArrayElements(overallGCP, &is_copy); -#ifdef DEBUG - assert(readBasesArray && "readBasesArray not initialized in JNI"); - assert(haplotypeBasesArray && "haplotypeBasesArray not initialized in JNI"); - assert(readQualsArray && "readQualsArray not initialized in JNI"); - assert(insertionGOPArray && "insertionGOP array not initialized in JNI"); - assert(deletionGOPArray && "deletionGOP array not initialized in JNI"); - assert(overallGCPArray && "OverallGCP array not initialized in JNI"); - assert(readLength < INT_STORE_ARRAY_SIZE); -#endif - for(unsigned i=0;i(&tc); - double result = log10(result_avxd) - log10(ldexp(1.0, 1020)); -#ifdef DEBUG - debug_dump("return_values_jni.txt",to_string(result),true); -#endif - - - env->ReleaseByteArrayElements(overallGCP, overallGCPArray, JNI_ABORT); - env->ReleaseByteArrayElements(deletionGOP, deletionGOPArray, JNI_ABORT); - env->ReleaseByteArrayElements(insertionGOP, insertionGOPArray, JNI_ABORT); - env->ReleaseByteArrayElements(readQuals, readQualsArray, JNI_ABORT); - env->ReleaseByteArrayElements(haplotypeBases, haplotypeBasesArray, JNI_ABORT); - env->ReleaseByteArrayElements(readBases, readBasesArray, JNI_ABORT); - - return 0.0; -} - JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializeProbabilities (JNIEnv* env, jclass thisObject, jobjectArray transition, jbyteArray insertionGOP, jbyteArray deletionGOP, jbyteArray overallGCP @@ -142,3 +84,342 @@ Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializePrior ) { return 0.0; } + +#define DIRECT_ACCESS_TO_JAVA_HEAP_MEMORY 1 + +#ifdef DIRECT_ACCESS_TO_JAVA_HEAP_MEMORY +//Gets direct access to Java arrays +#define GET_BYTE_ARRAY_ELEMENTS env->GetPrimitiveArrayCritical +#define RELEASE_BYTE_ARRAY_ELEMENTS env->ReleasePrimitiveArrayCritical +#define JNI_RELEASE_MODE JNI_ABORT +#define GET_DOUBLE_ARRAY_ELEMENTS env->GetPrimitiveArrayCritical +#define RELEASE_DOUBLE_ARRAY_ELEMENTS env->ReleasePrimitiveArrayCritical + +#else +//Likely makes copy of Java arrays to JNI C++ space +#define GET_BYTE_ARRAY_ELEMENTS env->GetByteArrayElements +#define RELEASE_BYTE_ARRAY_ELEMENTS env->ReleaseByteArrayElements +#define JNI_RELEASE_MODE JNI_ABORT +#define GET_DOUBLE_ARRAY_ELEMENTS env->GetDoubleArrayElements +#define RELEASE_DOUBLE_ARRAY_ELEMENTS env->ReleaseDoubleArrayElements + +#endif //ifdef DIRECT_ACCESS_TO_JAVA_HEAP_MEMORY + +JNIEXPORT jdouble JNICALL +Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniSubComputeReadLikelihoodGivenHaplotypeLog10( + JNIEnv* env, jobject thisObject, + jint readLength, jint haplotypeLength, + jbyteArray readBases, jbyteArray haplotypeBases, jbyteArray readQuals, + jbyteArray insertionGOP, jbyteArray deletionGOP, jbyteArray overallGCP, + jint hapStartIndex + ) +{ + jboolean is_copy = JNI_FALSE; + jbyte* readBasesArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(readBases, &is_copy); + jbyte* haplotypeBasesArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(haplotypeBases, &is_copy); + jbyte* readQualsArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(readQuals, &is_copy); + jbyte* insertionGOPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(insertionGOP, &is_copy); + jbyte* deletionGOPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(deletionGOP, &is_copy); + jbyte* overallGCPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(overallGCP, &is_copy); +#ifdef DEBUG + assert(readBasesArray && "readBasesArray not initialized in JNI"); + assert(haplotypeBasesArray && "haplotypeBasesArray not initialized in JNI"); + assert(readQualsArray && "readQualsArray not initialized in JNI"); + assert(insertionGOPArray && "insertionGOP array not initialized in JNI"); + assert(deletionGOPArray && "deletionGOP array not initialized in JNI"); + assert(overallGCPArray && "OverallGCP array not initialized in JNI"); + assert(readLength < MROWS); +#endif + testcase tc; + tc.rslen = readLength; + tc.haplen = haplotypeLength; + tc.rs = (char*)readBasesArray; + tc.hap = (char*)haplotypeBasesArray; + for(unsigned i=0;i(&tc); + double result = log10(result_avxd) - log10(ldexp(1.0, 1020)); +#ifdef DEBUG + debug_dump("return_values_jni.txt",to_string(result),true); +#endif + + + RELEASE_BYTE_ARRAY_ELEMENTS(overallGCP, overallGCPArray, JNI_RELEASE_MODE); + RELEASE_BYTE_ARRAY_ELEMENTS(deletionGOP, deletionGOPArray, JNI_RELEASE_MODE); + RELEASE_BYTE_ARRAY_ELEMENTS(insertionGOP, insertionGOPArray, JNI_RELEASE_MODE); + RELEASE_BYTE_ARRAY_ELEMENTS(readQuals, readQualsArray, JNI_RELEASE_MODE); + RELEASE_BYTE_ARRAY_ELEMENTS(haplotypeBases, haplotypeBasesArray, JNI_RELEASE_MODE); + RELEASE_BYTE_ARRAY_ELEMENTS(readBases, readBasesArray, JNI_RELEASE_MODE); + + return 0.0; +} + +//Should be called only once for the whole Java process - initializes field ids for the classes JNIReadDataHolderClass +//and JNIHaplotypeDataHolderClass +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniGlobalInit + (JNIEnv* env, jobject thisObject, jclass readDataHolderClass, jclass haplotypeDataHolderClass) +{ + assert(readDataHolderClass); + assert(haplotypeDataHolderClass); + jfieldID fid; + fid = env->GetFieldID(readDataHolderClass, "readBases", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for readBases"); + g_load_time_initializer.m_readBasesFID = fid; + fid = env->GetFieldID(readDataHolderClass, "readQuals", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for readQuals"); + g_load_time_initializer.m_readQualsFID = fid; + fid = env->GetFieldID(readDataHolderClass, "insertionGOP", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for insertionGOP"); + g_load_time_initializer.m_insertionGOPFID = fid; + fid = env->GetFieldID(readDataHolderClass, "deletionGOP", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for deletionGOP"); + g_load_time_initializer.m_deletionGOPFID = fid; + fid = env->GetFieldID(readDataHolderClass, "overallGCP", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for overallGCP"); + g_load_time_initializer.m_overallGCPFID = fid; + + fid = env->GetFieldID(haplotypeDataHolderClass, "haplotypeBases", "[B"); + assert(fid && "JNI pairHMM: Could not get FID for haplotypeBases"); + g_load_time_initializer.m_haplotypeBasesFID = fid; +} + + +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniComputeLikelihoods + (JNIEnv* env, jobject thisObject, jint numReads, jint numHaplotypes, + jobjectArray readDataArray, jobjectArray haplotypeDataArray, jdoubleArray likelihoodArray, jint maxNumThreadsToUse) +{ + jboolean is_copy = JNI_FALSE; + //To ensure, GET_BYTE_ARRAY_ELEMENTS is invoked only once for each haplotype, store bytearrays in a vector + vector > haplotypeBasesArrayVector; + haplotypeBasesArrayVector.clear(); + haplotypeBasesArrayVector.resize(numHaplotypes); +#ifdef DEBUG0_1 + cout << "JNI numReads "<GetObjectArrayElement(haplotypeDataArray, j); + jbyteArray haplotypeBases = (jbyteArray)env->GetObjectField(haplotypeObject, g_load_time_initializer.m_haplotypeBasesFID); +#ifdef DEBUG + assert(haplotypeBases && ("haplotypeBases is NULL at index : "+to_string(j)+"\n").c_str()); +#endif + jbyte* haplotypeBasesArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(haplotypeBases, &is_copy); +#ifdef DEBUG + assert(haplotypeBasesArray && "haplotypeBasesArray not initialized in JNI"); + assert(env->GetArrayLength(haplotypeBases) < MCOLS); +#ifdef DEBUG0_1 + cout << "JNI haplotype length "<GetArrayLength(haplotypeBases)<<"\n"; +#endif +#endif + haplotypeBasesArrayVector[j] = make_pair(haplotypeBases, haplotypeBasesArray); +#ifdef DEBUG3 + for(unsigned k=0;kGetArrayLength(haplotypeBases);++k) + debug_dump("haplotype_bases_jni.txt",to_string((int)haplotypeBasesArray[k]),true); +#endif + } + + unsigned numTestCases = numReads*numHaplotypes; + vector tc_array; + tc_array.clear(); + tc_array.resize(numTestCases); + unsigned tc_idx = 0; + //Store arrays for release later + vector > > readBasesArrayVector; + readBasesArrayVector.clear(); + readBasesArrayVector.resize(numReads); + for(unsigned i=0;iGetObjectArrayElement(readDataArray, i); + jbyteArray readBases = (jbyteArray)env->GetObjectField(readObject, g_load_time_initializer.m_readBasesFID); + jbyteArray insertionGOP = (jbyteArray)env->GetObjectField(readObject, g_load_time_initializer.m_insertionGOPFID); + jbyteArray deletionGOP = (jbyteArray)env->GetObjectField(readObject, g_load_time_initializer.m_deletionGOPFID); + jbyteArray overallGCP = (jbyteArray)env->GetObjectField(readObject, g_load_time_initializer.m_overallGCPFID); + jbyteArray readQuals = (jbyteArray)env->GetObjectField(readObject, g_load_time_initializer.m_readQualsFID); + +#ifdef DEBUG + assert(readBases && ("readBases is NULL at index : "+to_string(i)+"\n").c_str()); + assert(insertionGOP && ("insertionGOP is NULL at index : "+to_string(i)+"\n").c_str()); + assert(deletionGOP && ("deletionGOP is NULL at index : "+to_string(i)+"\n").c_str()); + assert(overallGCP && ("overallGCP is NULL at index : "+to_string(i)+"\n").c_str()); + assert(readQuals && ("readQuals is NULL at index : "+to_string(i)+"\n").c_str()); +#endif + jsize readLength = env->GetArrayLength(readBases); + + jbyte* readBasesArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(readBases, &is_copy); //order of GET-RELEASE is important + jbyte* readQualsArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(readQuals, &is_copy); + jbyte* insertionGOPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(insertionGOP, &is_copy); + jbyte* deletionGOPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(deletionGOP, &is_copy); + jbyte* overallGCPArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(overallGCP, &is_copy); +#ifdef DEBUG + assert(readBasesArray && "readBasesArray not initialized in JNI"); + assert(readQualsArray && "readQualsArray not initialized in JNI"); + assert(insertionGOPArray && "insertionGOP array not initialized in JNI"); + assert(deletionGOPArray && "deletionGOP array not initialized in JNI"); + assert(overallGCPArray && "overallGCP array not initialized in JNI"); + assert(readLength < MROWS); + assert(readLength == env->GetArrayLength(readQuals)); + assert(readLength == env->GetArrayLength(insertionGOP)); + assert(readLength == env->GetArrayLength(deletionGOP)); + assert(readLength == env->GetArrayLength(overallGCP)); +#ifdef DEBUG0_1 + cout << "JNI read length "<GetArrayLength(haplotypeBasesArrayVector[j].first); + jbyte* haplotypeBasesArray = haplotypeBasesArrayVector[j].second; + tc_array[tc_idx].rslen = (int)readLength; + tc_array[tc_idx].haplen = (int)haplotypeLength; + tc_array[tc_idx].rs = (char*)readBasesArray; + tc_array[tc_idx].hap = (char*)haplotypeBasesArray; + //Can be avoided + for(unsigned k=0;kGetArrayLength(likelihoodArray) == numTestCases); +#pragma omp parallel for schedule (guided,2) private(tc_idx) num_threads(maxNumThreadsToUse) + for(tc_idx=0;tc_idx(&(tc_array[tc_idx])); + double result = log10(result_avxd) - log10(ldexp(1.0, 1020)); + likelihoodDoubleArray[tc_idx] = result; + } +#ifdef DEBUG + for(tc_idx=0;tc_idx=0;--i)//note the order - reverse of GET + { + for(int j=readBasesArrayVector.size()-1;j>=0;--j) + RELEASE_BYTE_ARRAY_ELEMENTS(readBasesArrayVector[i][j].first, readBasesArrayVector[i][j].second, JNI_RELEASE_MODE); + readBasesArrayVector[i].clear(); + } + readBasesArrayVector.clear(); + + //Now release haplotype arrays + for(int j=numHaplotypes-1;j>=0;--j) //note the order - reverse of GET + RELEASE_BYTE_ARRAY_ELEMENTS(haplotypeBasesArrayVector[j].first, haplotypeBasesArrayVector[j].second, JNI_RELEASE_MODE); + haplotypeBasesArrayVector.clear(); + tc_array.clear(); +} + + + + + /* + protected final double computeReadLikelihoodGivenHaplotypeLog10( final byte[] haplotypeBases, + final byte[] readBases, + final byte[] readQuals, + final byte[] insertionGOP, + final byte[] deletionGOP, + final byte[] overallGCP, + final boolean recacheReadValues, + final byte[] nextHaploytpeBases) { + + paddedReadLength = readBases.length + 1; + paddedHaplotypeLength = haplotypeBases.length + 1; + double result = subComputeReadLikelihoodGivenHaplotypeLog10(haplotypeBases, readBases, readQuals, insertionGOP, deletionGOP, overallGCP, 0, true, 0); + if ( ! MathUtils.goodLog10Probability(result) ) + throw new IllegalStateException("PairHMM Log Probability cannot be greater than 0: " + String.format("haplotype: %s, read: %s, result: %f", Arrays.toString(haplotypeBases), Arrays.toString(readBases), result)); + // Warning: Careful if using the PairHMM in parallel! (this update has to be taken care of). + // Warning: This assumes no downstream modification of the haplotype bases (saves us from copying the array). It is okay for the haplotype caller and the Unified Genotyper. + // KG: seems pointless is not being used anywhere + previousHaplotypeBases = haplotypeBases; + return result; + return 0; + }*/ + +//public PerReadAlleleLikelihoodMap computeLikelihoods(final List reads, final Map alleleHaplotypeMap, final Map GCPArrayMap) + //{ + ////PerReadAlleleLikelihoodMap retValue = super.computeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap); + //// (re)initialize the pairHMM only if necessary + //final int readMaxLength = findMaxReadLength(reads); + //final int haplotypeMaxLength = findMaxHaplotypeLength(alleleHaplotypeMap); + //if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength) + //{ initialize(readMaxLength, haplotypeMaxLength); } + + //if ( ! initialized ) + //throw new IllegalStateException("Must call initialize before calling jniComputeLikelihoods"); + //final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap(); + //jniComputeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap, likelihoodMap); + //for(GATKSAMRecord read : reads){ + //final byte[] readBases = read.getReadBases(); + //final byte[] readQuals = read.getBaseQualities(); + //final byte[] readInsQuals = read.getBaseInsertionQualities(); + //final byte[] readDelQuals = read.getBaseDeletionQualities(); + //final byte[] overallGCP = GCPArrayMap.get(read); + + //// peak at the next haplotype in the list (necessary to get nextHaplotypeBases, which is required for caching in the array implementation) + //byte[] currentHaplotypeBases = null; + //boolean isFirstHaplotype = true; + //Allele currentAllele = null; + //double log10l; + //for (final Allele allele : alleleHaplotypeMap.keySet()){ + //final Haplotype haplotype = alleleHaplotypeMap.get(allele); + //final byte[] nextHaplotypeBases = haplotype.getBases(); + //if (currentHaplotypeBases != null) { + //log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, + //readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, nextHaplotypeBases); + //likelihoodMap.add(read, currentAllele, log10l); + //} + //// update the current haplotype + //currentHaplotypeBases = nextHaplotypeBases; + //currentAllele = allele; + //} + //// process the final haplotype + //if (currentHaplotypeBases != null) { + + //// there is no next haplotype, so pass null for nextHaplotypeBases. + //log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, + //readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, null); + //likelihoodMap.add(read, currentAllele, log10l); + //} + //} + //if(debug) + //debugClose(); + //return likelihoodMap; + //} diff --git a/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h b/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h index 5408d2f57..1fff9fa5a 100644 --- a/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h +++ b/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h @@ -54,6 +54,21 @@ JNIEXPORT jdouble JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILogless JNIEXPORT jdouble JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniSubComputeReadLikelihoodGivenHaplotypeLog10 (JNIEnv *, jobject, jint, jint, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jbyteArray, jint); +/* + * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM + * Method: jniGlobalInit + * Signature: (Ljava/lang/Class;Ljava/lang/Class;)V + */ +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniGlobalInit + (JNIEnv *, jobject, jclass, jclass); + +/* + * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM + * Method: jniComputeLikelihoods + * Signature: ([Lorg/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM/JNIReadDataHolderClass;[Lorg/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM/JNIHaplotypeDataHolderClass;[D)V + */ +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniComputeLikelihoods + (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jdoubleArray, jint); #ifdef __cplusplus } diff --git a/pairhmm-1-base.cc b/pairhmm-1-base.cc index dacd02536..90f95aeec 100644 --- a/pairhmm-1-base.cc +++ b/pairhmm-1-base.cc @@ -1,10 +1,6 @@ -#include -#include -#include -#include - -#include - +//#define DEBUG 1 +//#define DEBUG0_1 1 +//#define DEBUG3 1 #define MM 0 #define GapM 1 #define MX 2 @@ -12,123 +8,34 @@ #define MY 4 #define YY 5 -//#define DEBUG +#include +#include +#include +#include +#include +#include +#include "template.h" -/* - q: read quality - i: insertion penalty - d: deletion penalty - c: gap continuation penalty -*/ +//#include "define-float.h" +//#include "shift_template.c" +//#include "pairhmm-template-kernel.cc" -typedef struct +#include "define-double.h" +#include "shift_template.c" +#include "pairhmm-template-kernel.cc" + + +using namespace std; +class LoadTimeInitializer { - int rslen, haplen, *q, *i, *d, *c; - char *hap, *rs; -} testcase; - -int normalize(char c) -{ - return ((int) (c - 33)); -} - -int read_testcase(testcase *tc) -{ - char *q, *i, *d, *c, *line = NULL; - int _q, _i, _d, _c; - int x, size = 0; - ssize_t read; - - read = getline(&line, (size_t *) &size, stdin); - if (read == -1) - return -1; - - tc->hap = (char *) malloc(size); - tc->rs = (char *) malloc(size); - q = (char *) malloc(size); - i = (char *) malloc(size); - d = (char *) malloc(size); - c = (char *) malloc(size); - - if (sscanf(line, "%s %s %s %s %s %s\n", tc->hap, tc->rs, q, i, d, c) != 6) - return -1; - - tc->haplen = strlen(tc->hap); - tc->rslen = strlen(tc->rs); - tc->q = (int *) malloc(sizeof(int) * tc->rslen); - tc->i = (int *) malloc(sizeof(int) * tc->rslen); - tc->d = (int *) malloc(sizeof(int) * tc->rslen); - tc->c = (int *) malloc(sizeof(int) * tc->rslen); - - for (x = 0; x < tc->rslen; x++) - { - _q = normalize(q[x]); - _i = normalize(i[x]); - _d = normalize(d[x]); - _c = normalize(c[x]); - tc->q[x] = (_q < 6) ? 6 : _q; - tc->i[x] = _i; - tc->d[x] = _d; - tc->c[x] = _c; - } - - free(q); - free(i); - free(d); - free(c); - free(line); - - return 0; -} - -template -struct Context{}; - -template<> -struct Context -{ - Context() - { - for (int x = 0; x < 128; x++) - ph2pr[x] = pow(10.0, -((double)x) / 10.0); - - INITIAL_CONSTANT = ldexp(1.0, 1020.0); - LOG10_INITIAL_CONSTANT = log10(INITIAL_CONSTANT); - RESULT_THRESHOLD = 0.0; - } - - double LOG10(double v){ return log10(v); } - - static double _(double n){ return n; } - static double _(float n){ return ((double) n); } - double ph2pr[128]; - double INITIAL_CONSTANT; - double LOG10_INITIAL_CONSTANT; - double RESULT_THRESHOLD; + public: + LoadTimeInitializer() //will be called when library is loaded + { + ConvertChar::init(); + } }; +LoadTimeInitializer g_load_time_initializer; -template<> -struct Context -{ - Context() - { - for (int x = 0; x < 128; x++) - ph2pr[x] = powf(10.f, -((float)x) / 10.f); - - INITIAL_CONSTANT = ldexpf(1.f, 120.f); - LOG10_INITIAL_CONSTANT = log10f(INITIAL_CONSTANT); - RESULT_THRESHOLD = ldexpf(1.f, -110.f); - } - - float LOG10(float v){ return log10f(v); } - - static float _(double n){ return ((float) n); } - static float _(float n){ return n; } - float ph2pr[128]; - float INITIAL_CONSTANT; - float LOG10_INITIAL_CONSTANT; - float RESULT_THRESHOLD; -}; template NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) @@ -139,10 +46,10 @@ NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) Context ctx; - NUMBER M[ROWS][COLS]; - NUMBER X[ROWS][COLS]; - NUMBER Y[ROWS][COLS]; - NUMBER p[ROWS][6]; + NUMBER M[MROWS][MCOLS]; + NUMBER X[MROWS][MCOLS]; + NUMBER Y[MROWS][MCOLS]; + NUMBER p[MROWS][6]; p[0][MM] = ctx._(0.0); p[0][GapM] = ctx._(0.0); @@ -159,8 +66,10 @@ NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) p[r][GapM] = ctx._(1.0) - ctx.ph2pr[_c]; p[r][MX] = ctx.ph2pr[_i]; p[r][XX] = 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]; + p[r][MY] = ctx.ph2pr[_d]; + 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++) @@ -188,6 +97,8 @@ NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) 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]; @@ -201,76 +112,127 @@ NUMBER compute_full_prob(testcase *tc, NUMBER *before_last_log = NULL) if (before_last_log != NULL) *before_last_log = result; - return result; //ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT; + return ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT; } #define BATCH_SIZE 10000 #define RUN_HYBRID -int main() +int main(int argc, char** argv) { - testcase tc[BATCH_SIZE]; - float result[BATCH_SIZE], result_avxf; - double result_avxd; - struct timeval start, end; - long long aggregateTimeRead = 0L; - long long aggregateTimeCompute = 0L; - long long aggregateTimeWrite = 0L; + if(argc < 2) + { + cerr << "Needs path to input file as argument\n"; + exit(0); + } + bool use_old_read_testcase = false; + if(argc >= 3 && string(argv[2]) == "1") + use_old_read_testcase = true; - bool noMoreData = false; - int count =0; - while (!noMoreData) - { - int read_count = BATCH_SIZE; - gettimeofday(&start, NULL); - for (int b=0;b= 0) + { + double result_avxd = GEN_INTRINSIC(compute_full_prob_avx, d)(&tc); + double result = log10(result_avxd) - log10(ldexp(1.0, 1020)); - gettimeofday(&start, NULL); - for (int b=0;b(&tc[b]); + cout << std::scientific << compute_full_prob(&tc) << " "< tokens; + ifptr.open(argv[1]); + assert(ifptr.is_open()); + while(1) + { + tokens.clear(); + if(read_mod_testcase(ifptr, &tc, false) < 0) + break; + //double result = 0; + double result_avxd = GEN_INTRINSIC(compute_full_prob_avx, d)(&tc); + double result = log10(result_avxd) - log10(ldexp(1.0, 1020)); - #ifdef RUN_HYBRID - #define MIN_ACCEPTED 1e-28f - if (result_avxf < MIN_ACCEPTED) { - count++; - result_avxd = compute_full_prob(&tc[b]); - result[b] = log10(result_avxd) - log10(ldexp(1.0, 1020.f)); - } - else - result[b] = log10f(result_avxf) - log10f(ldexpf(1.f, 120.f)); - #endif + cout << std::scientific << compute_full_prob(&tc) << " "<(&tc[b]); - } +#ifdef RUN_HYBRID +#define MIN_ACCEPTED 1e-28f + if (result_avxf < MIN_ACCEPTED) { + count++; + result_avxd = compute_full_prob(&tc[b]); + result[b] = log10(result_avxd) - log10(ldexp(1.0, 1020.f)); + } + else + result[b] = log10f(result_avxf) - log10f(ldexpf(1.f, 120.f)); +#endif - printf("AVX Read Time: %ld\n", aggregateTimeRead); - printf("AVX Compute Time: %ld\n", aggregateTimeCompute); - printf("AVX Write Time: %ld\n", aggregateTimeWrite); - printf("AVX Total Time: %ld\n", aggregateTimeRead + aggregateTimeCompute + aggregateTimeWrite); - printf("# Double called: %d\n", count); +#ifndef RUN_HYBRID + result[b] = log10f(result_avxf) - log10f(ldexpf(1.f, 120.f)); +#endif - return 0; + } + gettimeofday(&end, NULL); + aggregateTimeCompute += ((end.tv_sec * 1000000 + end.tv_usec) - (start.tv_sec * 1000000 + start.tv_usec)); + + gettimeofday(&start, NULL); + for (int b=0;b #include +#include +#include +#include +#include + #define MROWS 500 #define MCOLS 1000 @@ -127,18 +132,20 @@ struct Context typedef struct { - int rslen, haplen, *q, *i, *d, *c; + int rslen, haplen; + /*int *q, *i, *d, *c;*/ + int q[MROWS], i[MROWS], d[MROWS], c[MROWS]; char *hap, *rs; - int *ihap; - int *irs; + int *ihap; + int *irs; } testcase; template std::string to_string(T obj) { - stringstream ss; - string ret_string; + std::stringstream ss; + std::string ret_string; ss.clear(); ss << std::scientific << obj; ss >> ret_string; @@ -148,66 +155,8 @@ std::string to_string(T obj) void debug_dump(std::string filename, std::string s, bool to_append, bool add_newline=true); int normalize(char c); - - -int read_testcase(testcase *tc) -{ - char *q, *i, *d, *c, *line = NULL; - int _q, _i, _d, _c; - int x, size = 0; - ssize_t read; - - read = getline(&line, (size_t *) &size, stdin); - if (read == -1) - return -1; - - tc->hap = (char *) malloc(size); - tc->rs = (char *) malloc(size); - q = (char *) malloc(size); - i = (char *) malloc(size); - d = (char *) malloc(size); - c = (char *) malloc(size); - tc->ihap = (int *) malloc(size*sizeof(int)); - tc->irs = (int *) malloc(size*sizeof(int)); - - - if (sscanf(line, "%s %s %s %s %s %s\n", tc->hap, tc->rs, q, i, d, c) != 6) - return -1; - - tc->haplen = strlen(tc->hap); - tc->rslen = strlen(tc->rs); - tc->q = (int *) malloc(sizeof(int) * tc->rslen); - tc->i = (int *) malloc(sizeof(int) * tc->rslen); - tc->d = (int *) malloc(sizeof(int) * tc->rslen); - tc->c = (int *) malloc(sizeof(int) * tc->rslen); - - tc->ihap = (int *) malloc(tc->haplen*sizeof(int)); - tc->irs = (int *) malloc(tc->rslen*sizeof(int)); - - for (x = 0; x < tc->rslen; x++) - { - _q = normalize(q[x]); - _i = normalize(i[x]); - _d = normalize(d[x]); - _c = normalize(c[x]); - tc->q[x] = (_q < 6) ? 6 : _q; - tc->i[x] = _i; - tc->d[x] = _d; - tc->c[x] = _c; - tc->irs[x] = tc->rs[x]; - } - - 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; -} +int read_testcase(testcase *tc, FILE* ifp); +int read_mod_testcase(std::ifstream& fptr, testcase* tc, bool reformat=false); #define NUM_DISTINCT_CHARS 5 #define AMBIG_CHAR 4