http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java deleted file mode 100644 index 4260606..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolatingFunction.java +++ /dev/null @@ -1,482 +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.TrivariateFunction; -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/Tricubic_interpolation";> - * tricubic spline interpolation</a>, as proposed in - * <quote> - * Tricubic interpolation in three dimensions<br/> - * F. Lekien and J. Marsden<br/> - * <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471 - * </quote> - * - * @since 2.2 - * @deprecated To be removed in 4.0 (see MATH-1166). - */ -@Deprecated -public class TricubicSplineInterpolatingFunction - implements TrivariateFunction { - /** - * 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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 }, - {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,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 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 }, - { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 }, - { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 }, - { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 }, - { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 }, - { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 }, - { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 }, - { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 }, - { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 }, - { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 }, - { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 }, - { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 }, - { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 }, - { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 } - }; - - /** Samples x-coordinates */ - private final double[] xval; - /** Samples y-coordinates */ - private final double[] yval; - /** Samples z-coordinates */ - private final double[] zval; - /** Set of cubic splines pacthing the whole data grid */ - private final TricubicSplineFunction[][][] splines; - - /** - * @param x Sample values of the x-coordinate, in increasing order. - * @param y Sample values of the y-coordinate, in increasing order. - * @param z 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 dFdZ Values of the partial derivative of function with respect to z on every grid point. - * @param d2FdXdY Values of the cross partial derivative of function on every grid point. - * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. - * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. - * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. - * @throws NoDataException if any of the arrays has zero length. - * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. - * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. - */ - public TricubicSplineInterpolatingFunction(double[] x, - double[] y, - double[] z, - double[][][] f, - double[][][] dFdX, - double[][][] dFdY, - double[][][] dFdZ, - double[][][] d2FdXdY, - double[][][] d2FdXdZ, - double[][][] d2FdYdZ, - double[][][] d3FdXdYdZ) - throws NoDataException, - DimensionMismatchException, - NonMonotonicSequenceException { - final int xLen = x.length; - final int yLen = y.length; - final int zLen = z.length; - - if (xLen == 0 || yLen == 0 || z.length == 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 != dFdZ.length) { - throw new DimensionMismatchException(xLen, dFdZ.length); - } - if (xLen != d2FdXdY.length) { - throw new DimensionMismatchException(xLen, d2FdXdY.length); - } - if (xLen != d2FdXdZ.length) { - throw new DimensionMismatchException(xLen, d2FdXdZ.length); - } - if (xLen != d2FdYdZ.length) { - throw new DimensionMismatchException(xLen, d2FdYdZ.length); - } - if (xLen != d3FdXdYdZ.length) { - throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); - } - - MathArrays.checkOrder(x); - MathArrays.checkOrder(y); - MathArrays.checkOrder(z); - - xval = x.clone(); - yval = y.clone(); - zval = z.clone(); - - final int lastI = xLen - 1; - final int lastJ = yLen - 1; - final int lastK = zLen - 1; - splines = new TricubicSplineFunction[lastI][lastJ][lastK]; - - 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 (dFdZ[i].length != yLen) { - throw new DimensionMismatchException(dFdZ[i].length, yLen); - } - if (d2FdXdY[i].length != yLen) { - throw new DimensionMismatchException(d2FdXdY[i].length, yLen); - } - if (d2FdXdZ[i].length != yLen) { - throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); - } - if (d2FdYdZ[i].length != yLen) { - throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); - } - if (d3FdXdYdZ[i].length != yLen) { - throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); - } - - final int ip1 = i + 1; - for (int j = 0; j < lastJ; j++) { - if (f[i][j].length != zLen) { - throw new DimensionMismatchException(f[i][j].length, zLen); - } - if (dFdX[i][j].length != zLen) { - throw new DimensionMismatchException(dFdX[i][j].length, zLen); - } - if (dFdY[i][j].length != zLen) { - throw new DimensionMismatchException(dFdY[i][j].length, zLen); - } - if (dFdZ[i][j].length != zLen) { - throw new DimensionMismatchException(dFdZ[i][j].length, zLen); - } - if (d2FdXdY[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen); - } - if (d2FdXdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen); - } - if (d2FdYdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen); - } - if (d3FdXdYdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen); - } - - final int jp1 = j + 1; - for (int k = 0; k < lastK; k++) { - final int kp1 = k + 1; - - final double[] beta = new double[] { - f[i][j][k], f[ip1][j][k], - f[i][jp1][k], f[ip1][jp1][k], - f[i][j][kp1], f[ip1][j][kp1], - f[i][jp1][kp1], f[ip1][jp1][kp1], - - dFdX[i][j][k], dFdX[ip1][j][k], - dFdX[i][jp1][k], dFdX[ip1][jp1][k], - dFdX[i][j][kp1], dFdX[ip1][j][kp1], - dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1], - - dFdY[i][j][k], dFdY[ip1][j][k], - dFdY[i][jp1][k], dFdY[ip1][jp1][k], - dFdY[i][j][kp1], dFdY[ip1][j][kp1], - dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1], - - dFdZ[i][j][k], dFdZ[ip1][j][k], - dFdZ[i][jp1][k], dFdZ[ip1][jp1][k], - dFdZ[i][j][kp1], dFdZ[ip1][j][kp1], - dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1], - - d2FdXdY[i][j][k], d2FdXdY[ip1][j][k], - d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k], - d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1], - d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1], - - d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k], - d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k], - d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1], - d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1], - - d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k], - d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k], - d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1], - d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1], - - d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k], - d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k], - d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1], - d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1], - }; - - splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta)); - } - } - } - } - - /** - * {@inheritDoc} - * - * @throws OutOfRangeException if any of the variables is outside its interpolation range. - */ - public double value(double x, double y, double z) - throws OutOfRangeException { - final int i = searchIndex(x, xval); - if (i == -1) { - throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); - } - final int j = searchIndex(y, yval); - if (j == -1) { - throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); - } - final int k = searchIndex(z, zval); - if (k == -1) { - throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); - } - - final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); - final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); - final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); - - return splines[i][j][k].value(xN, yN, zN); - } - - /** - * @param c Coordinate. - * @param val Coordinate samples. - * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} - * if {@code c} is out of the range defined by the end values of {@code val}. - */ - private int searchIndex(double c, double[] val) { - if (c < val[0]) { - return -1; - } - - final int max = val.length; - for (int i = 1; i < max; i++) { - if (c <= val[i]) { - return i - 1; - } - } - - return -1; - } - - /** - * 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,0)</li> - * <li>f(1,0,0)</li> - * <li>f(0,1,0)</li> - * <li>f(1,1,0)</li> - * <li>f(0,0,1)</li> - * <li>f(1,0,1)</li> - * <li>f(0,1,1)</li> - * <li>f(1,1,1)</li> - * - * <li>f<sub>x</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>x</sub>(1,1,1)</li> - * - * <li>f<sub>y</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>y</sub>(1,1,1)</li> - * - * <li>f<sub>z</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>z</sub>(1,1,1)</li> - * - * <li>f<sub>xy</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xy</sub>(1,1,1)</li> - * - * <li>f<sub>xz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xz</sub>(1,1,1)</li> - * - * <li>f<sub>yz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>yz</sub>(1,1,1)</li> - * - * <li>f<sub>xyz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xyz</sub>(1,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 int sz = 64; - final double[] a = new double[sz]; - - for (int i = 0; i < sz; i++) { - double result = 0; - final double[] row = AINV[i]; - for (int j = 0; j < sz; j++) { - result += row[j] * beta[j]; - } - a[i] = result; - } - - return a; - } -} - -/** - * 3D-spline function. - * - */ -class TricubicSplineFunction - implements TrivariateFunction { - /** Number of points. */ - private static final short N = 4; - /** Coefficients */ - private final double[][][] a = new double[N][N][N]; - - /** - * @param aV List of spline coefficients. - */ - public TricubicSplineFunction(double[] aV) { - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - a[i][j][k] = aV[i + N * (j + N * k)]; - } - } - } - } - - /** - * @param x x-coordinate of the interpolation point. - * @param y y-coordinate of the interpolation point. - * @param z z-coordinate of the interpolation point. - * @return the interpolated value. - * @throws OutOfRangeException if {@code x}, {@code y} or - * {@code z} are not in the interval {@code [0, 1]}. - */ - public double value(double x, double y, double z) - throws OutOfRangeException { - if (x < 0 || x > 1) { - throw new OutOfRangeException(x, 0, 1); - } - if (y < 0 || y > 1) { - throw new OutOfRangeException(y, 0, 1); - } - if (z < 0 || z > 1) { - throw new OutOfRangeException(z, 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 }; - - final double z2 = z * z; - final double z3 = z2 * z; - final double[] pZ = { 1, z, z2, z3 }; - - double result = 0; - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; - } - } - } - - return result; - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java deleted file mode 100644 index da19986..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicSplineInterpolator.java +++ /dev/null @@ -1,201 +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 tricubic interpolating function. - * - * @since 2.2 - * @deprecated To be removed in 4.0 (see MATH-1166). - */ -@Deprecated -public class TricubicSplineInterpolator - implements TrivariateGridInterpolator { - /** - * {@inheritDoc} - */ - public TricubicSplineInterpolatingFunction interpolate(final double[] xval, - final double[] yval, - final double[] zval, - final double[][][] fval) - throws NoDataException, NumberIsTooSmallException, - DimensionMismatchException, NonMonotonicSequenceException { - if (xval.length == 0 || yval.length == 0 || zval.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); - MathArrays.checkOrder(zval); - - final int xLen = xval.length; - final int yLen = yval.length; - final int zLen = zval.length; - - // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets - // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k]) - // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k]) - final double[][][] fvalXY = new double[zLen][xLen][yLen]; - final double[][][] fvalZX = new double[yLen][zLen][xLen]; - for (int i = 0; i < xLen; i++) { - if (fval[i].length != yLen) { - throw new DimensionMismatchException(fval[i].length, yLen); - } - - for (int j = 0; j < yLen; j++) { - if (fval[i][j].length != zLen) { - throw new DimensionMismatchException(fval[i][j].length, zLen); - } - - for (int k = 0; k < zLen; k++) { - final double v = fval[i][j][k]; - fvalXY[k][i][j] = v; - fvalZX[j][k][i] = v; - } - } - } - - final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true); - - // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z - final BicubicSplineInterpolatingFunction[] xSplineYZ - = new BicubicSplineInterpolatingFunction[xLen]; - for (int i = 0; i < xLen; i++) { - xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); - } - - // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x - final BicubicSplineInterpolatingFunction[] ySplineZX - = new BicubicSplineInterpolatingFunction[yLen]; - for (int j = 0; j < yLen; j++) { - ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); - } - - // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y - final BicubicSplineInterpolatingFunction[] zSplineXY - = new BicubicSplineInterpolatingFunction[zLen]; - for (int k = 0; k < zLen; k++) { - zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); - } - - // Partial derivatives wrt x and wrt y - final double[][][] dFdX = new double[xLen][yLen][zLen]; - final double[][][] dFdY = new double[xLen][yLen][zLen]; - final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; - for (int k = 0; k < zLen; k++) { - final BicubicSplineInterpolatingFunction f = zSplineXY[k]; - for (int i = 0; i < xLen; i++) { - final double x = xval[i]; - for (int j = 0; j < yLen; j++) { - final double y = yval[j]; - dFdX[i][j][k] = f.partialDerivativeX(x, y); - dFdY[i][j][k] = f.partialDerivativeY(x, y); - d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); - } - } - } - - // Partial derivatives wrt y and wrt z - final double[][][] dFdZ = new double[xLen][yLen][zLen]; - final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; - for (int i = 0; i < xLen; i++) { - final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; - for (int j = 0; j < yLen; j++) { - final double y = yval[j]; - for (int k = 0; k < zLen; k++) { - final double z = zval[k]; - dFdZ[i][j][k] = f.partialDerivativeY(y, z); - d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); - } - } - } - - // Partial derivatives wrt x and wrt z - final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; - for (int j = 0; j < yLen; j++) { - final BicubicSplineInterpolatingFunction f = ySplineZX[j]; - for (int k = 0; k < zLen; k++) { - final double z = zval[k]; - for (int i = 0; i < xLen; i++) { - final double x = xval[i]; - d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); - } - } - } - - // Third partial cross-derivatives - final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; - for (int i = 0; i < xLen ; i++) { - final int nI = nextIndex(i, xLen); - final int pI = previousIndex(i); - for (int j = 0; j < yLen; j++) { - final int nJ = nextIndex(j, yLen); - final int pJ = previousIndex(j); - for (int k = 0; k < zLen; k++) { - final int nK = nextIndex(k, zLen); - final int pK = previousIndex(k); - - // XXX Not sure about this formula - d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - - fval[pI][nJ][nK] + fval[pI][pJ][nK] - - fval[nI][nJ][pK] + fval[nI][pJ][pK] + - fval[pI][nJ][pK] - fval[pI][pJ][pK]) / - ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ; - } - } - } - - // Create the interpolating splines - return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, - dFdX, dFdY, dFdZ, - d2FdXdY, d2FdZdX, d2FdYdZ, - d3FdXdYdZ); - } - - /** - * Compute the next index of an array, clipping if necessary. - * It is assumed (but not checked) that {@code i} is larger than or equal to 0}. - * - * @param i Index - * @param max Upper limit of the array - * @return the next index - */ - private int nextIndex(int i, int max) { - final int index = i + 1; - return index < max ? index : index - 1; - } - /** - * Compute the previous index of an array, clipping if necessary. - * It is assumed (but not checked) that {@code i} is smaller than the size of the array. - * - * @param i Index - * @return the previous index - */ - private int previousIndex(int i) { - final int index = i - 1; - return index >= 0 ? index : 0; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java deleted file mode 100644 index ec69715..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TrivariateGridInterpolator.java +++ /dev/null @@ -1,54 +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.TrivariateFunction; -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; - -/** - * Interface representing a trivariate real interpolating function where the - * sample points must be specified on a regular grid. - * - * @since 2.2 - */ -public interface TrivariateGridInterpolator { - /** - * Compute an interpolating function for the dataset. - * - * @param xval All the x-coordinates of the interpolation points, sorted - * in increasing order. - * @param yval All the y-coordinates of the interpolation points, sorted - * in increasing order. - * @param zval All the z-coordinates of the interpolation points, sorted - * in increasing order. - * @param fval the values of the interpolation points on all the grid knots: - * {@code fval[i][j][k] = f(xval[i], yval[j], zval[k])}. - * @return a function that interpolates the data set. - * @throws NoDataException if any of the arrays has zero length. - * @throws DimensionMismatchException if the array lengths are inconsistent. - * @throws NonMonotonicSequenceException if arrays are not sorted - * @throws NumberIsTooSmallException if the number of points is too small for - * the order of the interpolation - */ - TrivariateFunction interpolate(double[] xval, double[] yval, double[] zval, - double[][][] fval) - throws NoDataException, NumberIsTooSmallException, - DimensionMismatchException, NonMonotonicSequenceException; -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java deleted file mode 100644 index f7a1bd1..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariateInterpolator.java +++ /dev/null @@ -1,41 +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.UnivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.MathIllegalArgumentException; - -/** - * Interface representing a univariate real interpolating function. - * - */ -public interface UnivariateInterpolator { - /** - * Compute an interpolating function for the dataset. - * - * @param xval Arguments for the interpolation points. - * @param yval Values for the interpolation points. - * @return a function which interpolates the dataset. - * @throws MathIllegalArgumentException - * if the arguments violate assumptions made by the interpolation - * algorithm. - * @throws DimensionMismatchException if arrays lengthes do not match - */ - UnivariateFunction interpolate(double xval[], double yval[]) - throws MathIllegalArgumentException, DimensionMismatchException; -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java deleted file mode 100644 index 6b788b1..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/UnivariatePeriodicInterpolator.java +++ /dev/null @@ -1,123 +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.UnivariateFunction; -import org.apache.commons.math3.util.MathUtils; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.exception.MathIllegalArgumentException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; - -/** - * Adapter for classes implementing the {@link UnivariateInterpolator} - * interface. - * The data to be interpolated is assumed to be periodic. Thus values that are - * outside of the range can be passed to the interpolation function: They will - * be wrapped into the initial range before being passed to the class that - * actually computes the interpolation. - * - */ -public class UnivariatePeriodicInterpolator - implements UnivariateInterpolator { - /** Default number of extension points of the samples array. */ - public static final int DEFAULT_EXTEND = 5; - /** Interpolator. */ - private final UnivariateInterpolator interpolator; - /** Period. */ - private final double period; - /** Number of extension points. */ - private final int extend; - - /** - * Builds an interpolator. - * - * @param interpolator Interpolator. - * @param period Period. - * @param extend Number of points to be appended at the beginning and - * end of the sample arrays in order to avoid interpolation failure at - * the (periodic) boundaries of the orginal interval. The value is the - * number of sample points which the original {@code interpolator} needs - * on each side of the interpolated point. - */ - public UnivariatePeriodicInterpolator(UnivariateInterpolator interpolator, - double period, - int extend) { - this.interpolator = interpolator; - this.period = period; - this.extend = extend; - } - - /** - * Builds an interpolator. - * Uses {@link #DEFAULT_EXTEND} as the number of extension points on each side - * of the original abscissae range. - * - * @param interpolator Interpolator. - * @param period Period. - */ - public UnivariatePeriodicInterpolator(UnivariateInterpolator interpolator, - double period) { - this(interpolator, period, DEFAULT_EXTEND); - } - - /** - * {@inheritDoc} - * - * @throws NumberIsTooSmallException if the number of extension points - * is larger than the size of {@code xval}. - */ - public UnivariateFunction interpolate(double[] xval, - double[] yval) - throws NumberIsTooSmallException, NonMonotonicSequenceException { - if (xval.length < extend) { - throw new NumberIsTooSmallException(xval.length, extend, true); - } - - MathArrays.checkOrder(xval); - final double offset = xval[0]; - - final int len = xval.length + extend * 2; - final double[] x = new double[len]; - final double[] y = new double[len]; - for (int i = 0; i < xval.length; i++) { - final int index = i + extend; - x[index] = MathUtils.reduce(xval[i], period, offset); - y[index] = yval[i]; - } - - // Wrap to enable interpolation at the boundaries. - for (int i = 0; i < extend; i++) { - int index = xval.length - extend + i; - x[i] = MathUtils.reduce(xval[index], period, offset) - period; - y[i] = yval[index]; - - index = len - extend + i; - x[index] = MathUtils.reduce(xval[i], period, offset) + period; - y[index] = yval[i]; - } - - MathArrays.sortInPlace(x, y); - - final UnivariateFunction f = interpolator.interpolate(x, y); - return new UnivariateFunction() { - public double value(final double x) throws MathIllegalArgumentException { - return f.value(MathUtils.reduce(x, period, offset)); - } - }; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/package-info.java deleted file mode 100644 index b4b25dd..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/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. - */ -/** - * - * Univariate real functions interpolation algorithms. - * - */ -package org.apache.commons.math3.analysis.interpolation; http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/package-info.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/package-info.java b/src/main/java/org/apache/commons/math3/analysis/package-info.java deleted file mode 100644 index 46e0477..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/package-info.java +++ /dev/null @@ -1,32 +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. - */ -/** - * - * <p> - * Parent package for common numerical analysis procedures, including root finding, - * function interpolation and integration. Note that optimization (i.e. minimization - * and maximization) is a separate top-level package. - * </p> - * <p> - * Function interfaces are intended to be implemented by user code to represent - * domain problems. The algorithms provided by the library operate on these - * functions to find their roots, or integrate them, or ... Functions can be multivariate - * or univariate, real vectorial or matrix-valued, and they can be differentiable or not. - * </p> - * - */ -package org.apache.commons.math3.analysis; http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunction.java b/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunction.java deleted file mode 100644 index d424a88..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunction.java +++ /dev/null @@ -1,412 +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.polynomials; - -import java.io.Serializable; -import java.util.Arrays; - -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.analysis.ParametricUnivariateFunction; -import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; -import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathUtils; - -/** - * Immutable representation of a real polynomial function with real coefficients. - * <p> - * <a href="http://mathworld.wolfram.com/HornersMethod.html";>Horner's Method</a> - * is used to evaluate the function.</p> - * - */ -public class PolynomialFunction implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction, Serializable { - /** - * Serialization identifier - */ - private static final long serialVersionUID = -7726511984200295583L; - /** - * The coefficients of the polynomial, ordered by degree -- i.e., - * coefficients[0] is the constant term and coefficients[n] is the - * coefficient of x^n where n is the degree of the polynomial. - */ - private final double coefficients[]; - - /** - * Construct a polynomial with the given coefficients. The first element - * of the coefficients array is the constant term. Higher degree - * coefficients follow in sequence. The degree of the resulting polynomial - * is the index of the last non-null element of the array, or 0 if all elements - * are null. - * <p> - * The constructor makes a copy of the input array and assigns the copy to - * the coefficients property.</p> - * - * @param c Polynomial coefficients. - * @throws NullArgumentException if {@code c} is {@code null}. - * @throws NoDataException if {@code c} is empty. - */ - public PolynomialFunction(double c[]) - throws NullArgumentException, NoDataException { - super(); - MathUtils.checkNotNull(c); - int n = c.length; - if (n == 0) { - throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); - } - while ((n > 1) && (c[n - 1] == 0)) { - --n; - } - this.coefficients = new double[n]; - System.arraycopy(c, 0, this.coefficients, 0, n); - } - - /** - * Compute the value of the function for the given argument. - * <p> - * The value returned is <br/> - * <code>coefficients[n] * x^n + ... + coefficients[1] * x + coefficients[0]</code> - * </p> - * - * @param x Argument for which the function value should be computed. - * @return the value of the polynomial at the given point. - * @see UnivariateFunction#value(double) - */ - public double value(double x) { - return evaluate(coefficients, x); - } - - /** - * Returns the degree of the polynomial. - * - * @return the degree of the polynomial. - */ - public int degree() { - return coefficients.length - 1; - } - - /** - * Returns a copy of the coefficients array. - * <p> - * Changes made to the returned copy will not affect the coefficients of - * the polynomial.</p> - * - * @return a fresh copy of the coefficients array. - */ - public double[] getCoefficients() { - return coefficients.clone(); - } - - /** - * Uses Horner's Method to evaluate the polynomial with the given coefficients at - * the argument. - * - * @param coefficients Coefficients of the polynomial to evaluate. - * @param argument Input value. - * @return the value of the polynomial. - * @throws NoDataException if {@code coefficients} is empty. - * @throws NullArgumentException if {@code coefficients} is {@code null}. - */ - protected static double evaluate(double[] coefficients, double argument) - throws NullArgumentException, NoDataException { - MathUtils.checkNotNull(coefficients); - int n = coefficients.length; - if (n == 0) { - throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); - } - double result = coefficients[n - 1]; - for (int j = n - 2; j >= 0; j--) { - result = argument * result + coefficients[j]; - } - return result; - } - - - /** {@inheritDoc} - * @since 3.1 - * @throws NoDataException if {@code coefficients} is empty. - * @throws NullArgumentException if {@code coefficients} is {@code null}. - */ - public DerivativeStructure value(final DerivativeStructure t) - throws NullArgumentException, NoDataException { - MathUtils.checkNotNull(coefficients); - int n = coefficients.length; - if (n == 0) { - throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); - } - DerivativeStructure result = - new DerivativeStructure(t.getFreeParameters(), t.getOrder(), coefficients[n - 1]); - for (int j = n - 2; j >= 0; j--) { - result = result.multiply(t).add(coefficients[j]); - } - return result; - } - - /** - * Add a polynomial to the instance. - * - * @param p Polynomial to add. - * @return a new polynomial which is the sum of the instance and {@code p}. - */ - public PolynomialFunction add(final PolynomialFunction p) { - // identify the lowest degree polynomial - final int lowLength = FastMath.min(coefficients.length, p.coefficients.length); - final int highLength = FastMath.max(coefficients.length, p.coefficients.length); - - // build the coefficients array - double[] newCoefficients = new double[highLength]; - for (int i = 0; i < lowLength; ++i) { - newCoefficients[i] = coefficients[i] + p.coefficients[i]; - } - System.arraycopy((coefficients.length < p.coefficients.length) ? - p.coefficients : coefficients, - lowLength, - newCoefficients, lowLength, - highLength - lowLength); - - return new PolynomialFunction(newCoefficients); - } - - /** - * Subtract a polynomial from the instance. - * - * @param p Polynomial to subtract. - * @return a new polynomial which is the difference the instance minus {@code p}. - */ - public PolynomialFunction subtract(final PolynomialFunction p) { - // identify the lowest degree polynomial - int lowLength = FastMath.min(coefficients.length, p.coefficients.length); - int highLength = FastMath.max(coefficients.length, p.coefficients.length); - - // build the coefficients array - double[] newCoefficients = new double[highLength]; - for (int i = 0; i < lowLength; ++i) { - newCoefficients[i] = coefficients[i] - p.coefficients[i]; - } - if (coefficients.length < p.coefficients.length) { - for (int i = lowLength; i < highLength; ++i) { - newCoefficients[i] = -p.coefficients[i]; - } - } else { - System.arraycopy(coefficients, lowLength, newCoefficients, lowLength, - highLength - lowLength); - } - - return new PolynomialFunction(newCoefficients); - } - - /** - * Negate the instance. - * - * @return a new polynomial. - */ - public PolynomialFunction negate() { - double[] newCoefficients = new double[coefficients.length]; - for (int i = 0; i < coefficients.length; ++i) { - newCoefficients[i] = -coefficients[i]; - } - return new PolynomialFunction(newCoefficients); - } - - /** - * Multiply the instance by a polynomial. - * - * @param p Polynomial to multiply by. - * @return a new polynomial. - */ - public PolynomialFunction multiply(final PolynomialFunction p) { - double[] newCoefficients = new double[coefficients.length + p.coefficients.length - 1]; - - for (int i = 0; i < newCoefficients.length; ++i) { - newCoefficients[i] = 0.0; - for (int j = FastMath.max(0, i + 1 - p.coefficients.length); - j < FastMath.min(coefficients.length, i + 1); - ++j) { - newCoefficients[i] += coefficients[j] * p.coefficients[i-j]; - } - } - - return new PolynomialFunction(newCoefficients); - } - - /** - * Returns the coefficients of the derivative of the polynomial with the given coefficients. - * - * @param coefficients Coefficients of the polynomial to differentiate. - * @return the coefficients of the derivative or {@code null} if coefficients has length 1. - * @throws NoDataException if {@code coefficients} is empty. - * @throws NullArgumentException if {@code coefficients} is {@code null}. - */ - protected static double[] differentiate(double[] coefficients) - throws NullArgumentException, NoDataException { - MathUtils.checkNotNull(coefficients); - int n = coefficients.length; - if (n == 0) { - throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); - } - if (n == 1) { - return new double[]{0}; - } - double[] result = new double[n - 1]; - for (int i = n - 1; i > 0; i--) { - result[i - 1] = i * coefficients[i]; - } - return result; - } - - /** - * Returns the derivative as a {@link PolynomialFunction}. - * - * @return the derivative polynomial. - */ - public PolynomialFunction polynomialDerivative() { - return new PolynomialFunction(differentiate(coefficients)); - } - - /** - * Returns the derivative as a {@link UnivariateFunction}. - * - * @return the derivative function. - */ - public UnivariateFunction derivative() { - return polynomialDerivative(); - } - - /** - * Returns a string representation of the polynomial. - * - * <p>The representation is user oriented. Terms are displayed lowest - * degrees first. The multiplications signs, coefficients equals to - * one and null terms are not displayed (except if the polynomial is 0, - * in which case the 0 constant term is displayed). Addition of terms - * with negative coefficients are replaced by subtraction of terms - * with positive coefficients except for the first displayed term - * (i.e. we display <code>-3</code> for a constant negative polynomial, - * but <code>1 - 3 x + x^2</code> if the negative coefficient is not - * the first one displayed).</p> - * - * @return a string representation of the polynomial. - */ - @Override - public String toString() { - StringBuilder s = new StringBuilder(); - if (coefficients[0] == 0.0) { - if (coefficients.length == 1) { - return "0"; - } - } else { - s.append(toString(coefficients[0])); - } - - for (int i = 1; i < coefficients.length; ++i) { - if (coefficients[i] != 0) { - if (s.length() > 0) { - if (coefficients[i] < 0) { - s.append(" - "); - } else { - s.append(" + "); - } - } else { - if (coefficients[i] < 0) { - s.append("-"); - } - } - - double absAi = FastMath.abs(coefficients[i]); - if ((absAi - 1) != 0) { - s.append(toString(absAi)); - s.append(' '); - } - - s.append("x"); - if (i > 1) { - s.append('^'); - s.append(Integer.toString(i)); - } - } - } - - return s.toString(); - } - - /** - * Creates a string representing a coefficient, removing ".0" endings. - * - * @param coeff Coefficient. - * @return a string representation of {@code coeff}. - */ - private static String toString(double coeff) { - final String c = Double.toString(coeff); - if (c.endsWith(".0")) { - return c.substring(0, c.length() - 2); - } else { - return c; - } - } - - /** {@inheritDoc} */ - @Override - public int hashCode() { - final int prime = 31; - int result = 1; - result = prime * result + Arrays.hashCode(coefficients); - return result; - } - - /** {@inheritDoc} */ - @Override - public boolean equals(Object obj) { - if (this == obj) { - return true; - } - if (!(obj instanceof PolynomialFunction)) { - return false; - } - PolynomialFunction other = (PolynomialFunction) obj; - if (!Arrays.equals(coefficients, other.coefficients)) { - return false; - } - return true; - } - - /** - * Dedicated parametric polynomial class. - * - * @since 3.0 - */ - public static class Parametric implements ParametricUnivariateFunction { - /** {@inheritDoc} */ - public double[] gradient(double x, double ... parameters) { - final double[] gradient = new double[parameters.length]; - double xn = 1.0; - for (int i = 0; i < parameters.length; ++i) { - gradient[i] = xn; - xn *= x; - } - return gradient; - } - - /** {@inheritDoc} */ - public double value(final double x, final double ... parameters) - throws NoDataException { - return PolynomialFunction.evaluate(parameters, x); - } - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionLagrangeForm.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionLagrangeForm.java b/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionLagrangeForm.java deleted file mode 100644 index 9d812df..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionLagrangeForm.java +++ /dev/null @@ -1,326 +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.polynomials; - -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.util.FastMath; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.util.LocalizedFormats; - -/** - * Implements the representation of a real polynomial function in - * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html";> - * Lagrange Form</a>. For reference, see <b>Introduction to Numerical - * Analysis</b>, ISBN 038795452X, chapter 2. - * <p> - * The approximated function should be smooth enough for Lagrange polynomial - * to work well. Otherwise, consider using splines instead.</p> - * - * @since 1.2 - */ -public class PolynomialFunctionLagrangeForm implements UnivariateFunction { - /** - * The coefficients of the polynomial, ordered by degree -- i.e. - * coefficients[0] is the constant term and coefficients[n] is the - * coefficient of x^n where n is the degree of the polynomial. - */ - private double coefficients[]; - /** - * Interpolating points (abscissas). - */ - private final double x[]; - /** - * Function values at interpolating points. - */ - private final double y[]; - /** - * Whether the polynomial coefficients are available. - */ - private boolean coefficientsComputed; - - /** - * Construct a Lagrange polynomial with the given abscissas and function - * values. The order of interpolating points are not important. - * <p> - * The constructor makes copy of the input arrays and assigns them.</p> - * - * @param x interpolating points - * @param y function values at interpolating points - * @throws DimensionMismatchException if the array lengths are different. - * @throws NumberIsTooSmallException if the number of points is less than 2. - * @throws NonMonotonicSequenceException - * if two abscissae have the same value. - */ - public PolynomialFunctionLagrangeForm(double x[], double y[]) - throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException { - this.x = new double[x.length]; - this.y = new double[y.length]; - System.arraycopy(x, 0, this.x, 0, x.length); - System.arraycopy(y, 0, this.y, 0, y.length); - coefficientsComputed = false; - - if (!verifyInterpolationArray(x, y, false)) { - MathArrays.sortInPlace(this.x, this.y); - // Second check in case some abscissa is duplicated. - verifyInterpolationArray(this.x, this.y, true); - } - } - - /** - * Calculate the function value at the given point. - * - * @param z Point at which the function value is to be computed. - * @return the function value. - * @throws DimensionMismatchException if {@code x} and {@code y} have - * different lengths. - * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException - * if {@code x} is not sorted in strictly increasing order. - * @throws NumberIsTooSmallException if the size of {@code x} is less - * than 2. - */ - public double value(double z) { - return evaluateInternal(x, y, z); - } - - /** - * Returns the degree of the polynomial. - * - * @return the degree of the polynomial - */ - public int degree() { - return x.length - 1; - } - - /** - * Returns a copy of the interpolating points array. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * - * @return a fresh copy of the interpolating points array - */ - public double[] getInterpolatingPoints() { - double[] out = new double[x.length]; - System.arraycopy(x, 0, out, 0, x.length); - return out; - } - - /** - * Returns a copy of the interpolating values array. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * - * @return a fresh copy of the interpolating values array - */ - public double[] getInterpolatingValues() { - double[] out = new double[y.length]; - System.arraycopy(y, 0, out, 0, y.length); - return out; - } - - /** - * Returns a copy of the coefficients array. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * <p> - * Note that coefficients computation can be ill-conditioned. Use with caution - * and only when it is necessary.</p> - * - * @return a fresh copy of the coefficients array - */ - public double[] getCoefficients() { - if (!coefficientsComputed) { - computeCoefficients(); - } - double[] out = new double[coefficients.length]; - System.arraycopy(coefficients, 0, out, 0, coefficients.length); - return out; - } - - /** - * Evaluate the Lagrange polynomial using - * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html";> - * Neville's Algorithm</a>. It takes O(n^2) time. - * - * @param x Interpolating points array. - * @param y Interpolating values array. - * @param z Point at which the function value is to be computed. - * @return the function value. - * @throws DimensionMismatchException if {@code x} and {@code y} have - * different lengths. - * @throws NonMonotonicSequenceException - * if {@code x} is not sorted in strictly increasing order. - * @throws NumberIsTooSmallException if the size of {@code x} is less - * than 2. - */ - public static double evaluate(double x[], double y[], double z) - throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException { - if (verifyInterpolationArray(x, y, false)) { - return evaluateInternal(x, y, z); - } - - // Array is not sorted. - final double[] xNew = new double[x.length]; - final double[] yNew = new double[y.length]; - System.arraycopy(x, 0, xNew, 0, x.length); - System.arraycopy(y, 0, yNew, 0, y.length); - - MathArrays.sortInPlace(xNew, yNew); - // Second check in case some abscissa is duplicated. - verifyInterpolationArray(xNew, yNew, true); - return evaluateInternal(xNew, yNew, z); - } - - /** - * Evaluate the Lagrange polynomial using - * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html";> - * Neville's Algorithm</a>. It takes O(n^2) time. - * - * @param x Interpolating points array. - * @param y Interpolating values array. - * @param z Point at which the function value is to be computed. - * @return the function value. - * @throws DimensionMismatchException if {@code x} and {@code y} have - * different lengths. - * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException - * if {@code x} is not sorted in strictly increasing order. - * @throws NumberIsTooSmallException if the size of {@code x} is less - * than 2. - */ - private static double evaluateInternal(double x[], double y[], double z) { - int nearest = 0; - final int n = x.length; - final double[] c = new double[n]; - final double[] d = new double[n]; - double min_dist = Double.POSITIVE_INFINITY; - for (int i = 0; i < n; i++) { - // initialize the difference arrays - c[i] = y[i]; - d[i] = y[i]; - // find out the abscissa closest to z - final double dist = FastMath.abs(z - x[i]); - if (dist < min_dist) { - nearest = i; - min_dist = dist; - } - } - - // initial approximation to the function value at z - double value = y[nearest]; - - for (int i = 1; i < n; i++) { - for (int j = 0; j < n-i; j++) { - final double tc = x[j] - z; - final double td = x[i+j] - z; - final double divider = x[j] - x[i+j]; - // update the difference arrays - final double w = (c[j+1] - d[j]) / divider; - c[j] = tc * w; - d[j] = td * w; - } - // sum up the difference terms to get the final value - if (nearest < 0.5*(n-i+1)) { - value += c[nearest]; // fork down - } else { - nearest--; - value += d[nearest]; // fork up - } - } - - return value; - } - - /** - * Calculate the coefficients of Lagrange polynomial from the - * interpolation data. It takes O(n^2) time. - * Note that this computation can be ill-conditioned: Use with caution - * and only when it is necessary. - */ - protected void computeCoefficients() { - final int n = degree() + 1; - coefficients = new double[n]; - for (int i = 0; i < n; i++) { - coefficients[i] = 0.0; - } - - // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1]) - final double[] c = new double[n+1]; - c[0] = 1.0; - for (int i = 0; i < n; i++) { - for (int j = i; j > 0; j--) { - c[j] = c[j-1] - c[j] * x[i]; - } - c[0] *= -x[i]; - c[i+1] = 1; - } - - final double[] tc = new double[n]; - for (int i = 0; i < n; i++) { - // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1]) - double d = 1; - for (int j = 0; j < n; j++) { - if (i != j) { - d *= x[i] - x[j]; - } - } - final double t = y[i] / d; - // Lagrange polynomial is the sum of n terms, each of which is a - // polynomial of degree n-1. tc[] are the coefficients of the i-th - // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]). - tc[n-1] = c[n]; // actually c[n] = 1 - coefficients[n-1] += t * tc[n-1]; - for (int j = n-2; j >= 0; j--) { - tc[j] = c[j+1] + tc[j+1] * x[i]; - coefficients[j] += t * tc[j]; - } - } - - coefficientsComputed = true; - } - - /** - * Check that the interpolation arrays are valid. - * The arrays features checked by this method are that both arrays have the - * same length and this length is at least 2. - * - * @param x Interpolating points array. - * @param y Interpolating values array. - * @param abort Whether to throw an exception if {@code x} is not sorted. - * @throws DimensionMismatchException if the array lengths are different. - * @throws NumberIsTooSmallException if the number of points is less than 2. - * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException - * if {@code x} is not sorted in strictly increasing order and {@code abort} - * is {@code true}. - * @return {@code false} if the {@code x} is not sorted in increasing order, - * {@code true} otherwise. - * @see #evaluate(double[], double[], double) - * @see #computeCoefficients() - */ - public static boolean verifyInterpolationArray(double x[], double y[], boolean abort) - throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException { - if (x.length != y.length) { - throw new DimensionMismatchException(x.length, y.length); - } - if (x.length < 2) { - throw new NumberIsTooSmallException(LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true); - } - - return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java b/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java deleted file mode 100644 index fc2f1fd..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/polynomials/PolynomialFunctionNewtonForm.java +++ /dev/null @@ -1,245 +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.polynomials; - -import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; -import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.util.MathUtils; - -/** - * Implements the representation of a real polynomial function in - * Newton Form. For reference, see <b>Elementary Numerical Analysis</b>, - * ISBN 0070124477, chapter 2. - * <p> - * The formula of polynomial in Newton form is - * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... + - * a[n](x-c[0])(x-c[1])...(x-c[n-1]) - * Note that the length of a[] is one more than the length of c[]</p> - * - * @since 1.2 - */ -public class PolynomialFunctionNewtonForm implements UnivariateDifferentiableFunction { - - /** - * The coefficients of the polynomial, ordered by degree -- i.e. - * coefficients[0] is the constant term and coefficients[n] is the - * coefficient of x^n where n is the degree of the polynomial. - */ - private double coefficients[]; - - /** - * Centers of the Newton polynomial. - */ - private final double c[]; - - /** - * When all c[i] = 0, a[] becomes normal polynomial coefficients, - * i.e. a[i] = coefficients[i]. - */ - private final double a[]; - - /** - * Whether the polynomial coefficients are available. - */ - private boolean coefficientsComputed; - - /** - * Construct a Newton polynomial with the given a[] and c[]. The order of - * centers are important in that if c[] shuffle, then values of a[] would - * completely change, not just a permutation of old a[]. - * <p> - * The constructor makes copy of the input arrays and assigns them.</p> - * - * @param a Coefficients in Newton form formula. - * @param c Centers. - * @throws NullArgumentException if any argument is {@code null}. - * @throws NoDataException if any array has zero length. - * @throws DimensionMismatchException if the size difference between - * {@code a} and {@code c} is not equal to 1. - */ - public PolynomialFunctionNewtonForm(double a[], double c[]) - throws NullArgumentException, NoDataException, DimensionMismatchException { - - verifyInputArray(a, c); - this.a = new double[a.length]; - this.c = new double[c.length]; - System.arraycopy(a, 0, this.a, 0, a.length); - System.arraycopy(c, 0, this.c, 0, c.length); - coefficientsComputed = false; - } - - /** - * Calculate the function value at the given point. - * - * @param z Point at which the function value is to be computed. - * @return the function value. - */ - public double value(double z) { - return evaluate(a, c, z); - } - - /** - * {@inheritDoc} - * @since 3.1 - */ - public DerivativeStructure value(final DerivativeStructure t) { - verifyInputArray(a, c); - - final int n = c.length; - DerivativeStructure value = new DerivativeStructure(t.getFreeParameters(), t.getOrder(), a[n]); - for (int i = n - 1; i >= 0; i--) { - value = t.subtract(c[i]).multiply(value).add(a[i]); - } - - return value; - - } - - /** - * Returns the degree of the polynomial. - * - * @return the degree of the polynomial - */ - public int degree() { - return c.length; - } - - /** - * Returns a copy of coefficients in Newton form formula. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * - * @return a fresh copy of coefficients in Newton form formula - */ - public double[] getNewtonCoefficients() { - double[] out = new double[a.length]; - System.arraycopy(a, 0, out, 0, a.length); - return out; - } - - /** - * Returns a copy of the centers array. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * - * @return a fresh copy of the centers array. - */ - public double[] getCenters() { - double[] out = new double[c.length]; - System.arraycopy(c, 0, out, 0, c.length); - return out; - } - - /** - * Returns a copy of the coefficients array. - * <p> - * Changes made to the returned copy will not affect the polynomial.</p> - * - * @return a fresh copy of the coefficients array. - */ - public double[] getCoefficients() { - if (!coefficientsComputed) { - computeCoefficients(); - } - double[] out = new double[coefficients.length]; - System.arraycopy(coefficients, 0, out, 0, coefficients.length); - return out; - } - - /** - * Evaluate the Newton polynomial using nested multiplication. It is - * also called <a href="http://mathworld.wolfram.com/HornersRule.html";> - * Horner's Rule</a> and takes O(N) time. - * - * @param a Coefficients in Newton form formula. - * @param c Centers. - * @param z Point at which the function value is to be computed. - * @return the function value. - * @throws NullArgumentException if any argument is {@code null}. - * @throws NoDataException if any array has zero length. - * @throws DimensionMismatchException if the size difference between - * {@code a} and {@code c} is not equal to 1. - */ - public static double evaluate(double a[], double c[], double z) - throws NullArgumentException, DimensionMismatchException, NoDataException { - verifyInputArray(a, c); - - final int n = c.length; - double value = a[n]; - for (int i = n - 1; i >= 0; i--) { - value = a[i] + (z - c[i]) * value; - } - - return value; - } - - /** - * Calculate the normal polynomial coefficients given the Newton form. - * It also uses nested multiplication but takes O(N^2) time. - */ - protected void computeCoefficients() { - final int n = degree(); - - coefficients = new double[n+1]; - for (int i = 0; i <= n; i++) { - coefficients[i] = 0.0; - } - - coefficients[0] = a[n]; - for (int i = n-1; i >= 0; i--) { - for (int j = n-i; j > 0; j--) { - coefficients[j] = coefficients[j-1] - c[i] * coefficients[j]; - } - coefficients[0] = a[i] - c[i] * coefficients[0]; - } - - coefficientsComputed = true; - } - - /** - * Verifies that the input arrays are valid. - * <p> - * The centers must be distinct for interpolation purposes, but not - * for general use. Thus it is not verified here.</p> - * - * @param a the coefficients in Newton form formula - * @param c the centers - * @throws NullArgumentException if any argument is {@code null}. - * @throws NoDataException if any array has zero length. - * @throws DimensionMismatchException if the size difference between - * {@code a} and {@code c} is not equal to 1. - * @see org.apache.commons.math3.analysis.interpolation.DividedDifferenceInterpolator#computeDividedDifference(double[], - * double[]) - */ - protected static void verifyInputArray(double a[], double c[]) - throws NullArgumentException, NoDataException, DimensionMismatchException { - MathUtils.checkNotNull(a); - MathUtils.checkNotNull(c); - if (a.length == 0 || c.length == 0) { - throw new NoDataException(LocalizedFormats.EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY); - } - if (a.length != c.length + 1) { - throw new DimensionMismatchException(LocalizedFormats.ARRAY_SIZES_SHOULD_HAVE_DIFFERENCE_1, - a.length, c.length); - } - } - -}
