Examples with Max org.apache.commons.math3.stat.descriptive.rank.Max used on opensource projects

Example 21 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project gatk-protected by broadinstitute.

the class RobustBrentSolver method doSolve.

@Override
protected double doSolve() throws TooManyEvaluationsException, NoBracketingException {
    final double min = getMin();
    final double max = getMax();
    final double[] xSearchGrid = createHybridSearchGrid(min, max, numBisections, depth);
    final double[] fSearchGrid = Arrays.stream(xSearchGrid).map(this::computeObjectiveValue).toArray();
    /* find bracketing intervals on the search grid */
    final List<Bracket> bracketsList = detectBrackets(xSearchGrid, fSearchGrid);
    if (bracketsList.isEmpty()) {
        throw new NoBracketingException(min, max, fSearchGrid[0], fSearchGrid[fSearchGrid.length - 1]);
    }
    final BrentSolver solver = new BrentSolver(getRelativeAccuracy(), getAbsoluteAccuracy(), getFunctionValueAccuracy());
    final List<Double> roots = bracketsList.stream().map(b -> solver.solve(getMaxEvaluations(), this::computeObjectiveValue, b.min, b.max, 0.5 * (b.min + b.max))).collect(Collectors.toList());
    if (roots.size() == 1 || meritFunc == null) {
        return roots.get(0);
    }
    final double[] merits = roots.stream().mapToDouble(meritFunc::value).toArray();
    final int bestRootIndex = IntStream.range(0, roots.size()).boxed().max((i, j) -> (int) (merits[i] - merits[j])).get();
    return roots.get(bestRootIndex);
}

Example 22 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project gatk-protected by broadinstitute.

the class CoverageModelParameters method generateRandomModel.

/**
     * Generates random coverage model parameters.
     *
     * @param targetList list of targets
     * @param numLatents number of latent variables
     * @param seed random seed
     * @param randomMeanLogBiasStandardDeviation std of mean log bias (mean is set to 0)
     * @param randomBiasCovariatesStandardDeviation std of bias covariates (mean is set to 0)
     * @param randomMaxUnexplainedVariance max value of unexplained variance (samples are taken from a uniform
     *                                     distribution [0, {@code randomMaxUnexplainedVariance}])
     * @param initialBiasCovariatesARDCoefficients initial row vector of ARD coefficients
     * @return an instance of {@link CoverageModelParameters}
     */
public static CoverageModelParameters generateRandomModel(final List<Target> targetList, final int numLatents, final long seed, final double randomMeanLogBiasStandardDeviation, final double randomBiasCovariatesStandardDeviation, final double randomMaxUnexplainedVariance, final INDArray initialBiasCovariatesARDCoefficients) {
    Utils.validateArg(numLatents >= 0, "Dimension of the bias space must be non-negative");
    Utils.validateArg(randomBiasCovariatesStandardDeviation >= 0, "Standard deviation of random bias covariates" + " must be non-negative");
    Utils.validateArg(randomMeanLogBiasStandardDeviation >= 0, "Standard deviation of random mean log bias" + " must be non-negative");
    Utils.validateArg(randomMaxUnexplainedVariance >= 0, "Max random unexplained variance must be non-negative");
    Utils.validateArg(initialBiasCovariatesARDCoefficients == null || numLatents > 0 && initialBiasCovariatesARDCoefficients.length() == numLatents, "If ARD is enabled, the dimension" + " of the bias latent space must be positive and match the length of ARD coeffecient vector");
    final boolean biasCovariatesEnabled = numLatents > 0;
    final int numTargets = targetList.size();
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(seed));
    /* Gaussian random for mean log bias */
    final INDArray initialMeanLogBias = Nd4j.create(getNormalRandomNumbers(numTargets, 0, randomMeanLogBiasStandardDeviation, rng), new int[] { 1, numTargets });
    /* Uniform random for unexplained variance */
    final INDArray initialUnexplainedVariance = Nd4j.create(getUniformRandomNumbers(numTargets, 0, randomMaxUnexplainedVariance, rng), new int[] { 1, numTargets });
    final INDArray initialMeanBiasCovariates;
    if (biasCovariatesEnabled) {
        /* Gaussian random for bias covariates */
        initialMeanBiasCovariates = Nd4j.create(getNormalRandomNumbers(numTargets * numLatents, 0, randomBiasCovariatesStandardDeviation, rng), new int[] { numTargets, numLatents });
    } else {
        initialMeanBiasCovariates = null;
    }
    return new CoverageModelParameters(targetList, initialMeanLogBias, initialUnexplainedVariance, initialMeanBiasCovariates, initialBiasCovariatesARDCoefficients);
}

Example 23 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project gatk by broadinstitute.

the class RobustBrentSolver method doSolve.

