diff --git a/PairHMM_JNI/.gitignore b/PairHMM_JNI/.gitignore index cb70b53a7..9722ad970 100644 --- a/PairHMM_JNI/.gitignore +++ b/PairHMM_JNI/.gitignore @@ -1,6 +1,6 @@ .svn *.o -*.so +#*.so tests .deps hmm_Mohammad diff --git a/PairHMM_JNI/JNI_README b/PairHMM_JNI/JNI_README index 2a02d41ce..9fc3aa8d1 100644 --- a/PairHMM_JNI/JNI_README +++ b/PairHMM_JNI/JNI_README @@ -1 +1,33 @@ -TEST +Implementation overview: +Created a new Java class called JNILoglessPairHMM which extends LoglessPairHMM and +overrides functions from both LoglessPairHMM and PairHMM. +1. Constructor: Call base class constructors. Then, load the JNI library located in this +directory and call a global init function in the library to determine fields ids for the +members of classes JNIReadDataHolder and JNIHaplotypeDataHolders +2. initialize(): Override and do nothing as all matrix manipulation is done in the native +library. +3. computeLikelihoods(): Copies array references for readBases/quals etc and +haplotypeBases to arrays of JNIReadDataHolder and JNIHaplotypeDataHolders classes. Invokes +the JNI function to perform the computation and updates the likelihoodMap. + +Note: Lots of debug code still left in the files, however, debug code is not executed +(hopefully, not even compiled into bytecode) because the debug flags are set to false. + +On the C++ side, the primary function called is +Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM_jniComputeLikelihoods. It +uses standard JNI calls to get and return data from/to the Java class JNILoglessPairHMM. +The last argument to the function is the maximum number of OpenMP threads to use while +computing PairHMM in C++. This option is set when the native function call is made from +JNILoglessPairHMM computeLikelihoods - currently it is set to 12 (no logical reason). + +Compiling: +Make sure you have icc (Intel C compiler) available. Currently, gcc does not seem to +support all AVX intrinsics. +Type 'make'. This should create a library called libJNILoglessPairHMM.so .Comment out +OMPCFLAGS if OpenMP is not desired. + +Running: +If libJNILoglessPairHMM.so is compiled using icc, make sure that the Intel Composer XE +libraries are in your LD_LIBRARY_PATH : +source /bin/compilervars.sh intel64 +See run.sh in this directory on how to invoke HaplotypeCaller with the native library diff --git a/PairHMM_JNI/libJNILoglessPairHMM.so b/PairHMM_JNI/libJNILoglessPairHMM.so index f7a3569bf..435034f78 100755 Binary files a/PairHMM_JNI/libJNILoglessPairHMM.so and b/PairHMM_JNI/libJNILoglessPairHMM.so differ diff --git a/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc index 852caa00d..59d00928e 100644 --- a/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc +++ b/PairHMM_JNI/org_broadinstitute_sting_utils_pairhmm_JNILoglessPairHMM.cc @@ -18,6 +18,8 @@ #include #include using namespace std; + +#define ENABLE_ASSERTIONS 1 //#define DEBUG 1 //#define DEBUG0_1 1 //#define DEBUG3 1 @@ -189,7 +191,11 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai g_load_time_initializer.m_haplotypeBasesFID = fid; } - +//JNI function to invoke compute_full_prob_avx +//readDataArray - array of JNIReadDataHolderClass objects which contain the readBases, readQuals etc +//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 (JNIEnv* env, jobject thisObject, jint numReads, jint numHaplotypes, jobjectArray readDataArray, jobjectArray haplotypeDataArray, jdoubleArray likelihoodArray, jint maxNumThreadsToUse) @@ -206,16 +212,16 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai { jobject haplotypeObject = env->GetObjectArrayElement(haplotypeDataArray, j); jbyteArray haplotypeBases = (jbyteArray)env->GetObjectField(haplotypeObject, g_load_time_initializer.m_haplotypeBasesFID); -#ifdef DEBUG +#ifdef ENABLE_ASSERTIONS 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 +#ifdef ENABLE_ASSERTIONS assert(haplotypeBasesArray && "haplotypeBasesArray not initialized in JNI"); assert(env->GetArrayLength(haplotypeBases) < MCOLS); +#endif #ifdef DEBUG0_1 cout << "JNI haplotype length "<GetArrayLength(haplotypeBases)<<"\n"; -#endif #endif haplotypeBasesArrayVector[j] = make_pair(haplotypeBases, haplotypeBasesArray); #ifdef DEBUG3 @@ -243,7 +249,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai 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 +#ifdef ENABLE_ASSERTIONS 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()); @@ -257,7 +263,7 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai 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 +#ifdef ENABLE_ASSERTIONS assert(readBasesArray && "readBasesArray not initialized in JNI"); assert(readQualsArray && "readQualsArray not initialized in JNI"); assert(insertionGOPArray && "insertionGOP array not initialized in JNI"); @@ -268,10 +274,10 @@ JNIEXPORT void JNICALL Java_org_broadinstitute_sting_utils_pairhmm_JNILoglessPai assert(readLength == env->GetArrayLength(insertionGOP)); assert(readLength == env->GetArrayLength(deletionGOP)); assert(readLength == env->GetArrayLength(overallGCP)); +#endif #ifdef DEBUG0_1 cout << "JNI read length "<GetArrayLength(likelihoodArray) == numTestCases); +#endif #pragma omp parallel for schedule (dynamic,10) private(tc_idx) num_threads(maxNumThreadsToUse) for(tc_idx=0;tc_idx readDataHolderClass, Class haplotypeDataHolderClass); + + private static boolean isJNILoglessPairHMMLibraryLoaded = false; + + public JNILoglessPairHMM() + { + super(); + if(!isJNILoglessPairHMMLibraryLoaded) + { + System.loadLibrary("JNILoglessPairHMM"); + isJNILoglessPairHMMLibraryLoaded = true; + jniGlobalInit(JNIReadDataHolderClass.class, JNIHaplotypeDataHolderClass.class); //need to do this only once } } - - private native void jniInitialize(final int readMaxLength, final int haplotypeMaxLength); - - private native static void jniInitializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP); + //Used to test parts of the compute kernel separately + private native void jniInitialize(final int readMaxLength, final int haplotypeMaxLength); + private native static void jniInitializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP); private native double jniInitializePriorsAndUpdateCells( boolean doInitialization, final int paddedReadLength, final int paddedHaplotypeLength, final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals, final int hapStartIndex); - private native double jniSubComputeReadLikelihoodGivenHaplotypeLog10( final int readLength, final int haplotypeLength, final byte[] readBases, final byte[] haplotypeBases, final byte[] readQuals, @@ -101,14 +119,15 @@ public class JNILoglessPairHMM extends LoglessPairHMM { final int hapStartIndex); - + //Used only when testing parts of the compute kernel /** * {@inheritDoc} */ @Override public void initialize(final int readMaxLength, final int haplotypeMaxLength) { - super.initialize(readMaxLength, haplotypeMaxLength); + if(debug) + super.initialize(readMaxLength, haplotypeMaxLength); if(debug3) { System.out.println("Java: alloc initialized readMaxLength : "+readMaxLength+" haplotypeMaxLength : "+haplotypeMaxLength); @@ -119,97 +138,142 @@ public class JNILoglessPairHMM extends LoglessPairHMM { jniInitialize(readMaxLength, haplotypeMaxLength); } + //Real compute kernel + private native void jniComputeLikelihoods(int numReads, int numHaplotypes, + JNIReadDataHolderClass[] readDataArray, JNIHaplotypeDataHolderClass[] haplotypeDataArray, + double[] likelihoodArray, int maxNumThreadsToUse); /** * {@inheritDoc} */ @Override 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 = debug ? findMaxReadLength(reads) : 0; + final int haplotypeMaxLength = debug ? findMaxHaplotypeLength(alleleHaplotypeMap) : 0; if(debug) + { + if (!initialized || readMaxLength > maxReadLength || haplotypeMaxLength > maxHaplotypeLength) + { initialize(readMaxLength, haplotypeMaxLength); } + if ( ! initialized ) + throw new IllegalStateException("Must call initialize before calling jniComputeLikelihoods in debug mode"); + } + int readListSize = reads.size(); + int alleleHaplotypeMapSize = alleleHaplotypeMap.size(); + if(debug0_1) + System.out.println("Java numReads "+readListSize+" numHaplotypes "+alleleHaplotypeMapSize); + JNIReadDataHolderClass[] readDataArray = new JNIReadDataHolderClass[readListSize]; + int idx = 0; + for(GATKSAMRecord read : reads) + { + readDataArray[idx] = new JNIReadDataHolderClass(); + readDataArray[idx].readBases = read.getReadBases(); + readDataArray[idx].readQuals = read.getBaseQualities(); + readDataArray[idx].insertionGOP = read.getBaseInsertionQualities(); + readDataArray[idx].deletionGOP = read.getBaseDeletionQualities(); + readDataArray[idx].overallGCP = GCPArrayMap.get(read); + + if(debug0_1) + System.out.println("Java read length "+readDataArray[idx].readBases.length); + if(debug3) + { + for(int i=0;i currEntry : alleleHaplotypeMap.entrySet()) //order is important - access in same order always + { + haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass(); + haplotypeDataArray[idx].haplotypeBases = currEntry.getValue().getBases(); + if(debug0_1) + System.out.println("Java haplotype length "+haplotypeDataArray[idx].haplotypeBases.length); + if(debug3) + { + for(int i=0;i currEntry : alleleHaplotypeMap.entrySet())//order is important - access in same order always + { + likelihoodMap.add(read, currEntry.getKey(), likelihoodArray[idx]); + ++idx; + } + if(debug) + { + //to compare values + likelihoodMap = super.computeLikelihoods(reads, alleleHaplotypeMap, GCPArrayMap); + //for floating point values, no exact equality + //check whether numbers are close in terms of abs_error or relative_error + //For very large values, relative_error is relevant + //For very small values, abs_error is relevant + boolean toDump = false; + for(int i=0;i 1e-5 && relative_error > 1e-5) + { + toDump = true; + break; + } + } + //if numbers are not close, then dump out the data that produced the inconsistency + if(toDump) + { + idx = 0; + System.out.println("Dump : Java numReads "+readListSize+" numHaplotypes "+alleleHaplotypeMapSize); + for(int i=0;i maxReadLength || haplotypeMaxLength > maxHaplotypeLength) - //{ initialize(readMaxLength, haplotypeMaxLength); } - - //if ( ! initialized ) - //throw new IllegalStateException("Must call initialize before calling computeReadLikelihoodGivenHaplotypeLog10"); - //final PerReadAlleleLikelihoodMap likelihoodMap = new PerReadAlleleLikelihoodMap(); - //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); - //} - //} - //return likelihoodMap; - //} - - ///** - //* {@inheritDoc} - //*/ - //@Override - //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) { - - //if ( haplotypeBases == null ) throw new IllegalArgumentException("haplotypeBases cannot be null"); - //if ( haplotypeBases.length > maxHaplotypeLength ) throw new IllegalArgumentException("Haplotype bases is too long, got " + haplotypeBases.length + " but max is " + maxHaplotypeLength); - //if ( readBases == null ) throw new IllegalArgumentException("readBases cannot be null"); - //if ( readBases.length > maxReadLength ) throw new IllegalArgumentException("readBases is too long, got " + readBases.length + " but max is " + maxReadLength); - //if ( readQuals.length != readBases.length ) throw new IllegalArgumentException("Read bases and read quals aren't the same size: " + readBases.length + " vs " + readQuals.length); - //if ( insertionGOP.length != readBases.length ) throw new IllegalArgumentException("Read bases and read insertion quals aren't the same size: " + readBases.length + " vs " + insertionGOP.length); - //if ( deletionGOP.length != readBases.length ) throw new IllegalArgumentException("Read bases and read deletion quals aren't the same size: " + readBases.length + " vs " + deletionGOP.length); - //if ( overallGCP.length != readBases.length ) throw new IllegalArgumentException("Read bases and overall GCP aren't the same size: " + readBases.length + " vs " + overallGCP.length); - - //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; - //} /** * {@inheritDoc} @@ -232,10 +296,11 @@ public class JNILoglessPairHMM extends LoglessPairHMM { //} //System.out.println("#### END STACK TRACE ####"); // + if(debug1) jniSubComputeReadLikelihoodGivenHaplotypeLog10(readBases.length, haplotypeBases.length, - readBases, haplotypeBases, readQuals, - insertionGOP, deletionGOP, overallGCP, - hapStartIndex); + readBases, haplotypeBases, readQuals, + insertionGOP, deletionGOP, overallGCP, + hapStartIndex); boolean doInitialization = (previousHaplotypeBases == null || previousHaplotypeBases.length != haplotypeBases.length); if (doInitialization) { @@ -310,6 +375,9 @@ public class JNILoglessPairHMM extends LoglessPairHMM { // initialize the pBaseReadLog10 matrix for all combinations of read x haplotype bases // Abusing the fact that java initializes arrays with 0.0, so no need to fill in rows and columns below 2. + if(debug3) + System.out.println("hapStartIndex "+startIndex); + for (int i = 0; i < readBases.length; i++) { final byte x = readBases[i]; final byte qual = readQuals[i]; 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 ff883c5ae..3b4498776 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java +++ b/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java @@ -57,6 +57,8 @@ public abstract class PairHMM { ORIGINAL, /* Optimized version of the PairHMM which caches per-read computations and operations in real space to avoid costly sums of log10'ed likelihoods */ LOGLESS_CACHING, + /* Optimized AVX implementation of LOGLESS_CACHING called through JNI */ + JNI_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 } @@ -70,6 +72,9 @@ public abstract class PairHMM { protected boolean doNotUseTristateCorrection = false; protected void doNotUseTristateCorrection() { doNotUseTristateCorrection = true; } + //debug array + protected double[] mLikelihoodArray; + /** * Initialize this PairHMM, making it suitable to run against a read and haplotype with given lengths * @@ -132,6 +137,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; for(GATKSAMRecord read : reads){ final byte[] readBases = read.getReadBases(); final byte[] readQuals = read.getBaseQualities(); @@ -144,12 +151,16 @@ public abstract class PairHMM { boolean isFirstHaplotype = true; Allele currentAllele = null; double log10l; - for (final Allele allele : alleleHaplotypeMap.keySet()){ - final Haplotype haplotype = alleleHaplotypeMap.get(allele); + //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 byte[] nextHaplotypeBases = haplotype.getBases(); if (currentHaplotypeBases != null) { log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, nextHaplotypeBases); + mLikelihoodArray[idx++] = log10l; likelihoodMap.add(read, currentAllele, log10l); } // update the current haplotype @@ -163,6 +174,7 @@ public abstract class PairHMM { log10l = computeReadLikelihoodGivenHaplotypeLog10(currentHaplotypeBases, readBases, readQuals, readInsQuals, readDelQuals, overallGCP, isFirstHaplotype, null); likelihoodMap.add(read, currentAllele, log10l); + mLikelihoodArray[idx++] = log10l; } } return likelihoodMap;