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 java.util.ArrayList;
020import java.util.List;
021
022import org.apache.commons.math3.FieldElement;
023import org.apache.commons.math3.exception.DimensionMismatchException;
024import org.apache.commons.math3.exception.MathArithmeticException;
025import org.apache.commons.math3.exception.NoDataException;
026import org.apache.commons.math3.exception.NullArgumentException;
027import org.apache.commons.math3.exception.ZeroException;
028import org.apache.commons.math3.exception.util.LocalizedFormats;
029import org.apache.commons.math3.util.MathArrays;
030import org.apache.commons.math3.util.MathUtils;
031
032/** Polynomial interpolator using both sample values and sample derivatives.
033 * <p>
034 * The interpolation polynomials match all sample points, including both values
035 * and provided derivatives. There is one polynomial for each component of
036 * the values vector. All polynomials have the same degree. The degree of the
037 * polynomials depends on the number of points and number of derivatives at each
038 * point. For example the interpolation polynomials for n sample points without
039 * any derivatives all have degree n-1. The interpolation polynomials for n
040 * sample points with the two extreme points having value and first derivative
041 * and the remaining points having value only all have degree n+1. The
042 * interpolation polynomial for n sample points with value, first and second
043 * derivative for all points all have degree 3n-1.
044 * </p>
045 *
046 * @param <T> Type of the field elements.
047 *
048 * @since 3.2
049 */
050public class FieldHermiteInterpolator<T extends FieldElement<T>> {
051
052    /** Sample abscissae. */
053    private final List<T> abscissae;
054
055    /** Top diagonal of the divided differences array. */
056    private final List<T[]> topDiagonal;
057
058    /** Bottom diagonal of the divided differences array. */
059    private final List<T[]> bottomDiagonal;
060
061    /** Create an empty interpolator.
062     */
063    public FieldHermiteInterpolator() {
064        this.abscissae      = new ArrayList<T>();
065        this.topDiagonal    = new ArrayList<T[]>();
066        this.bottomDiagonal = new ArrayList<T[]>();
067    }
068
069    /** Add a sample point.
070     * <p>
071     * This method must be called once for each sample point. It is allowed to
072     * mix some calls with values only with calls with values and first
073     * derivatives.
074     * </p>
075     * <p>
076     * The point abscissae for all calls <em>must</em> be different.
077     * </p>
078     * @param x abscissa of the sample point
079     * @param value value and derivatives of the sample point
080     * (if only one row is passed, it is the value, if two rows are
081     * passed the first one is the value and the second the derivative
082     * and so on)
083     * @exception ZeroException if the abscissa difference between added point
084     * and a previous point is zero (i.e. the two points are at same abscissa)
085     * @exception MathArithmeticException if the number of derivatives is larger
086     * than 20, which prevents computation of a factorial
087     * @throws DimensionMismatchException if derivative structures are inconsistent
088     * @throws NullArgumentException if x is null
089     */
090    public void addSamplePoint(final T x, final T[] ... value)
091        throws ZeroException, MathArithmeticException,
092               DimensionMismatchException, NullArgumentException {
093
094        MathUtils.checkNotNull(x);
095        T factorial = x.getField().getOne();
096        for (int i = 0; i < value.length; ++i) {
097
098            final T[] y = value[i].clone();
099            if (i > 1) {
100                factorial = factorial.multiply(i);
101                final T inv = factorial.reciprocal();
102                for (int j = 0; j < y.length; ++j) {
103                    y[j] = y[j].multiply(inv);
104                }
105            }
106
107            // update the bottom diagonal of the divided differences array
108            final int n = abscissae.size();
109            bottomDiagonal.add(n - i, y);
110            T[] bottom0 = y;
111            for (int j = i; j < n; ++j) {
112                final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
113                if (x.equals(abscissae.get(n - (j + 1)))) {
114                    throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
115                }
116                final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal();
117                for (int k = 0; k < y.length; ++k) {
118                    bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
119                }
120                bottom0 = bottom1;
121            }
122
123            // update the top diagonal of the divided differences array
124            topDiagonal.add(bottom0.clone());
125
126            // update the abscissae array
127            abscissae.add(x);
128
129        }
130
131    }
132
133    /** Interpolate value at a specified abscissa.
134     * @param x interpolation abscissa
135     * @return interpolated value
136     * @exception NoDataException if sample is empty
137     * @throws NullArgumentException if x is null
138     */
139    public T[] value(T x) throws NoDataException, NullArgumentException {
140
141        // safety check
142        MathUtils.checkNotNull(x);
143        if (abscissae.isEmpty()) {
144            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
145        }
146
147        final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
148        T valueCoeff = x.getField().getOne();
149        for (int i = 0; i < topDiagonal.size(); ++i) {
150            T[] dividedDifference = topDiagonal.get(i);
151            for (int k = 0; k < value.length; ++k) {
152                value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff));
153            }
154            final T deltaX = x.subtract(abscissae.get(i));
155            valueCoeff = valueCoeff.multiply(deltaX);
156        }
157
158        return value;
159
160    }
161
162    /** Interpolate value and first derivatives at a specified abscissa.
163     * @param x interpolation abscissa
164     * @param order maximum derivation order
165     * @return interpolated value and derivatives (value in row 0,
166     * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
167     * @exception NoDataException if sample is empty
168     * @throws NullArgumentException if x is null
169     */
170    public T[][] derivatives(T x, int order) throws NoDataException, NullArgumentException {
171
172        // safety check
173        MathUtils.checkNotNull(x);
174        if (abscissae.isEmpty()) {
175            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
176        }
177
178        final T zero = x.getField().getZero();
179        final T one  = x.getField().getOne();
180        final T[] tj = MathArrays.buildArray(x.getField(), order + 1);
181        tj[0] = zero;
182        for (int i = 0; i < order; ++i) {
183            tj[i + 1] = tj[i].add(one);
184        }
185
186        final T[][] derivatives =
187                MathArrays.buildArray(x.getField(), order + 1, topDiagonal.get(0).length);
188        final T[] valueCoeff = MathArrays.buildArray(x.getField(), order + 1);
189        valueCoeff[0] = x.getField().getOne();
190        for (int i = 0; i < topDiagonal.size(); ++i) {
191            T[] dividedDifference = topDiagonal.get(i);
192            final T deltaX = x.subtract(abscissae.get(i));
193            for (int j = order; j >= 0; --j) {
194                for (int k = 0; k < derivatives[j].length; ++k) {
195                    derivatives[j][k] =
196                            derivatives[j][k].add(dividedDifference[k].multiply(valueCoeff[j]));
197                }
198                valueCoeff[j] = valueCoeff[j].multiply(deltaX);
199                if (j > 0) {
200                    valueCoeff[j] = valueCoeff[j].add(tj[j].multiply(valueCoeff[j - 1]));
201                }
202            }
203        }
204
205        return derivatives;
206
207    }
208
209}