diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 4f74b143e..089ba266e 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -402,11 +402,12 @@ namespace storm { this->startMeasureProgress(); - - auto statusIters = storm::solver::helper::solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, + storm::solver::helper::OptimisticValueIterationHelper helper(*this->A); + auto statusIters = helper.solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, b, [&] (std::vector*& y, std::vector*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI) { this->showProgressIterative(i); - return performValueIteration(env, dir, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle); + auto result = performValueIteration(env, dir, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle); + return std::make_pair(result.iterations, result.status); }, [&] (std::vector* y, std::vector* yPrime, uint64_t const& i) { this->showProgressIterative(i); @@ -423,6 +424,7 @@ namespace storm { env.solver().minMax().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().minMax().getPrecision()), env.solver().minMax().getMaximalNumberOfIterations(), + dir, relevantValues); auto two = storm::utility::convertNumber(2.0); storm::utility::vector::applyPointwise(*lowerX, *upperX, x, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index a27988ba6..9ca9683fc 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -675,11 +675,12 @@ namespace storm { std::vector* auxVector = this->cachedRowVector2.get(); this->startMeasureProgress(); - - auto statusIters = storm::solver::helper::solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, + storm::solver::helper::OptimisticValueIterationHelper helper(*this->A); + auto statusIters = helper.solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, b, [&] (std::vector*& y, std::vector*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI) { this->showProgressIterative(i); - return performPowerIteration(env, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle); + auto result = performPowerIteration(env, y, yPrime, b, precision, relative, guarantee, i, maxI, multiplicationStyle); + return std::make_pair(result.iterations, result.status); }, [&] (std::vector* y, std::vector* yPrime, uint64_t const& i) { this->showProgressIterative(i); @@ -696,6 +697,7 @@ namespace storm { env.solver().native().getRelativeTerminationCriterion(), storm::utility::convertNumber(env.solver().native().getPrecision()), env.solver().native().getMaximalNumberOfIterations(), + boost::none, // No optimization dir relevantValues); auto two = storm::utility::convertNumber(2.0); storm::utility::vector::applyPointwise(*lowerX, *upperX, x, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); diff --git a/src/storm/solver/helper/OptimisticValueIterationHelper.cpp b/src/storm/solver/helper/OptimisticValueIterationHelper.cpp new file mode 100644 index 000000000..d52d41d37 --- /dev/null +++ b/src/storm/solver/helper/OptimisticValueIterationHelper.cpp @@ -0,0 +1,362 @@ +#include "OptimisticValueIterationHelper.h" + +#include "storm/utility/vector.h" +#include "storm/utility/SignalHandler.h" +#include "storm/environment/solver/OviSolverEnvironment.h" + +#include "storm/exceptions/NotSupportedException.h" + +#include "storm/utility/macros.h" + +namespace storm { + + namespace solver { + namespace helper { + namespace oviinternal { + + template + ValueType computeMaxAbsDiff(std::vector const& allOldValues, std::vector const& allNewValues, storm::storage::BitVector const& relevantValues) { + ValueType result = storm::utility::zero(); + for (auto value : relevantValues) { + result = storm::utility::max(result, storm::utility::abs(allNewValues[value] - allOldValues[value])); + } + return result; + } + + template + ValueType computeMaxAbsDiff(std::vector const& allOldValues, std::vector const& allNewValues) { + ValueType result = storm::utility::zero(); + for (uint64_t i = 0; i < allOldValues.size(); ++i) { + result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i])); + } + return result; + } + + template + ValueType computeMaxRelDiff(std::vector const& allOldValues, std::vector const& allNewValues, storm::storage::BitVector const& relevantValues) { + ValueType result = storm::utility::zero(); + for (auto const& i : relevantValues) { + STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector."); + if (!storm::utility::isZero(allNewValues[i])) { + result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i]) / allNewValues[i]); + } + } + return result; + } + + template + ValueType computeMaxRelDiff(std::vector const& allOldValues, std::vector const& allNewValues) { + ValueType result = storm::utility::zero(); + for (uint64_t i = 0; i < allOldValues.size(); ++i) { + STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector."); + if (!storm::utility::isZero(allNewValues[i])) { + result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i]) / allNewValues[i]); + } + } + return result; + } + + template + ValueType updateIterationPrecision(storm::Environment const& env, std::vector const& currentX, std::vector const& newX, bool const& relative, boost::optional const& relevantValues) { + auto factor = storm::utility::convertNumber(env.solver().ovi().getPrecisionUpdateFactor()); + bool useRelevant = relevantValues.is_initialized() && env.solver().ovi().useRelevantValuesForPrecisionUpdate(); + if (relative) { + return (useRelevant ? computeMaxRelDiff(newX, currentX, relevantValues.get()) : computeMaxRelDiff(newX, currentX)) * factor; + } else { + return (useRelevant ? computeMaxAbsDiff(newX, currentX, relevantValues.get()) : computeMaxAbsDiff(newX, currentX)) * factor; + } + } + + template + void guessUpperBoundRelative(std::vector const& x, std::vector &target, ValueType const& relativeBoundGuessingScaler) { + storm::utility::vector::applyPointwise(x, target, [&relativeBoundGuessingScaler] (ValueType const& argument) -> ValueType { return argument * relativeBoundGuessingScaler; }); + } + + template + void guessUpperBoundAbsolute(std::vector const& x, std::vector &target, ValueType const& precision) { + storm::utility::vector::applyPointwise(x, target, [&precision] (ValueType const& argument) -> ValueType { return argument + precision; }); + } + + template + UpperBoundIterator::UpperBoundIterator(storm::storage::SparseMatrix const& matrix) { + STORM_LOG_THROW(static_cast(std::numeric_limits::max()) > matrix.getRowCount() + 1, storm::exceptions::NotSupportedException, "Matrix dimensions too large."); + STORM_LOG_THROW(static_cast(std::numeric_limits::max()) > matrix.getEntryCount(), storm::exceptions::NotSupportedException, "Matrix dimensions too large."); + matrixValues.reserve(matrix.getNonzeroEntryCount()); + matrixColumns.reserve(matrix.getColumnCount()); + rowIndications.reserve(matrix.getRowCount() + 1); + rowIndications.push_back(0); + for (IndexType r = 0; r < static_cast(matrix.getRowCount()); ++r) { + for (auto const& entry : matrix.getRow(r)) { + matrixValues.push_back(entry.getValue()); + matrixColumns.push_back(entry.getColumn()); + } + rowIndications.push_back(matrixValues.size()); + } + if (!matrix.hasTrivialRowGrouping()) { + rowGroupIndices = &matrix.getRowGroupIndices(); + } + } + + template + typename UpperBoundIterator::IterateResult UpperBoundIterator::iterate(std::vector& x, std::vector const& b, bool takeMinOfOldAndNew) { + return iterateInternal(x, b, takeMinOfOldAndNew); + } + + template + typename UpperBoundIterator::IterateResult UpperBoundIterator::iterate(storm::solver::OptimizationDirection const& dir, std::vector& x, std::vector const& b, bool takeMinOfOldAndNew) { + if (minimize(dir)) { + return iterateInternal(x, b, takeMinOfOldAndNew); + } else { + return iterateInternal(x, b, takeMinOfOldAndNew); + } + } + + template + template + typename UpperBoundIterator::IterateResult UpperBoundIterator::iterateInternal(std::vector& x, std::vector const& b, bool takeMinOfOldAndNew) { + // For each row compare the new upper bound candidate with the old one + bool newUpperBoundAlwaysHigherEqual = true; + bool newUpperBoundAlwaysLowerEqual = true; + // Do a backwards gauss-seidel style iteration + for (IndexType i = x.size(); i > 0;) { + --i; + ValueType newXi = HasRowGroups ? multiplyRowGroup(i, b, x) : multiplyRow(i, b[i], x); + ValueType& oldXi = x[i]; + if (newXi > oldXi) { + newUpperBoundAlwaysLowerEqual = false; + if (!takeMinOfOldAndNew) { + oldXi = newXi; + } + } else if (newXi != oldXi) { + assert(newXi < oldXi); + newUpperBoundAlwaysHigherEqual = false; + oldXi = newXi; + } + } + // Return appropriate result + if (newUpperBoundAlwaysLowerEqual) { + if (newUpperBoundAlwaysHigherEqual) { + return IterateResult::Equal; + } else { + return IterateResult::AlwaysLowerOrEqual; + } + } else { + if (newUpperBoundAlwaysHigherEqual) { + return IterateResult::AlwaysHigherOrEqual; + } else { + return IterateResult::Incomparable; + } + } + } + + template + ValueType UpperBoundIterator::multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector const& x) { + assert(rowIndex < rowIndications.size()); + ValueType xRes = bi; + + 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]; + } + return xRes; + } + + template + template + ValueType UpperBoundIterator::multiplyRowGroup(IndexType const& rowGroupIndex, std::vector const& b, std::vector const& x) { + STORM_LOG_ASSERT(rowGroupIndices != nullptr, "No row group indices available."); + auto row = (*rowGroupIndices)[rowGroupIndex]; + auto const& groupEnd = (*rowGroupIndices)[rowGroupIndex + 1]; + STORM_LOG_ASSERT(row < groupEnd, "Empty row group not expected."); + ValueType xRes = multiplyRow(row, b[row], x); + for (++row; row < groupEnd; ++row) { + ValueType xCur = multiplyRow(row, b[row], x); + xRes = minimize(Dir) ? std::min(xRes, xCur) : std::max(xRes, xCur); + } + return xRes; + } + } + + template + OptimisticValueIterationHelper::OptimisticValueIterationHelper(storm::storage::SparseMatrix const& matrix) : upperBoundIterator(matrix) { + // Intentionally left empty. + } + + template + std::pair OptimisticValueIterationHelper::solveEquationsOptimisticValueIteration(Environment const& env, std::vector* lowerX, std::vector* upperX, std::vector* auxVector, std::vector const& b, ValueIterationCallBackType const& valueIterationCallback, SingleIterationCallBackType const& singleIterationCallback, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional dir, boost::optional relevantValues) { + STORM_LOG_ASSERT(lowerX->size() == upperX->size(), "Dimension missmatch."); + STORM_LOG_ASSERT(lowerX->size() == auxVector->size(), "Dimension missmatch."); + + // As we will shuffle pointers around, let's store the original positions here. + std::vector* initLowerX = lowerX; + std::vector* initUpperX = upperX; + std::vector* initAux = auxVector; + + uint64_t overallIterations = 0; + uint64_t lastValueIterationIterations = 0; + uint64_t currentVerificationIterations = 0; + uint64_t valueIterationInvocations = 0; + + // Get some parameters for the algorithm + // 2 + ValueType two = storm::utility::convertNumber(2.0); + // Use no termination guaranteed upper bound iteration method + bool noTerminationGuarantee = env.solver().ovi().useNoTerminationGuaranteeMinimumMethod(); + // Desired max difference between upperX and lowerX + ValueType doublePrecision = precision * two; + // Upper bound only iterations + uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations(); + ValueType relativeBoundGuessingScaler = (storm::utility::one() + storm::utility::convertNumber(env.solver().ovi().getUpperBoundGuessingFactor()) * precision); + // Initial precision for the value iteration calls + ValueType iterationPrecision = precision; + + SolverStatus status = SolverStatus::InProgress; + + while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { + + // Perform value iteration until convergence + ++valueIterationInvocations; + auto result = valueIterationCallback(lowerX, auxVector, iterationPrecision, relative, overallIterations, maxOverallIterations); + lastValueIterationIterations = result.first; + overallIterations += result.first; + + if (result.second != SolverStatus::Converged) { + status = result.second; + } else { + bool intervalIterationNeeded = false; + currentVerificationIterations = 0; + + if (relative) { + oviinternal::guessUpperBoundRelative(*lowerX, *upperX, relativeBoundGuessingScaler); + } else { + oviinternal::guessUpperBoundAbsolute(*lowerX, *upperX, precision); + } + + bool cancelGuess = false; + while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { + ++overallIterations; + ++currentVerificationIterations; + // Perform value iteration stepwise for lower bound and guessed upper bound + + // Upper bound iteration + auto upperBoundIterResult = dir ? upperBoundIterator.iterate(dir.get(), *upperX, b, !noTerminationGuarantee) : upperBoundIterator.iterate(*upperX, b, !noTerminationGuarantee); + + if (upperBoundIterResult == oviinternal::UpperBoundIterator::IterateResult::AlwaysHigherOrEqual) { + // All values moved up (and did not stay the same) + // That means the guess for an upper bound is actually a lower bound + iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *upperX, relative, relevantValues); + // We assume to have a single fixed point. We can thus safely set the new lower bound, to the wrongly guessed upper bound + // Set lowerX to the upper bound candidate + std::swap(lowerX, upperX); + break; + } else if (upperBoundIterResult == oviinternal::UpperBoundIterator::IterateResult::AlwaysLowerOrEqual) { + // All values moved down (and stayed not the same) + // This is a valid upper bound. We still need to check the precision. + // We can safely use twice the requested precision, as we calculate the center of both vectors + bool reachedPrecision; + if (relevantValues) { + reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, relevantValues.get(), doublePrecision, relative); + } else { + reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, doublePrecision, relative); + } + if (reachedPrecision) { + status = SolverStatus::Converged; + break; + } else { + // From now on, we keep updating both bounds + intervalIterationNeeded = true; + } + // The following case below covers that both vectors (old and new) are equal. + // Theoretically, this means that the precise fixpoint has been reached. However, numerical instabilities can be tricky and this detection might be incorrect (see the haddad-monmege model). + // We therefore disable it. It is very unlikely that we guessed the right fixpoint anyway. + //} else if (upperBoundIterResult == oviinternal::UpperBoundIterator::IterateResult::Equal) { + // In this case, the guessed upper bound is the precise fixpoint + // status = SolverStatus::Converged; + // std::swap(lowerX, auxVector); + // break; + } + + // Check whether we tried this guess for too long + ValueType scaledIterationCount = storm::utility::convertNumber(currentVerificationIterations) * storm::utility::convertNumber(env.solver().ovi().getMaxVerificationIterationFactor()); + if (!intervalIterationNeeded && scaledIterationCount >= storm::utility::convertNumber(lastValueIterationIterations)) { + cancelGuess = true; + // In this case we will make one more iteration on the lower bound (mainly to obtain a new iterationPrecision) + } + + // Lower bound iteration (only if needed) + if (cancelGuess || intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { + singleIterationCallback(lowerX, auxVector, overallIterations); + // At this point, auxVector contains the old values for the lower bound whereas lowerX contains the new ones. + + // Check whether the upper and lower bounds have crossed, i.e., the upper bound is smaller than the lower bound. + bool valuesCrossed = false; + for (uint64_t i = 0; i < lowerX->size(); ++i) { + if ((*upperX)[i] < (*lowerX)[i]) { + valuesCrossed = true; + break; + } + } + + if (cancelGuess || valuesCrossed) { + // A new guess is needed. + iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *lowerX, relative, relevantValues); + break; + } + } + } + if (storm::utility::resources::isTerminate()) { + status = SolverStatus::Aborted; + } + } + } // end while + + // Swap the results into the output vectors. + if (initLowerX == lowerX) { + // lowerX is already at the correct position. We still have to care for upperX + if (initUpperX != upperX) { + // UpperX is not at the correct position. It has to be at the auxVector + assert(initAux == upperX); + std::swap(*initUpperX, *initAux); + } + } else if (initUpperX == upperX) { + // UpperX is already at the correct position. + // We already know that lowerX is at the wrong position. It has to be at the auxVector + assert(initAux == lowerX); + std::swap(*initLowerX, *initAux); + } else if (initAux == auxVector) { + // We know that upperX and lowerX are swapped. + assert(initLowerX == upperX); + assert(initUpperX == lowerX); + std::swap(*initUpperX, *initLowerX); + } else { + // Now we know that all vectors are at the wrong position. There are only two possibilities left + if (initLowerX == upperX) { + assert(initUpperX == auxVector); + assert(initAux == lowerX); + std::swap(*initLowerX, *initAux); + std::swap(*initUpperX, *initAux); + } else { + assert(initLowerX == auxVector); + assert(initUpperX == lowerX); + assert (initAux == upperX); + std::swap(*initUpperX, *initAux); + std::swap(*initLowerX, *initAux); + } + } + + if (overallIterations > maxOverallIterations) { + status = SolverStatus::MaximalIterationsExceeded; + } + + return {status, overallIterations}; + } + + + template class OptimisticValueIterationHelper; + template class OptimisticValueIterationHelper; + } + } +} + diff --git a/src/storm/solver/helper/OptimisticValueIterationHelper.h b/src/storm/solver/helper/OptimisticValueIterationHelper.h index 75adafb09..abf67dd3c 100644 --- a/src/storm/solver/helper/OptimisticValueIterationHelper.h +++ b/src/storm/solver/helper/OptimisticValueIterationHelper.h @@ -3,312 +3,104 @@ #include #include +#include "storm/storage/SparseMatrix.h" + #include "storm/solver/OptimizationDirection.h" #include "storm/solver/SolverStatus.h" -#include "storm/utility/vector.h" -#include "storm/utility/ProgressMeasurement.h" -#include "storm/utility/SignalHandler.h" #include "storm/storage/BitVector.h" -#include "storm/environment/solver/MinMaxSolverEnvironment.h" -#include "storm/environment/solver/OviSolverEnvironment.h" -#include "storm/utility/macros.h" namespace storm { + class Environment; namespace solver { namespace helper { namespace oviinternal { - template - ValueType computeMaxAbsDiff(std::vector const& allOldValues, std::vector const& allNewValues, storm::storage::BitVector const& relevantValues) { - ValueType result = storm::utility::zero(); - for (auto value : relevantValues) { - result = storm::utility::max(result, storm::utility::abs(allNewValues[value] - allOldValues[value])); - } - return result; - } - - template - ValueType computeMaxAbsDiff(std::vector const& allOldValues, std::vector const& allNewValues) { - ValueType result = storm::utility::zero(); - for (uint64_t i = 0; i < allOldValues.size(); ++i) { - result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i])); - } - return result; - } - - template - ValueType computeMaxRelDiff(std::vector const& allOldValues, std::vector const& allNewValues, storm::storage::BitVector const& relevantValues) { - ValueType result = storm::utility::zero(); - for (auto const& i : relevantValues) { - STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector."); - if (!storm::utility::isZero(allNewValues[i])) { - result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i]) / allNewValues[i]); - } - } - return result; - } + template + class UpperBoundIterator { + public: + typedef uint32_t IndexType; + UpperBoundIterator(storm::storage::SparseMatrix const& matrix); + enum class IterateResult { + AlwaysHigherOrEqual, + AlwaysLowerOrEqual, + Equal, + Incomparable + }; + + IterateResult iterate(std::vector& x, std::vector const& b, bool takeMinOfOldAndNew); + IterateResult iterate(storm::solver::OptimizationDirection const& dir, std::vector& x, std::vector const& b, bool takeMinOfOldAndNew); - template - ValueType computeMaxRelDiff(std::vector const& allOldValues, std::vector const& allNewValues) { - ValueType result = storm::utility::zero(); - for (uint64_t i = 0; i < allOldValues.size(); ++i) { - STORM_LOG_ASSERT(!storm::utility::isZero(allNewValues[i]) || storm::utility::isZero(allOldValues[i]), "Unexpected entry in iteration vector."); - if (!storm::utility::isZero(allNewValues[i])) { - result = storm::utility::max(result, storm::utility::abs(allNewValues[i] - allOldValues[i]) / allNewValues[i]); - } - } - return result; - } - - template - ValueType updateIterationPrecision(storm::Environment const& env, std::vector const& currentX, std::vector const& newX, bool const& relative, boost::optional const& relevantValues) { - auto factor = storm::utility::convertNumber(env.solver().ovi().getPrecisionUpdateFactor()); - bool useRelevant = relevantValues.is_initialized() && env.solver().ovi().useRelevantValuesForPrecisionUpdate(); - if (relative) { - return (useRelevant ? computeMaxRelDiff(newX, currentX, relevantValues.get()) : computeMaxRelDiff(newX, currentX)) * factor; - } else { - return (useRelevant ? computeMaxAbsDiff(newX, currentX, relevantValues.get()) : computeMaxAbsDiff(newX, currentX)) * factor; - } - } - - template - void guessUpperBoundRelative(std::vector const& x, std::vector &target, ValueType const& relativeBoundGuessingScaler) { - storm::utility::vector::applyPointwise(x, target, [&relativeBoundGuessingScaler] (ValueType const& argument) -> ValueType { return argument * relativeBoundGuessingScaler; }); - } - - template - void guessUpperBoundAbsolute(std::vector const& x, std::vector &target, ValueType const& precision) { - storm::utility::vector::applyPointwise(x, target, [&precision] (ValueType const& argument) -> ValueType { return argument + precision; }); - } - + private: + + template + IterateResult iterateInternal(std::vector& x, std::vector const& b, bool takeMinOfOldAndNew); + ValueType multiplyRow(IndexType const& rowIndex, ValueType const& bi, std::vector const& x); + template + ValueType multiplyRowGroup(IndexType const& rowGroupIndex, std::vector const& b, std::vector const& x); + + std::vector matrixValues; + std::vector matrixColumns; + std::vector rowIndications; + std::vector const* rowGroupIndices; + }; } - /*! * Performs Optimistic value iteration. - * See https://arxiv.org/abs/1910.01100 for more information on this algorithm - * - * @tparam ValueType - * @tparam ValueType - * @param env - * @param lowerX Needs to be some arbitrary lower bound on the actual values initially - * @param upperX Does not need to be an upper bound initially - * @param auxVector auxiliary storage - * @param valueIterationCallback Function that should perform standard value iteration on the input vector - * @param singleIterationCallback Function that should perform a single value iteration step on the input vector e.g. ( x' = min/max(A*x + b)) - * @param relevantValues If given, we only check the precision at the states with the given indices. - * @return The status upon termination as well as the number of iterations Also, the maximum (relative/absolute) difference between lowerX and upperX will be 2*epsilon - * with precision parameters as given by the environment env. + * @see Hartmanns, Kaminski: Optimistic Value Iteration. https://doi.org/10.1007/978-3-030-53291-8\_26 */ - template - std::pair solveEquationsOptimisticValueIteration(Environment const& env, std::vector* lowerX, std::vector* upperX, std::vector* auxVector, ValueIterationCallback const& valueIterationCallback, SingleIterationCallback const& singleIterationCallback, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional relevantValues = boost::none) { - STORM_LOG_ASSERT(lowerX->size() == upperX->size(), "Dimension missmatch."); - STORM_LOG_ASSERT(lowerX->size() == auxVector->size(), "Dimension missmatch."); - - // As we will shuffle pointers around, let's store the original positions here. - std::vector* initLowerX = lowerX; - std::vector* initUpperX = upperX; - std::vector* initAux = auxVector; - - uint64_t overallIterations = 0; - uint64_t lastValueIterationIterations = 0; - uint64_t currentVerificationIterations = 0; - uint64_t valueIterationInvocations = 0; - - // Get some parameters for the algorithm - // 2 - ValueType two = storm::utility::convertNumber(2.0); - // Use no termination guaranteed upper bound iteration method - bool noTerminationGuarantee = env.solver().ovi().useNoTerminationGuaranteeMinimumMethod(); - // Desired max difference between upperX and lowerX - ValueType doublePrecision = precision * two; - // Upper bound only iterations - uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations(); - ValueType relativeBoundGuessingScaler = (storm::utility::one() + storm::utility::convertNumber(env.solver().ovi().getUpperBoundGuessingFactor()) * precision); - // Initial precision for the value iteration calls - ValueType iterationPrecision = precision; + template + class OptimisticValueIterationHelper { + public: + OptimisticValueIterationHelper(storm::storage::SparseMatrix const& matrix); - SolverStatus status = SolverStatus::InProgress; - - while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { - - // Perform value iteration until convergence - ++valueIterationInvocations; - auto result = valueIterationCallback(lowerX, auxVector, iterationPrecision, relative, overallIterations, maxOverallIterations); - lastValueIterationIterations = result.iterations; - overallIterations += result.iterations; - - if (result.status != SolverStatus::Converged) { - status = result.status; - } else { - bool intervalIterationNeeded = false; - currentVerificationIterations = 0; - - if (relative) { - oviinternal::guessUpperBoundRelative(*lowerX, *upperX, relativeBoundGuessingScaler); - } else { - oviinternal::guessUpperBoundAbsolute(*lowerX, *upperX, precision); - } - - bool cancelGuess = false; - while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { - ++overallIterations; - ++currentVerificationIterations; - // Perform value iteration stepwise for lower bound and guessed upper bound - - // Upper bound iteration - singleIterationCallback(upperX, auxVector, overallIterations); - // At this point, auxVector contains the old values for the upper bound whereas upperX contains the new ones. - - // Compare the new upper bound candidate with the old one - bool newUpperBoundAlwaysHigherEqual = true; - bool newUpperBoundAlwaysLowerEqual = true; - if (noTerminationGuarantee) { - bool cancelOuterScan = false; - for (uint64_t i = 0; i < upperX->size() && !cancelOuterScan; ++i) { - if ((*upperX)[i] < (*auxVector)[i]) { - newUpperBoundAlwaysHigherEqual = false; - for (++i; i < upperX->size(); ++i) { - if ((*upperX)[i] > (*auxVector)[i]) { - newUpperBoundAlwaysLowerEqual = false; - cancelOuterScan = true; - break; - } - } - } else if ((*upperX)[i] != (*auxVector)[i]) { - newUpperBoundAlwaysLowerEqual = false; - for (++i; i < upperX->size(); ++i) { - if ((*upperX)[i] < (*auxVector)[i]) { - newUpperBoundAlwaysHigherEqual = false; - cancelOuterScan = true; - break; - } - } - } - } - } else { - for (uint64_t i = 0; i < upperX->size(); ++i) { - if ((*upperX)[i] < (*auxVector)[i]) { - newUpperBoundAlwaysHigherEqual = false; - } else if ((*upperX)[i] != (*auxVector)[i]) { - newUpperBoundAlwaysLowerEqual = false; - std::swap((*upperX)[i], (*auxVector)[i]); - } - } - } - - if (newUpperBoundAlwaysHigherEqual && !newUpperBoundAlwaysLowerEqual) { - // All values moved up or stayed the same - // That means the guess for an upper bound is actually a lower bound - iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *upperX, relative, relevantValues); - // We assume to have a single fixed point. We can thus safely set the new lower bound, to the wrongly guessed upper bound - // Set lowerX to the upper bound candidate - std::swap(lowerX, upperX); - break; - } else if (newUpperBoundAlwaysLowerEqual && !newUpperBoundAlwaysHigherEqual) { - // All values moved down or stayed the same and we have a maximum difference of twice the requested precision - // We can safely use twice the requested precision, as we calculate the center of both vectors - bool reachedPrecision; - if (relevantValues) { - reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, relevantValues.get(), doublePrecision, relative); - } else { - reachedPrecision = storm::utility::vector::equalModuloPrecision(*lowerX, *upperX, doublePrecision, relative); - } - if (reachedPrecision) { - status = SolverStatus::Converged; - break; - } else { - // From now on, we keep updating both bounds - intervalIterationNeeded = true; - } - // The following case below covers that bouth vectors (old and new) are equal. - // Theoretically, this means that the precise fixpoint has been reached. However, numerical instabilities can be tricky and this detection might be incorrect (see the haddad-monmege model). - // We therefore disable it. It is very unlikely that we guessed the right fixpoint anyway. - //} else if (newUpperBoundAlwaysHigherEqual && newUpperBoundAlwaysLowerEqual) { - // In this case, the guessed upper bound is the precise fixpoint - // status = SolverStatus::Converged; - // std::swap(lowerX, auxVector); - // break; - } - // At this point, the old upper bounds (auxVector) are not needed anymore. - - // Check whether we tried this guess for too long - ValueType scaledIterationCount = storm::utility::convertNumber(currentVerificationIterations) * storm::utility::convertNumber(env.solver().ovi().getMaxVerificationIterationFactor()); - if (!intervalIterationNeeded && scaledIterationCount >= storm::utility::convertNumber(lastValueIterationIterations)) { - cancelGuess = true; - // In this case we will make one more iteration on the lower bound (mainly to obtain a new iterationPrecision) - } - - // Lower bound iteration (only if needed) - if (cancelGuess || intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { - singleIterationCallback(lowerX, auxVector, overallIterations); - // At this point, auxVector contains the old values for the lower bound whereas lowerX contains the new ones. - - // Check whether the upper and lower bounds have crossed, i.e., the upper bound is smaller than the lower bound. - bool valuesCrossed = false; - for (uint64_t i = 0; i < lowerX->size(); ++i) { - if ((*upperX)[i] < (*lowerX)[i]) { - valuesCrossed = true; - break; - } - } - - if (cancelGuess || valuesCrossed) { - // A new guess is needed. - iterationPrecision = oviinternal::updateIterationPrecision(env, *auxVector, *lowerX, relative, relevantValues); - break; - } - } - } - if (storm::utility::resources::isTerminate()) { - status = SolverStatus::Aborted; - } - } - } // end while + /*! + * Return type of value iteration callback. + * The first component shall be the number of performed iterations, the second component is the final convergence status. + */ + typedef std::pair ValueIterationReturnType; - // Swap the results into the output vectors. - if (initLowerX == lowerX) { - // lowerX is already at the correct position. We still have to care for upperX - if (initUpperX != upperX) { - // UpperX is not at the correct position. It has to be at the auxVector - assert(initAux == upperX); - std::swap(*initUpperX, *initAux); - } - } else if (initUpperX == upperX) { - // UpperX is already at the correct position. - // We already know that lowerX is at the wrong position. It has to be at the auxVector - assert(initAux == lowerX); - std::swap(*initLowerX, *initAux); - } else if (initAux == auxVector) { - // We know that upperX and lowerX are swapped. - assert(initLowerX == upperX); - assert(initUpperX == lowerX); - std::swap(*initUpperX, *initLowerX); - } else { - // Now we know that all vectors are at the wrong position. There are only two possibilities left - if (initLowerX == upperX) { - assert(initUpperX == auxVector); - assert(initAux == lowerX); - std::swap(*initLowerX, *initAux); - std::swap(*initUpperX, *initAux); - } else { - assert(initLowerX == auxVector); - assert(initUpperX == lowerX); - assert (initAux == upperX); - std::swap(*initUpperX, *initAux); - std::swap(*initLowerX, *initAux); - } - } + /*! + * Value iteration callback. Performs conventional value iteration for the given input. + * @pre y points to a vector that contains starting values + * @post y points to a vector that contains resulting values + * @param yPrime points to auxiliary storage + * @param precision is the target precision + * @param relative sets whether relative precision is considered + * @param maxI the maximal number of iterations to perform + */ + typedef std::function*& y, std::vector*& yPrime, ValueType const& precision, bool const& relative, uint64_t const& i, uint64_t const& maxI)> ValueIterationCallBackType; - if (overallIterations > maxOverallIterations) { - status = SolverStatus::MaximalIterationsExceeded; - } + /*! + * Should perform a single value iteration step (using conventional value iteration). + * @pre y points to a vector that contains starting values + * @post y points to a vector that contains resulting values + * @param yPrime points to auxiliary storage + * @param currI the current iteration (can be used to display progress within the callback) + */ + typedef std::function* y, std::vector* yPrime, uint64_t const& currI)> SingleIterationCallBackType; - return {status, overallIterations}; - } + /*! + * @param env + * @param lowerX Needs to be some arbitrary lower bound on the actual values initially + * @param upperX Does not need to be an upper bound initially + * @param auxVector auxiliary storage (same size as lowerX and upperX) + * @param b the values added to each matrix row (the b in A*x+b) + * @param valueIterationCallback Function that should perform standard value iteration on the input vector + * @param singleIterationCallback Function that should perform a single value iteration step on the input vector e.g. ( x' = min/max(A*x + b)) + * @param dir The optimization direction + * @param relevantValues If given, we only check the precision at the states with the given indices. + * @return The status upon termination as well as the number of iterations Also, the maximum (relative/absolute) difference between lowerX and upperX will be 2*epsilon + * with the provided precision parameters. + */ + std::pair solveEquationsOptimisticValueIteration(Environment const& env, std::vector* lowerX, std::vector* upperX, std::vector* auxVector, std::vector const& b, ValueIterationCallBackType const& valueIterationCallback, SingleIterationCallBackType const& singleIterationCallback, bool relative, ValueType precision, uint64_t maxOverallIterations, boost::optional dir, boost::optional relevantValues = boost::none); + + private: + oviinternal::UpperBoundIterator upperBoundIterator; + + }; } } } -