plugins/org.eclipse.app4mc.amalthea.model/src/org/eclipse/app4mc/amalthea/model/util/WeibullUtil.java - app4mc/org.eclipse.app4mc - Git at Google

 /**
  ********************************************************************************
  * 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);
 	}

 }
	/**
	********************************************************************************
	* 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);
	}

	}