(interval iteration) include/Measures.h: provide helper functions for interval iteration

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

prism/include/Measures.h

@@ -103,4 +103,86 @@ public:
 
 };
 
+/**
+ * Measure for determining the difference between the upper and lower values in
+ * an interval iteration.
+ * In relative mode, the difference is scaled by the mid-point between upper and lower value.
+ */
+class MeasureSupNormInterval {
+private:
+	// relative mode?
+	const bool relative;
+	// the current maximal value
+	double sup_norm;
+
+public:
+	/** Constructor, set relative flag */
+	MeasureSupNormInterval(bool relative) : relative(relative) {
+		reset();
+	}
+
+	/** Reset for new measurement */
+	void reset() {
+		sup_norm = 0.0;
+	}
+
+	/** Relative mode? */
+	bool isRelative() const {
+		return relative;
+	}
+
+	/**
+	 * Do the measurement for an upper and lower value pair.
+	 * For relative mode, the lower value is used as the divisor.
+	 */
+	inline void measure(double lower, double upper) {
+		double x;
+
+		// special case: one of the values is infinite (this can happen e.g. for non-converging
+		// iterations when the values grow extremely large and overflow to infinity).
+		if (isinf(lower) || isinf(upper)) {
+			x = std::numeric_limits<double>::infinity();
+		} else {
+			// compute difference
+			// we don't use fabs like for MeasureSupNorm, as we want x to become negative
+			// in situations where upper < lower (should only happen due to numerical inaccuracies / rounding)
+			x = upper - lower;
+			// we clamp to zero for negative values
+			if (x < 0)
+				x = 0;
+
+			if (relative && x != 0.0) {
+				// for relative mode: divide by lower
+				// taking lower ensures that if the actual value should happen to be
+				// the lower value, that then the relative precision is satisfied.
+				// We take the absolute value of the lower to ensure that
+				// x does not flip signs.
+				// Note: if lower is 0.0, then x will become +inf, as x!=0
+				x /= fabs(lower);
+			}
+		}
+
+		// sup_norm = max { x, sup_norm }
+		if (x > sup_norm) {
+			sup_norm = x;
+		}
+	}
+
+	/**
+	 * Do the measurement for a pair of arrays (lower and upper bounds) of size n.
+	 * For relative mode, the midpoints are used as the divisors.
+	 */
+	inline void measure(const double *lower, const double *upper, std::size_t n) {
+		for (std::size_t i = 0; i < n; i++) {
+			measure(lower[i], upper[i]);
+		}
+	}
+
+	/** Return the measured value */
+	double value() const {
+		return sup_norm;
+	}
+
+};
+
 #endif