From a24de86ce1ee1ce107de01a248f348e243d86046 Mon Sep 17 00:00:00 2001 From: TimQu Date: Tue, 13 Mar 2018 13:41:12 +0100 Subject: [PATCH] Avoided duplicated code for sound value iteration --- .../IterativeMinMaxLinearEquationSolver.cpp | 406 +---------------- .../solver/NativeLinearEquationSolver.cpp | 213 +-------- .../helper/SoundValueIterationHelper.cpp | 418 ++++++++++++++++++ .../solver/helper/SoundValueIterationHelper.h | 147 ++++++ 4 files changed, 572 insertions(+), 612 deletions(-) create mode 100644 src/storm/solver/helper/SoundValueIterationHelper.cpp create mode 100644 src/storm/solver/helper/SoundValueIterationHelper.h diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index decf30b34..4bb4cc766 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -10,6 +10,7 @@ #include "storm/utility/KwekMehlhorn.h" #include "storm/utility/NumberTraits.h" +#include "storm/solver/helper/SoundValueIterationHelper.h" #include "storm/utility/Stopwatch.h" #include "storm/utility/vector.h" #include "storm/utility/macros.h" @@ -604,395 +605,6 @@ namespace storm { return status == SolverStatus::Converged; } - template - class SoundValueIterationHelper { - public: - - typedef uint32_t IndexType; - - SoundValueIterationHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision), rowGroupIndices(matrix.getRowGroupIndices()) { - STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type."); - uint64_t sizeOfLargestRowGroup = matrix.getSizeOfLargestRowGroup(); - xTmp.resize(sizeOfLargestRowGroup); - yTmp.resize(sizeOfLargestRowGroup); - x.assign(x.size(), storm::utility::zero()); - y.assign(x.size(), storm::utility::one()); - hasDecisionValue = false; - decisionValueBlocks = false; - convergencePhase1 = true; - firstIndexViolatingConvergence = 0; - - numRows = matrix.getRowCount(); - matrixValues.clear(); - matrixColumns.clear(); - rowIndications.clear(); - matrixValues.reserve(matrix.getNonzeroEntryCount()); - matrixColumns.reserve(matrix.getColumnCount()); - rowIndications.reserve(numRows + 1); - rowIndications.push_back(0); - for (IndexType r = 0; r < numRows; ++r) { - for (auto const& entry : matrix.getRow(r)) { - matrixValues.push_back(entry.getValue()); - matrixColumns.push_back(entry.getColumn()); - } - rowIndications.push_back(matrixValues.size()); - } - } - - - inline void setLowerBound(ValueType const& value) { - hasLowerBound = true; - lowerBound = value; - } - - inline void setUpperBound(ValueType const& value) { - hasUpperBound = true; - upperBound = value; - } - - template - inline bool better(ValueType const& val1, ValueType const& val2) { - return maximize(dir) ? val1 > val2 : val1 < val2; - } - - template - inline ValueType& getPrimaryBound() { - return maximize(dir) ? upperBound : lowerBound; - } - - template - inline bool& hasPrimaryBound() { - return maximize(dir) ? hasUpperBound : hasLowerBound; - } - - template - inline ValueType& getSecondaryBound() { - return maximize(dir) ? lowerBound : upperBound; - } - - template - inline uint64_t& getPrimaryIndex() { - return maximize(dir) ? maxIndex : minIndex; - } - - template - inline uint64_t& getSecondaryIndex() { - return maximize(dir) ? minIndex : maxIndex; - } - - void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) { - assert(rowIndex < numRows); - ValueType xRes = bi; - ValueType yRes = storm::utility::zero(); - - auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; - auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; - auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; - for (; entryIt != entryItE; ++entryIt, ++colIt) { - xRes += *entryIt * x[*colIt]; - yRes += *entryIt * y[*colIt]; - } - xi = std::move(xRes); - yi = std::move(yRes); - } - - template - void performIterationStep(std::vector const& b) { - if (!decisionValueBlocks) { - performIterationStepUpdateDecisionValue(b); - } else { - assert(decisionValue == getPrimaryBound()); - auto xIt = x.rbegin(); - auto yIt = y.rbegin(); - auto groupStartIt = rowGroupIndices.rbegin(); - uint64_t groupEnd = *groupStartIt; - ++groupStartIt; - for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { - // Perform the iteration for the first row in the group - IndexType row = *groupStartIt; - ValueType xBest, yBest; - multiplyRow(row, b[row], xBest, yBest); - ++row; - // Only do more work if there are still rows in this row group - if (row != groupEnd) { - ValueType xi, yi; - ValueType bestValue = xBest + yBest * getPrimaryBound(); - for (;row < groupEnd; ++row) { - // Get the multiplication results - multiplyRow(row, b[row], xi, yi); - ValueType currentValue = xi + yi * getPrimaryBound(); - // Check if the current row is better then the previously found one - if (better(currentValue, bestValue)) { - xBest = std::move(xi); - yBest = std::move(yi); - bestValue = std::move(currentValue); - } else if (currentValue == bestValue && yBest > yi) { - // If the value for this row is not strictly better, it might still be equal and have a better y value - xBest = std::move(xi); - yBest = std::move(yi); - } - } - } - *xIt = std::move(xBest); - *yIt = std::move(yBest); - } - } - } - - template - void performIterationStepUpdateDecisionValue(std::vector const& b) { - auto xIt = x.rbegin(); - auto yIt = y.rbegin(); - auto groupStartIt = rowGroupIndices.rbegin(); - uint64_t groupEnd = *groupStartIt; - ++groupStartIt; - for (auto groupStartIte = rowGroupIndices.rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { - // Perform the iteration for the first row in the group - uint64_t row = *groupStartIt; - ValueType xBest, yBest; - multiplyRow(row, b[row], xBest, yBest); - ++row; - // Only do more work if there are still rows in this row group - if (row != groupEnd) { - ValueType xi, yi; - uint64_t xyTmpIndex = 0; - if (hasPrimaryBound()) { - ValueType bestValue = xBest + yBest * getPrimaryBound(); - for (;row < groupEnd; ++row) { - // Get the multiplication results - multiplyRow(row, b[row], xi, yi); - ValueType currentValue = xi + yi * getPrimaryBound(); - // Check if the current row is better then the previously found one - if (better(currentValue, bestValue)) { - if (yBest < yi) { - // We need to store the 'old' best value as it might be relevant for the decision value - xTmp[xyTmpIndex] = std::move(xBest); - yTmp[xyTmpIndex] = std::move(yBest); - ++xyTmpIndex; - } - xBest = std::move(xi); - yBest = std::move(yi); - bestValue = std::move(currentValue); - } else if (yBest > yi) { - // If the value for this row is not strictly better, it might still be equal and have a better y value - if (currentValue == bestValue) { - xBest = std::move(xi); - yBest = std::move(yi); - } else { - xTmp[xyTmpIndex] = std::move(xi); - yTmp[xyTmpIndex] = std::move(yi); - ++xyTmpIndex; - } - } - } - } else { - for (;row < groupEnd; ++row) { - multiplyRow(row, b[row], xi, yi); - // Update the best choice - if (yi > yBest || (yi == yBest && better(xi, xBest))) { - xTmp[xyTmpIndex] = std::move(xBest); - yTmp[xyTmpIndex] = std::move(yBest); - ++xyTmpIndex; - xBest = std::move(xi); - yBest = std::move(yi); - } else { - xTmp[xyTmpIndex] = std::move(xi); - yTmp[xyTmpIndex] = std::move(yi); - ++xyTmpIndex; - } - } - } - - // Update the decision value - for (uint64_t i = 0; i < xyTmpIndex; ++i) { - ValueType deltaY = yBest - yTmp[i]; - if (deltaY > storm::utility::zero()) { - ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY; - if (!hasDecisionValue || better(newDecisionValue, decisionValue)) { - decisionValue = std::move(newDecisionValue); - hasDecisionValue = true; - } - } - } - } - *xIt = std::move(xBest); - *yIt = std::move(yBest); - } - } - - template - bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr) { - - if (convergencePhase1) { - if (checkConvergencePhase1()) { - firstIndexViolatingConvergence = 0; - if (relevantValues != nullptr) { - firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); - } - } else { - return false; - } - } - STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point."); - - // Reaching this point means that we are in Phase 2: - // The difference between lower and upper bound has to be < precision at every (relevant) value - - // For efficiency reasons we first check whether it is worth to compute the actual bounds. We do so by considering possibly too tight bounds - ValueType lowerBoundCandidate, upperBoundCandidate; - if (preliminaryConvergenceCheck(lowerBoundCandidate, upperBoundCandidate)) { - updateLowerUpperBound(lowerBoundCandidate, upperBoundCandidate); - checkIfDecisionValueBlocks(); - return checkConvergencePhase2(relevantValues); - } - return false; - } - - void setSolutionVector() { - STORM_LOG_WARN_COND(hasLowerBound && hasUpperBound, "No lower or upper result bound could be computed within the given number of Iterations."); - - ValueType meanBound = (upperBound + lowerBound) / storm::utility::convertNumber(2.0); - storm::utility::vector::applyPointwise(x, y, x, [&meanBound] (ValueType const& xi, ValueType const& yi) { return xi + yi * meanBound; }); - - STORM_LOG_INFO("Sound Value Iteration terminated with lower value bound " - << (hasLowerBound ? lowerBound : storm::utility::zero()) << (hasLowerBound ? "" : "(none)") - << " and upper value bound " - << (hasUpperBound ? upperBound : storm::utility::zero()) << (hasUpperBound ? "" : "(none)") - << ". Decision value is " - << (hasDecisionValue ? decisionValue : storm::utility::zero()) << (hasDecisionValue ? "" : "(none)") - << "."); - - } - - private: - - bool checkConvergencePhase1() { - // Return true if y ('the probability to stay within the matrix') is < 1 at every entry - for (; firstIndexViolatingConvergence != y.size(); ++firstIndexViolatingConvergence) { - static_assert(NumberTraits::IsExact || std::is_same::value, "Considered ValueType not handled."); - if (NumberTraits::IsExact) { - if (storm::utility::isOne(y[firstIndexViolatingConvergence])) { - return false; - } - } else { - if (storm::utility::isAlmostOne(storm::utility::convertNumber(y[firstIndexViolatingConvergence]))) { - return false; - } - } - } - convergencePhase1 = false; - return true; - } - - - bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub) { - return yi * (ub - lb) <= storm::utility::abs((relative ? (precision * xi) : (precision * storm::utility::convertNumber(2.0)))); - } - - template - bool preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { - lowerBoundCandidate = x[minIndex] / (storm::utility::one() - y[minIndex]); - upperBoundCandidate = x[maxIndex] / (storm::utility::one() - y[maxIndex]); - // Make sure that these candidates are at least as tight as the already known bounds - if (hasLowerBound && lowerBoundCandidate < lowerBound) { - lowerBoundCandidate = lowerBound; - } - if (hasUpperBound && upperBoundCandidate > upperBound) { - upperBoundCandidate = upperBound; - } - if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) { - return true; - } - if (!decisionValueBlocks) { - return hasDecisionValue && better(decisionValue, getPrimaryBound()); - } - return false; - } - - template - void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { - auto xIt = x.begin(); - auto xIte = x.end(); - auto yIt = y.begin(); - for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) { - ValueType currentBound = *xIt / (storm::utility::one() - *yIt); - if (decisionValueBlocks) { - if (better(getSecondaryBound(), currentBound)) { - getSecondaryIndex() = index; - getSecondaryBound() = std::move(currentBound); - } - } else { - if (currentBound < lowerBoundCandidate) { - minIndex = index; - lowerBoundCandidate = std::move(currentBound); - } else if (currentBound > upperBoundCandidate) { - maxIndex = index; - upperBoundCandidate = std::move(currentBound); - } - } - } - if ((maximize(dir) || !decisionValueBlocks) && (!hasLowerBound || lowerBoundCandidate > lowerBound)) { - setLowerBound(lowerBoundCandidate); - } - if ((minimize(dir) || !decisionValueBlocks) && (!hasUpperBound || upperBoundCandidate < upperBound)) { - setUpperBound(upperBoundCandidate); - } - } - - template - void checkIfDecisionValueBlocks() { - // Check whether the decision value blocks now (i.e. further improvement of the primary bound would lead to a non-optimal scheduler). - if (!decisionValueBlocks && hasDecisionValue && better(decisionValue, getPrimaryBound())) { - getPrimaryBound() = decisionValue; - decisionValueBlocks = true; - } - } - - template - bool checkConvergencePhase2(storm::storage::BitVector const* relevantValues = nullptr) { - // Check whether the desired precision is reached - if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { - // The current index satisfies the desired bound. We now move to the next index that violates it - while (true) { - ++firstIndexViolatingConvergence; - if (relevantValues != nullptr) { - firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); - } - if (firstIndexViolatingConvergence == x.size()) { - // Converged! - return true; - } else { - if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { - // not converged yet - return false; - } - } - } - } - return false; - } - - std::vector& x; - std::vector& y; - std::vector xTmp, yTmp; - - ValueType lowerBound, upperBound, decisionValue; - bool hasLowerBound, hasUpperBound, hasDecisionValue; - uint64_t minIndex, maxIndex; - bool convergencePhase1; - bool decisionValueBlocks; - uint64_t firstIndexViolatingConvergence; - - bool relative; - ValueType precision; - - std::vector matrixValues; - std::vector matrixColumns; - std::vector rowIndications; - std::vector::index_type> const& rowGroupIndices; - IndexType numRows; - }; - template bool IterativeMinMaxLinearEquationSolver::solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { @@ -1003,7 +615,7 @@ namespace storm { } // TODO: implement caching for the helper - SoundValueIterationHelper helper(*this->A, x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision())); + storm::solver::helper::SoundValueIterationHelper helper(*this->A, x, *this->auxiliaryRowGroupVector, env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision())); // Prepare initial bounds for the solution (if given) if (this->hasLowerBound()) { @@ -1023,17 +635,9 @@ namespace storm { uint64_t iterations = 0; while (status == SolverStatus::InProgress && iterations < env.solver().minMax().getMaximalNumberOfIterations()) { - if (minimize(dir)) { - helper.template performIterationStep(b); - if (helper.template checkConvergenceUpdateBounds(relevantValuesPtr)) { - status = SolverStatus::Converged; - } - } else { - assert(maximize(dir)); - helper.template performIterationStep(b); - if (helper.template checkConvergenceUpdateBounds(relevantValuesPtr)) { - status = SolverStatus::Converged; - } + helper.performIterationStep(dir, b); + if (helper.checkConvergenceUpdateBounds(dir, relevantValuesPtr)) { + status = SolverStatus::Converged; } // Update environment variables. diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 897b6522a..dd67522a1 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -9,6 +9,7 @@ #include "storm/utility/NumberTraits.h" #include "storm/utility/constants.h" #include "storm/utility/vector.h" +#include "storm/solver/helper/SoundValueIterationHelper.h" #include "storm/solver/Multiplier.h" #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidEnvironmentException.h" @@ -563,216 +564,6 @@ namespace storm { return converged; } - template - class SoundPowerHelper { - public: - - typedef uint32_t IndexType; - - SoundPowerHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), minIndex(0), maxIndex(0), relative(relative), precision(precision) { - STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type."); - x.assign(x.size(), storm::utility::zero()); - y.assign(x.size(), storm::utility::one()); - convergencePhase1 = true; - firstIndexViolatingConvergence = 0; - - numRows = matrix.getRowCount(); - matrixValues.clear(); - matrixColumns.clear(); - rowIndications.clear(); - matrixValues.reserve(matrix.getNonzeroEntryCount()); - matrixColumns.reserve(matrix.getColumnCount()); - rowIndications.reserve(numRows + 1); - rowIndications.push_back(0); - for (IndexType r = 0; r < numRows; ++r) { - for (auto const& entry : matrix.getRow(r)) { - matrixValues.push_back(entry.getValue()); - matrixColumns.push_back(entry.getColumn()); - } - rowIndications.push_back(matrixValues.size()); - } - } - - inline void setLowerBound(ValueType const& value) { - hasLowerBound = true; - lowerBound = value; - } - - inline void setUpperBound(ValueType const& value) { - hasUpperBound = true; - upperBound = value; - } - - void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) { - assert(rowIndex < numRows); - ValueType xRes = bi; - ValueType yRes = storm::utility::zero(); - - auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; - auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; - auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; - for (; entryIt != entryItE; ++entryIt, ++colIt) { - xRes += *entryIt * x[*colIt]; - yRes += *entryIt * y[*colIt]; - } - xi = std::move(xRes); - yi = std::move(yRes); - } - - void performIterationStep(std::vector const& b) { - auto xIt = x.rbegin(); - auto yIt = y.rbegin(); - IndexType row = numRows; - while (row > 0) { - --row; - multiplyRow(row, b[row], *xIt, *yIt); - ++xIt; - ++yIt; - } - } - - bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr) { - - if (convergencePhase1) { - if (checkConvergencePhase1()) { - firstIndexViolatingConvergence = 0; - if (relevantValues != nullptr) { - firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); - } - } else { - return false; - } - } - STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point."); - - // Reaching this point means that we are in Phase 2: - // The difference between lower and upper bound has to be < precision at every (relevant) value - - // For efficiency reasons we first check whether it is worth to compute the actual bounds. We do so by considering possibly too tight bounds - ValueType lowerBoundCandidate, upperBoundCandidate; - if (preliminaryConvergenceCheck(lowerBoundCandidate, upperBoundCandidate)) { - updateLowerUpperBound(lowerBoundCandidate, upperBoundCandidate); - return checkConvergencePhase2(relevantValues); - } - return false; - } - - void setSolutionVector() { - STORM_LOG_WARN_COND(hasLowerBound && hasUpperBound, "No lower or upper result bound could be computed within the given number of Iterations."); - - ValueType meanBound = (upperBound + lowerBound) / storm::utility::convertNumber(2.0); - storm::utility::vector::applyPointwise(x, y, x, [&meanBound] (ValueType const& xi, ValueType const& yi) { return xi + yi * meanBound; }); - - STORM_LOG_INFO("Sound Power Iteration terminated with lower value bound " - << (hasLowerBound ? lowerBound : storm::utility::zero()) << (hasLowerBound ? "" : "(none)") - << " and upper value bound " - << (hasUpperBound ? upperBound : storm::utility::zero()) << (hasUpperBound ? "" : "(none)") - << "."); - } - - private: - - bool checkConvergencePhase1() { - // Return true if y ('the probability to stay within the 'maybestates') is < 1 at every entry - for (; firstIndexViolatingConvergence != y.size(); ++firstIndexViolatingConvergence) { - static_assert(NumberTraits::IsExact || std::is_same::value, "Considered ValueType not handled."); - if (NumberTraits::IsExact) { - if (storm::utility::isOne(y[firstIndexViolatingConvergence])) { - return false; - } - } else { - if (storm::utility::isAlmostOne(storm::utility::convertNumber(y[firstIndexViolatingConvergence]))) { - return false; - } - } - } - convergencePhase1 = false; - return true; - } - - - bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub) { - return yi * (ub - lb) <= storm::utility::abs((relative ? (precision * xi) : (precision * storm::utility::convertNumber(2.0)))); - } - - bool preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { - lowerBoundCandidate = x[minIndex] / (storm::utility::one() - y[minIndex]); - upperBoundCandidate = x[maxIndex] / (storm::utility::one() - y[maxIndex]); - // Make sure that these candidates are at least as tight as the already known bounds - if (hasLowerBound && lowerBoundCandidate < lowerBound) { - lowerBoundCandidate = lowerBound; - } - if (hasUpperBound && upperBoundCandidate > upperBound) { - upperBoundCandidate = upperBound; - } - if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) { - return true; - } - return false; - } - - void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { - auto xIt = x.begin(); - auto xIte = x.end(); - auto yIt = y.begin(); - for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) { - ValueType currentBound = *xIt / (storm::utility::one() - *yIt); - if (currentBound < lowerBoundCandidate) { - minIndex = index; - lowerBoundCandidate = std::move(currentBound); - } else if (currentBound > upperBoundCandidate) { - maxIndex = index; - upperBoundCandidate = std::move(currentBound); - } - } - if (!hasLowerBound || lowerBoundCandidate > lowerBound) { - setLowerBound(lowerBoundCandidate); - } - if (!hasUpperBound || upperBoundCandidate < upperBound) { - setUpperBound(upperBoundCandidate); - } - } - - bool checkConvergencePhase2(storm::storage::BitVector const* relevantValues = nullptr) { - // Check whether the desired precision is reached - if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { - // The current index satisfies the desired bound. We now move to the next index that violates it - while (true) { - ++firstIndexViolatingConvergence; - if (relevantValues != nullptr) { - firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); - } - if (firstIndexViolatingConvergence == x.size()) { - // Converged! - return true; - } else { - if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { - // not converged yet - return false; - } - } - } - } - return false; - } - - std::vector& x; - std::vector& y; - - ValueType lowerBound, upperBound, decisionValue; - bool hasLowerBound, hasUpperBound, hasDecisionValue; - uint64_t minIndex, maxIndex; - bool convergencePhase1; - uint64_t firstIndexViolatingConvergence; - - std::vector matrixValues; - std::vector matrixColumns; - std::vector rowIndications; - IndexType numRows; - - bool relative; - ValueType precision; - }; template bool NativeLinearEquationSolver::solveEquationsSoundPower(Environment const& env, std::vector& x, std::vector const& b) const { @@ -784,7 +575,7 @@ namespace storm { } // TODO: implement caching for the helper - SoundPowerHelper helper(*this->A, x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().native().getPrecision())); + storm::solver::helper::SoundValueIterationHelper helper(*this->A, x, *this->cachedRowVector, env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().native().getPrecision())); // Prepare initial bounds for the solution (if given) if (this->hasLowerBound()) { diff --git a/src/storm/solver/helper/SoundValueIterationHelper.cpp b/src/storm/solver/helper/SoundValueIterationHelper.cpp new file mode 100644 index 000000000..c208cea55 --- /dev/null +++ b/src/storm/solver/helper/SoundValueIterationHelper.cpp @@ -0,0 +1,418 @@ +#include "storm/solver/helper/SoundValueIterationHelper.h" + +#include "storm/storage/SparseMatrix.h" +#include "storm/storage/BitVector.h" +#include "storm/utility/vector.h" +#include "storm/utility/macros.h" +#include "storm/utility/NumberTraits.h" + +#include "storm/exceptions/NotSupportedException.h" + + +namespace storm { + namespace solver { + namespace helper { + + template + SoundValueIterationHelper::SoundValueIterationHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), hasLowerBound(false), hasUpperBound(false), hasDecisionValue(false), convergencePhase1(true), decisionValueBlocks(false), firstIndexViolatingConvergence(0), minIndex(0), maxIndex(0), relative(relative), precision(precision), rowGroupIndices(nullptr) { + STORM_LOG_THROW(matrix.getEntryCount() < std::numeric_limits::max(), storm::exceptions::NotSupportedException, "The number of matrix entries is too large for the selected index type."); + if (!matrix.hasTrivialRowGrouping()) { + rowGroupIndices = &matrix.getRowGroupIndices(); + uint64_t sizeOfLargestRowGroup = matrix.getSizeOfLargestRowGroup(); + xTmp.resize(sizeOfLargestRowGroup); + yTmp.resize(sizeOfLargestRowGroup); + } + x.assign(x.size(), storm::utility::zero()); + y.assign(x.size(), storm::utility::one()); + + numRows = matrix.getRowCount(); + matrixValues.clear(); + matrixColumns.clear(); + rowIndications.clear(); + matrixValues.reserve(matrix.getNonzeroEntryCount()); + matrixColumns.reserve(matrix.getColumnCount()); + rowIndications.reserve(numRows + 1); + rowIndications.push_back(0); + for (IndexType r = 0; r < numRows; ++r) { + for (auto const& entry : matrix.getRow(r)) { + matrixValues.push_back(entry.getValue()); + matrixColumns.push_back(entry.getColumn()); + } + rowIndications.push_back(matrixValues.size()); + } + } + + template + SoundValueIterationHelper::SoundValueIterationHelper(SoundValueIterationHelper&& oldHelper, std::vector& x, std::vector& y, bool relative, ValueType const& precision) : x(x), y(y), xTmp(std::move(oldHelper.xTmp)), yTmp(std::move(oldHelper.yTmp)), hasLowerBound(false), hasUpperBound(false), hasDecisionValue(false), convergencePhase1(true), decisionValueBlocks(false), firstIndexViolatingConvergence(0), minIndex(0), maxIndex(0), relative(relative), precision(precision), numRows(std::move(oldHelper.numRows)), matrixValues(std::move(oldHelper.matrixValues)), matrixColumns(std::move(oldHelper.matrixColumns)), rowIndications(std::move(oldHelper.rowIndications)), rowGroupIndices(oldHelper.rowGroupIndices) { + + x.assign(x.size(), storm::utility::zero()); + y.assign(x.size(), storm::utility::one()); + } + + + template + void SoundValueIterationHelper::setLowerBound(ValueType const& value) { + hasLowerBound = true; + lowerBound = value; + } + + template + void SoundValueIterationHelper::setUpperBound(ValueType const& value) { + hasUpperBound = true; + upperBound = value; + } + + template + void SoundValueIterationHelper::multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi) { + assert(rowIndex < numRows); + ValueType xRes = bi; + ValueType yRes = storm::utility::zero(); + + auto entryIt = matrixValues.begin() + rowIndications[rowIndex]; + auto entryItE = matrixValues.begin() + rowIndications[rowIndex + 1]; + auto colIt = matrixColumns.begin() + rowIndications[rowIndex]; + for (; entryIt != entryItE; ++entryIt, ++colIt) { + xRes += *entryIt * x[*colIt]; + yRes += *entryIt * y[*colIt]; + } + xi = std::move(xRes); + yi = std::move(yRes); + } + + template + void SoundValueIterationHelper::performIterationStep(OptimizationDirection const& dir, std::vector const& b) { + if (rowGroupIndices) { + if (minimize(dir)) { + performIterationStep(b); + } else { + performIterationStep(b); + } + } else { + performIterationStep(b); + } + } + + template + void SoundValueIterationHelper::performIterationStep(std::vector const& b) { + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + IndexType row = numRows; + while (row > 0) { + --row; + multiplyRow(row, b[row], *xIt, *yIt); + ++xIt; + ++yIt; + } + } + + template + template::InternalOptimizationDirection dir> + void SoundValueIterationHelper::performIterationStep(std::vector const& b) { + if (!decisionValueBlocks) { + performIterationStepUpdateDecisionValue(b); + } else { + assert(decisionValue == getPrimaryBound()); + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = rowGroupIndices->rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = rowGroupIndices->rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + IndexType row = *groupStartIt; + ValueType xBest, yBest; + multiplyRow(row, b[row], xBest, yBest); + ++row; + // Only do more work if there are still rows in this row group + if (row != groupEnd) { + ValueType xi, yi; + ValueType bestValue = xBest + yBest * getPrimaryBound(); + for (;row < groupEnd; ++row) { + // Get the multiplication results + multiplyRow(row, b[row], xi, yi); + ValueType currentValue = xi + yi * getPrimaryBound(); + // Check if the current row is better then the previously found one + if (better(currentValue, bestValue)) { + xBest = std::move(xi); + yBest = std::move(yi); + bestValue = std::move(currentValue); + } else if (currentValue == bestValue && yBest > yi) { + // If the value for this row is not strictly better, it might still be equal and have a better y value + xBest = std::move(xi); + yBest = std::move(yi); + } + } + } + *xIt = std::move(xBest); + *yIt = std::move(yBest); + } + } + } + + template + template::InternalOptimizationDirection dir> + void SoundValueIterationHelper::performIterationStepUpdateDecisionValue(std::vector const& b) { + auto xIt = x.rbegin(); + auto yIt = y.rbegin(); + auto groupStartIt = rowGroupIndices->rbegin(); + uint64_t groupEnd = *groupStartIt; + ++groupStartIt; + for (auto groupStartIte = rowGroupIndices->rend(); groupStartIt != groupStartIte; groupEnd = *(groupStartIt++), ++xIt, ++yIt) { + // Perform the iteration for the first row in the group + uint64_t row = *groupStartIt; + ValueType xBest, yBest; + multiplyRow(row, b[row], xBest, yBest); + ++row; + // Only do more work if there are still rows in this row group + if (row != groupEnd) { + ValueType xi, yi; + uint64_t xyTmpIndex = 0; + if (hasPrimaryBound()) { + ValueType bestValue = xBest + yBest * getPrimaryBound(); + for (;row < groupEnd; ++row) { + // Get the multiplication results + multiplyRow(row, b[row], xi, yi); + ValueType currentValue = xi + yi * getPrimaryBound(); + // Check if the current row is better then the previously found one + if (better(currentValue, bestValue)) { + if (yBest < yi) { + // We need to store the 'old' best value as it might be relevant for the decision value + xTmp[xyTmpIndex] = std::move(xBest); + yTmp[xyTmpIndex] = std::move(yBest); + ++xyTmpIndex; + } + xBest = std::move(xi); + yBest = std::move(yi); + bestValue = std::move(currentValue); + } else if (yBest > yi) { + // If the value for this row is not strictly better, it might still be equal and have a better y value + if (currentValue == bestValue) { + xBest = std::move(xi); + yBest = std::move(yi); + } else { + xTmp[xyTmpIndex] = std::move(xi); + yTmp[xyTmpIndex] = std::move(yi); + ++xyTmpIndex; + } + } + } + } else { + for (;row < groupEnd; ++row) { + multiplyRow(row, b[row], xi, yi); + // Update the best choice + if (yi > yBest || (yi == yBest && better(xi, xBest))) { + xTmp[xyTmpIndex] = std::move(xBest); + yTmp[xyTmpIndex] = std::move(yBest); + ++xyTmpIndex; + xBest = std::move(xi); + yBest = std::move(yi); + } else { + xTmp[xyTmpIndex] = std::move(xi); + yTmp[xyTmpIndex] = std::move(yi); + ++xyTmpIndex; + } + } + } + + // Update the decision value + for (uint64_t i = 0; i < xyTmpIndex; ++i) { + ValueType deltaY = yBest - yTmp[i]; + if (deltaY > storm::utility::zero()) { + ValueType newDecisionValue = (xTmp[i] - xBest) / deltaY; + if (!hasDecisionValue || better(newDecisionValue, decisionValue)) { + decisionValue = std::move(newDecisionValue); + hasDecisionValue = true; + } + } + } + } + *xIt = std::move(xBest); + *yIt = std::move(yBest); + } + } + + template + bool SoundValueIterationHelper::checkConvergenceUpdateBounds(OptimizationDirection const& dir, storm::storage::BitVector const* relevantValues) { + if (rowGroupIndices) { + if (minimize(dir)) { + return checkConvergenceUpdateBounds(relevantValues); + } else { + return checkConvergenceUpdateBounds(relevantValues); + } + } else { + return checkConvergenceUpdateBounds(relevantValues); + } + } + + template + bool SoundValueIterationHelper::checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues) { + return checkConvergenceUpdateBounds(relevantValues); + } + + template + template::InternalOptimizationDirection dir> + bool SoundValueIterationHelper::checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues) { + + if (convergencePhase1) { + if (checkConvergencePhase1()) { + firstIndexViolatingConvergence = 0; + if (relevantValues != nullptr) { + firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); + } + } else { + return false; + } + } + STORM_LOG_ASSERT(!std::any_of(y.begin(), y.end(), [](ValueType value){return storm::utility::isOne(value);}), "Did not expect staying-probability 1 at this point."); + + // Reaching this point means that we are in Phase 2: + // The difference between lower and upper bound has to be < precision at every (relevant) value + + // For efficiency reasons we first check whether it is worth to compute the actual bounds. We do so by considering possibly too tight bounds + ValueType lowerBoundCandidate, upperBoundCandidate; + if (preliminaryConvergenceCheck(lowerBoundCandidate, upperBoundCandidate)) { + updateLowerUpperBound(lowerBoundCandidate, upperBoundCandidate); + if (dir != InternalOptimizationDirection::None) { + checkIfDecisionValueBlocks(); + } + return checkConvergencePhase2(relevantValues); + } + return false; + } + + template + void SoundValueIterationHelper::setSolutionVector() { + STORM_LOG_WARN_COND(hasLowerBound && hasUpperBound, "No lower or upper result bound could be computed within the given number of Iterations."); + + ValueType meanBound = (upperBound + lowerBound) / storm::utility::convertNumber(2.0); + storm::utility::vector::applyPointwise(x, y, x, [&meanBound] (ValueType const& xi, ValueType const& yi) { return xi + yi * meanBound; }); + + STORM_LOG_INFO("Sound Value Iteration terminated with lower value bound " + << (hasLowerBound ? lowerBound : storm::utility::zero()) << (hasLowerBound ? "" : "(none)") + << " and upper value bound " + << (hasUpperBound ? upperBound : storm::utility::zero()) << (hasUpperBound ? "" : "(none)") + << ". Decision value is " + << (hasDecisionValue ? decisionValue : storm::utility::zero()) << (hasDecisionValue ? "" : "(none)") + << "."); + + } + + + template + bool SoundValueIterationHelper::checkConvergencePhase1() { + // Return true if y ('the probability to stay within the matrix') is < 1 at every entry + for (; firstIndexViolatingConvergence != y.size(); ++firstIndexViolatingConvergence) { + static_assert(NumberTraits::IsExact || std::is_same::value, "Considered ValueType not handled."); + if (NumberTraits::IsExact) { + if (storm::utility::isOne(y[firstIndexViolatingConvergence])) { + return false; + } + } else { + if (storm::utility::isAlmostOne(storm::utility::convertNumber(y[firstIndexViolatingConvergence]))) { + return false; + } + } + } + convergencePhase1 = false; + return true; + } + + + template + bool SoundValueIterationHelper::isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub) { + return yi * (ub - lb) <= storm::utility::abs((relative ? (precision * xi) : (precision * storm::utility::convertNumber(2.0)))); + } + + template + template::InternalOptimizationDirection dir> + bool SoundValueIterationHelper::preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { + lowerBoundCandidate = x[minIndex] / (storm::utility::one() - y[minIndex]); + upperBoundCandidate = x[maxIndex] / (storm::utility::one() - y[maxIndex]); + // Make sure that these candidates are at least as tight as the already known bounds + if (hasLowerBound && lowerBoundCandidate < lowerBound) { + lowerBoundCandidate = lowerBound; + } + if (hasUpperBound && upperBoundCandidate > upperBound) { + upperBoundCandidate = upperBound; + } + if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBoundCandidate, upperBoundCandidate)) { + return true; + } + if (dir != InternalOptimizationDirection::None && !decisionValueBlocks) { + return hasDecisionValue && better(decisionValue, getPrimaryBound()); + } + return false; + } + + template + template::InternalOptimizationDirection dir> + void SoundValueIterationHelper::updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate) { + auto xIt = x.begin(); + auto xIte = x.end(); + auto yIt = y.begin(); + for (uint64_t index = 0; xIt != xIte; ++xIt, ++yIt, ++index) { + ValueType currentBound = *xIt / (storm::utility::one() - *yIt); + if (dir != InternalOptimizationDirection::None && decisionValueBlocks) { + if (better(getSecondaryBound(), currentBound)) { + getSecondaryIndex() = index; + getSecondaryBound() = std::move(currentBound); + } + } else { + if (currentBound < lowerBoundCandidate) { + minIndex = index; + lowerBoundCandidate = std::move(currentBound); + } else if (currentBound > upperBoundCandidate) { + maxIndex = index; + upperBoundCandidate = std::move(currentBound); + } + } + } + if ((dir != InternalOptimizationDirection::Minimize || !decisionValueBlocks) && (!hasLowerBound || lowerBoundCandidate > lowerBound)) { + setLowerBound(lowerBoundCandidate); + } + if ((dir != InternalOptimizationDirection::Maximize || !decisionValueBlocks) && (!hasUpperBound || upperBoundCandidate < upperBound)) { + setUpperBound(upperBoundCandidate); + } + } + + template + template::InternalOptimizationDirection dir> + void SoundValueIterationHelper::checkIfDecisionValueBlocks() { + // Check whether the decision value blocks now (i.e. further improvement of the primary bound would lead to a non-optimal scheduler). + if (!decisionValueBlocks && hasDecisionValue && better(decisionValue, getPrimaryBound())) { + getPrimaryBound() = decisionValue; + decisionValueBlocks = true; + } + } + + template + bool SoundValueIterationHelper::checkConvergencePhase2(storm::storage::BitVector const* relevantValues) { + // Check whether the desired precision is reached + if (isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { + // The current index satisfies the desired bound. We now move to the next index that violates it + while (true) { + ++firstIndexViolatingConvergence; + if (relevantValues != nullptr) { + firstIndexViolatingConvergence = relevantValues->getNextSetIndex(firstIndexViolatingConvergence); + } + if (firstIndexViolatingConvergence == x.size()) { + // Converged! + return true; + } else { + if (!isPreciseEnough(x[firstIndexViolatingConvergence], y[firstIndexViolatingConvergence], lowerBound, upperBound)) { + // not converged yet + return false; + } + } + } + } + return false; + } + + template class SoundValueIterationHelper; + template class SoundValueIterationHelper; + } + + } +} + diff --git a/src/storm/solver/helper/SoundValueIterationHelper.h b/src/storm/solver/helper/SoundValueIterationHelper.h new file mode 100644 index 000000000..d9f7bc476 --- /dev/null +++ b/src/storm/solver/helper/SoundValueIterationHelper.h @@ -0,0 +1,147 @@ +#pragma once + +#include + +#include "storm/solver/OptimizationDirection.h" + +namespace storm { + + namespace storage { + template + class SparseMatrix; + + class BitVector; + } + + namespace solver { + namespace helper { + + + template + class SoundValueIterationHelper { + public: + + typedef uint32_t IndexType; + + /*! + * Creates a new helper from the given data + */ + SoundValueIterationHelper(storm::storage::SparseMatrix const& matrix, std::vector& x, std::vector& y, bool relative, ValueType const& precision); + + /*! + * Creates a helper from the given data, considering the same matrix as the given old helper + */ + SoundValueIterationHelper(SoundValueIterationHelper&& oldHelper, std::vector& x, std::vector& y, bool relative, ValueType const& precision); + + /*! + * Sets the currently known lower / upper bound + */ + void setLowerBound(ValueType const& value); + void setUpperBound(ValueType const& value); + + void setSolutionVector(); + + /*! + * Performs one iteration step with respect to the given optimization direction. + */ + void performIterationStep(OptimizationDirection const& dir, std::vector const& b); + + /*! + * Performs one iteration step, assuming that the row grouping of the initial matrix is trivial. + */ + void performIterationStep(std::vector const& b); + + /*! + * Checks for convergence and updates the known lower/upper bounds. + */ + bool checkConvergenceUpdateBounds(OptimizationDirection const& dir, storm::storage::BitVector const* relevantValues = nullptr); + + /*! + * Checks for convergence and updates the known lower/upper bounds, assuming that the row grouping of the initial matrix is trivial. + */ + bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr); + + private: + + enum class InternalOptimizationDirection { + None, Minimize, Maximize + }; + + template + void performIterationStep(std::vector const& b); + + template + void performIterationStepUpdateDecisionValue(std::vector const& b); + + void multiplyRow(IndexType const& rowIndex, ValueType const& bi, ValueType& xi, ValueType& yi); + + template + bool checkConvergenceUpdateBounds(storm::storage::BitVector const* relevantValues = nullptr); + + bool checkConvergencePhase1(); + bool checkConvergencePhase2(storm::storage::BitVector const* relevantValues = nullptr); + + bool isPreciseEnough(ValueType const& xi, ValueType const& yi, ValueType const& lb, ValueType const& ub); + + template + bool preliminaryConvergenceCheck(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate); + + template + void updateLowerUpperBound(ValueType& lowerBoundCandidate, ValueType& upperBoundCandidate); + + template + void checkIfDecisionValueBlocks(); + + + // Auxiliary helper functions to avoid case distinctions due to different optimization directions + template + inline bool better(ValueType const& val1, ValueType const& val2) { + return (dir == InternalOptimizationDirection::Maximize) ? val1 > val2 : val1 < val2; + } + template + inline ValueType& getPrimaryBound() { + return (dir == InternalOptimizationDirection::Maximize) ? upperBound : lowerBound; + } + template + inline bool& hasPrimaryBound() { + return (dir == InternalOptimizationDirection::Maximize) ? hasUpperBound : hasLowerBound; + } + template + inline ValueType& getSecondaryBound() { + return (dir == InternalOptimizationDirection::Maximize) ? lowerBound : upperBound; + } + template + inline uint64_t& getPrimaryIndex() { + return (dir == InternalOptimizationDirection::Maximize) ? maxIndex : minIndex; + } + template + inline uint64_t& getSecondaryIndex() { + return (dir == InternalOptimizationDirection::Maximize) ? minIndex : maxIndex; + } + + + std::vector& x; + std::vector& y; + std::vector xTmp, yTmp; + + ValueType lowerBound, upperBound, decisionValue; + bool hasLowerBound, hasUpperBound, hasDecisionValue; + bool convergencePhase1; + bool decisionValueBlocks; + uint64_t firstIndexViolatingConvergence; + uint64_t minIndex, maxIndex; + + bool relative; + ValueType precision; + + IndexType numRows; + std::vector matrixValues; + std::vector matrixColumns; + std::vector rowIndications; + std::vector const* rowGroupIndices; + }; + + } + } +} +