Skip to content

Commit

Permalink
add tests to math utils. 90% now
Browse files Browse the repository at this point in the history
  • Loading branch information
@akiezun
akiezun committed Jun 8, 2015
1 parent 6ba4787 commit 6a9cfc6
Show file tree
Hide file tree
Showing 5 changed files with 278 additions and 186 deletions.
48 changes: 32 additions & 16 deletions src/main/java/org/broadinstitute/hellbender/utils/MathUtils.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,22 @@ public static double log10BinomialCoefficient(final int n, final int k) {
return log10Factorial(n) - log10Factorial(k) - log10Factorial(n - k);
}

/**
* Computes a binomial probability. This is computed using the formula
* <p/>
* B(k; n; p) = [ n! / ( k! (n - k)! ) ] (p^k)( (1-p)^k )
* <p/>
* where n is the number of trials, k is the number of successes, and p is the probability of success
*
* @param n number of Bernoulli trials
* @param k number of successes
* @param p probability of success
* @return the binomial probability of the specified configuration. Computes values down to about 1e-237.
*/
public static double binomialProbability(final int n, final int k, final double p) {
return Math.pow(10, log10BinomialProbability(n, k, Math.log10(p)));
}

/**
* binomial Probability(int, int, double) with log10 applied to result
*/
Expand All @@ -253,21 +269,21 @@ public static double log10sumLog10(final double[] log10values) {


public static double log10sumLog10(final double[] log10p, final int start, final int finish) {

if (start >= finish)
if (start >= finish) {
return Double.NEGATIVE_INFINITY;
final int maxElementIndex = MathUtils.maxElementIndex(log10p, start, finish);
}
final int maxElementIndex = maxElementIndex(log10p, start, finish);
final double maxValue = log10p[maxElementIndex];
if(maxValue == Double.NEGATIVE_INFINITY)
if(maxValue == Double.NEGATIVE_INFINITY) {
return maxValue;
}
double sum = 1.0;
for (int i = start; i < finish; i++) {
double curVal = log10p[i];
double scaled_val = curVal - maxValue;
if (i == maxElementIndex || curVal == Double.NEGATIVE_INFINITY) {
continue;
}
else {
} else {
sum += Math.pow(10.0, scaled_val);
}
}
Expand All @@ -277,6 +293,16 @@ public static double log10sumLog10(final double[] log10p, final int start, final
return maxValue + (sum != 1.0 ? Math.log10(sum) : 0.0);
}

/**
* normalizes the log10-based array. ASSUMES THAT ALL ARRAY ENTRIES ARE <= 0 (<= 1 IN REAL-SPACE).
*
* @param array the array to be normalized
* @return a newly allocated array corresponding the normalized values in array
*/
public static double[] normalizeFromLog10(final double[] array) {
return normalizeFromLog10(array, false);
}

/**
* normalizes the log10-based array. ASSUMES THAT ALL ARRAY ENTRIES ARE <= 0 (<= 1 IN REAL-SPACE).
*
Expand Down Expand Up @@ -334,16 +360,6 @@ public static double[] normalizeFromLog10(final double[] array, final boolean ta
return normalized;
}

/**
* normalizes the log10-based array. ASSUMES THAT ALL ARRAY ENTRIES ARE <= 0 (<= 1 IN REAL-SPACE).
*
* @param array the array to be normalized
* @return a newly allocated array corresponding the normalized values in array
*/
public static double[] normalizeFromLog10(final double[] array) {
return normalizeFromLog10(array, false);
}

/**
* normalizes the real-space probability array.
*
Expand Down
22 changes: 1 addition & 21 deletions src/main/java/org/broadinstitute/hellbender/utils/pairhmm/LoglessPairHMM.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

import static org.broadinstitute.hellbender.utils.pairhmm.PairHMMModel.*;

public class LoglessPairHMM extends N2MemoryPairHMM {
public final class LoglessPairHMM extends N2MemoryPairHMM {
protected static final double INITIAL_CONDITION = Math.pow(2, 1020);
protected static final double INITIAL_CONDITION_LOG10 = Math.log10(INITIAL_CONDITION);

Expand Down Expand Up @@ -102,24 +102,4 @@ protected void initializePriors(final byte[] haplotypeBases, final byte[] readBa
protected static void initializeProbabilities(final double[][] transition, final byte[] insertionGOP, final byte[] deletionGOP, final byte[] overallGCP) {
PairHMMModel.qualToTransProbs(transition,insertionGOP,deletionGOP,overallGCP);
}

/**
* Updates a cell in the HMM matrix
*
* The read and haplotype indices are offset by one because the state arrays have an extra column to hold the
* initial conditions
* @param indI row index in the matrices to update
* @param indJ column index in the matrices to update
* @param prior the likelihood editing distance matrix for the read x haplotype
* @param transition an array with the six transition relevant to this location
*/
protected void updateCell( final int indI, final int indJ, final double prior, final double[] transition) {

matchMatrix[indI][indJ] = prior * ( matchMatrix[indI - 1][indJ - 1] * transition[matchToMatch] +
insertionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] +
deletionMatrix[indI - 1][indJ - 1] * transition[indelToMatch] );
insertionMatrix[indI][indJ] = matchMatrix[indI - 1][indJ] * transition[matchToInsertion] + insertionMatrix[indI - 1][indJ] * transition[insertionToInsertion];
deletionMatrix[indI][indJ] = matchMatrix[indI][indJ - 1] * transition[matchToDeletion] + deletionMatrix[indI][indJ - 1] * transition[deletionToDeletion];
}
}
16 changes: 0 additions & 16 deletions src/main/java/org/broadinstitute/hellbender/utils/pairhmm/PairHMM.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,17 +115,6 @@ static int findMaxReadLength(final List<SAMRecord> reads) {
return listMaxReadLength;
}

