diff --git a/PairHMM_JNI/libJNILoglessPairHMM.so b/PairHMM_JNI/libJNILoglessPairHMM.so index 053b901cd..e18552b6b 100755 Binary files a/PairHMM_JNI/libJNILoglessPairHMM.so and b/PairHMM_JNI/libJNILoglessPairHMM.so differ diff --git a/PairHMM_JNI/pairhmm-1-base.cc b/PairHMM_JNI/pairhmm-1-base.cc index 3941c92de..a2d5dfa88 100644 --- a/PairHMM_JNI/pairhmm-1-base.cc +++ b/PairHMM_JNI/pairhmm-1-base.cc @@ -58,47 +58,53 @@ int main(int argc, char** argv) vector tc_vector; tc_vector.clear(); testcase tc; + double total_time = 0; while(1) { int break_value = use_old_read_testcase ? read_testcase(&tc, fptr) : read_mod_testcase(ifptr,&tc,true); + if(break_value >= 0) + tc_vector.push_back(tc); + if(tc_vector.size() == BATCH_SIZE || (break_value < 0 && tc_vector.size() > 0)) + { + vector results_vec; + results_vec.clear(); + results_vec.resize(tc_vector.size()); + double start_time = getCurrClk(); +#pragma omp parallel for schedule(dynamic,chunk_size) num_threads(12) + for(unsigned i=0;i(&tc); + baseline_result = log10(baseline_result) - log10(ldexp(1.0, 1020.0)); + double abs_error = fabs(baseline_result-results_vec[i]); + double rel_error = (baseline_result != 0) ? fabs(abs_error/baseline_result) : 0; + if(abs_error > 1e-5 && rel_error > 1e-5) + cout << std::scientific << baseline_result << " "< results_vec; - results_vec.clear(); - results_vec.resize(tc_vector.size()); - double start_time = getCurrClk(); -#pragma omp parallel for schedule(dynamic,chunk_size) num_threads(12) - for(unsigned i=0;i(&tc); - baseline_result = log10(baseline_result) - log10(ldexp(1.0, 1020.0)); - double abs_error = fabs(baseline_result-results_vec[i]); - double rel_error = (baseline_result != 0) ? fabs(abs_error/baseline_result) : 0; - if(abs_error > 1e-5 && rel_error > 1e-5) - cout << std::scientific << baseline_result << " "< void GEN_INTRINSIC(GEN_INTRINSIC(initializeVectors,SIMD_TYPE), PRECISION)(int ROWS, int COLS, NUMBER* shiftOutM, NUMBER *shiftOutX, NUMBER *shiftOutY, Context ctx, testcase *tc, _256_TYPE *p_MM, _256_TYPE *p_GAPM, _256_TYPE *p_MX, _256_TYPE *p_XX, _256_TYPE *p_MY, _256_TYPE *p_YY, _256_TYPE *distm1D); -template void GEN_INTRINSIC(GEN_INTRINSIC(stripINITIALIZATION,SIMD_TYPE), PRECISION)( +inline void GEN_INTRINSIC(GEN_INTRINSIC(_vector_shift,SIMD_TYPE), PRECISION) (UNION_TYPE &x, MAIN_TYPE shiftIn, MAIN_TYPE &shiftOut); +inline void GEN_INTRINSIC(GEN_INTRINSIC(_vector_shift_last,SIMD_TYPE), PRECISION) (UNION_TYPE &x, MAIN_TYPE shiftIn); +inline void GEN_INTRINSIC(GEN_INTRINSIC(precompute_masks_,SIMD_TYPE), PRECISION)(const testcase& tc, int COLS, int numMaskVecs, MASK_TYPE (*maskArr)[NUM_DISTINCT_CHARS]); +inline void GEN_INTRINSIC(GEN_INTRINSIC(init_masks_for_row_,SIMD_TYPE), PRECISION)(const testcase& tc, char* rsArr, MASK_TYPE* lastMaskShiftOut, int beginRowIndex, int numRowsToProcess); +inline void GEN_INTRINSIC(GEN_INTRINSIC(update_masks_for_cols_,SIMD_TYPE), PRECISION)(int maskIndex, MASK_VEC& currMaskVecLow, MASK_VEC& currMaskVecHigh, MASK_TYPE (*maskArr) [NUM_DISTINCT_CHARS], char* rsArr, MASK_TYPE* lastMaskShiftOut, MASK_TYPE maskBitCnt); +inline void GEN_INTRINSIC(GEN_INTRINSIC(computeDistVec,SIMD_TYPE), PRECISION) (MASK_VEC& currMaskVecLow, MASK_VEC& currMaskVecHigh, _256_TYPE& distm, _256_TYPE& _1_distm, _256_TYPE& distmChosen); +template inline void GEN_INTRINSIC(GEN_INTRINSIC(initializeVectors,SIMD_TYPE), PRECISION)(int ROWS, int COLS, NUMBER* shiftOutM, NUMBER *shiftOutX, NUMBER *shiftOutY, Context ctx, testcase *tc, _256_TYPE *p_MM, _256_TYPE *p_GAPM, _256_TYPE *p_MX, _256_TYPE *p_XX, _256_TYPE *p_MY, _256_TYPE *p_YY, _256_TYPE *distm1D); +template inline void GEN_INTRINSIC(GEN_INTRINSIC(stripINITIALIZATION,SIMD_TYPE), PRECISION)( int stripIdx, Context ctx, testcase *tc, _256_TYPE &pGAPM, _256_TYPE &pMM, _256_TYPE &pMX, _256_TYPE &pXX, _256_TYPE &pMY, _256_TYPE &pYY, _256_TYPE &rs, UNION_TYPE &rsN, _256_TYPE &distm, _256_TYPE &_1_distm, _256_TYPE *distm1D, _256_TYPE N_packed256, _256_TYPE *p_MM , _256_TYPE *p_GAPM , _256_TYPE *p_MX, _256_TYPE *p_XX , _256_TYPE *p_MY, _256_TYPE *p_YY, UNION_TYPE &M_t_2, UNION_TYPE &X_t_2, UNION_TYPE &M_t_1, UNION_TYPE &X_t_1, UNION_TYPE &Y_t_2, UNION_TYPE &Y_t_1, UNION_TYPE &M_t_1_y, NUMBER* shiftOutX, NUMBER* shiftOutM); -_256_TYPE GEN_INTRINSIC(GEN_INTRINSIC(computeDISTM,SIMD_TYPE), PRECISION)(int d, int COLS, testcase * tc, HAP_TYPE &hap, _256_TYPE rs, UNION_TYPE rsN, _256_TYPE N_packed256, +inline _256_TYPE GEN_INTRINSIC(GEN_INTRINSIC(computeDISTM,SIMD_TYPE), PRECISION)(int d, int COLS, testcase * tc, HAP_TYPE &hap, _256_TYPE rs, UNION_TYPE rsN, _256_TYPE N_packed256, _256_TYPE distm, _256_TYPE _1_distm); -void GEN_INTRINSIC(GEN_INTRINSIC(computeMXY,SIMD_TYPE), PRECISION)(UNION_TYPE &M_t, UNION_TYPE &X_t, UNION_TYPE &Y_t, UNION_TYPE &M_t_y, +inline void GEN_INTRINSIC(GEN_INTRINSIC(computeMXY,SIMD_TYPE), PRECISION)(UNION_TYPE &M_t, UNION_TYPE &X_t, UNION_TYPE &Y_t, UNION_TYPE &M_t_y, UNION_TYPE M_t_2, UNION_TYPE X_t_2, UNION_TYPE Y_t_2, UNION_TYPE M_t_1, UNION_TYPE X_t_1, UNION_TYPE M_t_1_y, UNION_TYPE Y_t_1, _256_TYPE pMM, _256_TYPE pGAPM, _256_TYPE pMX, _256_TYPE pXX, _256_TYPE pMY, _256_TYPE pYY, _256_TYPE distmSel); template NUMBER GEN_INTRINSIC(GEN_INTRINSIC(compute_full_prob_,SIMD_TYPE), PRECISION) (testcase *tc, NUMBER *before_last_log = NULL); 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 5ab49b531..caa880b41 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 @@ -332,7 +332,7 @@ public class PairHMMLikelihoodCalculationEngine implements LikelihoodCalculation } // initialize arrays to hold the probabilities of being in the match, insertion and deletion cases - pairHMMThreadLocal.get().initialize(X_METRIC_LENGTH, Y_METRIC_LENGTH); + pairHMMThreadLocal.get().initialize(haplotypes, perSampleReadList, X_METRIC_LENGTH, Y_METRIC_LENGTH); } 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 ca3db5ed3..ed628f064 100644 --- a/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java +++ b/protected/java/src/org/broadinstitute/sting/utils/pairhmm/JNILoglessPairHMM.java @@ -57,6 +57,7 @@ import org.broadinstitute.variant.variantcontext.Allele; import java.util.List; import java.util.Map; +import java.util.HashMap; /** @@ -147,6 +148,52 @@ public class JNILoglessPairHMM extends LoglessPairHMM { jniInitialize(readMaxLength, haplotypeMaxLength); } + private native void jniInitialize(final int numHaplotypes, final int numTotalReads, JNIHaplotypeDataHolderClass[] haplotypeDataArray, + JNIReadDataHolderClass[] readDataArray); + HashMap sampleToIndex = new HashMap(); + //Used to transfer data to JNI + /** + * {@inheritDoc} + */ + @Override + public void initialize( final List haplotypes, final Map> perSampleReadList, final int readMaxLength, final int haplotypeMaxLength ) { + if(verify) + super.initialize(haplotypes, perSampleReadList, readMaxLength, haplotypeMaxLength); + /* + int numHaplotypes = haplotypes.size(); + JNIHaplotypeDataHolderClass[] haplotypeDataArray = new JNIHaplotypeDataHolderClass[numHaplotypes]; + int idx = 0; + for(final Haplotype currHaplotype : haplotypes) + { + haplotypeDataArray[idx] = new JNIHaplotypeDataHolderClass(); + haplotypeDataArray[idx].haplotypeBases = currHaplotype.getBases(); + ++idx; + } + sampleToIndex.clear(); + int totalNumReads = 0; + for(final Map.Entry> currEntry : perSampleReadList) + { + sampleToIndex.put(currEntry.getKey(), totalNumReads); + totalNumReads += currEntry.getValue().size(); + } + + JNIHaplotypeDataHolderClass[] readDataArray = new JNIReadDataHolderClass[totalNumReads]; + idx = 0; + for(final Map.Entry> currEntry : perSampleReadList) + { + for(GATKSAMRecord read : currEntry.getValue()) + { + 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); + ++idx; + } + } + jniInitialize(numHaplotypes, numTotalReads, haplotypeDataArray, readDataArray);*/ + } //Real compute kernel private native void jniComputeLikelihoods(int numReads, int numHaplotypes, 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 7b09fe047..08a49cfb7 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java +++ b/public/java/src/org/broadinstitute/sting/utils/pairhmm/PairHMM.java @@ -100,6 +100,18 @@ public abstract class PairHMM { initialized = true; } + /** + * Initialize this PairHMM, making it suitable to run against a read and haplotype with given lengths + * This function is used by the JNI implementations to transfer all data once to the native code + * @param haplotypes the list of haplotypes + * @param perSampleReadList map from sample name to list of reads + * @param haplotypeMaxLength the max length of haplotypes we want to use with this 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); + } + protected int findMaxReadLength(final List reads) { int listMaxReadLength = 0; for(GATKSAMRecord read : reads){