/**
 ********************************************************************************
 * Copyright (c) 2020 Robert Bosch GmbH and others.
 *
 * This program and the accompanying materials are made
 * available under the terms of the Eclipse Public License 2.0
 * which is available at https://www.eclipse.org/legal/epl-2.0/
 *
 * SPDX-License-Identifier: EPL-2.0
 *
 * Contributors:
 *     Robert Bosch GmbH - initial API and implementation
 ********************************************************************************
 */

package org.eclipse.app4mc.amalthea.model.util;

import org.apache.commons.math3.distribution.WeibullDistribution;
import org.apache.commons.math3.special.Gamma;

/**
 * @author mlh2sh
 *
 */
public class WeibullUtil {

	private static final String NO_VALID_WEIBULL_DISTRIBUTION_FOUND = "invalid argument: there is no valid weibull distribution for the used parameter";

	// Suppress default constructor
	private WeibullUtil() {
		throw new IllegalStateException("Utility class");
	}

	public static class Parameters {
		public double given_lowerBound = 0.0;
		public double given_upperBound = 0.0;
		public double requested_mean = 0.0;
		public double requested_pRemainPromille = 0.0;

		public double shape = 1.0;
		public double scale = 1.0;

		public String error = null;

		@Override
		public String toString() {
			String output = "WeibullParam"
					+ "\n  shape = " + shape
					+ "\n  scale = " + scale
					+ "\n  lowerBound      = " + given_lowerBound
					+ "\n  upperBound      = " + given_upperBound
					+ "\n  mean            = " + requested_mean
					+ "\n  pRemainPromille = " + requested_pRemainPromille;
			if (error != null) {
				output = output + "\n  *** error *** : " + error;
			}
			return output;
		}

		public String qualityReport() {
			// check with Apache Commons Math
			final WeibullDistribution mathFunction = new WeibullDistribution(null, shape, scale);
			double cdfAvg = mathFunction.cumulativeProbability(requested_mean - given_lowerBound);
			double cdfMax = mathFunction.cumulativeProbability(given_upperBound - given_lowerBound);
			double remPromilleC = (1.0 - cdfMax) * 1000.0;
			double averageC = WeibullUtil.computeAverage(shape, scale, given_lowerBound, given_upperBound);
			double medianC = WeibullUtil.computeMedian(shape, scale, given_lowerBound, given_upperBound);
			return "Weibull Parameter Quality"
				+ "\nInput: [" + given_lowerBound + ", " + requested_mean + ", " + given_upperBound + "], " + requested_pRemainPromille
				+ "\nEstimated parameters: shape = " + shape + "\tscale = " + scale
				+ "\n computed average: " + averageC
				+ "\n computed median: " + medianC
				+ "\n computed pRemainPromille: " + remPromilleC
				+ "\n computed cdf(lower interval) = " + cdfAvg
				+ "\n computed cdf(upper interval) = " + (cdfMax - cdfAvg);
		}
	}

	/**
	 * See Crénin, François, Truncated Weibull Distribution Functions and Moments (October 30, 2015).
	 * Available at SSRN: https://ssrn.com/abstract=2690255 or http://dx.doi.org/10.2139/ssrn.2690255
	 *
	 * @param shape
	 * @param scale
	 * @param location
	 * @param upperBound
	 * @return
	 */
	public static double computeAverage(double shape, double scale, double location, Double upperBound) {
		if (shape <= 0.0 || scale <= 0.0) return Double.NaN;

		if (upperBound == null || upperBound.isNaN() || upperBound.isInfinite()) {
			// Standard Weibull distribution
			return location + scale * Gamma.gamma(1 + 1 / shape);

		} else if (location < upperBound) {
			// Right truncated Weibull distribution
			final double range = upperBound - location;
			final double power = Math.pow(range / scale, shape);
			double factor1 = scale / (1 - Math.exp(-1.0 * power));
			double factor2 = lowerIncompleteGammaFunction(1 + 1 / shape, power);
			return location + factor1 * factor2;
		}

		return Double.NaN;
	}

	/**
	 * See https://en.wikipedia.org/wiki/Incomplete_gamma_function
	 *
	 * @param s
	 * @param x
	 * @return
	 */
	private static double lowerIncompleteGammaFunction(double s, double x) {
		return Gamma.regularizedGammaP(s, x) * Gamma.gamma(s);
	}

	public static double computeMedian(double shape, double scale, double location, Double upperBound) {
		double pRemainPromille;
		if (upperBound == null || upperBound.isNaN() || upperBound.isInfinite()) {
			pRemainPromille = 0.0;
		} else {
			pRemainPromille = computePRemainPromille(shape, scale, location, upperBound);
		}

		return computeMedianWithPRemainPromille(shape, scale, location, pRemainPromille);
	}

	public static double computeMedianWithPRemainPromille(double shape, double scale, double location, double pRemainPromille) {
		if (shape <= 0.0 || scale <= 0.0 || pRemainPromille < 0 || pRemainPromille >= 1000) return Double.NaN;

		double x = scale * Math.pow(-1.0 * Math.log(0.5 + pRemainPromille / 2000.0), 1 / shape);
		return location + x;
	}

