Examples with GaussIntegrator org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator used on opensource projects

Example 1 with GaussIntegrator

use of org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator in project gatk by broadinstitute.

the class IntegrationUtils method integrate2d.

public static double integrate2d(final ToDoubleBiFunction<Double, Double> getIntegrand, final double xLowerBound, final double xUpperBound, final int xNumPoints, final double yLowerBound, final double yUpperBound, final int yNumPoints) {
    final GaussIntegrator xIntegrator = integratorFactory.legendre(xNumPoints, xLowerBound, xUpperBound);
    final GaussIntegrator yIntegrator = integratorFactory.legendre(yNumPoints, yLowerBound, yUpperBound);
    final double[] xIntegrationWeights = new IndexRange(0, xNumPoints).mapToDouble(xIntegrator::getWeight);
    final double[] xAbscissas = new IndexRange(0, xNumPoints).mapToDouble(xIntegrator::getPoint);
    final double[] yIntegrationWeights = new IndexRange(0, yNumPoints).mapToDouble(yIntegrator::getWeight);
    final double[] yAbscissas = new IndexRange(0, yNumPoints).mapToDouble(yIntegrator::getPoint);
    double integral = 0;
    for (int i = 0; i < xNumPoints; i++) {
        final double x = xAbscissas[i];
        for (int j = 0; j < yNumPoints; j++) {
            final double y = yAbscissas[j];
            final double integrand = getIntegrand.applyAsDouble(x, y);
            integral += xIntegrationWeights[i] * yIntegrationWeights[j] * integrand;
        }
    }
    return integral;
}

Example 2 with GaussIntegrator

use of org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator in project gatk by broadinstitute.

the class IntegrationUtils method integrate.

public static double integrate(final DoubleFunction<Double> getIntegrand, final double lowerBound, final double upperBound, final int numPoints) {
    final GaussIntegrator integrator = integratorFactory.legendre(numPoints, lowerBound, upperBound);
    final double[] gaussIntegrationWeights = new IndexRange(0, numPoints).mapToDouble(integrator::getWeight);
    final double[] gaussIntegrationAbscissas = new IndexRange(0, numPoints).mapToDouble(integrator::getPoint);
    final double[] integrands = MathUtils.applyToArrayInPlace(gaussIntegrationAbscissas, getIntegrand::apply);
    return GATKProtectedMathUtils.dotProduct(gaussIntegrationWeights, integrands);
}

Example 3 with GaussIntegrator

use of org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator in project gatk-protected by broadinstitute.

the class HeterogeneousHeterozygousPileupPriorModel method initializeIntegrationQuadrature.

/**
     * Initilizes the quadrature for calculating allele ratio integrals in
     * {@link HeterogeneousHeterozygousPileupPriorModel#getHetLogLikelihood(List)}
     *
     * @param numIntegPoints  number of points in the quadrature
     */
private void initializeIntegrationQuadrature(final int numIntegPoints) {
    /* get Gauss-Legendre quadrature factory of order @numIntegPoints */
    final GaussIntegratorFactory integratorFactory = new GaussIntegratorFactory();
    final GaussIntegrator gaussIntegrator = integratorFactory.legendre(numIntegPoints, minHetAlleleFraction, 1.0 - minHetAlleleFraction);
    /* abscissas */
    gaussIntegrationAbscissas.clear();
    gaussIntegrationAbscissas.addAll(IntStream.range(0, numIntegPoints).mapToDouble(gaussIntegrator::getPoint).boxed().collect(Collectors.toList()));
    /* weights */
    gaussIntegrationWeights.clear();
    gaussIntegrationWeights.addAll(IntStream.range(0, numIntegPoints).mapToDouble(gaussIntegrator::getWeight).boxed().collect(Collectors.toList()));
    /* log of weights */
    gaussIntegrationLogWeights.clear();
    gaussIntegrationLogWeights.addAll(gaussIntegrationWeights.stream().mapToDouble(FastMath::log).boxed().collect(Collectors.toList()));
}

Example 4 with GaussIntegrator

use of org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator in project gatk by broadinstitute.

the class HeterogeneousHeterozygousPileupPriorModel method initializeIntegrationQuadrature.

/**
     * Initilizes the quadrature for calculating allele ratio integrals in
     * {@link HeterogeneousHeterozygousPileupPriorModel#getHetLogLikelihood(List)}
     *
     * @param numIntegPoints  number of points in the quadrature
     */
private void initializeIntegrationQuadrature(final int numIntegPoints) {
    /* get Gauss-Legendre quadrature factory of order @numIntegPoints */
    final GaussIntegratorFactory integratorFactory = new GaussIntegratorFactory();
    final GaussIntegrator gaussIntegrator = integratorFactory.legendre(numIntegPoints, minHetAlleleFraction, 1.0 - minHetAlleleFraction);
    /* abscissas */
    gaussIntegrationAbscissas.clear();
    gaussIntegrationAbscissas.addAll(IntStream.range(0, numIntegPoints).mapToDouble(gaussIntegrator::getPoint).boxed().collect(Collectors.toList()));
    /* weights */
    gaussIntegrationWeights.clear();
    gaussIntegrationWeights.addAll(IntStream.range(0, numIntegPoints).mapToDouble(gaussIntegrator::getWeight).boxed().collect(Collectors.toList()));
    /* log of weights */
    gaussIntegrationLogWeights.clear();
    gaussIntegrationLogWeights.addAll(gaussIntegrationWeights.stream().mapToDouble(FastMath::log).boxed().collect(Collectors.toList()));
}

Example 5 with GaussIntegrator

use of org.apache.commons.math3.analysis.integration.gauss.GaussIntegrator in project gatk-protected by broadinstitute.

the class IntegrationUtils method integrate2d.

public static double integrate2d(final ToDoubleBiFunction<Double, Double> getIntegrand, final double xLowerBound, final double xUpperBound, final int xNumPoints, final double yLowerBound, final double yUpperBound, final int yNumPoints) {
    final GaussIntegrator xIntegrator = integratorFactory.legendre(xNumPoints, xLowerBound, xUpperBound);
    final GaussIntegrator yIntegrator = integratorFactory.legendre(yNumPoints, yLowerBound, yUpperBound);
    final double[] xIntegrationWeights = new IndexRange(0, xNumPoints).mapToDouble(xIntegrator::getWeight);
    final double[] xAbscissas = new IndexRange(0, xNumPoints).mapToDouble(xIntegrator::getPoint);
    final double[] yIntegrationWeights = new IndexRange(0, yNumPoints).mapToDouble(yIntegrator::getWeight);
    final double[] yAbscissas = new IndexRange(0, yNumPoints).mapToDouble(yIntegrator::getPoint);
    double integral = 0;
    for (int i = 0; i < xNumPoints; i++) {
        final double x = xAbscissas[i];
        for (int j = 0; j < yNumPoints; j++) {
            final double y = yAbscissas[j];
            final double integrand = getIntegrand.applyAsDouble(x, y);
            integral += xIntegrationWeights[i] * yIntegrationWeights[j] * integrand;
        }
    }
    return integral;
}