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