(interval iteration) PrismUtils: add helper methods for interval iteration

git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@12135 bbc10eb1-c90d-0410-af57-cb519fbb1720

prism/src/prism/PrismUtils.java

@@ -205,6 +205,73 @@ public class PrismUtils
 		return (Math.abs(d1 - d2) / d1);
 	}
 
+	/**
+	 * Measure supremum norm, for all the entries given by the {@code indizes}
+	 * iterator, for an interval iteration.
+	 */
+	public static double measureSupNormInterval(double[] lower, double[] upper, boolean abs, PrimitiveIterator.OfInt indizes)
+	{
+		int n = lower.length;
+		assert( n== upper.length);
+
+		double value = 0;
+		while (indizes.hasNext()) {
+			int i = indizes.nextInt();
+			double diff = measureSupNormInterval(lower[i], upper[i], abs);
+			if (diff > value)
+				value = diff;
+		}
+		return value;
+	}
+
+	/**
+	 * Measure supremum norm, either absolute or relative, for an interval iteration,
+	 * return the maximum difference.
+	 */
+	public static double measureSupNormInterval(double[] lower, double[] upper, boolean abs)
+	{
+		int n = lower.length;
+		assert( n== upper.length);
+
+		double value = 0;
+		for (int i=0; i < n; i++) {
+			double diff = measureSupNormInterval(lower[i], upper[i], abs);
+			if (diff > value)
+				value = diff;
+		}
+		return value;
+	}
+
+	/**
+	 * Measure supremum norm for two values, for interval iteration.
+	 */
+	public static double measureSupNormInterval(double lower, double upper, boolean abs)
+	{
+		// Deal with infinite cases
+		if (Double.isInfinite(lower)) {
+			if (Double.isInfinite(upper) && (lower > 0) == (upper > 0)) {
+				return 0;
+			} else {
+				return Double.POSITIVE_INFINITY;
+			}
+		} else if (Double.isInfinite(upper)) {
+			return Double.POSITIVE_INFINITY;
+		}
+
+		if (lower == upper)
+			return 0.0;
+
+		// Compute/check error
+		lower = Math.abs(lower);
+		upper = Math.abs(upper);
+		double result = upper - lower;
+		result = Math.abs(result);
+		if (!abs) {
+			result = result / lower;
+		}
+		return result;
+	}
+
 	/**
 	 * See if two doubles are (nearly) equal.
 	 */
@@ -231,6 +298,100 @@ public class PrismUtils
 		return max_v;
 	}
 
+	/**
+	 * Ensure monotonicity from below for interval iteration solution vectors.
+	 * Compares the old and new values and overwrites the new value with the old
+	 * value if that is larger.
+	 * @param old_values old solution vector
+	 * @param new_values new solution vector
+	 */
+	public static void ensureMonotonicityFromBelow(double[] old_values, double[] new_values)
+	{
+		assert(old_values.length == new_values.length);
+
+		for (int i = 0, n = old_values.length; i < n; i++) {
+			double old_value = old_values[i];
+			double new_value = new_values[i];
+			// from below: do max
+			if (old_value > new_value) {
+				new_values[i] = old_value;
+			}
+		}
+	}
+
+	/**
+	 * Ensure monotonicity from above for interval iteration solution vectors.
+	 * Compares the old and new values and overwrites the new value with the old
+	 * value if that is smaller.
+	 * @param old_values old solution vector
+	 * @param new_values new solution vector
+	 */
+	public static void ensureMonotonicityFromAbove(double[] old_values, double[] new_values)
+	{
+		assert(old_values.length == new_values.length);
+
+		for (int i = 0, n = old_values.length; i < n; i++) {
+			double old_value = old_values[i];
+			double new_value = new_values[i];
+			// from above: do min
+			if (old_value < new_value) {
+				new_values[i] = old_value;
+			}
+		}
+	}
+
+	/**
+	 * Check for monotonicity: If the new_values are not element-wise less-than-equal the older values
+	 * (for from_above == true), then throws an exception. If from_above == false, the logic is reversed,
+	 * i.e., if the new_value is not greater-than-equal, an exception is thrown.
+	 * @param old_values the old values
+	 * @param new_values the new values
+	 * @param from_above the direction
+	 */
+	public static void checkMonotonicity(double[] old_values, double[] new_values, boolean from_above) throws PrismException
+	{
+		assert(old_values.length == new_values.length);
+
+		for (int i = 0, n = old_values.length; i < n; i++) {
+			double old_value = old_values[i];
+			double new_value = new_values[i];
+			if (from_above && old_value < new_value) {
+				throw new PrismException("Monotonicity violated (from above): old value " + old_value + " < new value " + new_value);
+			}
+			if (!from_above && old_value > new_value) {
+				throw new PrismException("Monotonicity violated (from below): old value " + old_value + " > new value " + new_value);
+			}
+		}
+	}
+
+	/**
+	 * Select midpoint from two interval iteration solution vectors.
+	 * Stores the result in soln_below.
+	 * @param soln_below solution vector from below
+	 * @param soln_above solution vector from above
+	 */
+	public static void selectMidpoint(double[] soln_below, double[] soln_above)
+	{
+		assert(soln_below.length == soln_above.length);
+
+		for (int i = 0, n = soln_below.length; i < n; i++) {
+			double below = soln_below[i];
+			double above = soln_above[i];
+
+			if (below != above) {
+				// use below + ( above - below ) / 2 instead of (below+above)/2 for better numerical
+				// stability
+				double d = below + (above - below)/2.0;
+				if (d >= below && d <= above) {
+					// only store result if between below and above
+					// guard against rounding problems,
+					// fallback is to simply return below as is
+					soln_below[i] = d;
+				}
+			}
+		}
+	}
+
 	/**
 	 * Format a large integer, represented by a double, as a string. 
 	 */