com.meliorbis.numerics.generic.primitives.impl.Interpolation Maven / Gradle / Ivy
Show all versions of Numerics Show documentation
/**
*
*/
package com.meliorbis.numerics.generic.primitives.impl;
import static java.lang.Math.signum;
import java.util.Arrays;
import java.util.NoSuchElementException;
import java.util.PrimitiveIterator;
import java.util.logging.Logger;
import org.apache.commons.lang.builder.ToStringBuilder;
import org.apache.commons.math3.util.Precision;
import com.meliorbis.numerics.DoubleNumerics;
import com.meliorbis.numerics.generic.MultiDimensionalArrayException;
import com.meliorbis.numerics.generic.primitives.DoubleArray;
import com.meliorbis.numerics.generic.primitives.DoubleIndexedReduction;
import com.meliorbis.numerics.generic.primitives.DoubleNaryOp;
import com.meliorbis.numerics.generic.primitives.DoubleSettableIndexedIterator;
import com.meliorbis.numerics.generic.primitives.DoubleSubspaceSplitIterator;
import com.meliorbis.numerics.index.impl.Index;
import com.meliorbis.numerics.index.impl.Index.SubSpaceSplitIterator;
import com.meliorbis.utils.Utils;
/**
* Static methods to peform various types of interpolation on IDoubleArray
* classes.
*
* Presently, all interpolation is linear. The focus of these methods is
* primarily on treating the arrays as if they were domains and ranges of
* mathematical functions approximated over a grid of numbers (the domain).
* Interpolating to a different set of points along a dimension of the function
* can then be done efficiently by iterating along the domain and filling in
* points as they are encountered, rather than searching for each target point
* across the whole domain
*
* @author Tobias Grasl
*
* @see DoubleArray
*/
public abstract class Interpolation
{
private static final Logger LOG = Logger.getLogger(Interpolation.class.getName());
/* private default constructor prevents inheritance */
private Interpolation()
{
}
private static DoubleBlockedArrayFactory _arrayFactory;
static
{
_arrayFactory = DoubleNumerics.instance().getArrayFactory();
}
/**
* Used when interpolating an array, this value class determines, for a
* given dimension, what the source values of the interpolation along that
* dimension of the array are, and what the target value is.
*/
public static final class Specification implements Cloneable
{
private final int dimension;
public DoubleArray gridVals;
public double target;
/**
* @param dimension_ The dimensions this specification applies to
* @param gridVals_ The grid values at the points along that dimension
* @param targetVal_ The target value
*/
public Specification(int dimension_, DoubleArray gridVals_, double targetVal_)
{
dimension = dimension_;
gridVals = gridVals_;
target = targetVal_;
}
/**
* Shallow (!) Copy constructor
*
* @param src_ The specification to make a copy of
*/
public Specification(Specification src_)
{
// Shallow copy!
dimension = src_.dimension;
gridVals = src_.gridVals;
target = src_.target;
}
@Override
public String toString()
{
return new ToStringBuilder(this).append("dim", dimension).append("target", target).append("axis", gridVals).toString();
}
@Override
protected Object clone() throws CloneNotSupportedException
{
return new Specification(this);
}
}
/**
* A callback function which can be used during interpolation to be notified
* as each target value is found
*/
@FunctionalInterface
public static interface Callback
{
void callback(double interpResult_, int[] index_, int lowerTargetIndex_, double lowerTargetProportion_);
/**
* No-op callback
*/
public static Callback NOOP = (double interpResult_, int[] index_, int lowerTargetIndex_, double lowerTargetProportion_) -> {
};
}
@FunctionalInterface
public static interface Discontinuities
{
public PrimitiveIterator.OfDouble getDiscontinuitiesAt(int[] index_);
public Discontinuities NONE = index -> new PrimitiveIterator.OfDouble()
{
@Override
public boolean hasNext()
{
return false;
}
@Override
public double nextDouble()
{
throw new NoSuchElementException("The empty iterator has no elements");
}
};
}
/**
* Parameter class for passing the flags and extra information for
* interpolation to the interpolation functions, NOT the actual data
*/
public static class Params
{
/**
* Indicates whether the function being interpolated is constrained.
* TODO: Explain more
*
* Defaults to {@code false}.
*/
boolean _constrained = false;
/**
* Provides the discontinuities
*/
private Discontinuities _discontinuities = Discontinuities.NONE;
/**
* Only relevant when discontinuities are set, this flag indicates
* whether the function is assumed to be monotonous and this
* monotonicity should be preserved when target values are extrapolated
* around discontinuities.
*
* Defaults to {@code false}.
*/
private boolean _preserveMonotonicity = false;
/**
* The callback to notify with interpolated values.
*
* Defaults to {@link Callback#NOOP}
*/
private Callback _callback = Callback.NOOP;
public Params constrained()
{
_constrained = true;
return this;
}
public Params constrained(boolean flag_)
{
_constrained = flag_;
return this;
}
public Params withDiscontinuities(Discontinuities discontinuities_)
{
if(discontinuities_ == null) {
LOG.warning("Setting discontinuties to null is discouraged. The default is no discontinuities, or pass Discontinuities.NONE");
discontinuities_ = Discontinuities.NONE;
}
_discontinuities = discontinuities_;
return this;
}
public Params preserveMonotonicity()
{
_preserveMonotonicity = true;
return this;
}
public Params preserveMonotonicity(boolean value_)
{
_preserveMonotonicity = value_;
return this;
}
public Params withCallback(Callback callback_)
{
if(callback_ == null) {
LOG.warning("Setting callback to null is discouraged. The default is no callback, or pass Callback.NOOP");
callback_ = Callback.NOOP;
}
_callback = callback_;
return this;
}
private static Params DEFAULT = new Params();
}
/**
* @return A Params instance with default values.
*/
public static Params params() {
return new Params();
}
/**
* Creates a Specification detailing the interpolation settings in a single dimension
*
* @param dimension_ The dimension this specification applies to
* @param gridVals_ The grid values along that dimension
* @param target_ The interpolation target in that dimension
*
* @return A {@link Specification} object with the provided values set
*/
public static Specification spec(int dimension_, DoubleArray gridVals_, double target_)
{
return new Specification(dimension_, gridVals_, target_);
}
/**
* Interpolates a grid across one of more dimensions, where the axis values
* along the dimension and the target value are provided
*
* @param Y_
* The grid to interpolate
* @param dimSpecs_
* The specification of dimension index, axis and target values
* for each dimension to interpolate along
* @return An Array of the same size as the non-interpolated dimensions
* @throws com.meliorbis.numerics.generic.MultiDimensionalArrayException
* In the event of error
*/
public static DoubleArray interp(DoubleArray Y_, Specification... dimSpecs_)
{
// No interpolation specified (short-circuit for dynamic callers)
if (dimSpecs_.length == 0)
{
return Y_;
}
int[] upperIndices = new int[dimSpecs_.length];
final double[] upperProportions = new double[dimSpecs_.length];
// Count the dimensions that don't need to be interpolated because they
// are constant, or
// the interpolation exactly hits a grid point
int noInterpDimensions = 0;
for (int dimIndex = dimSpecs_.length - 1; dimIndex >= 0; dimIndex--)
{
Specification currentSpec = dimSpecs_[dimIndex];
if (currentSpec.gridVals.size()[0] == 1)
{
upperIndices[dimIndex] = 0;
upperProportions[dimIndex] = 1d;
noInterpDimensions++;
continue;
}
int upperIndex = 1;
double lower = currentSpec.gridVals.get(0);
double upper = currentSpec.gridVals.get(1);
while (upper < currentSpec.target && upperIndex < (currentSpec.gridVals.numberOfElements() - 1))
{
upperIndex++;
lower = upper;
upper = currentSpec.gridVals.get(upperIndex);
}
upperIndices[dimIndex] = upperIndex;
double upperProportion = (currentSpec.target - lower) / (upper - lower);
upperProportions[dimIndex] = upperProportion;
if (upperProportion == 1.0d || upperProportion == 0d)
{
noInterpDimensions++;
}
}
DoubleArray[] slices = new DoubleArray[(int) Math.pow(2, dimSpecs_.length - noInterpDimensions)];
final double[] proportions = new double[slices.length];
int[] index = Utils.repeatArray(-1, Y_.numberOfDimensions());
addBoundingPoints(Y_, dimSpecs_, upperIndices, upperProportions, slices, proportions, index);
return slices[0].with(Arrays.copyOfRange(slices, 1, slices.length)).map((DoubleNaryOp) inputs_ -> {
double sum = 0d;
for (int i = 0; i < inputs_.length; i++)
{
sum += inputs_[i] * proportions[i];
}
return sum;
});
}
/**
* Assumes that the arrays {@code sourceX_} and {@code sourceY_} are input
* and output values respectively of a function along dimension
* {@code dimension_}. The arrays must be of equal size, and the other
* dimensions are implicitly assumed to form a rectangular grid of other
* input variables, though they are irrelevant for this operation.
*
* {@code targetX_} must be either a one dimensional array or an array with
* the same number of dimensions as the sources and the same size in non-
* {@code dimension_} dimensions. In the former case the same x-values are
* interpolated too for each orthogonal point, whereas in the latter case
* they may vary.
*
* The returned array is the same size as {@code sourceX_} except in
* dimension {@code dimension_}, where it takes the size of {@code targetX_}
* . It represents the output values of a function of
* {@code targetX_.numberOfDimensions()} variables defined on a rectangular
* grid with the same input grid as {@code sourceX_} and the points
* {@code targetX_} in the interpolated dimension.
*
* Default {@link Params} are used.
*
* @param sourceX_
* The array of inputs, of arbitrary dimension. The values
* must be monotonous along {@code dimension_}
* @param sourceY_
* The array of function values, must be of the same dimension as
* {@code sourceX_}
* @param dimension_
* The dimension along which to interpolate, which must be
* smaller than the number of dimensions in {@code sourceX_}
* @param targetX_
* The x values to determine the function values for. This must
* be either a one dimensional array, or one with the same size
* as {@code sourceX_} except in dimension {@code dimension_}. It
* must be monotonous along {@code dimension_} with the same
* direction of monotonicity as sourceX_
*
* @return A value class containing three arrays: interpolated Y values,
* indexes of the point below each target X, and the proportion to
* assign to that lower point when distributing
*/
public static DoubleArray interpolateFunction(final DoubleArray sourceX_, final DoubleArray sourceY_, final int dimension_,
final DoubleArray targetX_)
{
return interpolateFunction(sourceX_, sourceY_, dimension_, targetX_, Params.DEFAULT);
}
/**
* Assumes that the arrays {@code sourceX_} and {@code sourceY_} are input
* and output values respectively of a function along dimension
* {@code dimension_}. The arrays must be of equal size, and the other
* dimensions are implicitly assumed to form a rectangular grid of other
* input variables, though they are irrelevant for this operation.
*
* {@code targetX_} must be either a one dimensional array or an array with
* the same number of dimensions as the sources and the same size in non-
* {@code dimension_} dimensions. In the former case the same x-values are
* interpolated too for each orthogonal point, whereas in the latter case
* they may vary.
*
* The returned array is the same size as {@code sourceX_} except in
* dimension {@code dimension_}, where it takes the size of {@code targetX_}
* . It represents the output values of a function of
* {@code targetX_.numberOfDimensions()} variables defined on a rectangular
* grid with the same input grid as {@code sourceX_} and the points
* {@code targetX_} in the interpolated dimension.
*
* @param sourceX_
* The array of inputs, of arbitrary dimension. The values
* must be monotonous along {@code dimension_}
* @param sourceY_
* The array of function values, must be of the same dimension as
* {@code sourceX_}
* @param dimension_
* The dimension along which to interpolate, which must be
* smaller than the number of dimensions in {@code sourceX_}
* @param targetX_
* The x values to determine the function values for. This must
* be either a one dimensional array, or one with the same size
* as {@code sourceX_} except in dimension {@code dimension_}. It
* must be monotonous along {@code dimension_} with the same
* direction of monotonicity as sourceX_
* @param params_
* Interpolation parameters
*
* @return A value class containing three arrays: interpolated Y values,
* indexes of the point below each target X, and the proportion to
* assign to that lower point when distributing
*/
public static DoubleArray interpolateFunction(final DoubleArray sourceX_, final DoubleArray sourceY_, final int dimension_,
final DoubleArray targetX_, Params params_)
{
DoubleArray targetY;
DoubleArray targetXExpanded;
// If targetX is one-D but sourceX is not, then by assumption a
// rectangular grid of values is called for
if (targetX_.numberOfDimensions() == 1 && sourceX_.numberOfDimensions() != 1)
{
int[] size = sourceX_.size().clone();
size[dimension_] = targetX_.numberOfElements();
targetY = _arrayFactory.newArray(size);
targetXExpanded = _arrayFactory.newArray(size);
targetXExpanded.fillDimensions(targetX_, dimension_);
} else
{
// TargetY should be of the correct size already
targetY = _arrayFactory.newArray(targetX_.size());
targetXExpanded = targetX_;
}
interpFunctionAndFill(sourceX_, sourceY_, dimension_, targetXExpanded, targetY, params_);
return targetY;
}
/**
* Iterates over the {@code across_} dimensions in {@code targetX_} and, for
* each one, interpolates the function defined by {@code sourceX_} and
* {@code sourceY_} to the orthogonal sub-array at that point exactly as
* described in
* {@link #interpolateFunction(DoubleArray, DoubleArray, int, DoubleArray)}
*
* Default {@link Params} are used.
*
* @param sourceX_
* The input values in the interpolation dimension. They
* need therefore not form a rectangular grid
* @param sourceY_
* The output values
* @param dimension_
* The dimension to interpolate along
* @param targetX_
* The x values to interpolate too
* @param across_
* The dimensions in targetX which are in excess of those in
* sourceX, and hence are iterated over and handled independently
*
* @return An array the same size as {@code targetX_} with the interpolated
* y value corresponding to each x value
*
* @see #interpolateFunction(DoubleArray, DoubleArray, int, DoubleArray,
* Params)
*/
public static DoubleArray interpolateFunctionAcross(final DoubleArray sourceX_, final DoubleArray sourceY_, int dimension_,
final DoubleArray targetX_, int... across_)
{
return interpolateFunctionAcross(sourceX_, sourceY_, dimension_, targetX_, Params.DEFAULT, across_);
}
/**
* Iterates over the {@code across_} dimensions in {@code targetX_} and, for
* each one, interpolates the function defined by {@code sourceX_} and
* {@code sourceY_} to the orthogonal sub-array at that point exactly as
* described in
* {@link #interpolateFunction(DoubleArray, DoubleArray, int, DoubleArray, Params)}
*
* @param sourceX_
* The input values in the interpolation dimension. They
* need therefore not form a rectangular grid
* @param sourceY_
* The output values
* @param dimension_
* The dimension to interpolate along
* @param targetX_
* The x values to interpolate too
* @param params_
* Additional parameters
* @param across_
* The dimensions in targetX which are in excess of those in
* sourceX, and hence are iterated over and handled independently
*
* @return An array the same size as {@code targetX_} with the interpolated
* y value corresponding to each x value
*
* @see #interpolateFunction(DoubleArray, DoubleArray, int, DoubleArray,
* Params)
*/
public static DoubleArray interpolateFunctionAcross(final DoubleArray sourceX_, final DoubleArray sourceY_, int dimension_,
final DoubleArray targetX_, Params params_, int... across_)
{
final DoubleArray targetY = _arrayFactory.newArray(targetX_.size());
/*
* Iterate across the 'across' dimension in the target X values
*/
Index index = new Index(targetX_.size());
SubSpaceSplitIterator acrossDimsIter = index.iterator(across_);
while (acrossDimsIter.hasNext())
{
acrossDimsIter.nextInt();
int[] currentFullIndex = acrossDimsIter.getCurrentFullIndex();
/*
* Perform the interpolation on each subspace at the across
* dimensions
*/
interpFunctionAndFill(sourceX_, sourceY_, dimension_, targetX_.at(currentFullIndex), targetY.at(currentFullIndex), params_);
}
return targetY;
}
/**
* Interpolates a set of function values along one dimension.
*
* Note that both xIterTgt_ and xIterSrc_ must be monotonous, the latter
* strictly so, and in the same direction.
*
* @param targetY_toFill_
* The area to be filled with the resulting y values
* @param tgtXIter_
* An iterator over the target x values
* @param srcYIter_
* An iterator over the source y values
* @param srcXIter_
* An iterator over the source x values
* @param constrained_
* Indicates that the function is constrained, so no
* extrapolation is performed prior to the first x
* @param callBack_
* Callback function to notify with each point found
*/
protected static void interpSingle(final DoubleArray targetY_toFill_, DoubleSubspaceSplitIterator tgtXIter_,
DoubleSettableIndexedIterator srcYIter_, DoubleSettableIndexedIterator srcXIter_, final boolean constrained_,
final Callback callBack_)
{
// Get the first two source x values
double lowerX = srcXIter_.nextDouble();
double upperX = srcXIter_.nextDouble();
// If the x's are decreasing, this will be true
boolean decreasing = lowerX > upperX;
// First two source Y values
double lowerY = srcYIter_.nextDouble();
double upperY = srcYIter_.nextDouble();
DoubleSettableIndexedIterator yIterTgt = targetY_toFill_.iterator();
/*
* The loop iterates over the target X values and finds the appropriate
* Y value for each
*/
while (tgtXIter_.hasNext())
{
double currentXTarget = tgtXIter_.nextDouble();
yIterTgt.nextDouble();
// Iterate over the x values until we find one that is
// bigger than the target
while (((upperX <= currentXTarget) ^ decreasing) && srcXIter_.hasNext())
{
lowerX = upperX;
upperX = srcXIter_.nextDouble();
lowerY = upperY;
upperY = srcYIter_.nextDouble();
}
double yAtTarget = interpValue(currentXTarget, lowerX, upperX, lowerY, upperY, decreasing, constrained_);
// Store the interpolated y value
yIterTgt.set(yAtTarget);
// Semantically, this check is unnecessary since the NULL_CALLBACK
// anyway does nothing, but the
// supporting bits are a performance drain
// TODO: Refactor callback API to put the supporting bits inside, or
// otherwise wrap with another callback
// to avoid this check
if (callBack_ != Callback.NOOP)
{
int[] currentTargetIndex = tgtXIter_.getFullIndex();
int currentUpperIndex = srcXIter_.getIndex()[0];
double lowerPropForPoint = (upperY - yAtTarget) / (upperY - lowerY);
callBack_.callback(yAtTarget, currentTargetIndex, currentUpperIndex - 1, lowerPropForPoint);
}
}
}
protected static void interpFunctionAndFill(final DoubleArray sourceX_, final DoubleArray sourceY_, final int dimension_,
final DoubleArray targetX_, final DoubleArray targetY_, Params params_)
{
targetX_.across(dimension_).reduce(new DoubleIndexedReduction()
{
@Override
public double perform(DoubleSubspaceSplitIterator iterator_)
{
DoubleSubspaceSplitIterator targetIterator = iterator_;
int[] fullIndex = targetIterator.getFullIndex();
DoubleSubspaceSplitIterator interpDimYIter = sourceY_.iteratorAt(fullIndex);
/*
* Get an iterator over the target values at this point in the
* non-interp dimension (the interp dimension will be -1)
*
* NOTE: The length of the target values along the interp
* dimension may be different than that of x and y, but it must
* be the same in other dimensions
*/
DoubleSettableIndexedIterator interpDimXIter = sourceX_.numberOfDimensions() == 1 ? sourceX_.iterator() : sourceX_
.iteratorAtArray(fullIndex);
if (params_._discontinuities == Discontinuities.NONE)
{
interpSingle(targetY_.at(fullIndex), targetIterator, interpDimYIter, interpDimXIter, params_._constrained, params_._callback);
} else
{
interpSingleDiscontinuous(targetY_.at(fullIndex), targetIterator, interpDimYIter, interpDimXIter,
params_._discontinuities.getDiscontinuitiesAt(iterator_.getOtherIndex()), params_._preserveMonotonicity,
params_._constrained, params_._callback);
}
return Double.NaN;
}
});
}
/**
* Interpolates a set of function values along one dimension, allowing for
* discontinuities.
*
* Points lying before (after) a discontinuity with no intervening values in
* {@code xIterSrc_} are treated as though they were beyond (before) the
* range of the fn up to (after) the discontinuity, and are thus
* extrapolated from that range. By implication, if a value in
* {@code xIterSrc_} is on a discontinuity it's corresponding y-value will
* never be used.
*
* Note that xIterTgt_, xIterSrc_ and discontinuities_ must be monotonous,
* the latter two strictly so, and in the same direction.
*
* @param targetY_toFill_
* The area to be filled with the resulting y values
* @param tgtXIter_
* An iterator over the target x values
* @param srcYIter_
* An iterator over the source y values
* @param srcXIter_
* An iterator over the source x values
* @param discontinuities_
* An iterator over the x values of discontinuities
* @param constrained_
* Indicates that the function is constrained, so no
* extrapolation is performed prior to the first source x
*/
static void interpSingleDiscontinuous(final DoubleArray targetY_toFill_, DoubleSettableIndexedIterator tgtXIter_,
DoubleSettableIndexedIterator srcYIter_, DoubleSettableIndexedIterator srcXIter_, PrimitiveIterator.OfDouble discontinuities_,
boolean preserveMonotonicity_, boolean constrained_, Callback callBack_)
{
// Get the first two source x values
double lowerX = srcXIter_.nextDouble();
double upperX = srcXIter_.nextDouble();
boolean decreasing = lowerX > upperX;
// First two source Y values
double lowerY = srcYIter_.nextDouble();
double upperY = srcYIter_.nextDouble();
DoubleSettableIndexedIterator yIterTgt = targetY_toFill_.iterator();
double nextDiscontinuity = Double.NaN;
if (discontinuities_.hasNext())
{
nextDiscontinuity = discontinuities_.nextDouble();
}
double nextX = Double.NaN;
double nextY = Double.NaN;
double preDiscontinuityY = Double.NaN;
double yAtTarget = Double.NaN;
/*
* The loop iterates over the target X values and finds the appropriate
* Y value for each
*/
while (tgtXIter_.hasNext())
{
double currentXTarget = tgtXIter_.nextDouble();
yIterTgt.nextDouble();
// Iterate over the x values while...
while (
// .. upperX is below the target
(((upperX <= currentXTarget) ^ decreasing) ||
// ... OR lowerX is below the next discontinuity and upperX is above
// the next discontinuity
(((Precision.equals(lowerX, nextDiscontinuity, 1e-10) || (lowerX < nextDiscontinuity) ^ decreasing)) && ((currentXTarget >= nextDiscontinuity) ^ decreasing)))
&& srcXIter_.hasNext())
{
// If nextX is set then the discontinuous zone has already been
// reached
if (Double.isNaN(nextX))
{
// Peek at the next X
nextX = srcXIter_.nextDouble();
nextY = srcYIter_.nextDouble();
}
// If the next src X is beyond the discontinuity but the current
// target is not, then it is in the zone to be extrapolated
if (!Double.isNaN(nextDiscontinuity)
&& (Precision.equals(nextX, nextDiscontinuity, 1e-10) || ((nextX > nextDiscontinuity) ^ decreasing))
&& ((currentXTarget < nextDiscontinuity) ^ decreasing))
{
break;
}
// Once lower X is moved on the first time, the constraint has
// been passed
constrained_ = false;
lowerX = upperX;
upperX = nextX;
nextX = Double.NaN;
// Moving past a disconinuity
if (!Double.isNaN(nextDiscontinuity) && ((lowerX > nextDiscontinuity) ^ decreasing))
{
// Remember the last y achieved before the discontinuity
preDiscontinuityY = yAtTarget;
if (discontinuities_.hasNext())
{
nextDiscontinuity = discontinuities_.nextDouble();
} else
{
nextDiscontinuity = Double.NaN;
}
}
lowerY = upperY;
upperY = nextY;
nextY = Double.NaN;
}
yAtTarget = interpValue(currentXTarget, lowerX, upperX, lowerY, upperY, decreasing, constrained_);
if (preserveMonotonicity_)
{
// If there is a
if (!Double.isNaN(preDiscontinuityY))
{
// was the y extrapolated?
if ((currentXTarget < lowerX) ^ decreasing)
{
// If the change from the discontinuity-Y to the target
// Y is in the opposite direction as from lower to
// upper,
// the extrapolation was beyond the discontiuity and it
// should be bounded
if (signum(upperY - lowerY) != signum(yAtTarget - preDiscontinuityY))
{
// constrain the value
yAtTarget = preDiscontinuityY;
}
} else
{
// Can't reset it, the extrapolation zone has passed
preDiscontinuityY = Double.NaN;
}
}
// If there is a nextY, the value was extrapolated before the
// discontinuity.
// Make sure it does not exceed the next Y after the
// discontinuity.
if (!Double.isNaN(nextY))
{
// was the y extrapolated?
if ((currentXTarget > upperX) ^ decreasing)
{
// If the change from the discontinuity-Y to the target
// Y is in the opposite direction as from lower to
// upper,
// the extrapolation was beyond the discontiuity and it
// should be bounded
if (signum(upperY - lowerY) != signum(nextY - yAtTarget))
{
// constrain the value
yAtTarget = nextY;
}
}
}
}
// Store the interpolated y value
yIterTgt.set(yAtTarget);
// Semantically, this check is unnecessary since the NULL_CALLBACK
// anyway does nothing, but the
// supporting bits are a performance drain
// TODO: Refactor callback API to put the supporting bits inside, or
// otherwise wrap with another callback
// to avoid this check
if (callBack_ != Callback.NOOP)
{
int[] currentTargetIndex = (tgtXIter_ instanceof DoubleSubspaceSplitIterator) ? ((DoubleSubspaceSplitIterator) tgtXIter_)
.getOtherIndex() : tgtXIter_.getIndex();
int currentUpperIndex = srcXIter_.getIndex()[0];
double lowerPropForPoint = (upperY - yAtTarget) / (upperY - lowerY);
callBack_.callback(yAtTarget, currentTargetIndex, currentUpperIndex - 1, lowerPropForPoint);
}
}
}
/*
* Calculates the interpolated Y value as appropriate
*/
private static double interpValue(double targetX_, double lowerSrcX_, double upperSrcX_, double lowerSrcY_, double upperSrcY_,
boolean decreasing_, boolean constrained_)
{
double yAtTarget;
// if the lower source x is beyond the target x then the value will be
// either constrained to the
// lower y or extrapolated
if (constrained_ && ((lowerSrcX_ > targetX_) ^ decreasing_))
{
// The current lowerY should be used
return lowerSrcY_;
}
if (Double.isInfinite(upperSrcX_))
{
yAtTarget = lowerSrcY_;
} else if (upperSrcX_ == lowerSrcX_)
{
throw new MultiDimensionalArrayException("Source x values must be strictly monotonous, but identical successive values were found");
} else
{
if (Double.isInfinite(upperSrcY_))
{
yAtTarget = upperSrcY_;
} else if (Double.isInfinite(lowerSrcY_))
{
yAtTarget = lowerSrcY_;
} else
{
double lowerPropForPoint = (upperSrcX_ - targetX_) / (upperSrcX_ - lowerSrcX_);
yAtTarget = upperSrcY_ - (upperSrcY_ - lowerSrcY_) * lowerPropForPoint;
}
}
// TODO: When does this happen?
if (Double.isNaN(yAtTarget))
{
LOG.severe(String.format("NaN result, upperY: %s, lowerY: %s, upperX: %s, lowerX: %s,currentXTarget: %s", upperSrcY_, lowerSrcY_,
upperSrcX_, lowerSrcX_, targetX_));
throw new RuntimeException("WTF?!");
}
return yAtTarget;
}
private static void addBoundingPoints(DoubleArray Y_, Interpolation.Specification[] dimSpecs_, int[] upperIndices_,
double[] upperProportions_, DoubleArray[] slices_, double[] proportions_, int[] index_)
{
addBoundingPointsRecursive(Y_, dimSpecs_, upperIndices_, upperProportions_, slices_, proportions_, index_, 1d, 0, 0);
}
private static int addBoundingPointsRecursive(DoubleArray Y_, Interpolation.Specification[] dimSpecs_, int[] upperIndices_,
double[] upperProportions_, DoubleArray[] slices_, double[] proportions_, int[] index_, double currentProportion, int i_,
int propCount_)
{
if (upperProportions_[i_] != 1d)
{
// Set the lower bound for this dimension
index_[dimSpecs_[i_].dimension] = upperIndices_[i_] - 1;
// Multiply by the lower proportion for this dimension
double proportion = currentProportion * (1 - upperProportions_[i_]);
if (i_ == (dimSpecs_.length - 1))
{
slices_[propCount_] = Y_.at(index_);
proportions_[propCount_++] = proportion;
} else
{
// Otherwise call recursively up to the nextDouble dimension
propCount_ = addBoundingPointsRecursive(Y_, dimSpecs_, upperIndices_, upperProportions_, slices_, proportions_, index_, proportion,
i_ + 1, propCount_);
}
}
/*
* Now the exact same for the upper bound in this dimension
*/
if (upperProportions_[i_] != 0d)
{
// Set the upper bound for this dimension
index_[dimSpecs_[i_].dimension] = upperIndices_[i_];
// Multiply by the upper proportion for this dimension
double proportion = currentProportion * upperProportions_[i_];
if (i_ == (dimSpecs_.length - 1))
{
slices_[propCount_] = Y_.at(index_);
proportions_[propCount_++] = proportion;
} else
{
// Otherwise call recursively up to the nextDouble dimension
propCount_ = addBoundingPointsRecursive(Y_, dimSpecs_, upperIndices_, upperProportions_, slices_, proportions_, index_, proportion,
i_ + 1, propCount_);
}
}
return propCount_;
}
/**
* Interpolates array Y to a specific target value along a given dimension,
* assuming the array represents the output of a function along that
* dimension and the input values corresponding to points in the array are
* given by X
*
* @param Y_
* The y-values of the function
* @param X_
* The x-values of gridpoints along the given dimension
* @param target_
* The target value
* @param dimension_
* The dimension along which to interpolate
*
* @return An array with one fewer dimension than Y, where the dimension
* dimension_ has been removed, and for which each value corresponds
* to the interpolation result along that dimension.
*/
public static DoubleArray interpDimension(DoubleArray Y_, DoubleArray X_, double target_, int dimension_)
{
int upperXIndex = 1;
double lowerX = X_.get(0);
double upperX = X_.get(1);
double direction = Math.signum(upperX - lowerX);
while ((upperX - target_) * direction < 0 && upperXIndex < (X_.numberOfElements() - 1))
{
upperXIndex++;
lowerX = upperX;
upperX = X_.get(upperXIndex);
}
double upperProportion = (target_ - lowerX) / (upperX - lowerX);
double lowerProportion = 1 - upperProportion;
int[] lowerIndex = Utils.repeatArray(-1, Y_.numberOfDimensions());
lowerIndex[dimension_] = upperXIndex - 1;
int[] upperIndex = Arrays.copyOf(lowerIndex, lowerIndex.length);
upperIndex[dimension_]++;
return Y_.at(lowerIndex).multiply(lowerProportion).add(Y_.at(upperIndex).multiply(upperProportion));
}
}