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 }