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 }