/*
 *  Copyright 2013 Chris Pheby
 *
 *  Licensed 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
 *
 *      http://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.jadira.quant.coremath;

import org.jadira.quant.api.CompositeFunction;
import org.jadira.quant.api.QuantitativeFunction;
import org.jadira.quant.exception.IntervalBisectionOutOfRangeException;

/**
 * This class can be used to solve the equation f(x) = 0. It will determine the double value of x,
 * where f is a continuous function defined on an interval between [a, b] and f(a) and f(b) have opposite
 * signs. [a, b] are said to bracket the root since there must be at least one root between a and b that
 * yields f(x) = 0.
 *
 * For each iteration, the approach divides the interval in two by determining the midpoint of the interval, c,
 * such that c = (a+b)/2. The method then determines the interval that brackets the root, this being either
 * [a, c] or [c, b]. This is used as the interval for the next iteration. The process is terminated either
 * if the maximum number of iterations is reached, the gap between this result and the previous iteration
 * result is sufficiently small (if configured with {@link PrecisionStrategy#BETWEEN_RESULTS}) or the value is
 * within the desired distance of root (if configured with {@link PrecisionStrategy#TO_INTERVAL}.
 *
 * PrecisionStrategy defaults to {@link PrecisionStrategy#TO_INTERVAL}, the method defaults to 20 iterations,
 * and an accuracy distance of 1e-3 (0.001).
 *
 * Inspired by the example that appears in Java Methods for Financial Engineering (Philip Barker),
 * this implementation includes error handling, a functional programming model and two precision
 * strategies.
 */
public class IntervalBisection implements CompositeFunction<Double, QuantitativeFunction<Double, Double>> {

        private final double precision;
        private final int iterations;
        private final double lowerBound;
        private final double higherBound;
        private PrecisionStrategy precisionStrategy;

        public enum PrecisionStrategy {
                TO_INTERVAL, BETWEEN_RESULTS;
        }

        public IntervalBisection(double lowerBound, double higherBound) {
                this(lowerBound, higherBound, 20);
        }

        /**
         * Create a new instance with the given number of iterations and precision
         * @param lowerBound The minimum bound for the range to converge from
         * @param higherBound The maximum bound for the range to converge from
         * @param iterations The number of iterations to perform
         */
        public IntervalBisection(double lowerBound, double higherBound, int iterations) {
                // Defaults precision strategy to 1e-3
                this(lowerBound, higherBound, iterations, 0.001D, PrecisionStrategy.TO_INTERVAL);
        }

        /**
         * Create a new instance with the given number of iterations and precision
         * @param lowerBound The minimum bound for the range to converge from
         * @param higherBound The maximum bound for the range to converge from
         * @param iterations The number of iterations to perform
         * @param precision The requested precision
         * @param precisionStrategy Indicates approach for managing precision
         */
        public IntervalBisection(double lowerBound, double higherBound, int iterations, double precision, PrecisionStrategy precisionStrategy) {

                if (iterations < 1) {
                        throw new IllegalArgumentException("iterations must be greater than 1");
                }
                if (Double.compare(precision, 0.0D)  < 0) {
                        throw new IllegalArgumentException("precision must be a positive number");
                }
                if (precisionStrategy == null) {
                        throw new IllegalArgumentException("precisionStrategy may not be null");
                }

                this.lowerBound = lowerBound;
                this.higherBound = higherBound;
                this.iterations = iterations;
                this.precision = precision;
                this.precisionStrategy = precisionStrategy;
        }

        public double getLowerBound() {
                return lowerBound;
        }

        public double getHigherBound() {
                return higherBound;
        }

        public int getIterations() {
                return iterations;
        }

        public double getPrecision() {
                return precision;
        }

        @Override
        public Double apply(QuantitativeFunction<Double, Double> computeFunction) {

                final double resultA = computeFunction.apply(lowerBound);
                final double resultB = computeFunction.apply(higherBound);

                if (resultA != 0.0D && resultB != 0.0D
                                && (Double.compare(resultA, 0.0D) ^ Double.compare(0.0D, resultB)) != 0) {
                        throw new IntervalBisectionOutOfRangeException(lowerBound, higherBound);
                }

                double currentLower = lowerBound;
                double currentHigher = higherBound;

                double currentMiddle = Double.NaN;
                double midValueResult = Double.NaN;
                double preceedingMidValueResult = Double.NaN;

                for (int i = 0; i < iterations; i++) {

                        preceedingMidValueResult = midValueResult;

                        currentMiddle = currentLower + 0.5 * (currentHigher - currentLower);
                        midValueResult = computeFunction.apply(currentMiddle);

                        if (resultA * midValueResult < 0) {
                                currentHigher = currentMiddle;
                        } else if (resultA * midValueResult > 0) {
                                currentLower = currentMiddle;
                        }

                        if (PrecisionStrategy.TO_INTERVAL == precisionStrategy) {
                                if (Math.abs(midValueResult) <= precision) {
                                        break;
                                }
                        } else if ((PrecisionStrategy.BETWEEN_RESULTS == precisionStrategy || null == precisionStrategy)
                                        && (Math.abs(midValueResult - preceedingMidValueResult) <= precision)) {
                                        break;
                        }
                }
                return currentMiddle;
        }
}