From 99db5693caee957e56284ebba6dcbfa2d3a17c02 Mon Sep 17 00:00:00 2001
From: Jan Erik Karuc <jan.karuc@rwth-aachen.de>
Date: Thu, 27 Feb 2020 21:33:38 +0100
Subject: [PATCH] OVI: Implement upper bound only iterations

---
 .../IterativeMinMaxLinearEquationSolver.cpp   | 34 +++++++++++++------
 1 file changed, 24 insertions(+), 10 deletions(-)

diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
index 5b41909a1..1925390f2 100644
--- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
+++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
@@ -446,6 +446,8 @@ namespace storm {
             bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel;
             // Relative errors
             bool relative = env.solver().minMax().getRelativeTerminationCriterion();
+            // Upper bound only iterations
+            uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations();
 
             boost::optional<storm::storage::BitVector> relevantValues;
             if (this->hasRelevantValues()) {
@@ -480,6 +482,7 @@ namespace storm {
                 if (result.status != SolverStatus::Converged) {
                     status = result.status;
                 } else {
+                    bool intervalIterationNeeded = false;
                     currentVerificationIterations = 0;
 
                     if (relative) {
@@ -497,31 +500,39 @@ namespace storm {
                         if (useGaussSeidelMultiplication) {
                             // Copy over the current vectors so we can modify them in-place.
                             // This is necessary as we want to compare the new values with the current ones.
-                            *newX = *currentX;
                             newUpperBound = currentUpperBound;
-                            // Do the calculation
-                            multiplier.multiplyAndReduceGaussSeidel(env, dir, *newX, &b);
+                            // Do the calculation.
                             multiplier.multiplyAndReduceGaussSeidel(env, dir, newUpperBound, &b);
+                            if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) {
+                                // Now do interval iteration.
+                                *newX = *currentX;
+                                multiplier.multiplyAndReduceGaussSeidel(env, dir, *newX, &b);
+                            }
                         } else {
-                            multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX);
                             multiplier.multiplyAndReduce(env, dir, currentUpperBound, &b, newUpperBound);
+                            if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) {
+                                // Now do interval iteration.
+                                multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX);
+                            }
                         }
 
                         bool newUpperBoundAlwaysHigherEqual = true;
                         bool newUpperBoundAlwaysLowerEqual = true;
                         bool valuesCrossed = false;
-                        ValueType maxBoundDiff = storm::utility::zero<ValueType>();
                         for (uint64_t i = 0; i < x.size(); ++i) {
                             if (newUpperBound[i] < currentUpperBound[i]) {
                                 newUpperBoundAlwaysHigherEqual = false;
                             } else if (newUpperBound[i] != currentUpperBound[i]) {
                                 newUpperBoundAlwaysLowerEqual = false;
                             }
-                            if (newUpperBound[i] < (*newX)[i]) {
-                                valuesCrossed = true;
-                                break;
-                            } else {
-                                maxBoundDiff = std::max<ValueType>(maxBoundDiff, newUpperBound[i] - (*newX)[i]);
+                        }
+
+                        if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) {
+                            for (uint64_t i = 0; i < x.size(); ++i) {
+                                if (newUpperBound[i] < (*newX)[i]) {
+                                    valuesCrossed = true;
+                                    break;
+                                }
                             }
                         }
                         
@@ -569,6 +580,9 @@ namespace storm {
                                 storm::utility::vector::applyPointwise<ValueType, ValueType, ValueType>(*currentX, currentUpperBound, *currentX, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; });
                                 status = SolverStatus::Converged;
                             }
+                            else {
+                                intervalIterationNeeded = true;
+                            }
                         }
                         else {
                             ++no_case_occur;