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 * @version $Id: SmoothingPolynomialBicubicSplineInterpolator.java 1455194 2013-03-11 15:45:54Z luc $
037 * @since 2.2
038 */
039public class SmoothingPolynomialBicubicSplineInterpolator
040    extends BicubicSplineInterpolator {
041    /** Fitter for x. */
042    private final PolynomialFitter xFitter;
043    /** Degree of the fitting polynomial. */
044    private final int xDegree;
045    /** Fitter for y. */
046    private final PolynomialFitter yFitter;
047    /** Degree of the fitting polynomial. */
048    private final int yDegree;
049
050    /**
051     * Default constructor. The degree of the fitting polynomials is set to 3.
052     */
053    public SmoothingPolynomialBicubicSplineInterpolator() {
054        this(3);
055    }
056
057    /**
058     * @param degree Degree of the polynomial fitting functions.
059     * @exception NotPositiveException if degree is not positive
060     */
061    public SmoothingPolynomialBicubicSplineInterpolator(int degree)
062        throws NotPositiveException {
063        this(degree, degree);
064    }
065
066    /**
067     * @param xDegree Degree of the polynomial fitting functions along the
068     * x-dimension.
069     * @param yDegree Degree of the polynomial fitting functions along the
070     * y-dimension.
071     * @exception NotPositiveException if degrees are not positive
072     */
073    public SmoothingPolynomialBicubicSplineInterpolator(int xDegree, int yDegree)
074        throws NotPositiveException {
075        if (xDegree < 0) {
076            throw new NotPositiveException(xDegree);
077        }
078        if (yDegree < 0) {
079            throw new NotPositiveException(yDegree);
080        }
081        this.xDegree = xDegree;
082        this.yDegree = yDegree;
083
084        final double safeFactor = 1e2;
085        final SimpleVectorValueChecker checker
086            = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON,
087                                           safeFactor * Precision.SAFE_MIN);
088        xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
089        yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
090    }
091
092    /**
093     * {@inheritDoc}
094     */
095    @Override
096    public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
097                                                          final double[] yval,
098                                                          final double[][] fval)
099        throws NoDataException, NullArgumentException,
100               DimensionMismatchException, NonMonotonicSequenceException {
101        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
102            throw new NoDataException();
103        }
104        if (xval.length != fval.length) {
105            throw new DimensionMismatchException(xval.length, fval.length);
106        }
107
108        final int xLen = xval.length;
109        final int yLen = yval.length;
110
111        for (int i = 0; i < xLen; i++) {
112            if (fval[i].length != yLen) {
113                throw new DimensionMismatchException(fval[i].length, yLen);
114            }
115        }
116
117        MathArrays.checkOrder(xval);
118        MathArrays.checkOrder(yval);
119
120        // For each line y[j] (0 <= j < yLen), construct a polynomial, with
121        // respect to variable x, fitting array fval[][j]
122        final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
123        for (int j = 0; j < yLen; j++) {
124            xFitter.clearObservations();
125            for (int i = 0; i < xLen; i++) {
126                xFitter.addObservedPoint(1, xval[i], fval[i][j]);
127            }
128
129            // Initial guess for the fit is zero for each coefficients (of which
130            // there are "xDegree" + 1).
131            yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1]));
132        }
133
134        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
135        // values fval_1
136        final double[][] fval_1 = new double[xLen][yLen];
137        for (int j = 0; j < yLen; j++) {
138            final PolynomialFunction f = yPolyX[j];
139            for (int i = 0; i < xLen; i++) {
140                fval_1[i][j] = f.value(xval[i]);
141            }
142        }
143
144        // For each line x[i] (0 <= i < xLen), construct a polynomial, with
145        // respect to variable y, fitting array fval_1[i][]
146        final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
147        for (int i = 0; i < xLen; i++) {
148            yFitter.clearObservations();
149            for (int j = 0; j < yLen; j++) {
150                yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
151            }
152
153            // Initial guess for the fit is zero for each coefficients (of which
154            // there are "yDegree" + 1).
155            xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1]));
156        }
157
158        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
159        // values fval_2
160        final double[][] fval_2 = new double[xLen][yLen];
161        for (int i = 0; i < xLen; i++) {
162            final PolynomialFunction f = xPolyY[i];
163            for (int j = 0; j < yLen; j++) {
164                fval_2[i][j] = f.value(yval[j]);
165            }
166        }
167
168        return super.interpolate(xval, yval, fval_2);
169    }
170}