static int findMaxHaplotypeLength(final Collection<Haplotype> haplotypes) {
int listMaxHaplotypeLength = 0;
for(final Haplotype h : haplotypes) {
final int haplotypeLength = h.getBases().length;
if( haplotypeLength > listMaxHaplotypeLength ) {
listMaxHaplotypeLength = haplotypeLength;
}
}
return listMaxHaplotypeLength;
}

/**
* Given a list of reads and haplotypes, for every read compute the total probability of said read arising from
* each haplotype given base substitution, insertion, and deletion probabilities.
Expand Down Expand Up @@ -294,11 +283,6 @@ public static int findFirstPositionWhereHaplotypesDiffer(final byte[] haplotype1
return Math.min(haplotype1.length, haplotype2.length);
}

/**
* Return the results of the computeLikelihoods function
*/
public double[] getLikelihoodArray() { return mLikelihoodArray; }

/**
* Called at the end of the program to close files, print profiling information etc
*/
Expand Down
6 changes: 0 additions & 6 deletions src/main/java/org/broadinstitute/hellbender/utils/pairhmm/PairHMMModel.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,6 @@ public static void qualToTransProbs(final double[] dest, final byte insQual, fin
*
* @return never {@code null}. An array of length {@link #TRANS_PROB_ARRAY_LENGTH}.
*/
@SuppressWarnings("unused")
public static double[] qualToTransProbs(final byte insQual, final byte delQual, final byte gcp) {
final double[] dest = new double[TRANS_PROB_ARRAY_LENGTH];
qualToTransProbs(dest,insQual,delQual,gcp);
Expand All @@ -154,7 +153,6 @@ public static double[] qualToTransProbs(final byte insQual, final byte delQual,
* @throws ArrayIndexOutOfBoundsException if {@code dest} or any of its elements is not large enough to contain the
* transition matrix.
*/
@SuppressWarnings("unused")
public static void qualToTransProbs(final double[][] dest, final byte[] insQuals, final byte[] delQuals, final byte[] gcps) {
final int readLength = insQuals.length;
if (delQuals.length != readLength) throw new IllegalArgumentException("deletion quality array length does not match insert quality array length: " + readLength + " != " + delQuals.length);
Expand Down Expand Up @@ -185,7 +183,6 @@ public static void qualToTransProbs(final double[][] dest, final byte[] insQuals
*
* @return never {@code null}, an matrix of the dimensions explained above.
*/
@SuppressWarnings("unused")
public static double[][] qualToTransProbs(final byte[] insQuals, final byte[] delQuals, final byte[] gcps) {
final double[][] dest = createTransitionMatrix(insQuals.length);
qualToTransProbs(dest,insQuals,delQuals,gcps);
Expand Down Expand Up @@ -227,7 +224,6 @@ public static void qualToTransProbsLog10(final double[] dest, final byte insQual
*
* @return never {@code null}. An array of length {@link #TRANS_PROB_ARRAY_LENGTH}.
*/
@SuppressWarnings("unused")
public static double[] qualToTransProbsLog10(final byte insQual, final byte delQual, final byte gcp) {
final double[] dest = new double[TRANS_PROB_ARRAY_LENGTH];
qualToTransProbsLog10(dest,insQual,delQual,gcp);
Expand All @@ -253,7 +249,6 @@ public static double[] qualToTransProbsLog10(final byte insQual, final byte delQ
* @throws ArrayIndexOutOfBoundsException if {@code dest} or any of its elements is not large enough to contain the
* transition matrix.
*/
@SuppressWarnings("unused")
public static void qualToTransProbsLog10(final double[][] dest, final byte[] insQuals, final byte[] delQuals, final byte[] gcps) {
final int readLength = insQuals.length;
if (delQuals.length != readLength) throw new IllegalArgumentException("deletion quality array length does not match insert quality array length: " + readLength + " != " + delQuals.length);
Expand Down Expand Up @@ -284,7 +279,6 @@ public static void qualToTransProbsLog10(final double[][] dest, final byte[] ins
*
* @return never {@code null}, an matrix of the dimensions explained above.
*/
@SuppressWarnings("unused")
public static double[][] qualToTransProbsLog10(final byte[] insQuals, final byte[] delQuals, final byte[] gcps) {
final double[][] dest = createTransitionMatrix(insQuals.length);
qualToTransProbsLog10(dest,insQuals,delQuals,gcps);
Expand Down

0 comments on commit 6a9cfc6

Please sign in to comment.