diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 528de7d1c..9ea865138 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -2,6 +2,7 @@ #include #include "storm/solver/IterativeMinMaxLinearEquationSolver.h" +#include "storm/solver/helper/OptimisticValueIterationHelper.h" #include "storm/utility/ConstantsComparator.h" @@ -359,69 +360,9 @@ namespace storm { 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 bool IterativeMinMaxLinearEquationSolver::solveEquationsOptimisticValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { - uint64_t overallIterations = 0; - uint64_t maxOverallIterations = env.solver().minMax().getMaximalNumberOfIterations(); - uint64_t lastValueIterationIterations = 0; - uint64_t currentVerificationIterations = 0; - uint64_t valueIterationInvocations = 0; - if (!this->multiplierA) { this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); } @@ -440,10 +381,6 @@ namespace storm { // Allow aliased multiplications. storm::solver::MultiplicationStyle multiplicationStyle = env.solver().minMax().getMultiplicationStyle(); bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel; - // Relative errors - bool relative = env.solver().minMax().getRelativeTerminationCriterion(); - // Upper bound only iterations - uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations(); boost::optional relevantValues; if (this->hasRelevantValues()) { @@ -453,153 +390,34 @@ namespace storm { // x has to start with a lower bound. this->createLowerBoundsVector(x); - std::vector *currentX = &x; - std::vector *newX = auxiliaryRowGroupVector.get(); - std::vector *currentUpperBound = auxiliaryRowGroupVector2.get(); - std::vector newUpperBound(x.size()); - - ValueType two = storm::utility::convertNumber(2.0); - ValueType precision = storm::utility::convertNumber(env.solver().minMax().getPrecision()); - ValueType relativeBoundGuessingScaler = (storm::utility::one() + storm::utility::convertNumber(env.solver().ovi().getUpperBoundGuessingFactor()) * precision); - ValueType doublePrecision = precision * two; - ValueType iterationPrecision = precision; + std::vector* lowerX = &x; + std::vector* upperX = auxiliaryRowGroupVector.get(); + std::vector* auxVector = auxiliaryRowGroupVector2.get(); - SolverStatus status = SolverStatus::InProgress; this->startMeasureProgress(); - - while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations) { - - // Perform value iteration until convergence - ++valueIterationInvocations; - ValueIterationResult result = performValueIteration(env, dir, currentX, newX, b, iterationPrecision, relative, guarantee, overallIterations, env.solver().minMax().getMaximalNumberOfIterations(), multiplicationStyle); - lastValueIterationIterations = result.iterations; - overallIterations += result.iterations; - - if (result.status != SolverStatus::Converged) { - status = result.status; - } else { - bool intervalIterationNeeded = false; - currentVerificationIterations = 0; - - if (relative) { - guessUpperBoundRelative(*currentX, *currentUpperBound, relativeBoundGuessingScaler); - } else { - guessUpperBoundAbsolute(*currentX, *currentUpperBound, precision); - } - - bool cancelGuess = false; - while (status == SolverStatus::InProgress && overallIterations < maxOverallIterations && !cancelGuess) { - // Perform value iteration stepwise for lower bound and guessed upper bound - - // Lower and upper bound iteration - // Compute x' = min/max(A*x + b). + + + auto statusIters = storm::solver::helper::solveEquationsOptimisticValueIteration(env, lowerX, upperX, auxVector, + [&] (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); + }, + [&] (std::vector* y, std::vector* yPrime, uint64_t const& i) { + this->showProgressIterative(i); if (useGaussSeidelMultiplication) { // Copy over the current vectors so we can modify them in-place. // This is necessary as we want to compare the new values with the current ones. - newUpperBound = *currentUpperBound; - // Do the calculation. - multiplier.multiplyAndReduceGaussSeidel(env, dir, newUpperBound, &b); - if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { - // Now do interval iteration. - *newX = *currentX; - multiplier.multiplyAndReduceGaussSeidel(env, dir, *newX, &b); - } + *yPrime = *y; + multiplier.multiplyAndReduceGaussSeidel(env, dir, *y, &b); } else { - multiplier.multiplyAndReduce(env, dir, *currentUpperBound, &b, newUpperBound); - if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { - // Now do interval iteration. - multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX); - } - } - - bool newUpperBoundAlwaysHigherEqual = true; - bool newUpperBoundAlwaysLowerEqual = true; - bool valuesCrossed = false; - for (uint64_t i = 0; i < x.size(); ++i) { - if (newUpperBound[i] < (*currentUpperBound)[i]) { - newUpperBoundAlwaysHigherEqual = false; - } else if (newUpperBound[i] != (*currentUpperBound)[i]) { - newUpperBoundAlwaysLowerEqual = false; - } - } - - if (intervalIterationNeeded || currentVerificationIterations > upperBoundOnlyIterations) { - for (uint64_t i = 0; i < x.size(); ++i) { - if (newUpperBound[i] < (*newX)[i]) { - valuesCrossed = true; - break; - } - } - } - - // Update bounds - std::swap(currentX, newX); - std::swap(*currentUpperBound, newUpperBound); - - if (newUpperBoundAlwaysHigherEqual & ! newUpperBoundAlwaysLowerEqual) { - iterationPrecision = updateIterationPrecision(env, *currentX, *newX, relative, relevantValues); - // Not all values moved up or stayed the same - // If we have a single fixed point, we can safely set the new lower bound, to the wrongly guessed upper bound - if (this->hasUniqueSolution()) { - *currentX = *currentUpperBound; - } - break; - } else if (valuesCrossed) { - STORM_LOG_ASSERT(false, "Cross case occurred."); - iterationPrecision = updateIterationPrecision(env, *currentX, *newX, relative, relevantValues); - break; - } else if (newUpperBoundAlwaysLowerEqual) { - // 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 - // We can use max_if instead of computeMaxAbsDiff, as x is definitely a lower bound and ub is larger in all elements - // Recalculate terminationPrecision if relative error requested - bool reachedPrecision = true; - for (auto const& valueIndex : relevantValues ? relevantValues.get() : storm::storage::BitVector(x.size(), true)) { - ValueType absDiff = (*currentUpperBound)[valueIndex] - (*currentX)[valueIndex]; - if (relative) { - if (absDiff > doublePrecision * (*currentX)[valueIndex]) { - reachedPrecision = false; - break; - } - } else { - if (absDiff > doublePrecision) { - reachedPrecision = false; - break; - } - } - } - if (reachedPrecision) { - // Calculate the center of both vectors and store it in currentX - storm::utility::vector::applyPointwise(*currentX, *currentUpperBound, *currentX, [&two] (ValueType const& a, ValueType const& b) -> ValueType { return (a + b) / two; }); - status = SolverStatus::Converged; - } - else { - intervalIterationNeeded = true; - } + multiplier.multiplyAndReduce(env, dir, *y, &b, *yPrime); + std::swap(y, yPrime); } - - ValueType scaledIterationCount = storm::utility::convertNumber(currentVerificationIterations) * storm::utility::convertNumber(env.solver().ovi().getMaxVerificationIterationFactor()); - if (scaledIterationCount >= storm::utility::convertNumber(lastValueIterationIterations)) { - cancelGuess = true; - iterationPrecision = updateIterationPrecision(env, *currentX, *newX, relative, relevantValues); - } - - ++overallIterations; - ++currentVerificationIterations; - } - } - } + }, 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; }); - if (overallIterations > maxOverallIterations) { - status = SolverStatus::MaximalIterationsExceeded; - } - - // Swap the result into the output x. - if (currentX == auxiliaryRowGroupVector.get()) { - std::swap(x, *currentX); - } - - reportStatus(status, overallIterations); + reportStatus(statusIters.first, statusIters.second); // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { @@ -612,7 +430,7 @@ namespace storm { clearCache(); } - return status == SolverStatus::Converged || status == SolverStatus::TerminatedEarly; + return statusIters.first == SolverStatus::Converged || statusIters.first == SolverStatus::TerminatedEarly; } template diff --git a/src/storm/solver/helper/OptimisticValueIterationHelper.h b/src/storm/solver/helper/OptimisticValueIterationHelper.h new file mode 100644 index 000000000..1adb29b38 --- /dev/null +++ b/src/storm/solver/helper/OptimisticValueIterationHelper.h @@ -0,0 +1,279 @@ +#pragma once + +#include +#include + +#include "storm/solver/OptimizationDirection.h" +#include "storm/solver/SolverStatus.h" +#include "storm/utility/vector.h" +#include "storm/utility/ProgressMeasurement.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 { + + 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; }); + } + + } + + + /*! + * 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. + */ + template + std::pair solveEquationsOptimisticValueIteration(Environment const& env, std::vector* lowerX, std::vector* upperX, std::vector* auxVector, ValueIterationCallback const& valueIterationCallback, SingleIterationCallback const& singleIterationCallback, 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); + // Relative errors + bool relative = env.solver().minMax().getRelativeTerminationCriterion(); + // Goal precision + ValueType precision = storm::utility::convertNumber(env.solver().minMax().getPrecision()); + // Desired max difference between upperX and lowerX + ValueType doublePrecision = precision * two; + // Upper bound only iterations + uint64_t upperBoundOnlyIterations = env.solver().ovi().getUpperBoundOnlyIterations(); + // Maximum number of iterations done overall + uint64_t maxOverallIterations = env.solver().minMax().getMaximalNumberOfIterations(); + 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.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; + for (uint64_t i = 0; i < upperX->size(); ++i) { + if ((*auxVector)[i] > (*upperX)[i]) { + newUpperBoundAlwaysHigherEqual = false; + } else if ((*auxVector)[i] != (*upperX)[i]) { + newUpperBoundAlwaysLowerEqual = false; + } + } + if (newUpperBoundAlwaysHigherEqual & !newUpperBoundAlwaysLowerEqual) { + // All values moved up or stayed the same (but are not 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) { + // 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; + } + } + // 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; + } + } + } + } + } + + // 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}; + } + } + } +} +