@Override
protected double doSolve() throws TooManyEvaluationsException, NoBracketingException {
    final double min = getMin();
    final double max = getMax();
    final double[] xSearchGrid = createHybridSearchGrid(min, max, numBisections, depth);
    final double[] fSearchGrid = Arrays.stream(xSearchGrid).map(this::computeObjectiveValue).toArray();
    /* find bracketing intervals on the search grid */
    final List<Bracket> bracketsList = detectBrackets(xSearchGrid, fSearchGrid);
    if (bracketsList.isEmpty()) {
        throw new NoBracketingException(min, max, fSearchGrid[0], fSearchGrid[fSearchGrid.length - 1]);
    }
    final BrentSolver solver = new BrentSolver(getRelativeAccuracy(), getAbsoluteAccuracy(), getFunctionValueAccuracy());
    final List<Double> roots = bracketsList.stream().map(b -> solver.solve(getMaxEvaluations(), this::computeObjectiveValue, b.min, b.max, 0.5 * (b.min + b.max))).collect(Collectors.toList());
    if (roots.size() == 1 || meritFunc == null) {
        return roots.get(0);
    }
    final double[] merits = roots.stream().mapToDouble(meritFunc::value).toArray();
    final int bestRootIndex = IntStream.range(0, roots.size()).boxed().max((i, j) -> (int) (merits[i] - merits[j])).get();
    return roots.get(bestRootIndex);
}

Example 24 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project gatk-protected by broadinstitute.

the class HetPulldownCalculator method isPileupHetCompatible.

/**
     * Returns true if the distribution of major and other base-pair counts from a pileup at a locus is compatible with
     * allele fraction of 0.5.
     *
     * <p>
     *     Compatibility is defined by a p-value threshold.  That is, compute the two-sided p-value of observing
     *     a number of major read counts out of a total number of reads, assuming the given heterozygous
     *     allele fraction.  If the p-value is less than the given threshold, then reject the null hypothesis
     *     that the heterozygous allele fraction is 0.5 (i.e., SNP is likely to be homozygous) and return false,
     *     otherwise return true.
     * </p>
     * @param baseCounts        base-pair counts
     * @param totalBaseCount    total base-pair counts (excluding N, etc.)
     * @param pvalThreshold     p-value threshold for two-sided binomial test (should be in [0, 1], but no check is performed)
     * @return                  boolean compatibility with heterozygous allele fraction
     */
@VisibleForTesting
protected static boolean isPileupHetCompatible(final Nucleotide.Counter baseCounts, final int totalBaseCount, final double pvalThreshold) {
    final int majorReadCount = Arrays.stream(BASES).mapToInt(b -> (int) baseCounts.get(b)).max().getAsInt();
    if (majorReadCount == 0 || totalBaseCount - majorReadCount == 0) {
        return false;
    }
    final double pval = new BinomialTest().binomialTest(totalBaseCount, majorReadCount, HET_ALLELE_FRACTION, AlternativeHypothesis.TWO_SIDED);
    return pval >= pvalThreshold;
}

Example 25 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project gatk-protected by broadinstitute.

the class SomaticLikelihoodsEngineUnitTest method testEvidence.

@Test
public void testEvidence() {
    // one exact limit for the evidence is when the likelihoods of each read are so peaked (i.e. the most likely allele
    // of each read is much likelier than all other alleles) that the sum over latent read-to-allele assignments
    // (that is, over the indicator z in the notes) is dominated by the max-likelihood allele configuration
    // and thus the evidence reduces to exactly integrating out the Dirichlet allele fractions
    final double[] prior = new double[] { 1, 2 };
    final RealMatrix log10Likelihoods = new Array2DRowRealMatrix(2, 4);
    log10Likelihoods.setRow(0, new double[] { 0.1, 4.0, 3.0, -10 });
    log10Likelihoods.setRow(1, new double[] { -12, -9, -5.0, 0.5 });
    final double calculatedLog10Evidence = SomaticLikelihoodsEngine.log10Evidence(log10Likelihoods, prior);
    final double[] maxLikelihoodCounts = new double[] { 3, 1 };
    final double expectedLog10Evidence = SomaticLikelihoodsEngine.log10DirichletNormalization(prior) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior, maxLikelihoodCounts)) + new IndexRange(0, log10Likelihoods.getColumnDimension()).sum(read -> log10Likelihoods.getColumnVector(read).getMaxValue());
    Assert.assertEquals(calculatedLog10Evidence, expectedLog10Evidence, 1e-5);
    // when there's just one read we can calculate the likelihood exactly
    final double[] prior2 = new double[] { 1, 2 };
    final RealMatrix log10Likelihoods2 = new Array2DRowRealMatrix(2, 1);
    log10Likelihoods2.setRow(0, new double[] { 0.1 });
    log10Likelihoods2.setRow(1, new double[] { 0.5 });
    final double calculatedLog10Evidence2 = SomaticLikelihoodsEngine.log10Evidence(log10Likelihoods2, prior2);
    final double[] delta0 = new double[] { 1, 0 };
    final double[] delta1 = new double[] { 0, 1 };
    final double expectedLog10Evidence2 = MathUtils.log10SumLog10(log10Likelihoods2.getEntry(0, 0) + SomaticLikelihoodsEngine.log10DirichletNormalization(prior2) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior2, delta0)), +log10Likelihoods2.getEntry(1, 0) + SomaticLikelihoodsEngine.log10DirichletNormalization(prior2) - SomaticLikelihoodsEngine.log10DirichletNormalization(MathArrays.ebeAdd(prior2, delta1)));
    Assert.assertEquals(calculatedLog10Evidence2, expectedLog10Evidence2, 0.05);
}