View Javadoc
1   /*
2    * Licensed to the Hipparchus project 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 Hipparchus project 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  package org.hipparchus.analysis.interpolation;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.exception.LocalizedCoreFormats;
21  import org.hipparchus.exception.MathIllegalArgumentException;
22  import org.hipparchus.util.FastMath;
23  import org.hipparchus.util.MathArrays;
24  
25  import java.util.concurrent.atomic.AtomicInteger;
26  
27  /**
28   * Helper for finding interpolation nodes along one axis of grid data.
29   * <p>
30   * This class is intended to be used for interpolating inside grids.
31   * It works on any sorted data without duplication and size at least
32   * {@code n} where {@code n} is the number of points required for
33   * interpolation (i.e. 2 for linear interpolation, 3 for quadratic...)
34   * </p>
35   * <p>
36   * The method uses linear interpolation to select the nodes indices.
37   * It should be O(1) for sufficiently regular data, therefore much faster
38   * than bisection. It also features caching, which improves speed when
39   * interpolating several points in row in the close locations, i.e. when
40   * successive calls have a high probability to return the same interpolation
41   * nodes. This occurs for example when scanning with small steps a loose
42   * grid. The method also works on non-regular grids, but may be slower in
43   * this case.
44   * </p>
45   * <p>
46   * This class is thread-safe.
47   * </p>
48   * @param <T> Type of the field elements.
49   * @since 4.0
50   */
51  public class FieldGridAxis<T extends CalculusFieldElement<T>> {
52  
53      /** All the coordinates of the interpolation points, sorted in increasing order. */
54      private final T[] grid;
55  
56      /** Number of points required for interpolation. */
57      private final int n;
58  
59      /** Cached value of last x index. */
60      private final AtomicInteger cache;
61  
62      /** Simple constructor.
63       * @param grid coordinates of the interpolation points, sorted in increasing order
64       * @param n number of points required for interpolation, i.e. 2 for linear, 3
65       * for quadratic...
66       * @exception MathIllegalArgumentException if grid size is smaller than {@code n}
67       * or if the grid is not sorted in strict increasing order
68       */
69      public FieldGridAxis(final T[] grid, final int n)
70          throws MathIllegalArgumentException {
71  
72          // safety checks
73          if (grid.length < n) {
74              throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DIMENSION,
75                                                     grid.length, n);
76          }
77          MathArrays.checkOrder(grid);
78  
79          this.grid  = grid.clone();
80          this.n     = n;
81          this.cache = new AtomicInteger(0);
82  
83      }
84  
85      /** Get the number of points of the grid.
86       * @return number of points of the grid
87       */
88      public int size() {
89          return grid.length;
90      }
91  
92      /** Get the number of points required for interpolation.
93       * @return number of points required for interpolation
94       */
95      public int getN() {
96          return n;
97      }
98  
99      /** Get the interpolation node at specified index.
100      * @param index node index
101      * @return coordinate of the node at specified index
102      */
103     public T node(final int index) {
104         return grid[index];
105     }
106 
107     /** Get the index of the first interpolation node for some coordinate along the grid.
108      * <p>
109      * The index return is the one for the lowest interpolation node suitable for
110      * {@code t}. This means that if {@code i} is returned the nodes to use for
111      * interpolation at coordinate {@code t} are at indices {@code i}, {@code i+1},
112      * ..., {@code i+n-1}, where {@code n} is the number of points required for
113      * interpolation passed at construction.
114      * </p>
115      * <p>
116      * The index is selected in order to have the subset of nodes from {@code i} to
117      * {@code i+n-1} as balanced as possible around {@code t}:
118      * </p>
119      * <ul>
120      *   <li>
121      *     if {@code t} is inside the grid and sufficiently far from the endpoints
122      *     <ul>
123      *       <li>
124      *         if {@code n} is even, the returned nodes will be perfectly balanced:
125      *         there will be {@code n/2} nodes smaller than {@code t} and {@code n/2}
126      *         nodes larger than {@code t}
127      *       </li>
128      *       <li>
129      *         if {@code n} is odd, the returned nodes will be slightly unbalanced by
130      *         one point: there will be {@code (n+1)/2} nodes smaller than {@code t}
131      *         and {@code (n-1)/2} nodes larger than {@code t}
132      *       </li>
133      *     </ul>
134      *   </li>
135      *   <li>
136      *     if {@code t} is inside the grid and close to endpoints, the returned nodes
137      *     will be unbalanced: there will be less nodes on the endpoints side and
138      *     more nodes on the interior side
139      *   </li>
140      *   <li>
141      *     if {@code t} is outside of the grid, the returned nodes will completely
142      *     off balance: all nodes will be on the same side with respect to {@code t}
143      *   </li>
144      * </ul>
145      * <p>
146      * It is <em>not</em> an error to call this method with {@code t} outside of the grid,
147      * it simply implies that the interpolation will become an extrapolation and accuracy
148      * will decrease as {@code t} goes farther from the grid points. This is intended so
149      * interpolation does not fail near the end of the grid.
150      * </p>
151      * @param t coordinate of the point to interpolate
152      * @return index {@code i} such {@link #node(int) node(i)}, {@link #node(int) node(i+1)},
153      * ... {@link #node(int) node(i+n-1)} can be used for interpolating a value at
154      * coordinate {@code t}
155      * @since 1.4
156      */
157     public int interpolationIndex(final T t) {
158 
159         final int middleOffset = (n - 1) / 2;
160         int iInf = middleOffset;
161         int iSup = grid.length - (n - 1) + middleOffset;
162 
163         // first try to simply reuse the cached index,
164         // for faster return in a common case
165         final int    cached = cache.get();
166         final int    middle = cached + middleOffset;
167         final T aMid0  = grid[middle];
168         final T aMid1  = grid[middle + 1];
169         if (t.getReal() < aMid0.getReal()) {
170             if (middle == iInf) {
171                 // we are in the unbalanced low area
172                 return cached;
173             }
174         } else if (t.getReal() < aMid1.getReal()) {
175             // we are in the balanced middle area
176             return cached;
177         } else {
178             if (middle == iSup - 1) {
179                 // we are in the unbalanced high area
180                 return cached;
181             }
182         }
183 
184         // we need to find a new index
185         T aInf = grid[iInf];
186         T aSup = grid[iSup];
187         while (iSup - iInf > 1) {
188             final int iInterp = (int) ((iInf * (aSup.getReal() - t.getReal()) + iSup * (t.getReal() - aInf.getReal())) / (aSup.getReal() - aInf.getReal()));
189             final int iMed    = FastMath.max(iInf + 1, FastMath.min(iInterp, iSup - 1));
190             if (t.getReal() < grid[iMed].getReal()) {
191                 // keeps looking in the lower part of the grid
192                 iSup = iMed;
193                 aSup = grid[iSup];
194             } else {
195                 // keeps looking in the upper part of the grid
196                 iInf = iMed;
197                 aInf = grid[iInf];
198             }
199         }
200 
201        final int newCached = iInf - middleOffset;
202        cache.compareAndSet(cached, newCached);
203        return newCached;
204 
205     }
206 
207 }