GridAxis.java

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      https://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.hipparchus.analysis.interpolation;

import java.io.Serializable;
import java.util.concurrent.atomic.AtomicInteger;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;

/**
 * Helper for finding interpolation nodes along one axis of grid data.
 * <p>
 * This class is intended to be used for interpolating inside grids.
 * It works on any sorted data without duplication and size at least
 * {@code n} where {@code n} is the number of points required for
 * interpolation (i.e. 2 for linear interpolation, 3 for quadratic...)
 * </p>
 * <p>
 * The method uses linear interpolation to select the nodes indices.
 * It should be O(1) for sufficiently regular data, therefore much faster
 * than bisection. It also features caching, which improves speed when
 * interpolating several points in row in the close locations, i.e. when
 * successive calls have a high probability to return the same interpolation
 * nodes. This occurs for example when scanning with small steps a loose
 * grid. The method also works on non-regular grids, but may be slower in
 * this case.
 * </p>
 * <p>
 * This class is thread-safe.
 * </p>
 * @since 1.4
 */
public class GridAxis implements Serializable {

    /** Serializable UID. */
    private static final long serialVersionUID = 20180926L;

    /** All the coordinates of the interpolation points, sorted in increasing order. */
    private final double[] grid;

    /** Number of points required for interpolation. */
    private final int n;

    /** Cached value of last x index. */
    private final AtomicInteger cache;

    /** Simple constructor.
     * @param grid coordinates of the interpolation points, sorted in increasing order
     * @param n number of points required for interpolation, i.e. 2 for linear, 3
     * for quadratic...
     * @exception MathIllegalArgumentException if grid size is smaller than {@code n}
     * or if the grid is not sorted in strict increasing order
     */
    public GridAxis(final double[] grid, final int n)
        throws MathIllegalArgumentException {

        // safety checks
        if (grid.length < n) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DIMENSION,
                                                   grid.length, n);
        }
        MathArrays.checkOrder(grid);

        this.grid  = grid.clone();
        this.n     = n;
        this.cache = new AtomicInteger(0);

    }

    /** Get the number of points of the grid.
     * @return number of points of the grid
     */
    public int size() {
        return grid.length;
    }

    /** Get the number of points required for interpolation.
     * @return number of points required for interpolation
     */
    public int getN() {
        return n;
    }

    /** Get the interpolation node at specified index.
     * @param index node index
     * @return coordinate of the node at specified index
     */
    public double node(final int index) {
        return grid[index];
    }

    /** Get the index of the first interpolation node for some coordinate along the grid.
     * <p>
     * The index return is the one for the lowest interpolation node suitable for
     * {@code t}. This means that if {@code i} is returned the nodes to use for
     * interpolation at coordinate {@code t} are at indices {@code i}, {@code i+1},
     * ..., {@code i+n-1}, where {@code n} is the number of points required for
     * interpolation passed at construction.
     * </p>
     * <p>
     * The index is selected in order to have the subset of nodes from {@code i} to
     * {@code i+n-1} as balanced as possible around {@code t}:
     * </p>
     * <ul>
     *   <li>
     *     if {@code t} is inside the grid and sufficiently far from the endpoints
     *     <ul>
     *       <li>
     *         if {@code n} is even, the returned nodes will be perfectly balanced:
     *         there will be {@code n/2} nodes smaller than {@code t} and {@code n/2}
     *         nodes larger than {@code t}
     *       </li>
     *       <li>
     *         if {@code n} is odd, the returned nodes will be slightly unbalanced by
     *         one point: there will be {@code (n+1)/2} nodes smaller than {@code t}
     *         and {@code (n-1)/2} nodes larger than {@code t}
     *       </li>
     *     </ul>
     *   </li>
     *   <li>
     *     if {@code t} is inside the grid and close to endpoints, the returned nodes
     *     will be unbalanced: there will be less nodes on the endpoints side and
     *     more nodes on the interior side
     *   </li>
     *   <li>
     *     if {@code t} is outside of the grid, the returned nodes will completely
     *     off balance: all nodes will be on the same side with respect to {@code t}
     *   </li>
     * </ul>
     * <p>
     * It is <em>not</em> an error to call this method with {@code t} outside of the grid,
     * it simply implies that the interpolation will become an extrapolation and accuracy
     * will decrease as {@code t} goes farther from the grid points. This is intended so
     * interpolation does not fail near the end of the grid.
     * </p>
     * @param t coordinate of the point to interpolate
     * @return index {@code i} such {@link #node(int) node(i)}, {@link #node(int) node(i+1)},
     * ... {@link #node(int) node(i+n-1)} can be used for interpolating a value at
     * coordinate {@code t}
     * @since 1.4
     */
    public int interpolationIndex(final double t) {

        final int middleOffset = (n - 1) / 2;
        int iInf = middleOffset;
        int iSup = grid.length - (n - 1) + middleOffset;

        // first try to simply reuse the cached index,
        // for faster return in a common case
        final int    cached = cache.get();
        final int    middle = cached + middleOffset;
        final double aMid0  = grid[middle];
        final double aMid1  = grid[middle + 1];
        if (t < aMid0) {
            if (middle == iInf) {
                // we are in the unbalanced low area
                return cached;
            }
        } else if (t < aMid1) {
            // we are in the balanced middle area
            return cached;
        } else {
            if (middle == iSup - 1) {
                // we are in the unbalanced high area
                return cached;
            }
        }

        // we need to find a new index
        double aInf = grid[iInf];
        double aSup = grid[iSup];
        while (iSup - iInf > 1) {
            final int iInterp = (int) ((iInf * (aSup - t) + iSup * (t - aInf)) / (aSup - aInf));
            final int iMed    = FastMath.max(iInf + 1, FastMath.min(iInterp, iSup - 1));
            if (t < grid[iMed]) {
                // keeps looking in the lower part of the grid
                iSup = iMed;
                aSup = grid[iSup];
            } else {
                // keeps looking in the upper part of the grid
                iInf = iMed;
                aInf = grid[iInf];
            }
        }

       final int newCached = iInf - middleOffset;
       cache.compareAndSet(cached, newCached);
       return newCached;

    }

}