View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.interpolation;
23  
24  import java.util.Arrays;
25  
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.analysis.BivariateFunction;
28  import org.hipparchus.analysis.FieldBivariateFunction;
29  import org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction;
30  import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
31  import org.hipparchus.exception.LocalizedCoreFormats;
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.exception.NullArgumentException;
34  import org.hipparchus.util.MathArrays;
35  
36  /**
37   * Function that implements the
38   * <a href="http://www.paulinternet.nl/?page=bicubic">bicubic spline</a>
39   * interpolation.
40   * This implementation currently uses {@link AkimaSplineInterpolator} as the
41   * underlying one-dimensional interpolator, which requires 5 sample points;
42   * insufficient data will raise an exception when the
43   * {@link #value(double,double) value} method is called.
44   *
45   */
46  public class PiecewiseBicubicSplineInterpolatingFunction
47      implements BivariateFunction, FieldBivariateFunction {
48  
49      /** The minimum number of points that are needed to compute the function. */
50      private static final int MIN_NUM_POINTS = 5;
51      /** Samples x-coordinates */
52      private final double[] xval;
53      /** Samples y-coordinates */
54      private final double[] yval;
55      /** Set of cubic splines patching the whole data grid */
56      private final double[][] fval;
57  
58      /** Simple constructor.
59       * @param x Sample values of the x-coordinate, in increasing order.
60       * @param y Sample values of the y-coordinate, in increasing order.
61       * @param f Values of the function on every grid point. the expected number
62       *        of elements.
63       * @throws MathIllegalArgumentException if {@code x} or {@code y} are not
64       *         strictly increasing.
65       * @throws NullArgumentException if any of the arguments are null
66       * @throws MathIllegalArgumentException if any of the arrays has zero length.
67       * @throws MathIllegalArgumentException if the length of x and y don't match the row, column
68       *         height of f
69       */
70      public PiecewiseBicubicSplineInterpolatingFunction(double[] x,
71                                                         double[] y,
72                                                         double[][] f)
73          throws MathIllegalArgumentException, NullArgumentException {
74          if (x == null ||
75              y == null ||
76              f == null ||
77              f[0] == null) {
78              throw new NullArgumentException();
79          }
80  
81          final int xLen = x.length;
82          final int yLen = y.length;
83  
84          if (xLen == 0 ||
85              yLen == 0 ||
86              f.length == 0 ||
87              f[0].length == 0) {
88              throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
89          }
90  
91          if (xLen < MIN_NUM_POINTS ||
92              yLen < MIN_NUM_POINTS ||
93              f.length < MIN_NUM_POINTS ||
94              f[0].length < MIN_NUM_POINTS) {
95              throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DATA);
96          }
97  
98          if (xLen != f.length) {
99              throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
100                                                    xLen, f.length);
101         }
102 
103         if (yLen != f[0].length) {
104             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
105                                                    yLen, f[0].length);
106         }
107 
108         MathArrays.checkOrder(x);
109         MathArrays.checkOrder(y);
110 
111         xval = x.clone();
112         yval = y.clone();
113         fval = f.clone();
114     }
115 
116     /**
117      * {@inheritDoc}
118      */
119     @Override
120     public double value(double x,
121                         double y)
122         throws MathIllegalArgumentException {
123         final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
124         final int offset = 2;
125         final int count = offset + 3;
126         final int i = searchIndex(x, xval, offset, count);
127         final int j = searchIndex(y, yval, offset, count);
128 
129         final double xArray[] = new double[count];
130         final double yArray[] = new double[count];
131         final double zArray[] = new double[count];
132         final double interpArray[] = new double[count];
133 
134         for (int index = 0; index < count; index++) {
135             xArray[index] = xval[i + index];
136             yArray[index] = yval[j + index];
137         }
138 
139         for (int zIndex = 0; zIndex < count; zIndex++) {
140             for (int index = 0; index < count; index++) {
141                 zArray[index] = fval[i + index][j + zIndex];
142             }
143             final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
144             interpArray[zIndex] = spline.value(x);
145         }
146 
147         final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray);
148 
149         return spline.value(y);
150 
151     }
152 
153     /**
154      * {@inheritDoc}
155      * @since 1.5
156      */
157     @Override
158     public <T extends CalculusFieldElement<T>> T value(final T x, final T y)
159         throws MathIllegalArgumentException {
160         final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
161         final int offset = 2;
162         final int count = offset + 3;
163         final int i = searchIndex(x.getReal(), xval, offset, count);
164         final int j = searchIndex(y.getReal(), yval, offset, count);
165 
166         final double xArray[] = new double[count];
167         final T yArray[]      = MathArrays.buildArray(x.getField(), count);
168         final double zArray[] = new double[count];
169         final T interpArray[] = MathArrays.buildArray(x.getField(), count);
170 
171         final T zero = x.getField().getZero();
172         for (int index = 0; index < count; index++) {
173             xArray[index] = xval[i + index];
174             yArray[index] = zero.add(yval[j + index]);
175         }
176 
177         for (int zIndex = 0; zIndex < count; zIndex++) {
178             for (int index = 0; index < count; index++) {
179                 zArray[index] = fval[i + index][j + zIndex];
180             }
181             final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
182             interpArray[zIndex] = spline.value(x);
183         }
184 
185         final FieldPolynomialSplineFunction<T> spline = interpolator.interpolate(yArray, interpArray);
186 
187         return spline.value(y);
188 
189     }
190 
191     /**
192      * Indicates whether a point is within the interpolation range.
193      *
194      * @param x First coordinate.
195      * @param y Second coordinate.
196      * @return {@code true} if (x, y) is a valid point.
197      */
198     public boolean isValidPoint(double x,
199                                 double y) {
200         if (x < xval[0] ||
201             x > xval[xval.length - 1] ||
202             y < yval[0] ||
203             y > yval[yval.length - 1]) {
204             return false;
205         } else {
206             return true;
207         }
208     }
209 
210     /**
211      * @param c Coordinate.
212      * @param val Coordinate samples.
213      * @param offset how far back from found value to offset for querying
214      * @param count total number of elements forward from beginning that will be
215      *        queried
216      * @return the index in {@code val} corresponding to the interval containing
217      *         {@code c}.
218      * @throws MathIllegalArgumentException if {@code c} is out of the range defined by
219      *         the boundary values of {@code val}.
220      */
221     private int searchIndex(double c,
222                             double[] val,
223                             int offset,
224                             int count) {
225         int r = Arrays.binarySearch(val, c);
226 
227         if (r == -1 || r == -val.length - 1) {
228             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
229                                                    c, val[0], val[val.length - 1]);
230         }
231 
232         if (r < 0) {
233             // "c" in within an interpolation sub-interval, which returns
234             // negative
235             // need to remove the negative sign for consistency
236             r = -r - offset - 1;
237         } else {
238             r -= offset;
239         }
240 
241         if (r < 0) {
242             r = 0;
243         }
244 
245         if ((r + count) >= val.length) {
246             // "c" is the last sample of the range: Return the index
247             // of the sample at the lower end of the last sub-interval.
248             r = val.length - count;
249         }
250 
251         return r;
252     }
253 }