diff --git a/PairHMM_JNI/Makefile b/PairHMM_JNI/Makefile index 9b461a833..e207dd4cc 100644 --- a/PairHMM_JNI/Makefile +++ b/PairHMM_JNI/Makefile @@ -13,7 +13,7 @@ CXX=icc LDFLAGS=-lm -lrt $(OMPLDFLAGS) -BIN=libJNILoglessPairHMM.so pairhmm-template-main checker +BIN=libVectorLoglessPairHMM.so pairhmm-template-main checker #BIN=checker DEPDIR=.deps @@ -22,14 +22,14 @@ DF=$(DEPDIR)/$(*).d #Common across libJNI and sandbox COMMON_SOURCES=utils.cc avx_function_instantiations.cc baseline.cc sse_function_instantiations.cc LoadTimeInitializer.cc #Part of libJNI -LIBSOURCES=org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc $(COMMON_SOURCES) +LIBSOURCES=org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc $(COMMON_SOURCES) SOURCES=$(LIBSOURCES) pairhmm-template-main.cc pairhmm-1-base.cc LIBOBJECTS=$(LIBSOURCES:.cc=.o) COMMON_OBJECTS=$(COMMON_SOURCES:.cc=.o) #No vectorization for these files -NO_VECTOR_SOURCES=org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc pairhmm-template-main.cc pairhmm-1-base.cc utils.cc baseline.cc LoadTimeInitializer.cc +NO_VECTOR_SOURCES=org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc pairhmm-template-main.cc pairhmm-1-base.cc utils.cc baseline.cc LoadTimeInitializer.cc #Use -xAVX for these files AVX_SOURCES=avx_function_instantiations.cc #Use -xSSE4.2 for these files @@ -53,7 +53,7 @@ checker: pairhmm-1-base.o $(COMMON_OBJECTS) pairhmm-template-main: pairhmm-template-main.o $(COMMON_OBJECTS) $(CXX) $(OMPLFLAGS) -o $@ $^ $(LDFLAGS) -libJNILoglessPairHMM.so: $(LIBOBJECTS) +libVectorLoglessPairHMM.so: $(LIBOBJECTS) $(CXX) $(OMPLFLAGS) -shared -o $@ $(LIBOBJECTS) ${LDFLAGS} diff --git a/PairHMM_JNI/jni_common.h b/PairHMM_JNI/jni_common.h new file mode 100644 index 000000000..d1ce86e9c --- /dev/null +++ b/PairHMM_JNI/jni_common.h @@ -0,0 +1,32 @@ +#ifndef JNI_COMMON_H +#define JNI_COMMON_H + +#include +#define ENABLE_ASSERTIONS 1 +#define DO_PROFILING 1 +//#define DEBUG 1 +//#define DEBUG0_1 1 +//#define DEBUG3 1 + + +#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_RO_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_RO_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 + +#endif //ifndef JNI_COMMON_H diff --git a/PairHMM_JNI/jnidebug.h b/PairHMM_JNI/jnidebug.h new file mode 100644 index 000000000..7002b3637 --- /dev/null +++ b/PairHMM_JNI/jnidebug.h @@ -0,0 +1,166 @@ +#ifndef JNI_DEBUG_H +#define JNI_DEBUG_H + +template +class DataHolder +{ +#define INIT_MATRIX(X) \ + X = new NUMBER*[m_paddedMaxReadLength]; \ + for(int i=0;i ctx; + for (int r = 1; r <= length;r++) //in original code, r < ROWS (where ROWS = paddedReadLength) + { + int _i = insertionGOP[r-1]; //insertionGOP + int _d = deletionGOP[r-1]; //deletionGOP + int _c = overallGCP[r-1]; //overallGCP + m_transition[r][MM] = ctx._(1.0) - ctx.ph2pr[(_i + _d) & 127]; //lines 161-162 + m_transition[r][GapM] = ctx._(1.0) - ctx.ph2pr[_c]; //line 163 + m_transition[r][MX] = ctx.ph2pr[_i]; //164 + m_transition[r][XX] = ctx.ph2pr[_c]; //165 + m_transition[r][MY] = ctx.ph2pr[_d];//last row seems different, compared to line 166 + m_transition[r][YY] = ctx.ph2pr[_c];//same as above for line 167 + //m_transition[r][MY] = (r == length) ? ctx._(1.0) : ctx.ph2pr[_d];//last row seems different, compared to line 166 + //m_transition[r][YY] = (r == length) ? ctx._(1.0) : ctx.ph2pr[_c];//same as above for line 167 +#ifdef DEBUG3 + for(int j=0;j<6;++j) + debug_dump("transitions_jni.txt", to_string(m_transition[r][j]),true); +#endif + } + ++g_num_prob_init; + } + bool m_is_initialized; + int m_readMaxLength; + int m_haplotypeMaxLength; + int m_paddedMaxReadLength; + int m_paddedMaxHaplotypeLength; + NUMBER** m_matchMatrix; + NUMBER** m_insertionMatrix; + NUMBER** m_deletionMatrix; + NUMBER** m_prior; + NUMBER (*m_transition)[6]; +}; +extern DataHolder g_double_dataholder; + +template +NUMBER compute_full_prob(testcase *tc, NUMBER** M, NUMBER** X, NUMBER** Y, NUMBER (*p)[6], + bool do_initialization, jint hapStartIndex, NUMBER *before_last_log = NULL) +{ + int r, c; + int ROWS = tc->rslen + 1; //ROWS = paddedReadLength + int COLS = tc->haplen + 1; //COLS = paddedHaplotypeLength + + Context ctx; + //////NOTES + ////ctx.ph2pr[quality]; //This quantity is QualityUtils.qualToErrorProb(quality) + ////1-ctx.ph2pr[quality]; //This corresponds to QualityUtils.qualToProb(quality); + + //Initialization + if(do_initialization) + { + 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); //code from 87-90 in LoglessPairHMM + } + + for (r = 1; r < ROWS; r++) + { + M[r][0] = ctx._(0.0); + //deletionMatrix row 0 in above nest is initialized in the Java code + //However, insertionMatrix column 0 is not initialized in Java code, could it be that + //values are re-used from a previous iteration? + //Why even do this, X[0][0] = 0 from above loop nest, X[idx][0] = 0 from this computation + X[r][0] = X[r-1][0] * p[r][XX]; + Y[r][0] = ctx._(0.0); + } + } + + for (r = 1; r < ROWS; r++) + for (c = hapStartIndex+1; c < COLS; c++) + { + //The following lines correspond to initializePriors() + char _rs = tc->rs[r-1]; //line 137 + char _hap = tc->hap[c-1]; //line 140 + //int _q = tc->q[r-1] & 127; //line 138 - q is the quality (qual), should be byte hence int ANDed with 0xFF + int _q = tc->q[r-1]; //line 138 - q is the quality (qual), should be byte hence int ANDed with 0xFF + NUMBER distm = ctx.ph2pr[_q]; //This quantity is QualityUtils.qualToErrorProb(_q) + //The assumption here is that doNotUseTristateCorrection is true + //TOASK + if (_rs == _hap || _rs == 'N' || _hap == 'N') + distm = ctx._(1.0) - distm; //This is the quantity QualityUtils.qualToProb(qual) + else + distm = distm/3; +#ifdef DEBUG3 + debug_dump("priors_jni.txt",to_string(distm),true); +#endif + + //Computation inside updateCell + 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]; +#ifdef DEBUG3 + debug_dump("matrices_jni.txt",to_string(M[r][c]),true); + debug_dump("matrices_jni.txt",to_string(X[r][c]),true); + debug_dump("matrices_jni.txt",to_string(Y[r][c]),true); +#endif + } + + NUMBER result = ctx._(0.0); + for (c = 0; c < COLS; c++) + result += M[ROWS-1][c] + X[ROWS-1][c]; + + if (before_last_log != NULL) + *before_last_log = result; + +#ifdef DEBUG + debug_dump("return_values_jni.txt",to_string(ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT),true); +#endif + return ctx.LOG10(result) - ctx.LOG10_INITIAL_CONSTANT; +} + +#endif diff --git a/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc new file mode 100644 index 000000000..30f118c31 --- /dev/null +++ b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.cc @@ -0,0 +1,151 @@ +#include "headers.h" +#include "jni_common.h" +#include "org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM.h" +#include "template.h" +#include "utils.h" +#include "LoadTimeInitializer.h" +#include "jnidebug.h" +DataHolder g_double_dataholder; + +using namespace std; + +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniInitialize +(JNIEnv* env, jobject thisObject, + jint readMaxLength, jint haplotypeMaxLength) +{ + static int g_num_init_calls = 0; +#ifdef DEBUG3 + cout << "Entered alloc initialized .. readMaxLength "<GetArrayLength(insertionGOP); +#ifdef DEBUG3 + cout << "Entered initializeProbabilities .. length "<GetByteArrayElements(insertionGOP, &is_copy); + jbyte* deletionGOPArray = (env)->GetByteArrayElements(deletionGOP, &is_copy); + jbyte* overallGCPArray = (env)->GetByteArrayElements(overallGCP, &is_copy); +#ifdef DEBUG + if(insertionGOPArray == 0) + cerr << "insertionGOP array not initialized in JNI\n"; + ////assert(insertionGOPArray && "insertionGOP array not initialized in JNI"); + if(deletionGOPArray == 0) + cerr << "deletionGOP array not initialized in JNI\n"; + ////assert(deletionGOPArray && "deletionGOP array not initialized in JNI"); + assert(overallGCPArray && "OverallGCP array not initialized in JNI"); +#endif + + g_double_dataholder.initializeProbabilities(length, insertionGOPArray, deletionGOPArray, overallGCPArray); + + env->ReleaseByteArrayElements(overallGCP, overallGCPArray, JNI_ABORT); + env->ReleaseByteArrayElements(deletionGOP, deletionGOPArray, JNI_ABORT); + env->ReleaseByteArrayElements(insertionGOP, insertionGOPArray, JNI_ABORT); +} + +JNIEXPORT jdouble JNICALL +Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_jniInitializePriorsAndUpdateCells( + JNIEnv* env, jobject thisObject, + jboolean doInitialization, jint paddedReadLength, jint paddedHaplotypeLength, + jbyteArray readBases, jbyteArray haplotypeBases, jbyteArray readQuals, + jint hapStartIndex + ) +{ +#ifdef DEBUG3 + cout << "Entered mainCompute .. doInitialization "<<(doInitialization == JNI_TRUE)<<" hapStartIndex "<GetByteArrayElements(readBases, &is_copy); + jbyte* haplotypeBasesArray = (env)->GetByteArrayElements(haplotypeBases, &is_copy); + jbyte* readQualsArray = (env)->GetByteArrayElements(readQuals, &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"); +#endif + testcase tc; + + tc.rslen = paddedReadLength-1; + tc.haplen = paddedHaplotypeLength-1; + + tc.rs = (char*)readBasesArray; + tc.hap = (char*)haplotypeBasesArray; + tc.q = (char*)readQualsArray; //TOASK - q is now char* + + compute_full_prob(&tc, g_double_dataholder.m_matchMatrix, g_double_dataholder.m_insertionMatrix, + g_double_dataholder.m_deletionMatrix, g_double_dataholder.m_transition, + doInitialization == JNI_TRUE, hapStartIndex, NULL); + + env->ReleaseByteArrayElements(readBases, readBasesArray, JNI_ABORT); + env->ReleaseByteArrayElements(haplotypeBases, haplotypeBasesArray, JNI_ABORT); + env->ReleaseByteArrayElements(readQuals, readQualsArray, JNI_ABORT); + return 0.0; +} + +JNIEXPORT jdouble JNICALL +Java_org_broadinstitute_sting_utils_pairhmm_DebugJNILoglessPairHMM_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 -/* Header for class org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM */ - -#ifndef _Included_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM -#define _Included_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM -#ifdef __cplusplus -extern "C" { -#endif -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_TRISTATE_CORRECTION -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_TRISTATE_CORRECTION 3.0 -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToMatch -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToMatch 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_indelToMatch -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_indelToMatch 1L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToInsertion -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToInsertion 2L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_insertionToInsertion -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_insertionToInsertion 3L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToDeletion -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_matchToDeletion 4L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_deletionToDeletion -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_deletionToDeletion 5L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_verify -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_verify 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug0_1 -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug0_1 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug1 -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug1 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug2 -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug2 0L -#undef org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug3 -#define org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_debug3 0L -/* - * 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: jniClose - * Signature: ()V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniClose - (JNIEnv *, jobject); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniInitialize - * Signature: (II)V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitialize - (JNIEnv *, jobject, jint, jint); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniInitializeProbabilities - * Signature: ([[D[B[B[B)V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializeProbabilities - (JNIEnv *, jclass, jobjectArray, jbyteArray, jbyteArray, jbyteArray); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniInitializePriorsAndUpdateCells - * Signature: (ZII[B[B[BI)D - */ -JNIEXPORT jdouble JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializePriorsAndUpdateCells - (JNIEnv *, jobject, jboolean, jint, jint, jbyteArray, jbyteArray, jbyteArray, jint); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniSubComputeReadLikelihoodGivenHaplotypeLog10 - * Signature: (II[B[B[B[B[B[BI)D - */ -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: jniInitializeHaplotypes - * Signature: (I[Lorg/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM/JNIHaplotypeDataHolderClass;)V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializeHaplotypes - (JNIEnv *, jobject, jint, jobjectArray); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniFinalizeRegion - * Signature: ()V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniFinalizeRegion - (JNIEnv *, jobject); - -/* - * Class: org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM - * Method: jniComputeLikelihoods - * Signature: (II[Lorg/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM/JNIReadDataHolderClass;[Lorg/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM/JNIHaplotypeDataHolderClass;[DI)V - */ -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniComputeLikelihoods - (JNIEnv *, jobject, jint, jint, jobjectArray, jobjectArray, jdoubleArray, jint); - -#ifdef __cplusplus -} -#endif -#endif diff --git a/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc similarity index 72% rename from PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc rename to PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc index a4cd63bc4..54df4265f 100644 --- a/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc +++ b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.cc @@ -1,120 +1,16 @@ #include "headers.h" -#include -#include "org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.h" - - -#define ENABLE_ASSERTIONS 1 -#define DO_PROFILING 1 -//#define DEBUG 1 -//#define DEBUG0_1 1 -//#define DEBUG3 1 - +#include "jni_common.h" +#include "org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM.h" #include "template.h" #include "utils.h" #include "LoadTimeInitializer.h" using namespace std; -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializeProbabilities -(JNIEnv* env, jclass thisObject, - jobjectArray transition, jbyteArray insertionGOP, jbyteArray deletionGOP, jbyteArray overallGCP - ) -{} - - JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitialize -(JNIEnv* env, jobject thisObject, - jint readMaxLength, jint haplotypeMaxLength) -{} - -JNIEXPORT jdouble JNICALL -Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializePriorsAndUpdateCells( - JNIEnv* env, jobject thisObject, - jboolean doInitialization, jint paddedReadLength, jint paddedHaplotypeLength, - jbyteArray readBases, jbyteArray haplotypeBases, jbyteArray readQuals, - jint hapStartIndex - ) - -{ 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_RO_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_RO_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 > g_haplotypeBasesArrayVector; -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniInitializeHaplotypes +vector g_haplotypeBasesLengths; +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniInitializeHaplotypes (JNIEnv * env, jobject thisObject, jint numHaplotypes, jobjectArray haplotypeDataArray) { 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 = g_haplotypeBasesArrayVector; haplotypeBasesArrayVector.clear(); + g_haplotypeBasesLengths.clear(); haplotypeBasesArrayVector.resize(numHaplotypes); + g_haplotypeBasesLengths.resize(numHaplotypes); + jsize haplotypeBasesLength = 0; for(unsigned j=0;jGetObjectArrayElement(haplotypeDataArray, j); @@ -165,20 +65,22 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai #endif env->DeleteLocalRef(haplotypeBases); //free the local reference jbyte* haplotypeBasesArray = (jbyte*)GET_BYTE_ARRAY_ELEMENTS(haplotypeBasesGlobalRef, &is_copy); + haplotypeBasesLength = env->GetArrayLength(haplotypeBasesGlobalRef); #ifdef ENABLE_ASSERTIONS assert(haplotypeBasesArray && "haplotypeBasesArray not initialized in JNI"); - assert(env->GetArrayLength(haplotypeBasesGlobalRef) < MCOLS); + assert(haplotypeBasesLength < MCOLS); #endif #ifdef DEBUG0_1 - cout << "JNI haplotype length "<GetArrayLength(haplotypeBasesGlobalRef)<<"\n"; + cout << "JNI haplotype length "<GetArrayLength(haplotypeBases);++k) + for(unsigned k=0;kGetArrayLength(haplotypeBasesGlobalRef); + g_load_time_initializer.m_sumHaplotypeLengths += haplotypeBasesLength; #endif } } @@ -188,7 +90,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai //haplotypeDataArray - array of JNIHaplotypeDataHolderClass objects which contain the haplotypeBases //likelihoodArray - array of doubles to return results back to Java. Memory allocated by Java prior to JNI call //maxNumThreadsToUse - Max number of threads that OpenMP can use for the HMM computation -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniComputeLikelihoods +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniComputeLikelihoods (JNIEnv* env, jobject thisObject, jint numReads, jint numHaplotypes, jobjectArray readDataArray, jobjectArray haplotypeDataArray, jdoubleArray likelihoodArray, jint maxNumThreadsToUse) { @@ -265,7 +167,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai for(unsigned j=0;jGetArrayLength(haplotypeBasesArrayVector[j].first); + jsize haplotypeLength = (jsize)g_haplotypeBasesLengths[j]; jbyte* haplotypeBasesArray = haplotypeBasesArrayVector[j].second; tc_array[tc_idx].rslen = (int)readLength; tc_array[tc_idx].haplen = (int)haplotypeLength; @@ -363,7 +265,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai } //Release haplotypes at the end of a region -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniFinalizeRegion +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniFinalizeRegion (JNIEnv * env, jobject thisObject) { vector >& haplotypeBasesArrayVector = g_haplotypeBasesArrayVector; @@ -373,11 +275,12 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai RELEASE_BYTE_ARRAY_ELEMENTS(haplotypeBasesArrayVector[j].first, haplotypeBasesArrayVector[j].second, JNI_RO_RELEASE_MODE); env->DeleteGlobalRef(haplotypeBasesArrayVector[j].first); //free the global reference } - haplotypeBasesArrayVector.clear(); + haplotypeBasesArrayVector.clear(); + g_haplotypeBasesLengths.clear(); } -JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniClose +JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_VectorLoglessPairHMM_jniClose (JNIEnv* env, jobject thisObject) { #ifdef DO_PROFILING diff --git a/PairHMM_JNI/run.sh b/PairHMM_JNI/run.sh index 65a53289a..a1eacdefa 100755 --- a/PairHMM_JNI/run.sh +++ b/PairHMM_JNI/run.sh @@ -1,14 +1,20 @@ #!/bin/bash rm -f *.txt *.log GSA_ROOT_DIR=/home/karthikg/broad/gsa-unstable +pair_hmm_implementation="VECTOR_LOGLESS_CACHING"; +if [ "$#" -ge 1 ]; +then + pair_hmm_implementation=$1; +fi + #-Djava.library.path is needed if you are using JNI_LOGLESS_CACHING, else not needed java -Djava.library.path=${GSA_ROOT_DIR}/PairHMM_JNI/ -jar $GSA_ROOT_DIR/dist/GenomeAnalysisTK.jar -T HaplotypeCaller \ +--dbsnp /data/broad/samples/joint_variant_calling/dbSNP/00-All.vcf \ -R /data/broad/samples/joint_variant_calling/broad_reference/human_g1k_v37_decoy.fasta \ -I /data/broad/samples/joint_variant_calling/NA12878_low_coverage_alignment/NA12878.chrom11.ILLUMINA.bwa.CEU.low_coverage.20121211.bam \ ---dbsnp /data/broad/samples/joint_variant_calling/dbSNP/00-All.vcf \ -stand_call_conf 50.0 \ -stand_emit_conf 10.0 \ ---pair_hmm_implementation JNI_LOGLESS_CACHING \ +--pair_hmm_implementation $pair_hmm_implementation \ -o output.raw.snps.indels.vcf #--pair_hmm_implementation JNI_LOGLESS_CACHING \ diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java index 66d8ee33c..c534764a5 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/PairHMMLikelihoodCalculationEngine.java @@ -93,8 +93,8 @@ public class PairHMMLikelihoodCalculationEngine implements LikelihoodCalculation return new CnyPairHMM(); case VECTOR_LOGLESS_CACHING: return new VectorLoglessPairHMM(); - case DEBUG_JNI_LOGLESS_CACHING: - return new DebugJNILoglessPairHMM(hmmType); + case DEBUG_VECTOR_LOGLESS_CACHING: + return new DebugJNILoglessPairHMM(PairHMM.HMM_IMPLEMENTATION.VECTOR_LOGLESS_CACHING); case ARRAY_LOGLESS: if (noFpga || !CnyPairHMM.isAvailable()) return new ArrayLoglessPairHMM(); diff --git a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java index 294c5c451..55491bcd8 100644 --- a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java +++ b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/DebugJNILoglessPairHMM.java @@ -76,7 +76,7 @@ import java.io.IOException; public class DebugJNILoglessPairHMM extends LoglessPairHMM { private static final boolean debug = false; //simulates ifdef - private static final boolean verify = debug || false; //simulates ifdef + private static final boolean verify = debug || true; //simulates ifdef private static final boolean debug0_1 = false; //simulates ifdef private static final boolean debug1 = false; //simulates ifdef private static final boolean debug2 = false; diff --git a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java index 1c7adf4c3..9b0772941 100644 --- a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java +++ b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java @@ -58,4 +58,5 @@ import java.util.HashMap; */ public abstract class JNILoglessPairHMM extends LoglessPairHMM { public abstract HashMap getHaplotypeToHaplotypeListIdxMap(); + protected long setupTime = 0; } diff --git a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java index 9d6812c0c..6207c9679 100644 --- a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java +++ b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/VectorLoglessPairHMM.java @@ -70,6 +70,7 @@ public class VectorLoglessPairHMM extends JNILoglessPairHMM { //For machine capabilities public static final long sse42Mask = 1; public static final long avxMask = 2; + public static final long enableAll = 0xFFFFFFFFFFFFFFFFl; //Used to copy references to byteArrays to JNI from reads protected class JNIReadDataHolderClass { @@ -95,13 +96,15 @@ public class VectorLoglessPairHMM extends JNILoglessPairHMM { if(!isVectorLoglessPairHMMLibraryLoaded) { System.loadLibrary("VectorLoglessPairHMM"); isVectorLoglessPairHMMLibraryLoaded = true; - jniGlobalInit(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, 0xFFFFFFFFFFFFFFFFl); //need to do this only once + jniGlobalInit(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class, enableAll); //need to do this only once } } @Override public void close() { + System.out.println("Time spent in setup for JNI call : "+setupTime); + super.close(); jniClose(); } @@ -154,6 +157,8 @@ public class VectorLoglessPairHMM extends JNILoglessPairHMM { */ @Override public PerReadAlleleLikelihoodMap computeLikelihoods( final List reads, final Map alleleHaplotypeMap, final Map GCPArrayMap ) { + if(doProfiling) + startTime = System.nanoTime(); int readListSize = reads.size(); int numHaplotypes = alleleHaplotypeMap.size(); int numTestcases = readListSize*numHaplotypes; @@ -170,12 +175,13 @@ public class VectorLoglessPairHMM extends JNILoglessPairHMM { ++idx; } - JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes]; mLikelihoodArray = new double[readListSize*numHaplotypes]; //to store results + if(doProfiling) + setupTime += (System.nanoTime() - startTime); //for(reads) // for(haplotypes) // compute_full_prob() - jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, haplotypeDataArray, mLikelihoodArray, 12); + jniComputeLikelihoods(readListSize, numHaplotypes, readDataArray, null, mLikelihoodArray, 12); final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap(); idx = 0; @@ -193,6 +199,8 @@ public class VectorLoglessPairHMM extends JNILoglessPairHMM { } readIdx += numHaplotypes; } + if(doProfiling) + computeTime += (System.nanoTime() - startTime); return likelihoodMap; } } diff --git a/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java b/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java index 8b2123751..4ab17d5f3 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java +++ b/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java @@ -59,8 +59,8 @@ public abstract class PairHMM { LOGLESS_CACHING, /* Optimized AVX implementation of LOGLESS_CACHING called through JNI */ VECTOR_LOGLESS_CACHING, - /* Debugging for any JNI implementation of LOGLESS_CACHING */ - DEBUG_JNI_LOGLESS_CACHING, + /* Debugging for vector implementation of LOGLESS_CACHING */ + DEBUG_VECTOR_LOGLESS_CACHING, /* Logless caching PairHMM that stores computations in 1D arrays instead of matrices, and which proceeds diagonally over the (read x haplotype) intersection matrix */ ARRAY_LOGLESS } @@ -77,6 +77,11 @@ public abstract class PairHMM { //debug array protected double[] mLikelihoodArray; + //profiling information + protected static final boolean doProfiling = true; + protected long computeTime = 0; + protected long startTime = 0; + /** * Initialize this PairHMM, making it suitable to run against a read and haplotype with given lengths * @@ -107,7 +112,7 @@ public abstract class PairHMM { */ public void finalizeRegion() { - ; + ; } /** @@ -119,7 +124,7 @@ public abstract class PairHMM { * @param readMaxLength the max length of reads we want to use with this PairHMM */ public void initialize( final List haplotypes, final Map> perSampleReadList, final int readMaxLength, final int haplotypeMaxLength ) { - initialize(readMaxLength, haplotypeMaxLength); + initialize(readMaxLength, haplotypeMaxLength); } protected int findMaxReadLength(final List reads) { @@ -152,6 +157,8 @@ public abstract class PairHMM { * said read coming from the said haplotype under the provided error model */ public PerReadAlleleLikelihoodMap computeLikelihoods(final List reads, final Map alleleHaplotypeMap, final Map GCPArrayMap) { + if(doProfiling) + startTime = System.nanoTime(); // (re)initialize the pairHMM only if necessary final int readMaxLength = findMaxReadLength(reads); @@ -159,8 +166,8 @@ public abstract class PairHMM { if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength) { initialize(readMaxLength, haplotypeMaxLength); } final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap(); - mLikelihoodArray = new double[reads.size()*alleleHaplotypeMap.size()]; - int idx = 0; + mLikelihoodArray = new double[reads.size()*alleleHaplotypeMap.size()]; + int idx = 0; for(GATKSAMRecord read : reads){ final byte[] readBases = read.getReadBases(); final byte[] readQuals = read.getBaseQualities(); @@ -173,16 +180,16 @@ public abstract class PairHMM { boolean isFirstHaplotype = true; Allele currentAllele = null; double log10l; - //for (final Allele allele : alleleHaplotypeMap.keySet()){ + //for (final Allele allele : alleleHaplotypeMap.keySet()){ for (Map.Entry currEntry : alleleHaplotypeMap.entrySet()){ - //final Haplotype haplotype = alleleHaplotypeMap.get(allele); - final Allele allele = currEntry.getKey(); - final Haplotype haplotype = currEntry.getValue(); + //final Haplotype haplotype = alleleHaplotypeMap.get(allele); + final Allele allele = currEntry.getKey(); + final Haplotype haplotype = currEntry.getValue(); final byte[] nextHaplotypeBases = haplotype.getBases(); if (currentHaplotypeBases != null) { log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, nextHaplotypeBases); - mLikelihoodArray[idx++] = log10l; + mLikelihoodArray[idx++] = log10l; likelihoodMap.add(read, currentAllele, log10l); } // update the current haplotype @@ -196,9 +203,11 @@ public abstract class PairHMM { log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, null); likelihoodMap.add(read, currentAllele, log10l); - mLikelihoodArray[idx++] = log10l; + mLikelihoodArray[idx++] = log10l; } } + if(doProfiling) + computeTime += (System.nanoTime() - startTime); return likelihoodMap; } @@ -305,5 +314,9 @@ public abstract class PairHMM { public double[] getLikelihoodArray() { return mLikelihoodArray; } //Called at the end of all HC calls - public void close() { ; } + public void close() + { + if(doProfiling) + System.out.println("Total compute time in PairHMM computeLikelihoods() : "+computeTime); + } }