com.meliorbis.numerics.generic.primitives.impl.DoubleArrayFunctions Maven / Gradle / Ivy
Show all versions of Numerics Show documentation
/**
*
*/
package com.meliorbis.numerics.generic.primitives.impl;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.lang.ArrayUtils;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import com.meliorbis.numerics.DoubleNumerics;
import com.meliorbis.numerics.Numerics;
import com.meliorbis.numerics.NumericsException;
import com.meliorbis.numerics.generic.MultiDimensionalArray;
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.DoubleReduction;
import com.meliorbis.numerics.generic.primitives.DoubleSettableIndexedIterator;
import com.meliorbis.numerics.generic.primitives.DoubleSettableIterator;
import com.meliorbis.numerics.generic.primitives.DoubleSubspaceSplitIterator;
import com.meliorbis.numerics.generic.primitives.DoubleUnaryOp;
import com.meliorbis.utils.Pair;
import com.meliorbis.utils.Utils;
/**
* Some functions that could be in the array class but are too much of a pain to
* implement generically
*
* @author Tobias Grasl
*/
public final class DoubleArrayFunctions
{
private static DoubleBlockedArrayFactory _arrayFactory;
private static Numerics _numerics;
static
{
_numerics = DoubleNumerics.instance();
_arrayFactory = _numerics.getArrayFactory();
}
private DoubleArrayFunctions()
{
}
public static DoubleArray diff(final DoubleArray Y_, int dimension_)
{
// Create an identically sizes array as this to hold the result
final DoubleArray result = _arrayFactory.newArray(Y_.size());
Y_.across(dimension_).reduce((DoubleIndexedReduction) iterator_ -> {
DoubleSubspaceSplitIterator subSpaceIterator = (DoubleSubspaceSplitIterator) iterator_;
double y1 = subSpaceIterator.nextDouble();
int[] priorIndex = subSpaceIterator.getFullIndex();
while (subSpaceIterator.hasNext())
{
double y2 = subSpaceIterator.nextDouble();
double diff = y2 - y1;
result.set(diff, priorIndex);
priorIndex = subSpaceIterator.getFullIndex();
y1 = y2;
}
// Need to set the final one too!
result.set(0d, priorIndex);
return 0d;
});
return result;
}
/**
* Calculates the gradient of the grid along the given dimension, assuming
* that the points on the grid are y-values and the points given in points_
* are x-values.
*
* The grad is calculated as the weighted mean of the upper and lower
* gradients except at the points on the boundary, where it is only the
* upper or lower as appropriate. The weighting as in Fritsch & Carlson,
* 1980. The gradient is the difference between two successive points
* divided by the difference between the corresponding x values.
*
* @param Y_
* The function values
* @param X_
* The x-values assumed at each grid point in the given
* dimension. If X_ is one dimensional, they are constant across
* other dimensions; if not, then X_ must be the same size as Y_
* and each point can have a value for X_.
* @param dimension_
* The dimension along the specified dimension
* @return An array of the same dimension as this one with the gradient at
* each point
*/
public static DoubleArray gradWeigthedHarmonic(final DoubleArray Y_, final DoubleArray X_, int dimension_)
{
final boolean commonX = (X_.numberOfDimensions() == 1 || (X_.numberOfDimensions() == 2 && ArrayUtils.contains(X_.size(), 1)));
// Create an identically sizes array as this to hold the result
final DoubleArray result = _arrayFactory.newArray(Y_.size());
Y_.across(new int[]
{ dimension_ }).reduce(new DoubleIndexedReduction()
{
@Override
public double perform(DoubleSubspaceSplitIterator iterator_) throws RuntimeException
{
DoubleSubspaceSplitIterator subSpaceIterator_ = (DoubleSubspaceSplitIterator) iterator_;
DoubleSettableIterator xIter;
if (commonX)
{
xIter = (DoubleSettableIterator) X_.iterator();
} else
{
xIter = (DoubleSettableIterator) X_.iteratorAt(subSpaceIterator_.getFullIndex());
}
if (!subSpaceIterator_.hasNext())
{
return Double.NaN;
}
double y1 = subSpaceIterator_.nextDouble();
double x1 = xIter.nextDouble();
double priorXDiff = Double.NaN;
int[] priorIndex = subSpaceIterator_.getFullIndex();
double priorGrad = Double.NaN;
while (subSpaceIterator_.hasNext())
{
double y2 = subSpaceIterator_.nextDouble();
double x2 = xIter.nextDouble();
double xDiff = x2 - x1;
double currentGrad = (y2 - y1) / xDiff;
// Use the weighted Harmonic mean a la Fritsch & Carlson,
// except at the first point where it is just the current
// grad
double grad = Double.isNaN(priorGrad) ? currentGrad :
// If the signs are different or one is 0, the outcome
// should be 0
((priorGrad * currentGrad <= 0) ? 0 : 3d * (xDiff + priorXDiff)
/ ((2d * xDiff + priorXDiff) / priorGrad + (xDiff + 2d * priorXDiff) / currentGrad));
result.set(grad, priorIndex);
priorIndex = subSpaceIterator_.getFullIndex();
priorGrad = currentGrad;
priorXDiff = xDiff;
y1 = y2;
x1 = x2;
}
// Need to set the final one too! - this only considers the
// final
// segment, so that extrapolation is linear
result.set(priorGrad, priorIndex);
return Double.NaN;
}
});
return result;
}
/**
* Calculates the gradient of the grid along the given dimension, assuming
* that the points on the grid are y-values and the points given in points_
* are x-values.
*
* The grad is calculated as the weighted mean of the upper and lower
* gradients except at the points on the boundary, where it is only the
* upper or lower as appropriate. The weighting is linear by distance. The
* gradient is the difference between two successive points divided by the
* difference between the corresponding x values.
*
* @param Y_
* The function values
* @param X_
* The x-values assumed at each grid point in the given
* dimension. If X_ is one dimensional, they are constant across
* other dimensions; if not, then X_ must be the same size as Y_
* and each point can have a value for X_.
* @param dimension_
* The dimension along the specified dimension
*
* @return An array of the same dimension as this one with the gradient at
* each point
*/
public static DoubleArray grad(final DoubleArray Y_, final DoubleArray X_, int dimension_)
{
// Create an identically sizes array as this to hold the result
final DoubleArray result = _arrayFactory.newArray(Y_.size());
Y_.across(dimension_).reduce(new DoubleIndexedReduction()
{
@Override
public double perform(DoubleSubspaceSplitIterator iterator_) throws RuntimeException
{
DoubleSettableIterator xIter = (DoubleSettableIterator) X_.iterator();
if (!iterator_.hasNext())
{
return 0d;
}
double y1 = iterator_.nextDouble();
double x1 = xIter.nextDouble();
double priorXDiff = Double.NaN;
int[] priorIndex = iterator_.getFullIndex();
double priorGrad = Double.NaN;
while (iterator_.hasNext())
{
double y2 = iterator_.nextDouble();
double x2 = xIter.nextDouble();
double xDiff = x2 - x1;
double currentGrad = (y2 - y1) / xDiff;
// The proportion to assign to the prior - note this is the
// opposite of what it would
// be if we were calculating the average gradient
double priorProp = xDiff / (xDiff + priorXDiff);
double grad = Double.isNaN(priorGrad) ? currentGrad : (priorGrad * priorProp + currentGrad * (1 - priorProp));
result.set(grad, priorIndex);
priorIndex = iterator_.getFullIndex();
priorGrad = currentGrad;
priorXDiff = xDiff;
y1 = y2;
x1 = x2;
}
// Need to set the final one too!
result.set(priorGrad, priorIndex);
return 0d;
}
});
return result;
}
public static DoubleArray interpolateFixedPoint(DoubleArray fn_, final DoubleArray levels_, int dimension_)
{
DoubleReduction fpOp = new DoubleReduction()
{
@Override
public double perform(DoubleSettableIterator valuesIterator_) throws RuntimeException
{
DoubleSettableIndexedIterator levelIter;
levelIter = levels_.iterator();
double lastX = Double.NEGATIVE_INFINITY;
double lastY = Double.NEGATIVE_INFINITY;
while (valuesIterator_.hasNext())
{
double yVal = valuesIterator_.nextDouble();
double xVal = levelIter.nextDouble();
// Did we pass the fixed point?
if (xVal - yVal > 1e-5)
{
// Linearly interpolate the fixed point
return (lastX * (yVal - lastY) - lastY * (xVal - lastX)) / (yVal + lastX - lastY - xVal);
} else
{
lastX = xVal;
lastY = yVal;
}
}
return Double.NaN;
}
};
return fn_.across(dimension_).reduce(fpOp);
}
/**
* Finds the unique internal fixed point of the mapping from the regular
* grid defined by the provided x values to the provided y values.
*
* If there is no internal fixed point, or more than one, an
* MultiDimensionalArrayException is thrown
*
* @param y_
* The y values of the mapping
* @param xCoords_
* The x values along the dimensions of y that define the regular
* grid on which the mapping is defined. The length of xCoords_
* must be equal to the number of dimensions of x_ minus 1,
* unless y_ is one dimensional and xCoords_ has just one element
* (since no additional dimension is required for the one-
* dimensional mapping target)
*
* @return A one dimensional array of length ndim(y_) - 1
*
* @throws com.meliorbis.numerics.generic.MultiDimensionalArrayException
* If there are 0 or more than 1 internal fixed points, or
* another error occurs
*/
public static DoubleArray interpolateFixedPoint(DoubleArray y_, DoubleArray[] xCoords_) throws NumericsException
{
DoubleArray yAdjusted = y_.copy();
for (int i = 0; i < xCoords_.length; i++)
{
// Subtract the appropriate inputs from each output to turn the
// problem into a rootfinding one
(y_.numberOfDimensions() == 1 ? yAdjusted : yAdjusted.lastDimSlice(i)).modifying().across(i).subtract(xCoords_[i]);
}
try
{
return findRoot(yAdjusted, xCoords_);
} catch (NumericsException e_)
{
throw new NumericsException("No fixed point found", e_);
}
}
/**
* Finds the unique internal root of the mapping from the regular grid
* defined by the provided x values to the provided y values.
*
* If there is no root, or more than one, an NumericsException is thrown
*
* @param y_
* The y values of the mapping
* @param xCoords_
* The x values along the dimensions of y that define the regular
* grid on which the mapping is defined. The length of xCoords_
* must be equal to the number of dimensions of x_ minus 1,
* unless y_ is one dimensional and xCoords_ has just one element
* (since no additional dimension is required for the one-
* dimensional mapping target)
*
* @return A one dimensional array of length ndim(y_) - 1
*
* @throws NumericsException
* If there are 0 or more than 1 internal roots, or another
* error occurs
*/
public static DoubleArray findRoot(DoubleArray y_, DoubleArray[] xCoords_) throws NumericsException
{
final int nDims = xCoords_.length;
// Short-cut - which also assumes monotonicity
if (nDims == 1)
{
return Interpolation.interpDimension(xCoords_[0], y_, 0, 0);
}
double[] x_ = new double[nDims];
RealMatrix jacobean = MatrixUtils.createRealMatrix(nDims, nDims);
double relError;
int count = 0;
if (nDims == 1)
{
// Start half way between points either side of 0
for (int i = 0; i < y_.numberOfElements() - 1; i++)
{
if (y_.get(i) * y_.get(i + 1) < 0)
{
x_[0] = (xCoords_[0].get(i) + xCoords_[0].get(i + 1)) / 2d;
break;
}
}
if (x_[0] == 0.0d)
{
throw new NumericsException("y values do not cross x axis");
}
} else
{
for (int i = 0; i < nDims; i++)
{
// Subtract the appropriate inputs from each output to turn the
// problem into a rootfinding one
x_[i] = xCoords_[i].sum() / ((double) xCoords_[i].numberOfElements()) * 0.999;
}
}
do
{
final Pair, DoubleArray> jacAndY = jacobeanAtPoint(y_, xCoords_, x_);
final DoubleArray negError = jacAndY.getRight().multiply(-1d);
// Set the jacobean rows
for (int dim = 0; dim < nDims; dim++)
{
jacAndY.getLeft().at(dim, dim).modifying();
jacobean.setColumn(dim, jacAndY.getLeft().at(-1, dim).toArray());
}
RealMatrix solved = new LUDecomposition(jacobean).getSolver().solve(MatrixUtils.createColumnRealMatrix(negError.toArray()));
final double[] lastX = x_.clone();
relError = 0d;
solved = solved.scalarMultiply(0.1);
for (int dim = 0; dim < nDims; dim++)
{
final double adjustment = solved.getEntry(dim, 0);
x_[dim] += adjustment;
relError = Math.max(relError, Math.abs(lastX[dim] == 0d ? adjustment : adjustment / lastX[dim]));
}
} while (relError > 1e-10 && count++ < 10000);
if (relError > 1e-10)
{
throw new NumericsException("No root found!");
}
return _arrayFactory.newArray(x_);
}
/**
* Calculates the jacobean of the mapping from the regular grid defined by
* xCoords to y at the point x. The mapping is assumed to be piecewise
* linear, and x must be internal.
*
* @param y_
* The y values of the mapping
* @param xCoords_
* Coordinates along each dimension that define the regular grid
* on which the mapping is specified
* @param x_
* The point at which to get the jacobean. It must be internal to
* the grid.
* @return The jacobean at the specified point
*/
static Pair, DoubleArray> jacobeanAtPoint(DoubleArray y_, DoubleArray[] xCoords_, double[] x_)
{
/*
* TODO: HIGHLY INEFFICIENT! FOR SPEED, CHANGE!
*/
int ndims = xCoords_.length;
// Create interpolation specs for each dimension to the provided point
// x_
List dimSpecs = new ArrayList(ndims);
for (int dim = 0; dim < ndims; dim++)
{
dimSpecs.add(Interpolation.spec(dim, xCoords_[dim], x_[dim]));
}
final DoubleArray yAtX = Interpolation.interp(y_, dimSpecs.toArray(new Interpolation.Specification[ndims]));
final DoubleArray result = _numerics.getArrayFactory().newArray(ndims, ndims);
for (int dim = 0; dim < ndims; dim++)
{
// Don't interpolate along the dimension the grad is to be
// calculated along
final Interpolation.Specification removed = dimSpecs.remove(dim);
final DoubleArray interpolated = Interpolation.interp(y_, dimSpecs.toArray(new Interpolation.Specification[ndims - 1]));
final DoubleArray gradForDim = grad(interpolated, xCoords_[dim],/*
* already
* interpolated
* to
* here
*/0);
final DoubleArray gradAtX = Interpolation.interpDimension(gradForDim, xCoords_[dim], x_[dim], 0);
result.fillAt(gradAtX, -1, dim);
// Put back the removed one for the nextDouble iteration
dimSpecs.add(dim, removed);
}
return new Pair, DoubleArray>(result, yAtX);
}
/**
* Interpolates a fixed point in a function (mapping) from one set of values
* along a given dimension.
*
* The method assumes that the input values are monotonic along the
* dimension of interpolation
*
* @param x_
* The array of inputs, of arbitrary dimension
* @param y_
* The array of function values, must be of the same dimension as
* {@code x_}
* @param dimension_
* The dimension along which to interpolate
*
* @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 interpolateFixedPoints(DoubleArray x_, DoubleArray y_, int dimension_)
{
return interpolateFixedPoints(x_, y_, Double.NaN, dimension_);
}
/**
* Interpolates a fixed point function (mapping) from one set of values
* along a given dimension.
*
* The method assumes that the input values are monotonic along the
* dimension of interpolation
*
* @param x_
* The array of inputs, of arbitrary dimension
* @param y_
* The array of function values, must be of the same dimension as
* {@code x_}
*
* @param anchor_
* In the case that there are multiple fixed points, returns the
* one closed to anchor
* @param dimension_
* The dimension along which to interpolate
*
* @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 interpolateFixedPoints(DoubleArray x_, DoubleArray y_, double anchor_, int dimension_)
{
// Get an iterator across the dimensions not being interpolated
DoubleSubspaceSplitIterator nonInterpXIter = x_.iteratorAcross(new int[]
{ dimension_ }).getOrthogonalIterator();
DoubleSubspaceSplitIterator nonInterpYIter = y_.iteratorAcross(new int[]
{ dimension_ }).getOrthogonalIterator();
DoubleArray results = _arrayFactory.newArray(nonInterpXIter.getSizes());
DoubleSettableIndexedIterator resultIter = results.iterator();
while (nonInterpXIter.hasNext())
{
nonInterpXIter.nextDouble();
nonInterpYIter.nextDouble();
resultIter.nextDouble();
// Get an iterator along the dimension we are to interpolate
DoubleSubspaceSplitIterator interpDimXIter = nonInterpXIter.getOrthogonalIterator();
DoubleSubspaceSplitIterator interpDimYIter = nonInterpYIter.getOrthogonalIterator();
// Initial upper and lower X values; the upper X will definitely be
// set at least once in the loop
int[] currentIndex = interpDimXIter.getFullIndex();
currentIndex[dimension_]++;
double priorX = Double.NaN;
double priorY = Double.NaN;
double priorDiff = Double.NaN;
double minX = Double.NaN;
double maxX = Double.NaN;
List fixedPoints = new ArrayList(1);
while (interpDimXIter.hasNext())
{
double currentX = interpDimXIter.nextDouble();
double currentY = interpDimYIter.nextDouble();
double currentDiff = currentY - currentX;
if (Double.isNaN(minX))
{
minX = currentX;
}
// Are the signs of the two diffs the same?
if (priorDiff * currentDiff < 0)
{
// The signs are different, i.e. one Y is lower than the X
// and the other higher, hence they must cross-over
// inbetween
// Set the value and continue;
fixedPoints.add((priorX * (currentY - priorY) - priorY * (currentX - priorX)) / (currentY - priorY - (currentX - priorX)));
}
// Move to the nextDouble point
priorX = currentX;
priorY = currentY;
priorDiff = currentDiff;
maxX = currentX;
}
if (fixedPoints.size() == 0)
{
resultIter.set(Double.NaN);
} else if (fixedPoints.size() == 1)
{
resultIter.set(fixedPoints.get(0).doubleValue());
} else
{
double minDistanceFromEdge = Double.NEGATIVE_INFINITY;
double currentMatch = Double.NaN;
for (double fixedPoint : fixedPoints)
{
double distForPoint;
if (Double.isNaN(anchor_))
{
distForPoint = Math.min(fixedPoint - minX, maxX - fixedPoint);
} else
{
distForPoint = -Math.abs(fixedPoint - anchor_);
}
if (distForPoint >= minDistanceFromEdge)
{
minDistanceFromEdge = distForPoint;
currentMatch = fixedPoint;
}
}
resultIter.set(currentMatch);
}
}
return results;
}
public static DoubleUnaryOp cutToBounds(double min_, double max_)
{
return new CutToBounds(min_, max_);
}
/**
* Returns the input if it is within the bounds, otherwise, the bound it is
* beyond.
*
* @author toby
*/
public static final class CutToBounds implements DoubleUnaryOp
{
private final double _max;
private final double _min;
private CutToBounds(double min_, double max_)
{
_max = max_;
_min = min_;
}
@Override
public double perform(double input_)
{
return (input_ < _min) ? _min : (input_ > _max) ? _max : input_;
}
}
public static double maximumRelativeDifferenceSpecial(DoubleArray a_, DoubleArray b_)
{
double criterion;
DoubleSettableIndexedIterator bIter = b_.iterator();
DoubleSettableIndexedIterator aIter = a_.iterator();
criterion = 0d;
while (bIter.hasNext())
{
double newVal = aIter.nextDouble();
double oldVal = bIter.nextDouble();
if (Math.abs(oldVal) > 1e-5 && Math.abs(newVal) > 1e-5)
{
criterion = Math.max(criterion, Math.abs((newVal - oldVal) / (1 + Math.abs(oldVal))));
} else
{
criterion = Math.max(criterion, Math.abs(newVal - oldVal));
}
}
return criterion;
}
public static double maximumRelativeDifference(DoubleArray aggKapitalTransitionRule, DoubleArray capitalTransition)
{
double criterion;
DoubleSettableIndexedIterator newTransIter = capitalTransition.iterator();
DoubleSettableIndexedIterator oldTransIter = aggKapitalTransitionRule.iterator();
criterion = 0d;
while (newTransIter.hasNext())
{
double oldVal = oldTransIter.nextDouble();
double newVal = newTransIter.nextDouble();
if (Math.abs(oldVal) > 1e-5 && Math.abs(newVal) > 1e-5)
{
criterion = Math.max(criterion, Math.abs(newVal / oldVal - 1));
} else
{
criterion = Math.max(criterion, Math.abs(newVal - oldVal));
}
}
return criterion;
}
/**
* Given two multi-dimensional transition matrices, combines them as though
* they were independent.
*
* @param firstTrans_
* The first transition matrix. Must have as many source as
* target states.
* @param secondTrans_
* The second transition matrix. Must have as many source as
* target states.
*
* @return A transition matrix between the combined sets of states
*/
public static DoubleArray combineTransitions(DoubleArray firstTrans_, DoubleArray secondTrans_)
{
int[] firstSize = firstTrans_.size();
int[] secondSize = secondTrans_.size();
int[] resultSize = new int[firstSize.length + secondSize.length];
int firstSizeElements = firstSize.length / 2;
System.arraycopy(firstSize, 0, resultSize, 0, firstSizeElements);
int secondSizeElements = secondSize.length / 2;
System.arraycopy(secondSize, 0, resultSize, firstSizeElements, secondSizeElements);
System.arraycopy(firstSize, firstSizeElements, resultSize, firstSizeElements + secondSizeElements, firstSizeElements);
System.arraycopy(secondSize, secondSizeElements, resultSize, firstSize.length + secondSizeElements, secondSizeElements);
DoubleArray result = _arrayFactory.newArray(resultSize);
result.fillDimensions(
firstTrans_,
ArrayUtils.addAll(Utils.sequence(0, firstSizeElements),
Utils.sequence(firstSizeElements + secondSizeElements, 2 * firstSizeElements + secondSizeElements)));
return result.across(
ArrayUtils.addAll(Utils.sequence(firstSizeElements, firstSizeElements + secondSizeElements),
Utils.sequence(2 * firstSizeElements + secondSizeElements, 2 * firstSizeElements + 2 * secondSizeElements))).multiply(
secondTrans_);
}
public static final DoubleUnaryOp pow(final double power_)
{
return x -> Math.pow(x, power_);
}
public static final DoubleUnaryOp abs = x -> Math.abs(x);
public static final DoubleUnaryOp exp = x -> Math.exp(x);
public static final DoubleUnaryOp log = x -> Math.log(x);
public static final DoubleNaryOp weightedMean(final double... weights_)
{
double finalWeightCalc = 1d;
for (int i = 0; i < weights_.length; i++)
{
finalWeightCalc -= weights_[i];
}
final double finalWeight = finalWeightCalc;
return inputs_ -> {
double sum = 0d;
for (int i = 0; i < weights_.length; i++)
{
sum += inputs_[i] * weights_[i];
}
sum += finalWeight * inputs_[weights_.length];
return sum;
};
}
/**
* Gets a sequence of values that are numLogs-log distributed, i.e.
* e^numLogs times of the values will be evenly distributed
*
* @param maxValue_
* The maximum value
* @param numValues_
* The number of values
* @param numLogs_
* How many times to apply log
*
* @return A one dimensional array of values as described
*/
static public DoubleArray nLogSequence(double maxValue_, int numValues_, int numLogs_)
{
int logCount = 0;
// Take the log a number of times specific, adding 1 at each juncture
// (because log 1 = 0, but we actually want the base sequence to have a
// 0)
while (logCount++ < numLogs_)
{
maxValue_ = Math.log(maxValue_ + 1);
}
// Create an equally spaced sequence between 0 and the modified maxValue
double[] capLevels = Utils.sequence(0d, maxValue_, numValues_);
// Create an array with the values
DoubleArray values = _numerics.getArrayFactory().newArray(capLevels);
// Now reverse the previous step on each value
while (--logCount > 0)
{
values.modifying().map((DoubleUnaryOp) input_ -> Math.exp(input_) - 1);
}
return values;
}
/**
* Creates an array where values are more tightly clustered around the mean
*
* @param values_
* The array to pull
* @param maxMeanWeight_
* The maximum proportional weight to give to the mean when
* calculating the new value at each point. Points furthest from
* the mean are not pulled in
* @return The pulled array
*/
static public DoubleArray pullToMean(final DoubleArray values_, final double maxMeanWeight_)
{
final double meanValue = values_.sum() / ((double) values_.numberOfElements());
final double distToMax = values_.max() - meanValue;
final double distToMin = meanValue - values_.min();
final double highestDist = Math.max(distToMax, distToMin);
DoubleArray individualTransition = values_.map((DoubleUnaryOp) input_ -> {
double meanWeight = (1 - Math.abs(meanValue - input_) / highestDist) * maxMeanWeight_;
return input_ * (1d - meanWeight) + meanValue * meanWeight;
});
return individualTransition;
}
/**
* Creates an array where each value is pulled closer to the mean by the
* same proportion
*
* @param values_
* The array to squeeze
*
* @param meanWeight_
* The proportional weight to give to the mean when calculating
* the new value at each point
*
* @return The squeezed array
*/
static public DoubleArray squeezeToMean(final DoubleArray values_, final double meanWeight_)
{
return squeezeTo(values_, values_.sum() / ((double) values_.numberOfElements()), meanWeight_);
}
/**
* Creates an array where each value is pulled closer to the target by the
* same proportion
*
* @param values_
* The array to squeeze
*
* @param target_
* The target to squeeze toward
*
* @param targetWeight_
* The proportional weight to give to the mean when calculating
* the new value at each point
*
* @return The squeezed array
*/
static public DoubleArray squeezeTo(final DoubleArray values_, final double target_, final double targetWeight_)
{
final double weightedTarget = targetWeight_ * target_;
return values_.map((DoubleUnaryOp) input_ -> input_ * (1d - targetWeight_) + weightedTarget);
}
public static void ensureAllEqualSize(MultiDimensionalArray[] arrays_)
{
int[] size = arrays_[0].size();
ensureAllOfSize(size, arrays_, 1);
}
public static void ensureAllEqualSize(MultiDimensionalArray reference_, MultiDimensionalArray[] arrays_)
{
int[] size = reference_.size();
ensureAllOfSize(size, arrays_, 0);
}
private static void ensureAllOfSize(int[] size_, MultiDimensionalArray[] arrays_, int startAt_)
{
for (int i = startAt_; i < arrays_.length; i++)
{
if (!equalSize(size_, arrays_[i].size()))
{
throw new MultiDimensionalArrayException("Arrays must be of same size");
}
}
}
static private boolean equalSize(int[] a_, int[] b_)
{
int aIndex = 0;
int bIndex = 0;
// Iterate over the first size array
while (aIndex < a_.length)
{
// If we are done with the second one already, the first one must
// only have 1s
if (bIndex == b_.length)
{
if (a_[aIndex] == 1)
{
aIndex++;
continue;
}
return false;
}
// Standard case: they are the same size
if (a_[aIndex] == b_[bIndex])
{
aIndex++;
bIndex++;
continue;
}
// Otherwise, one of them may be of size 1, so a 'skipping'
// dimension
if (a_[aIndex] == 1)
{
aIndex++;
continue;
}
if (b_[bIndex] == 1)
{
bIndex++;
continue;
}
// Oops, not the same
return false;
}
// If there are still elements of b left they must be ones
while (bIndex < b_.length)
{
if (b_[bIndex++] != 1)
{
return false;
}
}
return true;
}
}