From e5e4381eb8c059171d0dfd082294b8851894ae74 Mon Sep 17 00:00:00 2001 From: Jan Erik Karuc Date: Thu, 20 Feb 2020 17:44:54 +0100 Subject: [PATCH] Basic unfinished implementation, reference in header --- .../IterativeMinMaxLinearEquationSolver.cpp | 183 ++++++++++++++++-- .../IterativeMinMaxLinearEquationSolver.h | 3 + src/storm/solver/SolverSelectionOptions.h | 2 +- 3 files changed, 166 insertions(+), 22 deletions(-) diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index df53fb538..a9482a2e8 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -68,6 +68,9 @@ namespace storm { case MinMaxMethod::ValueIteration: result = solveEquationsValueIteration(env, dir, x, b); break; + case MinMaxMethod::OptimisticValueIteration: + result = solveEquationsOptimisticValueIteration(env, dir, x, b); + break; case MinMaxMethod::PolicyIteration: result = solveEquationsPolicyIteration(env, dir, x, b); break; @@ -226,6 +229,8 @@ namespace storm { } } + + template MinMaxLinearEquationSolverRequirements IterativeMinMaxLinearEquationSolver::getRequirements(Environment const& env, boost::optional const& direction, bool const& hasInitialScheduler) const { auto method = getMethod(env, storm::NumberTraits::IsExact); @@ -327,7 +332,163 @@ namespace storm { return ValueIterationResult(iterations - currentIterations, status); } - + + template + ValueType computeMaxAbsDiff(std::vector const& allValues, storm::storage::BitVector const& relevantValues, std::vector const& oldValues) { + ValueType result = storm::utility::zero(); + auto oldValueIt = oldValues.begin(); + for (auto value : relevantValues) { + result = storm::utility::max(result, storm::utility::abs(allValues[value] - *oldValueIt)); + ++oldValueIt; + } + return result; + } + + 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 + bool IterativeMinMaxLinearEquationSolver::solveEquationsOptimisticValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { + + uint64_t overallIterations = 0; + uint64_t valueIterationInvocations = 0; + + if (!this->multiplierA) { + this->multiplierA = storm::solver::MultiplierFactory().create(env, *this->A); + } + + if (!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); + } + + // By default, we can not provide any guarantee + SolverGuarantee guarantee = SolverGuarantee::None; + + // Get handle to multiplier. + storm::solver::Multiplier const &multiplier = *this->multiplierA; + + std::vector *newX = auxiliaryRowGroupVector.get(); + std::vector *currentX = &x; + std::vector ubDiffV(newX->size()); + std::vector boundsDiffV(currentX->size()); + + ValueType precision = storm::utility::convertNumber(env.solver().minMax().getPrecision()); + ValueType iterationPrecision = precision; + + SolverStatus status = SolverStatus::InProgress; + this->startMeasureProgress(); + + while (status == SolverStatus::InProgress && + overallIterations < env.solver().minMax().getMaximalNumberOfIterations()) { + + storm::solver::MultiplicationStyle multiplicationStyle = env.solver().minMax().getMultiplicationStyle(); + + // Perform value iteration until convergence + ++valueIterationInvocations; + ValueIterationResult result = performValueIteration(env, dir, currentX, newX, b, iterationPrecision, + env.solver().minMax().getRelativeTerminationCriterion(), + guarantee, overallIterations, + env.solver().minMax().getMaximalNumberOfIterations(), + multiplicationStyle); + + overallIterations += result.iterations; + + if (result.status != SolverStatus::Converged) { + status = result.status; + } else { + std::vector upperBound = guessUpperBound(env, *currentX, precision); + std::vector newUpperBound; + while (status == SolverStatus::InProgress && + overallIterations < env.solver().minMax().getMaximalNumberOfIterations()) { + // Perform value iteration stepwise for lower bound and guessed upper bound + + // Allow aliased multiplications. + bool useGaussSeidelMultiplication = multiplicationStyle == storm::solver::MultiplicationStyle::GaussSeidel; + + // Lower and upper bound iteration + // Compute x' = min/max(A*x + b). + if (useGaussSeidelMultiplication) { + // Copy over the current vectors so we can modify them in-place. + *newX = *currentX; + newUpperBound = upperBound; + // Do the calculation + multiplier.multiplyAndReduceGaussSeidel(env, dir, *newX, &b); + multiplier.multiplyAndReduceGaussSeidel(env, dir, newUpperBound, &b); + } else { + multiplier.multiplyAndReduce(env, dir, *currentX, &b, *newX); + multiplier.multiplyAndReduce(env, dir, upperBound, &b, newUpperBound); + } + + // Set iteration precision beforehand + // Calculate the maximum difference between the old and new iteration + // This method adapts utility::vector::equalModuloPrecision from performValueIteration() + iterationPrecision = computeMaxAbsDiff(*newX, *currentX, this->getRelevantValues()); + + bool up, down = true; + bool cross = false; + + // Calculate difference vector of new and old upper bound + storm::utility::vector::subtractVectors(upperBound, newUpperBound, ubDiffV); + if (storm::utility::vector::hasPositiveEntry(ubDiffV)) { + // up = false (Not all values moved up) + // TODO: Optimisation (if(this->hasUniqueSolution()) ?) + up = false; + break; + } + if (storm::utility::vector::hasNegativeEntry(ubDiffV)) { + down = false; + } + + // Calculate difference vector of upper bound and x + storm::utility::vector::subtractVectors(newUpperBound, *newX, boundsDiffV); + if (storm::utility::vector::hasNegativeEntry(boundsDiffV)) { + cross = true; + break; + } + + // TODO: else if down und s_i precision okay, calculate middle of both vectors + + ++overallIterations; + } + } + + iterationPrecision = iterationPrecision / storm::utility::convertNumber(2.0); + } + + // Swap the result into the output x. + if (currentX == auxiliaryRowGroupVector.get()) { + std::swap(x, *currentX); + } + + reportStatus(status, overallIterations); + + // If requested, we store the scheduler for retrieval. + if (this->isTrackSchedulerSet()) { + this->schedulerChoices = std::vector(this->A->getRowGroupCount()); + this->multiplierA->multiplyAndReduce(env, dir, x, &b, *auxiliaryRowGroupVector.get(), + &this->schedulerChoices.get()); + } + + if (!this->isCachingEnabled()) { + clearCache(); + } + + return status == SolverStatus::Converged || status == SolverStatus::TerminatedEarly; + } + + template + std::vector IterativeMinMaxLinearEquationSolver::guessUpperBound(Environment const& env, std::vector const& x, ValueType const& precision) const { + std::vector ub = x; + storm::utility::vector::scaleVectorInPlace(ub, storm::utility::one() + precision); + return x; + } + template bool IterativeMinMaxLinearEquationSolver::solveEquationsValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { if (!this->multiplierA) { @@ -418,26 +579,6 @@ namespace storm { storm::utility::vector::selectVectorValues(oldValues, relevantValues, allValues); } - template - ValueType computeMaxAbsDiff(std::vector const& allValues, storm::storage::BitVector const& relevantValues, std::vector const& oldValues) { - ValueType result = storm::utility::zero(); - auto oldValueIt = oldValues.begin(); - for (auto value : relevantValues) { - result = storm::utility::max(result, storm::utility::abs(allValues[value] - *oldValueIt)); - ++oldValueIt; - } - return result; - } - - 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; - } - /*! * This version of value iteration is sound, because it approaches the solution from below and above. This * technique is due to Haddad and Monmege (Interval iteration algorithm for MDPs and IMDPs, TCS 2017) and was diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h index 4166a47aa..32ba7d32e 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h @@ -40,10 +40,13 @@ namespace storm { bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; bool solveEquationsValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; + bool solveEquationsOptimisticValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsIntervalIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsSoundValueIteration(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsViToPi(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; + std::vector guessUpperBound(Environment const& env, std::vector const& x, ValueType const& precision) const; + bool solveEquationsRationalSearch(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const; template diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h index d21262099..8c1a39e51 100644 --- a/src/storm/solver/SolverSelectionOptions.h +++ b/src/storm/solver/SolverSelectionOptions.h @@ -6,7 +6,7 @@ namespace storm { namespace solver { - ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, TopologicalCuda, ViToPi) + ExtendEnumsWithSelectionField(MinMaxMethod, ValueIteration, OptimisticValueIteration, PolicyIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, TopologicalCuda, ViToPi) ExtendEnumsWithSelectionField(MultiplierType, Native, Gmmxx) ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration) ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration, GainBiasEquations, LraDistributionEquations)