001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import org.apache.commons.math3.exception.DimensionMismatchException;
020import org.apache.commons.math3.exception.NoDataException;
021import org.apache.commons.math3.exception.NonMonotonicSequenceException;
022import org.apache.commons.math3.exception.NotPositiveException;
023import org.apache.commons.math3.exception.NullArgumentException;
024import org.apache.commons.math3.util.MathArrays;
025import org.apache.commons.math3.util.Precision;
026import org.apache.commons.math3.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer;
027import org.apache.commons.math3.fitting.PolynomialFitter;
028import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
029import org.apache.commons.math3.optim.SimpleVectorValueChecker;
030
031/**
032 * Generates a bicubic interpolation function.
033 * Prior to generating the interpolating function, the input is smoothed using
034 * polynomial fitting.
035 *
036 * @since 2.2
037 * @deprecated To be removed in 4.0 (see MATH-1166).
038 */
039@Deprecated
040public class SmoothingPolynomialBicubicSplineInterpolator
041    extends BicubicSplineInterpolator {
042    /** Fitter for x. */
043    private final PolynomialFitter xFitter;
044    /** Degree of the fitting polynomial. */
045    private final int xDegree;
046    /** Fitter for y. */
047    private final PolynomialFitter yFitter;
048    /** Degree of the fitting polynomial. */
049    private final int yDegree;
050
051    /**
052     * Default constructor. The degree of the fitting polynomials is set to 3.
053     */
054    public SmoothingPolynomialBicubicSplineInterpolator() {
055        this(3);
056    }
057
058    /**
059     * @param degree Degree of the polynomial fitting functions.
060     * @exception NotPositiveException if degree is not positive
061     */
062    public SmoothingPolynomialBicubicSplineInterpolator(int degree)
063        throws NotPositiveException {
064        this(degree, degree);
065    }
066
067    /**
068     * @param xDegree Degree of the polynomial fitting functions along the
069     * x-dimension.
070     * @param yDegree Degree of the polynomial fitting functions along the
071     * y-dimension.
072     * @exception NotPositiveException if degrees are not positive
073     */
074    public SmoothingPolynomialBicubicSplineInterpolator(int xDegree, int yDegree)
075        throws NotPositiveException {
076        if (xDegree < 0) {
077            throw new NotPositiveException(xDegree);
078        }
079        if (yDegree < 0) {
080            throw new NotPositiveException(yDegree);
081        }
082        this.xDegree = xDegree;
083        this.yDegree = yDegree;
084
085        final double safeFactor = 1e2;
086        final SimpleVectorValueChecker checker
087            = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON,
088                                           safeFactor * Precision.SAFE_MIN);
089        xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
090        yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
091    }
092
093    /**
094     * {@inheritDoc}
095     */
096    @Override
097    public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
098                                                          final double[] yval,
099                                                          final double[][] fval)
100        throws NoDataException, NullArgumentException,
101               DimensionMismatchException, NonMonotonicSequenceException {
102        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
103            throw new NoDataException();
104        }
105        if (xval.length != fval.length) {
106            throw new DimensionMismatchException(xval.length, fval.length);
107        }
108
109        final int xLen = xval.length;
110        final int yLen = yval.length;
111
112        for (int i = 0; i < xLen; i++) {
113            if (fval[i].length != yLen) {
114                throw new DimensionMismatchException(fval[i].length, yLen);
115            }
116        }
117
118        MathArrays.checkOrder(xval);
119        MathArrays.checkOrder(yval);
120
121        // For each line y[j] (0 <= j < yLen), construct a polynomial, with
122        // respect to variable x, fitting array fval[][j]
123        final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
124        for (int j = 0; j < yLen; j++) {
125            xFitter.clearObservations();
126            for (int i = 0; i < xLen; i++) {
127                xFitter.addObservedPoint(1, xval[i], fval[i][j]);
128            }
129
130            // Initial guess for the fit is zero for each coefficients (of which
131            // there are "xDegree" + 1).
132            yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1]));
133        }
134
135        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
136        // values fval_1
137        final double[][] fval_1 = new double[xLen][yLen];
138        for (int j = 0; j < yLen; j++) {
139            final PolynomialFunction f = yPolyX[j];
140            for (int i = 0; i < xLen; i++) {
141                fval_1[i][j] = f.value(xval[i]);
142            }
143        }
144
145        // For each line x[i] (0 <= i < xLen), construct a polynomial, with
146        // respect to variable y, fitting array fval_1[i][]
147        final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
148        for (int i = 0; i < xLen; i++) {
149            yFitter.clearObservations();
150            for (int j = 0; j < yLen; j++) {
151                yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
152            }
153
154            // Initial guess for the fit is zero for each coefficients (of which
155            // there are "yDegree" + 1).
156            xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1]));
157        }
158
159        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
160        // values fval_2
161        final double[][] fval_2 = new double[xLen][yLen];
162        for (int i = 0; i < xLen; i++) {
163            final PolynomialFunction f = xPolyY[i];
164            for (int j = 0; j < yLen; j++) {
165                fval_2[i][j] = f.value(yval[j]);
166            }
167        }
168
169        return super.interpolate(xval, yval, fval_2);
170    }
171}