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.math4.analysis.interpolation;
018
019import java.util.ArrayList;
020import java.util.List;
021
022import org.apache.commons.math4.FieldElement;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.MathArithmeticException;
025import org.apache.commons.math4.exception.NoDataException;
026import org.apache.commons.math4.exception.NullArgumentException;
027import org.apache.commons.math4.exception.ZeroException;
028import org.apache.commons.math4.exception.util.LocalizedFormats;
029import org.apache.commons.math4.util.MathArrays;
030import org.apache.commons.math4.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<>();
065        this.topDiagonal    = new ArrayList<>();
066        this.bottomDiagonal = new ArrayList<>();
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    @SafeVarargs
091    public final void addSamplePoint(final T x, final T[] ... value)
092        throws ZeroException, MathArithmeticException,
093               DimensionMismatchException, NullArgumentException {
094
095        MathUtils.checkNotNull(x);
096        T factorial = x.getField().getOne();
097        for (int i = 0; i < value.length; ++i) {
098
099            final T[] y = value[i].clone();
100            if (i > 1) {
101                factorial = factorial.multiply(i);
102                final T inv = factorial.reciprocal();
103                for (int j = 0; j < y.length; ++j) {
104                    y[j] = y[j].multiply(inv);
105                }
106            }
107
108            // update the bottom diagonal of the divided differences array
109            final int n = abscissae.size();
110            bottomDiagonal.add(n - i, y);
111            T[] bottom0 = y;
112            for (int j = i; j < n; ++j) {
113                final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
114                if (x.equals(abscissae.get(n - (j + 1)))) {
115                    throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
116                }
117                final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal();
118                for (int k = 0; k < y.length; ++k) {
119                    bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
120                }
121                bottom0 = bottom1;
122            }
123
124            // update the top diagonal of the divided differences array
125            topDiagonal.add(bottom0.clone());
126
127            // update the abscissae array
128            abscissae.add(x);
129
130        }
131
132    }
133
134    /** Interpolate value at a specified abscissa.
135     * @param x interpolation abscissa
136     * @return interpolated value
137     * @exception NoDataException if sample is empty
138     * @throws NullArgumentException if x is null
139     */
140    public T[] value(T x) throws NoDataException, NullArgumentException {
141
142        // safety check
143        MathUtils.checkNotNull(x);
144        if (abscissae.isEmpty()) {
145            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
146        }
147
148        final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
149        T valueCoeff = x.getField().getOne();
150        for (int i = 0; i < topDiagonal.size(); ++i) {
151            T[] dividedDifference = topDiagonal.get(i);
152            for (int k = 0; k < value.length; ++k) {
153                value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff));
154            }
155            final T deltaX = x.subtract(abscissae.get(i));
156            valueCoeff = valueCoeff.multiply(deltaX);
157        }
158
159        return value;
160
161    }
162
163    /** Interpolate value and first derivatives at a specified abscissa.
164     * @param x interpolation abscissa
165     * @param order maximum derivation order
166     * @return interpolated value and derivatives (value in row 0,
167     * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
168     * @exception NoDataException if sample is empty
169     * @throws NullArgumentException if x is null
170     */
171    public T[][] derivatives(T x, int order) throws NoDataException, NullArgumentException {
172
173        // safety check
174        MathUtils.checkNotNull(x);
175        if (abscissae.isEmpty()) {
176            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
177        }
178
179        final T zero = x.getField().getZero();
180        final T one  = x.getField().getOne();
181        final T[] tj = MathArrays.buildArray(x.getField(), order + 1);
182        tj[0] = zero;
183        for (int i = 0; i < order; ++i) {
184            tj[i + 1] = tj[i].add(one);
185        }
186
187        final T[][] derivatives =
188                MathArrays.buildArray(x.getField(), order + 1, topDiagonal.get(0).length);
189        final T[] valueCoeff = MathArrays.buildArray(x.getField(), order + 1);
190        valueCoeff[0] = x.getField().getOne();
191        for (int i = 0; i < topDiagonal.size(); ++i) {
192            T[] dividedDifference = topDiagonal.get(i);
193            final T deltaX = x.subtract(abscissae.get(i));
194            for (int j = order; j >= 0; --j) {
195                for (int k = 0; k < derivatives[j].length; ++k) {
196                    derivatives[j][k] =
197                            derivatives[j][k].add(dividedDifference[k].multiply(valueCoeff[j]));
198                }
199                valueCoeff[j] = valueCoeff[j].multiply(deltaX);
200                if (j > 0) {
201                    valueCoeff[j] = valueCoeff[j].add(tj[j].multiply(valueCoeff[j - 1]));
202                }
203            }
204        }
205
206        return derivatives;
207
208    }
209
210}