	public static double computePRemainPromille(double shape, double scale, double location, double upperBound) {
		if (shape <= 0.0 || scale <= 0.0) return Double.NaN;

		if (location < upperBound) {
			// Compute CDF at upper bound; pRemainPromille = (1 - CDF) * 1000
			final double range = upperBound - location;
			final double power = Math.pow(range / scale, shape);
			return Math.exp(-1.0 * power) * 1000.0;
		}

		return Double.NaN;
	}

	public static double computeUpperBound(double shape, double scale, double location, double pRemainPromille) {
		if (shape <= 0.0 || scale <= 0.0) return Double.NaN;

		double pRemain = pRemainPromille / 1000.0;
		if (pRemain > 0 && pRemain < 1) {
			// Compute CDF at upper bound; pRemain = (1 - CDF)
			return location + Math.pow(-1.0 * Math.log(pRemain), 1.0 / shape) * scale;
		}

		return Double.NaN;
	}

	/**
	 * This is a simple attempt to get a better parameter estimation for big pRemainPromille values.
	 *
	 * As an alternative to the regular estimation another estimation based on the median value is calculated.
	 * The better result is selected according to the relative error of both requested values.
	 *
	 * Better solutions are welcome !
	 *
	 * @param lowerBound
	 * @param average
	 * @param upperBound
	 * @param pRemainPromille
	 * @return
	 */
	public static Parameters findParameters(double lowerBound, double average, double upperBound, double pRemainPromille) {
		Parameters p1 = findParametersForAverage(lowerBound, average, upperBound, pRemainPromille);
		Parameters p2 = findParametersForMedian(lowerBound, average, upperBound, pRemainPromille);

		if (p2.error != null) return p1;

		// compute resulting values for estimated parameters

		double p1_average = computeAverage(p1.shape, p1.scale, lowerBound, upperBound);
		double p1_pRemainPromille = computePRemainPromille(p1.shape, p1.scale, lowerBound, upperBound);
		double p2_average = computeAverage(p2.shape, p2.scale, lowerBound, upperBound);
		double p2_pRemainPromille = computePRemainPromille(p2.shape, p2.scale, lowerBound, upperBound);

		// find the better one

		double error1 = Math.abs(1.0 - p1_average / average) + Math.abs(1.0 - p1_pRemainPromille / pRemainPromille);
		double error2 = Math.abs(1.0 - p2_average / average) + Math.abs(1.0 - p2_pRemainPromille / pRemainPromille);
		if (error2 < error1) return p2;

		return p1;
	}

	/**
	 * The approximation from the given parameters boil down to the following problem (latex formulas):
	 * <p>
	 * CDF_{weibull}(upperBound) = 1 - pRemainPromille/1000
	 * <p>
	 * Using the Weibull CDF (1) and the Expectation Value (average) (2)
	 * we can derive a formula for the scale dependent on shape and known values from (2)
	 * <p>
	 * (1) CDF_{weibull}(x) = 1 - e^{- (\frac{x}{scale})^{shape}}
	 * <p>
	 * (2) avg = scale \cdot \Gamma(1 + \frac{1}{shape})
	 * <p>
	 * (3) scale = \frac{avg}{\Gamma(1 + (1 / shape))}
	 * <p>
	 * with (1), (3) and pRemain = pRemainPromille/1000 follows
	 * <p>
	 * (4) pRemain = e^{-(\frac{upper limit}{E_{avg}} \cdot \Gamma(1 + \frac{1}{shape}))^{shape}}
	 * <p>
	 * from this we form a zero finding problem
	 * <p>
	 * (5) 0 = e^{-(\frac{upper limit}{E_{avg}}\cdot \Gamma(1 + \frac{1}{shape}))^{shape}} - pRemain
	 * <p>
	 * The work "Robust Scheduling of Real-Time Applications on Efficient Embedded Multicore Systems"
	 * (https://mediatum.ub.tum.de/download/1063381/1063381.pdf) proposes an algorithm from (4) which
	 * is basically the bisection method.
	 */
	public static Parameters findParametersForAverage(double lowerBound, double average, double upperBound, double pRemainPromille) {
		Parameters result = new Parameters();
		result.shape = 1;
		result.scale = 1;
		result.requested_pRemainPromille = pRemainPromille;

		result.given_lowerBound = lowerBound;
		result.requested_mean = average;
		result.given_upperBound = upperBound;

		// shifted values for standard 2 parameter Weibull distribution
		double avg = average - lowerBound;
		double max = upperBound - lowerBound;
		double pRemain = pRemainPromille / 1000;

		if ((max / avg) > 10000) {
			result.error = "finding: upper bound and mean have no valid value for parameter";
			return result;
		}

		if (pRemain > 1) {
			result.error = "invalid_argument: pRemainPromille needs to be smaller than 1";
			return result;
		}

		/* find zero interval */
		double right = 0.2;
		double left = 0.1;
		double right_f = paramZeroFunctionforAverage(right, avg, max, pRemain);
		double left_f = paramZeroFunctionforAverage(left, avg, max, pRemain);
		while ((right_f > 0 || left_f < 0) && right < 1000 && (pRemain + right_f) > 0) {
			left = right;
			right = right + Math.max((0.1 * right), 0.1);
			right_f = paramZeroFunctionforAverage(right, avg, max, pRemain);
			left_f = paramZeroFunctionforAverage(left, avg, max, pRemain);
		}
		if (right_f > 0 && left_f < 0) {
			result.error = NO_VALID_WEIBULL_DISTRIBUTION_FOUND;
			return result;
		}

		// if we found the exact solution on one of the interval borders, return it
		if (left_f == 0) {
			right = left;
			right_f = left_f;
		}

		if (right_f == 0) {
			result.shape = right;
			result.scale = computeScaleForAverage(result.shape, avg);
			return result;
		}

		while (true) {
			double half = left + (right - left) / 2;
			double half_f = paramZeroFunctionforAverage(half, avg, max, pRemain);
			// found solution
			if (half_f == 0.0 || left == right || half == right || half == left) {
				result.shape = half;
				result.scale = computeScaleForAverage(result.shape, avg);
				return result;
			} else {
				if (Math.signum(1.0) == Math.signum(half_f)) {
					left = half;
				} else {
					right = half;
				}
			}
			if (!Double.isFinite(left) || !Double.isFinite(right) || !Double.isFinite(half_f)) {
				result.error = NO_VALID_WEIBULL_DISTRIBUTION_FOUND;
				return result;
			}
		}
	}

