com.meliorbis.numerics.fixedpoint.UnivariateBoundedNewtonFixedPointFinder 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.fixedpoint;

import java.util.logging.Logger;

import com.meliorbis.numerics.NumericsException;
import com.meliorbis.numerics.function.MultiVariateVectorFunction;

/**
 * Does what it says on the tin.
 */
public class UnivariateBoundedNewtonFixedPointFinder implements FixedPointFinder
{
    private static final Logger LOG = Logger.getLogger(UnivariateBoundedNewtonFixedPointFinder.class.getName());

    private Double _min = Double.NEGATIVE_INFINITY;
    private Double _max = Double.POSITIVE_INFINITY;
    private double _damping;

    private final double _precision;
    private final double _gradDiffSize;

    public UnivariateBoundedNewtonFixedPointFinder(double precision_, double gradDiffSize_)
    {
        _precision = precision_;
        _gradDiffSize = gradDiffSize_;
        _damping = 1d;
    }

    @Override
    public double[] findFixedPoint(MultiVariateVectorFunction function_, double[] initialVals_)
    {
        if(initialVals_.length != 1)
        {
            throw new NumericsException("This fixed point finder is for univariate methods only!");
        }

        double lastInput = Double.NaN;
        double lastOutput = Double. NaN;

        double currentInput = initialVals_[0];

        double error = 1d;

        while(true)
        {
            // Calculate the implied value for the current input
            final Double currentOutput = function_.call(currentInput)[0];

            // Calculate the error
            error = Math.abs(Math.log(currentOutput/currentInput));

            LOG.fine("Error:  "+error);

            // If the error is smaller than the required precision - done!
            if(error < _precision)
            {
                LOG.fine("Fixed point found to required precision");
                return new double[]{currentInput};
            }

            // If lastInput is not set, we need to apply a small delta to calculate the grad
            if(Double.isNaN(lastInput))
            {
                lastInput = currentInput;
                lastOutput = currentOutput;

                if(currentInput == 0d) 
                {
                	currentInput = _gradDiffSize;
                }
                else 
                {
                	currentInput = currentInput*Math.exp(_gradDiffSize);
                }
            }
            else
            {
                // We have the two points for the grad, adjusting for fixed-point-ness
                double grad = (currentOutput - lastOutput)/(currentInput - lastInput)-1;

                // The last input is the new upper or lower bound, depending on the sign of the grad
                final double delta = -(lastOutput - lastInput) / grad;
                if(delta > 0)
                {
                    _min = lastInput;
                }
                else
                {
                    _max = lastInput;
                }

                // Adjust the input to remove the error
                double newInput = lastInput + delta *_damping;

                if(newInput > _max)
                {
                    LOG.fine("Computed point exceeds max, going half way");
                    newInput = (lastInput + _max)/2;

                }
                else if(newInput < _min)
                {
                    LOG.fine("Computed point exceeds min, going half way");
                    newInput = (lastInput + _min)/2;
                }

                LOG.fine(String.format("Min: %s\tMax: %s\tCurrent: %s",_min,_max,newInput));

                lastInput = Double.NaN;
                currentInput = newInput;
            }
        }
    }

    @Override
    public void setDamping(double damping_)
    {
        _damping = damping_;
    }
}