From c70b6baf8194a9974547d66b4b708428b0a9fa26 Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Thu, 12 Mar 2020 09:49:15 +0100 Subject: [PATCH] Abort unif+ also in inner iterations. Store the best known solution after each completed step. --- .../helper/SparseMarkovAutomatonCslHelper.cpp | 65 ++++++++++++++----- 1 file changed, 48 insertions(+), 17 deletions(-) diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp index 438181fe9..6647d83a6 100644 --- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp @@ -76,24 +76,31 @@ namespace storm { storm::storage::BitVector probabilisticMaybeStates = ~markovianStates & maybeStates; storm::storage::BitVector markovianStatesModMaybeStates = markovianMaybeStates % maybeStates; storm::storage::BitVector probabilisticStatesModMaybeStates = probabilisticMaybeStates % maybeStates; + // Catch the case where this query can be solved by solving the untimed variant instead. + // This is the case if there is no Markovian maybe state (e.g. if the initial state is already a psi state) of if the time bound is infinity. + if (markovianMaybeStates.empty() || storm::utility::isInfinity(upperTimeBound)) { + return SparseMarkovAutomatonCslHelper::computeUntilProbabilities(env, dir, transitionMatrix, transitionMatrix.transpose(true), phiStates, psiStates, false, false).values; + } + boost::optional relevantMaybeStates; if (relevantStates) { relevantMaybeStates = relevantStates.get() % maybeStates; } - // Catch the case where this is query can be solved by solving the untimed variant instead. - // This is the case if there is no Markovian maybe state (e.g. if the initial state is already a psi state) of if the time bound is infinity. - if (markovianMaybeStates.empty() || storm::utility::isInfinity(upperTimeBound)) { - return SparseMarkovAutomatonCslHelper::computeUntilProbabilities(env, dir, transitionMatrix, transitionMatrix.transpose(true), phiStates, psiStates, false, false).values; + // Store the best solution known so far (useful in cases where the computation gets aborted) + std::vector bestKnownSolution; + if (relevantMaybeStates) { + bestKnownSolution.resize(relevantStates->size()); } // Get the exit rates restricted to only markovian maybe states. std::vector markovianExitRates = storm::utility::vector::filterVector(exitRateVector, markovianMaybeStates); // Obtain parameters of the algorithm + auto two = storm::utility::convertNumber(2.0); // Truncation error ValueType kappa = storm::utility::convertNumber(env.solver().timeBounded().getUnifPlusKappa()); // Precision to be achieved - ValueType epsilon = storm::utility::convertNumber(2.0) * storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); + ValueType epsilon = two * storm::utility::convertNumber(env.solver().timeBounded().getPrecision()); bool relativePrecision = env.solver().timeBounded().getRelativeTerminationCriterion(); // Uniformization rate ValueType lambda = *std::max_element(markovianExitRates.begin(), markovianExitRates.end()); @@ -132,7 +139,7 @@ namespace storm { uint64_t iteration = 0; progressIterations.startNewMeasurement(iteration); bool converged = false; - + bool abortedInnerIterations = false; while (!converged) { // Maximal step size uint64_t N = storm::utility::ceil(lambda * upperTimeBound * std::exp(2) - storm::utility::log(kappa * epsilon)); @@ -155,8 +162,9 @@ namespace storm { progressSteps.setMaxCount(N); progressSteps.startNewMeasurement(0); bool firstIteration = true; // The first iterations can be irrelevant, because they will only produce zeroes anyway. + int64_t k = N; // Iteration k = N is always non-relevant - for (int64_t k = N - 1; k >= 0; --k) { + for (--k; k >= 0; --k) { // Check whether the iteration is relevant, that is, whether it will contribute non-zero values to the overall result if (computeLowerBound) { @@ -222,20 +230,36 @@ namespace storm { } progressSteps.updateProgress(N-k); + if (storm::utility::resources::isTerminate()) { + abortedInnerIterations = true; + break; + } } + if (computeLowerBound) { storm::utility::vector::scaleVectorInPlace(maybeStatesValuesLower, storm::utility::one() / foxGlynnResult.totalWeight); } else { storm::utility::vector::scaleVectorInPlace(maybeStatesValuesUpper, storm::utility::one() / foxGlynnResult.totalWeight); } - + + if (abortedInnerIterations || storm::utility::resources::isTerminate()) { + break; + } + // Check if the lower and upper bound are sufficiently close to each other converged = checkConvergence(maybeStatesValuesLower, maybeStatesValuesUpper, relevantMaybeStates, epsilon, relativePrecision, kappa); if (converged) { break; } - if (storm::utility::resources::isTerminate()) { - break; + + // Store the best solution we have found so far. + if (relevantMaybeStates) { + auto currentSolIt = bestKnownSolution.begin(); + for (auto const& state : relevantMaybeStates.get()) { + // We take the average of the lower and upper bounds + *currentSolIt = (maybeStatesValuesLower[state] + maybeStatesValuesUpper[state]) / two; + ++currentSolIt; + } } } @@ -244,7 +268,7 @@ namespace storm { // Double lambda. ValueType oldLambda = lambda; - lambda *= storm::utility::convertNumber(2.0); + lambda *= two; STORM_LOG_DEBUG("Increased lambda to " << lambda << "."); if (relativePrecision) { @@ -268,17 +292,24 @@ namespace storm { } progressIterations.updateProgress(++iteration); if (storm::utility::resources::isTerminate()) { + STORM_LOG_WARN("Aborted unif+ in iteration " << iteration << "."); break; } } - - // We take the average of the lower and upper bounds - auto two = storm::utility::convertNumber(2.0); - storm::utility::vector::applyPointwise(maybeStatesValuesLower, maybeStatesValuesUpper, maybeStatesValuesLower, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); - + + // Prepare the result vector std::vector result(transitionMatrix.getRowGroupCount(), storm::utility::zero()); storm::utility::vector::setVectorValues(result, psiStates, storm::utility::one()); - storm::utility::vector::setVectorValues(result, maybeStates, maybeStatesValuesLower); + + if (abortedInnerIterations && iteration > 1 && relevantMaybeStates && relevantStates) { + // We should take the stored solution instead of the current (probably more incorrect) lower/upper values + storm::utility::vector::setVectorValues(result, maybeStates & relevantStates.get(), bestKnownSolution); + } else { + // We take the average of the lower and upper bounds + storm::utility::vector::applyPointwise(maybeStatesValuesLower, maybeStatesValuesUpper, maybeStatesValuesLower, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); + + storm::utility::vector::setVectorValues(result, maybeStates, maybeStatesValuesLower); + } return result; }