http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java deleted file mode 100644 index 5c7b37f..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegrator.java +++ /dev/null @@ -1,129 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.util.Pair; - -/** - * Class that implements the Gaussian rule for - * {@link #integrate(UnivariateFunction) integrating} a weighted - * function. - * - * @since 3.1 - */ -public class GaussIntegrator { - /** Nodes. */ - private final double[] points; - /** Nodes weights. */ - private final double[] weights; - - /** - * Creates an integrator from the given {@code points} and {@code weights}. - * The integration interval is defined by the first and last value of - * {@code points} which must be sorted in increasing order. - * - * @param points Integration points. - * @param weights Weights of the corresponding integration nodes. - * @throws NonMonotonicSequenceException if the {@code points} are not - * sorted in increasing order. - * @throws DimensionMismatchException if points and weights don't have the same length - */ - public GaussIntegrator(double[] points, - double[] weights) - throws NonMonotonicSequenceException, DimensionMismatchException { - if (points.length != weights.length) { - throw new DimensionMismatchException(points.length, - weights.length); - } - - MathArrays.checkOrder(points, MathArrays.OrderDirection.INCREASING, true, true); - - this.points = points.clone(); - this.weights = weights.clone(); - } - - /** - * Creates an integrator from the given pair of points (first element of - * the pair) and weights (second element of the pair. - * - * @param pointsAndWeights Integration points and corresponding weights. - * @throws NonMonotonicSequenceException if the {@code points} are not - * sorted in increasing order. - * - * @see #GaussIntegrator(double[], double[]) - */ - public GaussIntegrator(Pair<double[], double[]> pointsAndWeights) - throws NonMonotonicSequenceException { - this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); - } - - /** - * Returns an estimate of the integral of {@code f(x) * w(x)}, - * where {@code w} is a weight function that depends on the actual - * flavor of the Gauss integration scheme. - * The algorithm uses the points and associated weights, as passed - * to the {@link #GaussIntegrator(double[],double[]) constructor}. - * - * @param f Function to integrate. - * @return the integral of the weighted function. - */ - public double integrate(UnivariateFunction f) { - double s = 0; - double c = 0; - for (int i = 0; i < points.length; i++) { - final double x = points[i]; - final double w = weights[i]; - final double y = w * f.value(x) - c; - final double t = s + y; - c = (t - s) - y; - s = t; - } - return s; - } - - /** - * @return the order of the integration rule (the number of integration - * points). - */ - public int getNumberOfPoints() { - return points.length; - } - - /** - * Gets the integration point at the given index. - * The index must be in the valid range but no check is performed. - * @param index index of the integration point - * @return the integration point. - */ - public double getPoint(int index) { - return points[index]; - } - - /** - * Gets the weight of the integration point at the given index. - * The index must be in the valid range but no check is performed. - * @param index index of the integration point - * @return the weight. - */ - public double getWeight(int index) { - return weights[index]; - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java deleted file mode 100644 index ebe9a5b..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/GaussIntegratorFactory.java +++ /dev/null @@ -1,169 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import java.math.BigDecimal; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.util.Pair; - -/** - * Class that provides different ways to compute the nodes and weights to be - * used by the {@link GaussIntegrator Gaussian integration rule}. - * - * @since 3.1 - */ -public class GaussIntegratorFactory { - /** Generator of Gauss-Legendre integrators. */ - private final BaseRuleFactory<Double> legendre = new LegendreRuleFactory(); - /** Generator of Gauss-Legendre integrators. */ - private final BaseRuleFactory<BigDecimal> legendreHighPrecision = new LegendreHighPrecisionRuleFactory(); - /** Generator of Gauss-Hermite integrators. */ - private final BaseRuleFactory<Double> hermite = new HermiteRuleFactory(); - - /** - * Creates a Gauss-Legendre integrator of the given order. - * The call to the - * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method will perform an integration on the natural interval - * {@code [-1 , 1]}. - * - * @param numberOfPoints Order of the integration rule. - * @return a Gauss-Legendre integrator. - */ - public GaussIntegrator legendre(int numberOfPoints) { - return new GaussIntegrator(getRule(legendre, numberOfPoints)); - } - - /** - * Creates a Gauss-Legendre integrator of the given order. - * The call to the - * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method will perform an integration on the given interval. - * - * @param numberOfPoints Order of the integration rule. - * @param lowerBound Lower bound of the integration interval. - * @param upperBound Upper bound of the integration interval. - * @return a Gauss-Legendre integrator. - * @throws NotStrictlyPositiveException if number of points is not positive - */ - public GaussIntegrator legendre(int numberOfPoints, - double lowerBound, - double upperBound) - throws NotStrictlyPositiveException { - return new GaussIntegrator(transform(getRule(legendre, numberOfPoints), - lowerBound, upperBound)); - } - - /** - * Creates a Gauss-Legendre integrator of the given order. - * The call to the - * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method will perform an integration on the natural interval - * {@code [-1 , 1]}. - * - * @param numberOfPoints Order of the integration rule. - * @return a Gauss-Legendre integrator. - * @throws NotStrictlyPositiveException if number of points is not positive - */ - public GaussIntegrator legendreHighPrecision(int numberOfPoints) - throws NotStrictlyPositiveException { - return new GaussIntegrator(getRule(legendreHighPrecision, numberOfPoints)); - } - - /** - * Creates an integrator of the given order, and whose call to the - * {@link GaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method will perform an integration on the given interval. - * - * @param numberOfPoints Order of the integration rule. - * @param lowerBound Lower bound of the integration interval. - * @param upperBound Upper bound of the integration interval. - * @return a Gauss-Legendre integrator. - * @throws NotStrictlyPositiveException if number of points is not positive - */ - public GaussIntegrator legendreHighPrecision(int numberOfPoints, - double lowerBound, - double upperBound) - throws NotStrictlyPositiveException { - return new GaussIntegrator(transform(getRule(legendreHighPrecision, numberOfPoints), - lowerBound, upperBound)); - } - - /** - * Creates a Gauss-Hermite integrator of the given order. - * The call to the - * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method will perform a weighted integration on the interval - * {@code [-&inf;, +&inf;]}: the computed value is the improper integral of - * <code> - * e<sup>-x<sup>2</sup></sup> f(x) - * </code> - * where {@code f(x)} is the function passed to the - * {@link SymmetricGaussIntegrator#integrate(org.apache.commons.math3.analysis.UnivariateFunction) - * integrate} method. - * - * @param numberOfPoints Order of the integration rule. - * @return a Gauss-Hermite integrator. - */ - public SymmetricGaussIntegrator hermite(int numberOfPoints) { - return new SymmetricGaussIntegrator(getRule(hermite, numberOfPoints)); - } - - /** - * @param factory Integration rule factory. - * @param numberOfPoints Order of the integration rule. - * @return the integration nodes and weights. - * @throws NotStrictlyPositiveException if number of points is not positive - * @throws DimensionMismatchException if the elements of the rule pair do not - * have the same length. - */ - private static Pair<double[], double[]> getRule(BaseRuleFactory<? extends Number> factory, - int numberOfPoints) - throws NotStrictlyPositiveException, DimensionMismatchException { - return factory.getRule(numberOfPoints); - } - - /** - * Performs a change of variable so that the integration can be performed - * on an arbitrary interval {@code [a, b]}. - * It is assumed that the natural interval is {@code [-1, 1]}. - * - * @param rule Original points and weights. - * @param a Lower bound of the integration interval. - * @param b Lower bound of the integration interval. - * @return the points and weights adapted to the new interval. - */ - private static Pair<double[], double[]> transform(Pair<double[], double[]> rule, - double a, - double b) { - final double[] points = rule.getFirst(); - final double[] weights = rule.getSecond(); - - // Scaling - final double scale = (b - a) / 2; - final double shift = a + scale; - - for (int i = 0; i < points.length; i++) { - points[i] = points[i] * scale + shift; - weights[i] *= scale; - } - - return new Pair<double[], double[]>(points, weights); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java deleted file mode 100644 index 3d873ab..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/HermiteRuleFactory.java +++ /dev/null @@ -1,178 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.util.Pair; -import org.apache.commons.math3.util.FastMath; - -/** - * Factory that creates a - * <a href="http://en.wikipedia.org/wiki/Gauss-Hermite_quadrature";> - * Gauss-type quadrature rule using Hermite polynomials</a> - * of the first kind. - * Such a quadrature rule allows the calculation of improper integrals - * of a function - * <code> - * f(x) e<sup>-x<sup>2</sup></sup> - * </code> - * <br/> - * Recurrence relation and weights computation follow - * <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun";> - * Abramowitz and Stegun, 1964</a>. - * <br/> - * The coefficients of the standard Hermite polynomials grow very rapidly; - * in order to avoid overflows, each Hermite polynomial is normalized with - * respect to the underlying scalar product. - * The initial interval for the application of the bisection method is - * based on the roots of the previous Hermite polynomial (interlacing). - * Upper and lower bounds of these roots are provided by - * <blockquote> - * I. Krasikov,<br> - * <em>Nonnegative quadratic forms and bounds on orthogonal polynomials</em>,<br> - * Journal of Approximation theory <b>111</b>, 31-49<br> - * </blockquote> - * - * @since 3.3 - */ -public class HermiteRuleFactory extends BaseRuleFactory<Double> { - /** π<sup>1/2</sup> */ - private static final double SQRT_PI = 1.77245385090551602729; - /** π<sup>-1/4</sup> */ - private static final double H0 = 7.5112554446494248286e-1; - /** π<sup>-1/4</sup> √2 */ - private static final double H1 = 1.0622519320271969145; - - /** {@inheritDoc} */ - @Override - protected Pair<Double[], Double[]> computeRule(int numberOfPoints) - throws DimensionMismatchException { - - if (numberOfPoints == 1) { - // Break recursion. - return new Pair<Double[], Double[]>(new Double[] { 0d }, - new Double[] { SQRT_PI }); - } - - // Get previous rule. - // If it has not been computed yet it will trigger a recursive call - // to this method. - final int lastNumPoints = numberOfPoints - 1; - final Double[] previousPoints = getRuleInternal(lastNumPoints).getFirst(); - - // Compute next rule. - final Double[] points = new Double[numberOfPoints]; - final Double[] weights = new Double[numberOfPoints]; - - final double sqrtTwoTimesLastNumPoints = FastMath.sqrt(2 * lastNumPoints); - final double sqrtTwoTimesNumPoints = FastMath.sqrt(2 * numberOfPoints); - - // Find i-th root of H[n+1] by bracketing. - final int iMax = numberOfPoints / 2; - for (int i = 0; i < iMax; i++) { - // Lower-bound of the interval. - double a = (i == 0) ? -sqrtTwoTimesLastNumPoints : previousPoints[i - 1].doubleValue(); - // Upper-bound of the interval. - double b = (iMax == 1) ? -0.5 : previousPoints[i].doubleValue(); - - // H[j-1](a) - double hma = H0; - // H[j](a) - double ha = H1 * a; - // H[j-1](b) - double hmb = H0; - // H[j](b) - double hb = H1 * b; - for (int j = 1; j < numberOfPoints; j++) { - // Compute H[j+1](a) and H[j+1](b) - final double jp1 = j + 1; - final double s = FastMath.sqrt(2 / jp1); - final double sm = FastMath.sqrt(j / jp1); - final double hpa = s * a * ha - sm * hma; - final double hpb = s * b * hb - sm * hmb; - hma = ha; - ha = hpa; - hmb = hb; - hb = hpb; - } - - // Now ha = H[n+1](a), and hma = H[n](a) (same holds for b). - // Middle of the interval. - double c = 0.5 * (a + b); - // P[j-1](c) - double hmc = H0; - // P[j](c) - double hc = H1 * c; - boolean done = false; - while (!done) { - done = b - a <= Math.ulp(c); - hmc = H0; - hc = H1 * c; - for (int j = 1; j < numberOfPoints; j++) { - // Compute H[j+1](c) - final double jp1 = j + 1; - final double s = FastMath.sqrt(2 / jp1); - final double sm = FastMath.sqrt(j / jp1); - final double hpc = s * c * hc - sm * hmc; - hmc = hc; - hc = hpc; - } - // Now h = H[n+1](c) and hm = H[n](c). - if (!done) { - if (ha * hc < 0) { - b = c; - hmb = hmc; - hb = hc; - } else { - a = c; - hma = hmc; - ha = hc; - } - c = 0.5 * (a + b); - } - } - final double d = sqrtTwoTimesNumPoints * hmc; - final double w = 2 / (d * d); - - points[i] = c; - weights[i] = w; - - final int idx = lastNumPoints - i; - points[idx] = -c; - weights[idx] = w; - } - - // If "numberOfPoints" is odd, 0 is a root. - // Note: as written, the test for oddness will work for negative - // integers too (although it is not necessary here), preventing - // a FindBugs warning. - if (numberOfPoints % 2 != 0) { - double hm = H0; - for (int j = 1; j < numberOfPoints; j += 2) { - final double jp1 = j + 1; - hm = -FastMath.sqrt(j / jp1) * hm; - } - final double d = sqrtTwoTimesNumPoints * hm; - final double w = 2 / (d * d); - - points[iMax] = 0d; - weights[iMax] = w; - } - - return new Pair<Double[], Double[]>(points, weights); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java deleted file mode 100644 index 93e1738..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreHighPrecisionRuleFactory.java +++ /dev/null @@ -1,215 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import java.math.BigDecimal; -import java.math.MathContext; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.util.Pair; - -/** - * Factory that creates Gauss-type quadrature rule using Legendre polynomials. - * In this implementation, the lower and upper bounds of the natural interval - * of integration are -1 and 1, respectively. - * The Legendre polynomials are evaluated using the recurrence relation - * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"; - * Abramowitz and Stegun, 1964</a>. - * - * @since 3.1 - */ -public class LegendreHighPrecisionRuleFactory extends BaseRuleFactory<BigDecimal> { - /** Settings for enhanced precision computations. */ - private final MathContext mContext; - /** The number {@code 2}. */ - private final BigDecimal two; - /** The number {@code -1}. */ - private final BigDecimal minusOne; - /** The number {@code 0.5}. */ - private final BigDecimal oneHalf; - - /** - * Default precision is {@link MathContext#DECIMAL128 DECIMAL128}. - */ - public LegendreHighPrecisionRuleFactory() { - this(MathContext.DECIMAL128); - } - - /** - * @param mContext Precision setting for computing the quadrature rules. - */ - public LegendreHighPrecisionRuleFactory(MathContext mContext) { - this.mContext = mContext; - two = new BigDecimal("2", mContext); - minusOne = new BigDecimal("-1", mContext); - oneHalf = new BigDecimal("0.5", mContext); - } - - /** {@inheritDoc} */ - @Override - protected Pair<BigDecimal[], BigDecimal[]> computeRule(int numberOfPoints) - throws DimensionMismatchException { - - if (numberOfPoints == 1) { - // Break recursion. - return new Pair<BigDecimal[], BigDecimal[]>(new BigDecimal[] { BigDecimal.ZERO }, - new BigDecimal[] { two }); - } - - // Get previous rule. - // If it has not been computed yet it will trigger a recursive call - // to this method. - final BigDecimal[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); - - // Compute next rule. - final BigDecimal[] points = new BigDecimal[numberOfPoints]; - final BigDecimal[] weights = new BigDecimal[numberOfPoints]; - - // Find i-th root of P[n+1] by bracketing. - final int iMax = numberOfPoints / 2; - for (int i = 0; i < iMax; i++) { - // Lower-bound of the interval. - BigDecimal a = (i == 0) ? minusOne : previousPoints[i - 1]; - // Upper-bound of the interval. - BigDecimal b = (iMax == 1) ? BigDecimal.ONE : previousPoints[i]; - // P[j-1](a) - BigDecimal pma = BigDecimal.ONE; - // P[j](a) - BigDecimal pa = a; - // P[j-1](b) - BigDecimal pmb = BigDecimal.ONE; - // P[j](b) - BigDecimal pb = b; - for (int j = 1; j < numberOfPoints; j++) { - final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); - final BigDecimal b_j = new BigDecimal(j, mContext); - final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); - - // Compute P[j+1](a) - // ppa = ((2 * j + 1) * a * pa - j * pma) / (j + 1); - - BigDecimal tmp1 = a.multiply(b_two_j_p_1, mContext); - tmp1 = pa.multiply(tmp1, mContext); - BigDecimal tmp2 = pma.multiply(b_j, mContext); - // P[j+1](a) - BigDecimal ppa = tmp1.subtract(tmp2, mContext); - ppa = ppa.divide(b_j_p_1, mContext); - - // Compute P[j+1](b) - // ppb = ((2 * j + 1) * b * pb - j * pmb) / (j + 1); - - tmp1 = b.multiply(b_two_j_p_1, mContext); - tmp1 = pb.multiply(tmp1, mContext); - tmp2 = pmb.multiply(b_j, mContext); - // P[j+1](b) - BigDecimal ppb = tmp1.subtract(tmp2, mContext); - ppb = ppb.divide(b_j_p_1, mContext); - - pma = pa; - pa = ppa; - pmb = pb; - pb = ppb; - } - // Now pa = P[n+1](a), and pma = P[n](a). Same holds for b. - // Middle of the interval. - BigDecimal c = a.add(b, mContext).multiply(oneHalf, mContext); - // P[j-1](c) - BigDecimal pmc = BigDecimal.ONE; - // P[j](c) - BigDecimal pc = c; - boolean done = false; - while (!done) { - BigDecimal tmp1 = b.subtract(a, mContext); - BigDecimal tmp2 = c.ulp().multiply(BigDecimal.TEN, mContext); - done = tmp1.compareTo(tmp2) <= 0; - pmc = BigDecimal.ONE; - pc = c; - for (int j = 1; j < numberOfPoints; j++) { - final BigDecimal b_two_j_p_1 = new BigDecimal(2 * j + 1, mContext); - final BigDecimal b_j = new BigDecimal(j, mContext); - final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); - - // Compute P[j+1](c) - tmp1 = c.multiply(b_two_j_p_1, mContext); - tmp1 = pc.multiply(tmp1, mContext); - tmp2 = pmc.multiply(b_j, mContext); - // P[j+1](c) - BigDecimal ppc = tmp1.subtract(tmp2, mContext); - ppc = ppc.divide(b_j_p_1, mContext); - - pmc = pc; - pc = ppc; - } - // Now pc = P[n+1](c) and pmc = P[n](c). - if (!done) { - if (pa.signum() * pc.signum() <= 0) { - b = c; - pmb = pmc; - pb = pc; - } else { - a = c; - pma = pmc; - pa = pc; - } - c = a.add(b, mContext).multiply(oneHalf, mContext); - } - } - final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); - BigDecimal tmp1 = pmc.subtract(c.multiply(pc, mContext), mContext); - tmp1 = tmp1.multiply(nP); - tmp1 = tmp1.pow(2, mContext); - BigDecimal tmp2 = c.pow(2, mContext); - tmp2 = BigDecimal.ONE.subtract(tmp2, mContext); - tmp2 = tmp2.multiply(two, mContext); - tmp2 = tmp2.divide(tmp1, mContext); - - points[i] = c; - weights[i] = tmp2; - - final int idx = numberOfPoints - i - 1; - points[idx] = c.negate(mContext); - weights[idx] = tmp2; - } - // If "numberOfPoints" is odd, 0 is a root. - // Note: as written, the test for oddness will work for negative - // integers too (although it is not necessary here), preventing - // a FindBugs warning. - if (numberOfPoints % 2 != 0) { - BigDecimal pmc = BigDecimal.ONE; - for (int j = 1; j < numberOfPoints; j += 2) { - final BigDecimal b_j = new BigDecimal(j, mContext); - final BigDecimal b_j_p_1 = new BigDecimal(j + 1, mContext); - - // pmc = -j * pmc / (j + 1); - pmc = pmc.multiply(b_j, mContext); - pmc = pmc.divide(b_j_p_1, mContext); - pmc = pmc.negate(mContext); - } - - // 2 / pow(numberOfPoints * pmc, 2); - final BigDecimal nP = new BigDecimal(numberOfPoints, mContext); - BigDecimal tmp1 = pmc.multiply(nP, mContext); - tmp1 = tmp1.pow(2, mContext); - BigDecimal tmp2 = two.divide(tmp1, mContext); - - points[iMax] = BigDecimal.ZERO; - weights[iMax] = tmp2; - } - - return new Pair<BigDecimal[], BigDecimal[]>(points, weights); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java deleted file mode 100644 index 225fa01..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/LegendreRuleFactory.java +++ /dev/null @@ -1,140 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.util.Pair; - -/** - * Factory that creates Gauss-type quadrature rule using Legendre polynomials. - * In this implementation, the lower and upper bounds of the natural interval - * of integration are -1 and 1, respectively. - * The Legendre polynomials are evaluated using the recurrence relation - * presented in <a href="http://en.wikipedia.org/wiki/Abramowitz_and_Stegun"; - * Abramowitz and Stegun, 1964</a>. - * - * @since 3.1 - */ -public class LegendreRuleFactory extends BaseRuleFactory<Double> { - /** {@inheritDoc} */ - @Override - protected Pair<Double[], Double[]> computeRule(int numberOfPoints) - throws DimensionMismatchException { - - if (numberOfPoints == 1) { - // Break recursion. - return new Pair<Double[], Double[]>(new Double[] { 0d }, - new Double[] { 2d }); - } - - // Get previous rule. - // If it has not been computed yet it will trigger a recursive call - // to this method. - final Double[] previousPoints = getRuleInternal(numberOfPoints - 1).getFirst(); - - // Compute next rule. - final Double[] points = new Double[numberOfPoints]; - final Double[] weights = new Double[numberOfPoints]; - - // Find i-th root of P[n+1] by bracketing. - final int iMax = numberOfPoints / 2; - for (int i = 0; i < iMax; i++) { - // Lower-bound of the interval. - double a = (i == 0) ? -1 : previousPoints[i - 1].doubleValue(); - // Upper-bound of the interval. - double b = (iMax == 1) ? 1 : previousPoints[i].doubleValue(); - // P[j-1](a) - double pma = 1; - // P[j](a) - double pa = a; - // P[j-1](b) - double pmb = 1; - // P[j](b) - double pb = b; - for (int j = 1; j < numberOfPoints; j++) { - final int two_j_p_1 = 2 * j + 1; - final int j_p_1 = j + 1; - // P[j+1](a) - final double ppa = (two_j_p_1 * a * pa - j * pma) / j_p_1; - // P[j+1](b) - final double ppb = (two_j_p_1 * b * pb - j * pmb) / j_p_1; - pma = pa; - pa = ppa; - pmb = pb; - pb = ppb; - } - // Now pa = P[n+1](a), and pma = P[n](a) (same holds for b). - // Middle of the interval. - double c = 0.5 * (a + b); - // P[j-1](c) - double pmc = 1; - // P[j](c) - double pc = c; - boolean done = false; - while (!done) { - done = b - a <= Math.ulp(c); - pmc = 1; - pc = c; - for (int j = 1; j < numberOfPoints; j++) { - // P[j+1](c) - final double ppc = ((2 * j + 1) * c * pc - j * pmc) / (j + 1); - pmc = pc; - pc = ppc; - } - // Now pc = P[n+1](c) and pmc = P[n](c). - if (!done) { - if (pa * pc <= 0) { - b = c; - pmb = pmc; - pb = pc; - } else { - a = c; - pma = pmc; - pa = pc; - } - c = 0.5 * (a + b); - } - } - final double d = numberOfPoints * (pmc - c * pc); - final double w = 2 * (1 - c * c) / (d * d); - - points[i] = c; - weights[i] = w; - - final int idx = numberOfPoints - i - 1; - points[idx] = -c; - weights[idx] = w; - } - // If "numberOfPoints" is odd, 0 is a root. - // Note: as written, the test for oddness will work for negative - // integers too (although it is not necessary here), preventing - // a FindBugs warning. - if (numberOfPoints % 2 != 0) { - double pmc = 1; - for (int j = 1; j < numberOfPoints; j += 2) { - pmc = -j * pmc / (j + 1); - } - final double d = numberOfPoints * pmc; - final double w = 2 / (d * d); - - points[iMax] = 0d; - weights[iMax] = w; - } - - return new Pair<Double[], Double[]>(points, weights); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java deleted file mode 100644 index 7fa4884..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/SymmetricGaussIntegrator.java +++ /dev/null @@ -1,103 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.integration.gauss; - -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.util.Pair; - -/** - * This class's implements {@link #integrate(UnivariateFunction) integrate} - * method assuming that the integral is symmetric about 0. - * This allows to reduce numerical errors. - * - * @since 3.3 - */ -public class SymmetricGaussIntegrator extends GaussIntegrator { - /** - * Creates an integrator from the given {@code points} and {@code weights}. - * The integration interval is defined by the first and last value of - * {@code points} which must be sorted in increasing order. - * - * @param points Integration points. - * @param weights Weights of the corresponding integration nodes. - * @throws NonMonotonicSequenceException if the {@code points} are not - * sorted in increasing order. - * @throws DimensionMismatchException if points and weights don't have the same length - */ - public SymmetricGaussIntegrator(double[] points, - double[] weights) - throws NonMonotonicSequenceException, DimensionMismatchException { - super(points, weights); - } - - /** - * Creates an integrator from the given pair of points (first element of - * the pair) and weights (second element of the pair. - * - * @param pointsAndWeights Integration points and corresponding weights. - * @throws NonMonotonicSequenceException if the {@code points} are not - * sorted in increasing order. - * - * @see #SymmetricGaussIntegrator(double[], double[]) - */ - public SymmetricGaussIntegrator(Pair<double[], double[]> pointsAndWeights) - throws NonMonotonicSequenceException { - this(pointsAndWeights.getFirst(), pointsAndWeights.getSecond()); - } - - /** - * {@inheritDoc} - */ - @Override - public double integrate(UnivariateFunction f) { - final int ruleLength = getNumberOfPoints(); - - if (ruleLength == 1) { - return getWeight(0) * f.value(0d); - } - - final int iMax = ruleLength / 2; - double s = 0; - double c = 0; - for (int i = 0; i < iMax; i++) { - final double p = getPoint(i); - final double w = getWeight(i); - - final double f1 = f.value(p); - final double f2 = f.value(-p); - - final double y = w * (f1 + f2) - c; - final double t = s + y; - - c = (t - s) - y; - s = t; - } - - if (ruleLength % 2 != 0) { - final double w = getWeight(iMax); - - final double y = w * f.value(0d) - c; - final double t = s + y; - - s = t; - } - - return s; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java b/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java deleted file mode 100644 index 066da30..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/gauss/package-info.java +++ /dev/null @@ -1,22 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -/** - * - * Gauss family of quadrature schemes. - * - */ -package org.apache.commons.math3.analysis.integration.gauss; http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java b/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java deleted file mode 100644 index f4df3ba..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/integration/package-info.java +++ /dev/null @@ -1,22 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -/** - * - * Numerical integration (quadrature) algorithms for univariate real functions. - * - */ -package org.apache.commons.math3.analysis.integration; http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java deleted file mode 100644 index 5b89dfe..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/AkimaSplineInterpolator.java +++ /dev/null @@ -1,215 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.util.Precision; - -/** - * Computes a cubic spline interpolation for the data set using the Akima - * algorithm, as originally formulated by Hiroshi Akima in his 1970 paper - * "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures." - * J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609 - * http://doi.acm.org/10.1145/321607.321609 - * <p> - * This implementation is based on the Akima implementation in the CubicSpline - * class in the Math.NET Numerics library. The method referenced is - * CubicSpline.InterpolateAkimaSorted - * </p> - * <p> - * The {@link #interpolate(double[], double[]) interpolate} method returns a - * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined - * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}. - * The Akima algorithm requires that {@code n >= 5}. - * </p> - */ -public class AkimaSplineInterpolator - implements UnivariateInterpolator { - /** The minimum number of points that are needed to compute the function. */ - private static final int MINIMUM_NUMBER_POINTS = 5; - - /** - * Computes an interpolating function for the data set. - * - * @param xvals the arguments for the interpolation points - * @param yvals the values for the interpolation points - * @return a function which interpolates the data set - * @throws DimensionMismatchException if {@code xvals} and {@code yvals} have - * different sizes. - * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in - * strict increasing order. - * @throws NumberIsTooSmallException if the size of {@code xvals} is smaller - * than 5. - */ - public PolynomialSplineFunction interpolate(double[] xvals, - double[] yvals) - throws DimensionMismatchException, - NumberIsTooSmallException, - NonMonotonicSequenceException { - if (xvals == null || - yvals == null) { - throw new NullArgumentException(); - } - - if (xvals.length != yvals.length) { - throw new DimensionMismatchException(xvals.length, yvals.length); - } - - if (xvals.length < MINIMUM_NUMBER_POINTS) { - throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, - xvals.length, - MINIMUM_NUMBER_POINTS, true); - } - - MathArrays.checkOrder(xvals); - - final int numberOfDiffAndWeightElements = xvals.length - 1; - - final double[] differences = new double[numberOfDiffAndWeightElements]; - final double[] weights = new double[numberOfDiffAndWeightElements]; - - for (int i = 0; i < differences.length; i++) { - differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]); - } - - for (int i = 1; i < weights.length; i++) { - weights[i] = FastMath.abs(differences[i] - differences[i - 1]); - } - - // Prepare Hermite interpolation scheme. - final double[] firstDerivatives = new double[xvals.length]; - - for (int i = 2; i < firstDerivatives.length - 2; i++) { - final double wP = weights[i + 1]; - final double wM = weights[i - 1]; - if (Precision.equals(wP, 0.0) && - Precision.equals(wM, 0.0)) { - final double xv = xvals[i]; - final double xvP = xvals[i + 1]; - final double xvM = xvals[i - 1]; - firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM); - } else { - firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM); - } - } - - firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2); - firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2); - firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2, - xvals.length - 3, xvals.length - 2, - xvals.length - 1); - firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1, - xvals.length - 3, xvals.length - 2, - xvals.length - 1); - - return interpolateHermiteSorted(xvals, yvals, firstDerivatives); - } - - /** - * Three point differentiation helper, modeled off of the same method in the - * Math.NET CubicSpline class. This is used by both the Apache Math and the - * Math.NET Akima Cubic Spline algorithms - * - * @param xvals x values to calculate the numerical derivative with - * @param yvals y values to calculate the numerical derivative with - * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around - * @param indexOfFirstSample index of the first element to sample for the three point method - * @param indexOfSecondsample index of the second element to sample for the three point method - * @param indexOfThirdSample index of the third element to sample for the three point method - * @return the derivative - */ - private double differentiateThreePoint(double[] xvals, double[] yvals, - int indexOfDifferentiation, - int indexOfFirstSample, - int indexOfSecondsample, - int indexOfThirdSample) { - final double x0 = yvals[indexOfFirstSample]; - final double x1 = yvals[indexOfSecondsample]; - final double x2 = yvals[indexOfThirdSample]; - - final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample]; - final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample]; - final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample]; - - final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2); - final double b = (x1 - x0 - a * t1 * t1) / t1; - - return (2 * a * t) + b; - } - - /** - * Creates a Hermite cubic spline interpolation from the set of (x,y) value - * pairs and their derivatives. This is modeled off of the - * InterpolateHermiteSorted method in the Math.NET CubicSpline class. - * - * @param xvals x values for interpolation - * @param yvals y values for interpolation - * @param firstDerivatives first derivative values of the function - * @return polynomial that fits the function - */ - private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals, - double[] yvals, - double[] firstDerivatives) { - if (xvals.length != yvals.length) { - throw new DimensionMismatchException(xvals.length, yvals.length); - } - - if (xvals.length != firstDerivatives.length) { - throw new DimensionMismatchException(xvals.length, - firstDerivatives.length); - } - - final int minimumLength = 2; - if (xvals.length < minimumLength) { - throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, - xvals.length, minimumLength, - true); - } - - final int size = xvals.length - 1; - final PolynomialFunction[] polynomials = new PolynomialFunction[size]; - final double[] coefficients = new double[4]; - - for (int i = 0; i < polynomials.length; i++) { - final double w = xvals[i + 1] - xvals[i]; - final double w2 = w * w; - - final double yv = yvals[i]; - final double yvP = yvals[i + 1]; - - final double fd = firstDerivatives[i]; - final double fdP = firstDerivatives[i + 1]; - - coefficients[0] = yv; - coefficients[1] = firstDerivatives[i]; - coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w; - coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2; - polynomials[i] = new PolynomialFunction(coefficients); - } - - return new PolynomialSplineFunction(xvals, polynomials); - - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java deleted file mode 100644 index 7fe947f..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolatingFunction.java +++ /dev/null @@ -1,325 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import java.util.Arrays; -import org.apache.commons.math3.analysis.BivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Function that implements the - * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation";> - * bicubic spline interpolation</a>. - * - * @since 3.4 - */ -public class BicubicInterpolatingFunction - implements BivariateFunction { - /** Number of coefficients. */ - private static final int NUM_COEFF = 16; - /** - * Matrix to compute the spline coefficients from the function values - * and function derivatives values - */ - private static final double[][] AINV = { - { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, - { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, - { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, - { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, - { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, - { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, - { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, - { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, - { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, - { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, - { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, - { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, - { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } - }; - - /** Samples x-coordinates */ - private final double[] xval; - /** Samples y-coordinates */ - private final double[] yval; - /** Set of cubic splines patching the whole data grid */ - private final BicubicFunction[][] splines; - - /** - * @param x Sample values of the x-coordinate, in increasing order. - * @param y Sample values of the y-coordinate, in increasing order. - * @param f Values of the function on every grid point. - * @param dFdX Values of the partial derivative of function with respect - * to x on every grid point. - * @param dFdY Values of the partial derivative of function with respect - * to y on every grid point. - * @param d2FdXdY Values of the cross partial derivative of function on - * every grid point. - * @throws DimensionMismatchException if the various arrays do not contain - * the expected number of elements. - * @throws NonMonotonicSequenceException if {@code x} or {@code y} are - * not strictly increasing. - * @throws NoDataException if any of the arrays has zero length. - */ - public BicubicInterpolatingFunction(double[] x, - double[] y, - double[][] f, - double[][] dFdX, - double[][] dFdY, - double[][] d2FdXdY) - throws DimensionMismatchException, - NoDataException, - NonMonotonicSequenceException { - final int xLen = x.length; - final int yLen = y.length; - - if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { - throw new NoDataException(); - } - if (xLen != f.length) { - throw new DimensionMismatchException(xLen, f.length); - } - if (xLen != dFdX.length) { - throw new DimensionMismatchException(xLen, dFdX.length); - } - if (xLen != dFdY.length) { - throw new DimensionMismatchException(xLen, dFdY.length); - } - if (xLen != d2FdXdY.length) { - throw new DimensionMismatchException(xLen, d2FdXdY.length); - } - - MathArrays.checkOrder(x); - MathArrays.checkOrder(y); - - xval = x.clone(); - yval = y.clone(); - - final int lastI = xLen - 1; - final int lastJ = yLen - 1; - splines = new BicubicFunction[lastI][lastJ]; - - for (int i = 0; i < lastI; i++) { - if (f[i].length != yLen) { - throw new DimensionMismatchException(f[i].length, yLen); - } - if (dFdX[i].length != yLen) { - throw new DimensionMismatchException(dFdX[i].length, yLen); - } - if (dFdY[i].length != yLen) { - throw new DimensionMismatchException(dFdY[i].length, yLen); - } - if (d2FdXdY[i].length != yLen) { - throw new DimensionMismatchException(d2FdXdY[i].length, yLen); - } - final int ip1 = i + 1; - final double xR = xval[ip1] - xval[i]; - for (int j = 0; j < lastJ; j++) { - final int jp1 = j + 1; - final double yR = yval[jp1] - yval[j]; - final double xRyR = xR * yR; - final double[] beta = new double[] { - f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], - dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR, - dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR, - d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR - }; - - splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta)); - } - } - } - - /** - * {@inheritDoc} - */ - public double value(double x, double y) - throws OutOfRangeException { - final int i = searchIndex(x, xval); - final int j = searchIndex(y, yval); - - final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); - final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); - - return splines[i][j].value(xN, yN); - } - - /** - * Indicates whether a point is within the interpolation range. - * - * @param x First coordinate. - * @param y Second coordinate. - * @return {@code true} if (x, y) is a valid point. - */ - public boolean isValidPoint(double x, double y) { - if (x < xval[0] || - x > xval[xval.length - 1] || - y < yval[0] || - y > yval[yval.length - 1]) { - return false; - } else { - return true; - } - } - - /** - * @param c Coordinate. - * @param val Coordinate samples. - * @return the index in {@code val} corresponding to the interval - * containing {@code c}. - * @throws OutOfRangeException if {@code c} is out of the - * range defined by the boundary values of {@code val}. - */ - private int searchIndex(double c, double[] val) { - final int r = Arrays.binarySearch(val, c); - - if (r == -1 || - r == -val.length - 1) { - throw new OutOfRangeException(c, val[0], val[val.length - 1]); - } - - if (r < 0) { - // "c" in within an interpolation sub-interval: Return the - // index of the sample at the lower end of the sub-interval. - return -r - 2; - } - final int last = val.length - 1; - if (r == last) { - // "c" is the last sample of the range: Return the index - // of the sample at the lower end of the last sub-interval. - return last - 1; - } - - // "c" is another sample point. - return r; - } - - /** - * Compute the spline coefficients from the list of function values and - * function partial derivatives values at the four corners of a grid - * element. They must be specified in the following order: - * <ul> - * <li>f(0,0)</li> - * <li>f(1,0)</li> - * <li>f(0,1)</li> - * <li>f(1,1)</li> - * <li>f<sub>x</sub>(0,0)</li> - * <li>f<sub>x</sub>(1,0)</li> - * <li>f<sub>x</sub>(0,1)</li> - * <li>f<sub>x</sub>(1,1)</li> - * <li>f<sub>y</sub>(0,0)</li> - * <li>f<sub>y</sub>(1,0)</li> - * <li>f<sub>y</sub>(0,1)</li> - * <li>f<sub>y</sub>(1,1)</li> - * <li>f<sub>xy</sub>(0,0)</li> - * <li>f<sub>xy</sub>(1,0)</li> - * <li>f<sub>xy</sub>(0,1)</li> - * <li>f<sub>xy</sub>(1,1)</li> - * </ul> - * where the subscripts indicate the partial derivative with respect to - * the corresponding variable(s). - * - * @param beta List of function values and function partial derivatives - * values. - * @return the spline coefficients. - */ - private double[] computeSplineCoefficients(double[] beta) { - final double[] a = new double[NUM_COEFF]; - - for (int i = 0; i < NUM_COEFF; i++) { - double result = 0; - final double[] row = AINV[i]; - for (int j = 0; j < NUM_COEFF; j++) { - result += row[j] * beta[j]; - } - a[i] = result; - } - - return a; - } -} - -/** - * Bicubic function. - */ -class BicubicFunction implements BivariateFunction { - /** Number of points. */ - private static final short N = 4; - /** Coefficients */ - private final double[][] a; - - /** - * Simple constructor. - * - * @param coeff Spline coefficients. - */ - public BicubicFunction(double[] coeff) { - a = new double[N][N]; - for (int j = 0; j < N; j++) { - final double[] aJ = a[j]; - for (int i = 0; i < N; i++) { - aJ[i] = coeff[i * N + j]; - } - } - } - - /** - * {@inheritDoc} - */ - public double value(double x, double y) { - if (x < 0 || x > 1) { - throw new OutOfRangeException(x, 0, 1); - } - if (y < 0 || y > 1) { - throw new OutOfRangeException(y, 0, 1); - } - - final double x2 = x * x; - final double x3 = x2 * x; - final double[] pX = {1, x, x2, x3}; - - final double y2 = y * y; - final double y3 = y2 * y; - final double[] pY = {1, y, y2, y3}; - - return apply(pX, pY, a); - } - - /** - * Compute the value of the bicubic polynomial. - * - * @param pX Powers of the x-coordinate. - * @param pY Powers of the y-coordinate. - * @param coeff Spline coefficients. - * @return the interpolated value. - */ - private double apply(double[] pX, double[] pY, double[][] coeff) { - double result = 0; - for (int i = 0; i < N; i++) { - final double r = MathArrays.linearCombination(coeff[i], pY); - result += r * pX[i]; - } - - return result; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java deleted file mode 100644 index 82ab44e..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/BicubicInterpolator.java +++ /dev/null @@ -1,112 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Generates a {@link BicubicInterpolatingFunction bicubic interpolating - * function}. - * <p> - * Caveat: Because the interpolation scheme requires that derivatives be - * specified at the sample points, those are approximated with finite - * differences (using the 2-points symmetric formulae). - * Since their values are undefined at the borders of the provided - * interpolation ranges, the interpolated values will be wrong at the - * edges of the patch. - * The {@code interpolate} method will return a function that overrides - * {@link BicubicInterpolatingFunction#isValidPoint(double,double)} to - * indicate points where the interpolation will be inaccurate. - * </p> - * - * @since 3.4 - */ -public class BicubicInterpolator - implements BivariateGridInterpolator { - /** - * {@inheritDoc} - */ - public BicubicInterpolatingFunction interpolate(final double[] xval, - final double[] yval, - final double[][] fval) - throws NoDataException, DimensionMismatchException, - NonMonotonicSequenceException, NumberIsTooSmallException { - if (xval.length == 0 || yval.length == 0 || fval.length == 0) { - throw new NoDataException(); - } - if (xval.length != fval.length) { - throw new DimensionMismatchException(xval.length, fval.length); - } - - MathArrays.checkOrder(xval); - MathArrays.checkOrder(yval); - - final int xLen = xval.length; - final int yLen = yval.length; - - // Approximation to the partial derivatives using finite differences. - final double[][] dFdX = new double[xLen][yLen]; - final double[][] dFdY = new double[xLen][yLen]; - final double[][] d2FdXdY = new double[xLen][yLen]; - for (int i = 1; i < xLen - 1; i++) { - final int nI = i + 1; - final int pI = i - 1; - - final double nX = xval[nI]; - final double pX = xval[pI]; - - final double deltaX = nX - pX; - - for (int j = 1; j < yLen - 1; j++) { - final int nJ = j + 1; - final int pJ = j - 1; - - final double nY = yval[nJ]; - final double pY = yval[pJ]; - - final double deltaY = nY - pY; - - dFdX[i][j] = (fval[nI][j] - fval[pI][j]) / deltaX; - dFdY[i][j] = (fval[i][nJ] - fval[i][pJ]) / deltaY; - - final double deltaXY = deltaX * deltaY; - - d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / deltaXY; - } - } - - // Create the interpolating function. - return new BicubicInterpolatingFunction(xval, yval, fval, - dFdX, dFdY, d2FdXdY) { - @Override - public boolean isValidPoint(double x, double y) { - if (x < xval[1] || - x > xval[xval.length - 2] || - y < yval[1] || - y > yval[yval.length - 2]) { - return false; - } else { - return true; - } - } - }; - } -}
