com.meliorbis.numerics.function.primitives.DoubleGridFunction Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of Numerics Show documentation
A library for working with large multi-dimensional arrays and the functions they represent
There is a newer version: 1.2
package com.meliorbis.numerics.function.primitives;

import java.util.ArrayList;
import java.util.List;

import com.meliorbis.numerics.function.FunctionException;
import com.meliorbis.numerics.function.grid.RectangularGridFunction;
import com.meliorbis.numerics.generic.primitives.DoubleArray;
import com.meliorbis.numerics.generic.primitives.impl.DoubleArrayFunctions;
import com.meliorbis.numerics.generic.primitives.impl.Interpolation;
import com.meliorbis.utils.Pair;
import com.meliorbis.utils.Utils;

/**
 * A function defined on a grid of doubles that is evaluated by interpolation.
 * 
 * This object is immutable. It can be used to evaluate a function but the function can not be changed.
 *
 * @author Tobias Grasl
 */
public class DoubleGridFunction implements MultiValuedDoubleFunction, DifferentiableDoubleFunction, RectangularGridFunction, DoubleArray>
{
	private final DoubleGridFunctionFactory _factory;
	
    private final DoubleRectangularGridDomain _domain;
    private final DoubleArray _values;
    private final int _dimensionality;
    
    DoubleGridFunction(DoubleRectangularGridDomain domain_, DoubleArray values_, DoubleGridFunctionFactory factory_)
    {
    	_domain = domain_;
    	_values = values_;
    	_factory = factory_;
    	
        _dimensionality = checkDimensions();
    }

	/**
	 * Checks that the dimensions of the arrays are consistent, throwing a runtime exception if not
	 * 
	 * @return An integer indicating the output dimensions.
 	 */
	protected int checkDimensions() 
	{
		if(getDomain().getNumberOfDimensions() > _values.numberOfDimensions())
        {
            throw new FunctionException("The number of dimensions of the value array must equal "
            		+ "or exceed by one the number of grid dimensions");
        }

        final int[] valueSize = _values.size();

        for(int dimension = 0; dimension < getDomain().getNumberOfDimensions(); dimension++)
        {
            if(getDomain().getAxis(dimension).numberOfDimensions() != 1)
            {
                throw new FunctionException("The gridpoint arrays must be one dimensional");
            }

            if(getDomain().getAxis(dimension).numberOfElements() != valueSize[dimension])
            {
                throw new FunctionException("The number of gridpoints in dimension "+dimension+
                    " must equal the size of the value array in that dimension");
            }
        }
        
        // If the values are the same size as the domain, the function is 1-dimensional
        if(getDomain().getNumberOfDimensions() == valueSize.length) 
        {
        	return 1;
        }
        else if(getDomain().getNumberOfDimensions() == valueSize.length-1)
        {
        	return valueSize[valueSize.length-1];
        }
        else
        {
        	 throw new FunctionException("The number of dimensions of the value array must equal "
             		+ "or exceed by one the number of grid dimensions");
        }
	}

    /**
     * @return The values at the grid coordinates that this function is defined on
     */
    public DoubleArray getValues()
    {
        return _values;
    }

    private Pair createDimSpecs(double[] values_){

        List dimSpecs = new ArrayList(getDomain().getNumberOfDimensions());
        List> remainingDims = new ArrayList>(getDomain().getNumberOfDimensions());

        // Create a specification to inerpolate each dimension, and collect the coordinates for the remaining dimensions
        for(int dimension = 0; dimension < getDomain().getNumberOfDimensions(); dimension++)
        {
            if(dimension >= values_.length || isInvalid(values_[dimension]))
            {
                remainingDims.add(getDomain().getAxis(dimension));
            }
            else
            {
                dimSpecs.add(new Interpolation.Specification(dimension, getDomain().getAxis(dimension), values_[dimension]));
            }
        }

        return new Pair(dimSpecs.toArray(new Interpolation.Specification[dimSpecs.size()]), 
        		_factory.createDomain(remainingDims));
    }

    private boolean isInvalid(double value_)
    {
        return Double.isNaN(value_);
    }

    @Override
    public DoubleArray callWithDouble(double... inputs_)
    {
        if(inputs_.length != getDomain().getNumberOfDimensions())
        {
            throw new FunctionException(
                    "The number of inputs passed must equal the number of dimensions of the function");
        }

        return Interpolation.interp(_values, createDimSpecs(inputs_).getLeft());
    }

    @Override
    public DoubleGridFunction restrict(double... partialInputs_)
    {
        final DoubleArray values = _values;

        return restrictArray(values, partialInputs_);
    }

    private DoubleGridFunction restrictArray(DoubleArray values_, double[] partialInputs_)
    {
        final Pair dimSpecs =
                createDimSpecs(partialInputs_);


        return _factory.createFunction(dimSpecs.getRight(), Interpolation.interp(values_, dimSpecs.getLeft()));
    }

    @Override
    public DoubleGridFunction partialDerivative(int dimension_)
    {
        return _factory.createFunction(getDomain(), pdArray(_values, dimension_));
    }

    private DoubleArray pdArray(DoubleArray values_, int dimension_)
    {
        return DoubleArrayFunctions.gradWeigthedHarmonic(values_, getDomain().getAxis(dimension_), dimension_);
    }

    @Override
    public DoubleGridFunction partialDerivative(int dimension_, double... at_)
    {
        /* First, restrict along dimensions that are not being differentiated. This saves unnecessary differentiations
         */
        // Copy the restrictions to an array that is as long as there are dimensions
        final double[] atCopied = Utils.repeatArray(Double.NaN, getDomain().getNumberOfDimensions());
        
        System.arraycopy(at_,0, atCopied, 0, at_.length);

        // Remember whether the dimension being differentiated is also restricted, but unset the restriction for now
        final double differentialDimensionRestrict = atCopied[dimension_];
        atCopied[dimension_] = Double.NaN;

        // Restrict the array according to the other restriction
        DoubleGridFunction restricted = restrictArray(_values, atCopied);

        /* Then, perform the differentiation
         */
        // Determine how many dimensions prior to the dimension being differentiated were restricted, and have thus
        // dropped out
        int i = 0;
        int lowerRestrictedDims = 0;

        while(i < dimension_)
        {
            if(!Double.isNaN(atCopied[i++]))
            {
                lowerRestrictedDims++;
            }
        }

        // Differentiate along the dimension requested, but taking into account the now-missing dimensions.
        // NOTE: This is assigning the values back to avoid creating a new object
        DoubleGridFunction derivative = _factory.createFunction(restricted.getDomain(), restricted.pdArray(restricted._values, dimension_ - lowerRestrictedDims));

        /* Finally, if necessary, restrict the differentiated dimension
         */
        // If the differentiated dimension was not restricted there is nothing left to do
        if(isInvalid(differentialDimensionRestrict))
        {
            return derivative;
        }

        // Otherwise, restrict that dimension also and return the result
        final double[] furtherRestrict = Utils.repeatArray(Double.NaN, dimension_ - lowerRestrictedDims + 1);
        furtherRestrict[dimension_ - lowerRestrictedDims] = differentialDimensionRestrict;

        return derivative.restrict(furtherRestrict);
    }

	public DoubleRectangularGridDomain getDomain()
	{
		return _domain;
	}

	public int getDimensionality()
	{
		return _dimensionality;
	}
}