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