	private static double paramZeroFunctionforAverage(double shape, double avg, double max, double pRemain) {
		return 1 / Math.exp(Math.pow((max / avg) * Gamma.gamma(1 + 1 / shape), shape)) - pRemain;
	}

	private static double computeScaleForAverage(double shape, double avg) {
		return avg / (Gamma.gamma(1 + (1 / shape)));
	}

	public static Parameters findParametersForMedian(double lowerBound, double median, double upperBound, double pRemainPromille) {
		Parameters result = new Parameters();
		result.shape = 1;
		result.scale = 1;
		result.requested_pRemainPromille = pRemainPromille;

		result.given_lowerBound = lowerBound;
		result.requested_mean = median;
		result.given_upperBound = upperBound;

		// shifted values for standard 2 parameter Weibull distribution
		double med = median - lowerBound;
		double max = upperBound - lowerBound;
		double pRemain = pRemainPromille / 1000;

		if ((max / med) > 10000) {
			result.error = "finding: upper bound and mean have no valid value for parameter";
			return result;
		}

		if (pRemain > 1) {
			result.error = "invalid_argument: pRemainPromille needs to be smaller than 1";
			return result;
		}

		/* find zero interval */
		double right = 0.2;
		double left = 0.1;
		double right_f = paramZeroFunctionForMedian(right, med, max, pRemain);
		double left_f = paramZeroFunctionForMedian(left, med, max, pRemain);
		while ((right_f > 0 || left_f < 0) && right < 1000 && (pRemain + right_f) > 0) {
			left = right;
			right = right + Math.max((0.1 * right), 0.1);
			right_f = paramZeroFunctionForMedian(right, med, max, pRemain);
			left_f = paramZeroFunctionForMedian(left, med, max, pRemain);
		}
		if (right_f > 0 && left_f < 0) {
			result.error = NO_VALID_WEIBULL_DISTRIBUTION_FOUND;
			return result;
		}

		// if we found the exact solution on one of the interval borders, return it
		if (left_f == 0) {
			right = left;
			right_f = left_f;
		}

		if (right_f == 0) {
			result.shape = right;
			result.scale = computeScaleForMedian(result.shape, med, pRemain);
			return result;
		}

		while (true) {
			double half = left + (right - left) / 2;
			double half_f = paramZeroFunctionForMedian(half, med, max, pRemain);
			// found solution
			if (half_f == 0.0 || left == right || half == right || half == left) {
				result.shape = half;
				result.scale = computeScaleForMedian(result.shape, med, pRemain);
				return result;
			} else {
				if (Math.signum(1.0) == Math.signum(half_f)) {
					left = half;
				} else {
					right = half;
				}
			}
			if (!Double.isFinite(left) || !Double.isFinite(right) || !Double.isFinite(half_f)) {
				result.error = NO_VALID_WEIBULL_DISTRIBUTION_FOUND;
				return result;
			}
		}
	}

	private static double paramZeroFunctionForMedian(double shape, double median, double max, double pRemain) {
		double factor1 = max / median;
		double factor2 = Math.pow(-1.0 * Math.log(0.5 + pRemain / 2.0), 1 / shape);
		return 1.0 / Math.exp(Math.pow(factor1 * factor2, shape)) - pRemain;
	}

	private static double computeScaleForMedian(double shape, double median, double pRemain) {
		return median / Math.pow(-1.0 * Math.log(0.5 + pRemain / 2.0), 1 / shape);
	}

}