From 51f9ec0a5c35335485c7b127c5f480839db16493 Mon Sep 17 00:00:00 2001 From: ebanks Date: Sun, 18 Oct 2009 05:20:15 +0000 Subject: [PATCH] subtract largest posterior value from all values; this hopefully solves any precision issues git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1870 348d0f76-0448-11de-a6fe-93d51630548a --- .../gatk/walkers/genotyper/GenotypeLikelihoods.java | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoods.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoods.java index b0d64ef43..5e24128ac 100644 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoods.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoods.java @@ -6,6 +6,7 @@ import org.broadinstitute.sting.utils.genotype.DiploidGenotype; import static java.lang.Math.log10; import static java.lang.Math.pow; +import java.util.Arrays; /** * Stable, error checking version of the Bayesian genotyper. Useful for calculating the likelihoods, priors, @@ -158,9 +159,17 @@ public abstract class GenotypeLikelihoods implements Cloneable { double[] normalized = new double[posteriors.length]; double sum = 0.0; + // for precision purposes, we need to add (or really subtract, since everything is negative) + // the largest posterior value from all entries so that numbers don't get too small + double maxValue = posteriors[0]; + for (int i = 1; i < posteriors.length; i++) { + if ( maxValue < posteriors[i] ) + maxValue = posteriors[i]; + } + // collect the posteriors for ( DiploidGenotype g : DiploidGenotype.values() ) { - double posterior = Math.pow(10, getPosterior(g)); + double posterior = Math.pow(10, getPosterior(g) - maxValue); normalized[g.ordinal()] = posterior; sum += posterior; }