From b55e92bef7c696fbb591364ce8358e4da309a807 Mon Sep 17 00:00:00 2001 From: TimQu Date: Mon, 4 Dec 2017 11:46:27 +0100 Subject: [PATCH] Make quick power iteration respect the relevant Values --- .../solver/NativeLinearEquationSolver.cpp | 94 +++++++++++++------ 1 file changed, 67 insertions(+), 27 deletions(-) diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 70d80e55d..a1e2c10ea 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -555,8 +555,6 @@ namespace storm { return converged; } - - template bool NativeLinearEquationSolver::solveEquationsQuickSoundPower(Environment const& env, std::vector& x, std::vector const& b) const { STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (QuickPower)"); @@ -577,53 +575,94 @@ namespace storm { std::vector* tmp = &x; - bool converged = false; - bool terminate = false; - uint64_t iterations = 0; ValueType precision = storm::utility::convertNumber(env.solver().native().getPrecision()); - ValueType lowerValueBound, upperValueBound; bool relative = env.solver().native().getRelativeTerminationCriterion(); if (!relative) { precision *= storm::utility::convertNumber(2.0); } uint64_t maxIter = env.solver().native().getMaximalNumberOfIterations(); + + //std::cout << *this->A << std::endl; + //std::cout << storm::utility::vector::toString(b) << std::endl; + + //std::cout << "solving eq sys.. " << std::endl; + + uint64_t iterations = 0; + bool converged = false; + bool terminate = false; + ValueType minValueBound, maxValueBound; + uint64_t minIndex(0), maxIndex(0); + uint64_t firstStayProb1Index = 0; + uint64_t firstIndexViolatingConvergence = this->hasRelevantValues() ? this->getRelevantValues().getNextSetIndex(0) : 0; this->startMeasureProgress(); - uint64_t firstProb1Entry = 0; while (!converged && !terminate && iterations < maxIter) { + + // Apply step this->multiplier.multAdd(*this->A, *stepBoundedValues, &b, *tmp); std::swap(tmp, stepBoundedValues); this->multiplier.multAdd(*this->A, *stepBoundedStayProbabilities, nullptr, *tmp); std::swap(tmp, stepBoundedStayProbabilities); - for (; firstProb1Entry != stepBoundedStayProbabilities->size(); ++firstProb1Entry) { + + //std::cout << "Iteration " << iterations << std::endl; + //std::cout << "x: " << storm::utility::vector::toString(*stepBoundedValues) << std::endl; + //std::cout << "y: " << storm::utility::vector::toString(*stepBoundedStayProbabilities) << std::endl; + + // Check for convergence + // Phase 1: the probability to 'stay within the matrix' has to be < 1 at every state + for (; firstStayProb1Index != stepBoundedStayProbabilities->size(); ++firstStayProb1Index) { static_assert(NumberTraits::IsExact || std::is_same::value, "Considered ValueType not handled."); if (NumberTraits::IsExact) { - if (storm::utility::isOne(stepBoundedStayProbabilities->at(firstProb1Entry))) { + if (storm::utility::isOne(stepBoundedStayProbabilities->at(firstStayProb1Index))) { break; } } else { - if (storm::utility::isAlmostOne(storm::utility::convertNumber(stepBoundedStayProbabilities->at(firstProb1Entry)))) { + if (storm::utility::isAlmostOne(storm::utility::convertNumber(stepBoundedStayProbabilities->at(firstStayProb1Index)))) { break; + // std::cout << "In Phase 1" << std::endl; } } } - if (firstProb1Entry == stepBoundedStayProbabilities->size()) { - auto valIt = stepBoundedValues->begin(); - auto valIte = stepBoundedValues->end(); - auto probIt = stepBoundedStayProbabilities->begin(); - STORM_LOG_ASSERT(!storm::utility::isOne(*probIt), "Did not expect staying-probability 1 at this point."); - lowerValueBound = *valIt / (storm::utility::one() - *probIt); - upperValueBound = lowerValueBound; - ValueType largestStayProb = *probIt; - for (; valIt != valIte; ++valIt, ++probIt) { - STORM_LOG_ASSERT(!storm::utility::isOne(*probIt), "Did not expect staying-probability 1 at this point."); - ValueType currentBound = *valIt / (storm::utility::one() - *probIt); - lowerValueBound = std::min(lowerValueBound, currentBound); - upperValueBound = std::max(upperValueBound, currentBound); - largestStayProb = std::max(largestStayProb, *probIt); + if (firstStayProb1Index == stepBoundedStayProbabilities->size()) { + STORM_LOG_ASSERT(!std::any_of(stepBoundedStayProbabilities->begin(), stepBoundedStayProbabilities->end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point."); + // Phase 2: the difference between lower and upper bound has to be < precision at every (relevant) value + // std::cout << "In Phase 2" << std::endl; + // First check with (possibly too tight) bounds from a previous iteration. Only compute the actual bounds if this first check passes. + minValueBound = stepBoundedValues->at(minIndex) / (storm::utility::one() - stepBoundedStayProbabilities->at(minIndex)); + maxValueBound = stepBoundedValues->at(maxIndex) / (storm::utility::one() - stepBoundedStayProbabilities->at(maxIndex)); + ValueType const& stayProb = stepBoundedStayProbabilities->at(firstIndexViolatingConvergence); + if (stayProb * (maxValueBound - minValueBound) < precision) { + // Compute the actual bounds now + auto valIt = stepBoundedValues->begin(); + auto valIte = stepBoundedValues->end(); + auto probIt = stepBoundedStayProbabilities->begin(); + for (uint64_t index = 0; valIt != valIte; ++valIt, ++probIt, ++index) { + ValueType currentBound = *valIt / (storm::utility::one() - *probIt); + if (currentBound < minValueBound) { + minIndex = index; + minValueBound = std::move(currentBound); + } else if (currentBound > maxValueBound) { + maxIndex = index; + maxValueBound = std::move(currentBound); + } + } + if (stayProb * (maxValueBound - minValueBound) < precision) { + // The current index satisfies the desired bound. We now move to the next index that violates it + do { + ++firstIndexViolatingConvergence; + if (this->hasRelevantValues()) { + firstIndexViolatingConvergence = this->getRelevantValues().getNextSetIndex(firstIndexViolatingConvergence); + } + if (firstIndexViolatingConvergence == stepBoundedStayProbabilities->size()) { + converged = true; + break; + } else if (stepBoundedStayProbabilities->at(firstIndexViolatingConvergence) * (maxValueBound - minValueBound) >= precision) { + // not converged yet + break; + } + } while (true); + } } STORM_LOG_ASSERT(!relative, "Relative termination criterion not implemented currently."); - converged = largestStayProb * (upperValueBound - lowerValueBound) < precision; - STORM_LOG_INFO_COND(!converged, "Lower value bound: " << lowerValueBound << " Upper value bound: " << upperValueBound << " Largest stay prob.: " << largestStayProb); } // Potentially show progress. @@ -635,7 +674,7 @@ namespace storm { } // Finally set up the solution vector - ValueType meanBound = (upperValueBound - lowerValueBound) / storm::utility::convertNumber(2.0); + ValueType meanBound = (maxValueBound + minValueBound) / storm::utility::convertNumber(2.0); storm::utility::vector::applyPointwise(*stepBoundedValues, *stepBoundedStayProbabilities, x, [&meanBound] (ValueType const& v, ValueType const& p) { return v + p * meanBound; }); if (!this->isCachingEnabled()) { @@ -643,6 +682,7 @@ namespace storm { } this->logIterations(converged, terminate, iterations); + STORM_LOG_INFO("Quick Power Iteration terminated with lower value bound " << minValueBound << " and upper value bound " << maxValueBound << "."); return converged;