From 254bc05e94a526fec50ed58bf74f2602c0e64481 Mon Sep 17 00:00:00 2001 From: dehnert Date: Wed, 27 Sep 2017 21:45:40 +0200 Subject: [PATCH] more work on RationalSearch --- .../IterativeMinMaxLinearEquationSolver.cpp | 234 ++++++++++++------ .../IterativeMinMaxLinearEquationSolver.h | 12 +- src/storm/utility/constants.cpp | 16 ++ src/storm/utility/constants.h | 3 + 4 files changed, 182 insertions(+), 83 deletions(-) diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index 63bcd1a54..f181cb1d0 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -305,6 +305,8 @@ namespace storm { // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. uint64_t iterations = currentIterations; + std::vector* originalX = currentX; + Status status = Status::InProgress; while (status == Status::InProgress) { // Compute x' = min/max(A*x + b). @@ -330,7 +332,12 @@ namespace storm { status = updateStatusIfNotConverged(status, *currentX, iterations, guarantee); } - return ValueIterationResult(iterations, status); + // Swap the pointers so that the output is always in currentX. + if (originalX == newX) { + std::swap(currentX, newX); + } + + return ValueIterationResult(iterations - currentIterations, status); } template @@ -386,12 +393,6 @@ namespace storm { reportStatus(result.status, result.iterations); - // If we performed an odd number of iterations, we need to swap the x and currentX, because the newest result - // is currently stored in currentX, but x is the output vector. - if (currentX == auxiliaryRowGroupVector.get()) { - std::swap(x, *currentX); - } - // If requested, we store the scheduler for retrieval. if (this->isTrackSchedulerSet()) { this->schedulerChoices = std::vector(this->A->getRowGroupCount()); @@ -656,38 +657,47 @@ namespace storm { } } - template - ValueType truncateToRational(double const& value, uint64_t precision) { + template + PreciseValueType truncateToRational(ImpreciseValueType const& value, uint64_t precision) { + STORM_LOG_ASSERT(precision < 16, "Truncating via double is not sound."); + PreciseValueType powerOfTen = storm::utility::pow(storm::utility::convertNumber(10), precision); + PreciseValueType truncated = storm::utility::trunc(value * powerOfTen); + return truncated / powerOfTen; + } + + template + PreciseValueType truncateToRational(double const& value, uint64_t precision) { STORM_LOG_ASSERT(precision < 16, "Truncating via double is not sound."); auto powerOfTen = std::pow(10, precision); double truncated = std::trunc(value * powerOfTen); - return storm::utility::convertNumber(truncated) / storm::utility::convertNumber(powerOfTen); + return storm::utility::convertNumber(truncated) / storm::utility::convertNumber(powerOfTen); } - - template - ValueType findRational(uint64_t precision, double const& value) { - ValueType truncatedFraction = truncateToRational(value, precision); - ValueType ten = storm::utility::convertNumber(10); - ValueType powerOfTen = storm::utility::pow(ten, precision); - ValueType multiplier = powerOfTen / storm::utility::denominator(truncatedFraction); - ValueType numerator = storm::utility::numerator(truncatedFraction) * multiplier; + template + PreciseValueType findRational(uint64_t precision, ImpreciseValueType const& value) { + PreciseValueType truncatedFraction = truncateToRational(value, precision); + + PreciseValueType ten = storm::utility::convertNumber(10); + PreciseValueType powerOfTen = storm::utility::pow(ten, precision); + PreciseValueType multiplier = powerOfTen / storm::utility::denominator(truncatedFraction); + PreciseValueType numerator = storm::utility::numerator(truncatedFraction) * multiplier; - std::pair result = findRational(numerator, powerOfTen, numerator + storm::utility::one(), powerOfTen); + std::pair result = findRational(numerator, powerOfTen, numerator + storm::utility::one(), powerOfTen); return result.first / result.second; } template - bool IterativeMinMaxLinearEquationSolver::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp) { + template + bool IterativeMinMaxLinearEquationSolver::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp) { for (uint64_t p = 1; p < precision; ++p) { for (uint64_t state = 0; state < x.size(); ++state) { - storm::RationalNumber rationalX = storm::utility::convertNumber(x[state]); - storm::RationalNumber integer = storm::utility::floor(rationalX); - storm::RationalNumber fraction = rationalX - integer; + PreciseValueType rationalX = storm::utility::convertNumber(x[state]); + PreciseValueType integer = storm::utility::floor(rationalX); + PreciseValueType fraction = rationalX - integer; - tmp[state] = integer + findRational(p, storm::utility::convertNumber(fraction)); + tmp[state] = integer + findRational(p, storm::utility::convertNumber(fraction)); } - if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { + if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { return true; } } @@ -701,71 +711,56 @@ namespace storm { template template - std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + typename std::enable_if::value && !NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - } - - template - template - std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Create a rational representation of the input so we can check for a proper solution later. + storm::storage::SparseMatrix rationalA = this->A->template toValueType(); + std::vector rationalX(x.size()); + std::vector rationalB = storm::utility::vector::convertNumericVector(b); - } - - template - bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - if (this->getSettings().getForceSoundness()) { - solveEquationsRationalSearchHelper(dir, x, b); - } else { - solveEquationsRationalSearchHelper(dir, x, b); + if (!this->linEqSolverA) { + this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A); + this->linEqSolverA->setCachingEnabled(true); } - // Create an imprecise solver with caching enabled. - IterativeMinMaxLinearEquationSolver impreciseSolver(std::make_unique>()); - impreciseSolver.setMatrix(this->A->template toValueType()); - impreciseSolver.createLinearEquationSolver(); - impreciseSolver.setCachingEnabled(true); - - std::vector impreciseX(x.size()); - { - std::vector tmp(x.size()); - this->createLowerBoundsVector(tmp); - auto targetIt = impreciseX.begin(); - for (auto sourceIt = tmp.begin(); targetIt != impreciseX.end(); ++targetIt, ++sourceIt) { - *targetIt = storm::utility::convertNumber(*sourceIt); - } - } - std::vector impreciseTmpX(x.size()); - std::vector impreciseB(b.size()); - auto targetIt = impreciseB.begin(); - for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) { - *targetIt = storm::utility::convertNumber(*sourceIt); + if (!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - - std::vector* currentX = &impreciseX; - std::vector* newX = &impreciseTmpX; - + + std::vector* newX = auxiliaryRowGroupVector.get(); + std::vector* currentX = &x; + Status status = Status::InProgress; uint64_t overallIterations = 0; ValueType precision = this->getSettings().getPrecision(); - impreciseSolver.startMeasureProgress(); + this->startMeasureProgress(); while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { // Perform value iteration with the current precision. - IterativeMinMaxLinearEquationSolver::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, impreciseB, storm::utility::convertNumber(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations); + ValueIterationResult result = performValueIteration(dir, currentX, newX, b, precision, this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations); // Count the iterations. overallIterations += result.iterations; // Try to sharpen the result. - uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); - bool foundSolution = sharpen(dir, p, *this->A, *currentX, b, x); + uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / storm::utility::convertNumber(precision)))); + bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, rationalX); if (foundSolution) { status = Status::Converged; + + auto rationalIt = rationalX.begin(); + for (auto it = x.begin(), ite = x.end(); it != ite; ++it, ++rationalIt) { + *it = storm::utility::convertNumber(*rationalIt); + } } else { precision = precision / 10; } } + if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { + status = Status::MaximalIterationsExceeded; + } + reportStatus(status, overallIterations); if (!this->isCachingEnabled()) { @@ -775,28 +770,25 @@ namespace storm { return status == Status::Converged || status == Status::TerminatedEarly; } - template<> - bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const { - // Create a rational representation of the input so we can check for a proper solution later. - storm::storage::SparseMatrix rationalA = this->A->toValueType(); - std::vector rationalX(x.size()); - std::vector rationalB = storm::utility::vector::convertNumericVector(b); - + template + template + typename std::enable_if::value && NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + if (!this->linEqSolverA) { this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A); this->linEqSolverA->setCachingEnabled(true); } if (!auxiliaryRowGroupVector) { - auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); + auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - std::vector* newX = auxiliaryRowGroupVector.get(); - std::vector* currentX = &x; + std::vector* newX = auxiliaryRowGroupVector.get(); + std::vector* currentX = &x; Status status = Status::InProgress; uint64_t overallIterations = 0; - double precision = this->getSettings().getPrecision(); + ValueType precision = this->getSettings().getPrecision(); this->startMeasureProgress(); while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { // Perform value iteration with the current precision. @@ -806,8 +798,76 @@ namespace storm { overallIterations += result.iterations; // Try to sharpen the result. - uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / storm::utility::convertNumber(precision)))); - bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, rationalX); + uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); + bool foundSolution = sharpen(dir, p, *this->A, *currentX, b, *newX); + + if (foundSolution) { + status = Status::Converged; + + // Swap the result in place. + if (&x == currentX) { + std::swap(*currentX, *newX); + } + } else { + precision = precision / 10; + } + } + + if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { + status = Status::MaximalIterationsExceeded; + } + + reportStatus(status, overallIterations); + + if (!this->isCachingEnabled()) { + clearCache(); + } + + return status == Status::Converged || status == Status::TerminatedEarly; + } + + template + template + typename std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Create an imprecise solver with caching enabled. + IterativeMinMaxLinearEquationSolver impreciseSolver(std::make_unique>()); + impreciseSolver.setMatrix(this->A->template toValueType()); + impreciseSolver.createLinearEquationSolver(); + impreciseSolver.setCachingEnabled(true); + + std::vector impreciseX(x.size()); + { + std::vector tmp(x.size()); + this->createLowerBoundsVector(tmp); + auto targetIt = impreciseX.begin(); + for (auto sourceIt = tmp.begin(); targetIt != impreciseX.end(); ++targetIt, ++sourceIt) { + *targetIt = storm::utility::convertNumber(*sourceIt); + } + } + std::vector impreciseTmpX(x.size()); + std::vector impreciseB(b.size()); + auto targetIt = impreciseB.begin(); + for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) { + *targetIt = storm::utility::convertNumber(*sourceIt); + } + + std::vector* currentX = &impreciseX; + std::vector* newX = &impreciseTmpX; + + Status status = Status::InProgress; + uint64_t overallIterations = 0; + ValueType precision = this->getSettings().getPrecision(); + impreciseSolver.startMeasureProgress(); + while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { + // Perform value iteration with the current precision. + typename IterativeMinMaxLinearEquationSolver::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, impreciseB, storm::utility::convertNumber(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations); + + // Count the iterations. + overallIterations += result.iterations; + + // Try to sharpen the result. + uint64_t p = storm::utility::convertNumber(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); + bool foundSolution = sharpen(dir, p, *this->A, *currentX, b, x); if (foundSolution) { status = Status::Converged; @@ -816,6 +876,11 @@ namespace storm { } } + if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { + status = Status::MaximalIterationsExceeded; + } + + reportStatus(status, overallIterations); if (!this->isCachingEnabled()) { @@ -824,6 +889,15 @@ namespace storm { return status == Status::Converged || status == Status::TerminatedEarly; } + + template + bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + if (this->getSettings().getForceSoundness()) { + return solveEquationsRationalSearchHelper(dir, x, b); + } else { + return solveEquationsRationalSearchHelper(dir, x, b); + } + } template bool IterativeMinMaxLinearEquationSolver::solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const { diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h index 1b1d40f60..974faa018 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h @@ -2,6 +2,8 @@ #include "storm/solver/MultiplicationStyle.h" +#include "storm/utility/NumberTraits.h" + #include "storm/solver/LinearEquationSolver.h" #include "storm/solver/StandardMinMaxLinearEquationSolver.h" @@ -70,12 +72,16 @@ namespace storm { bool solveEquationsSoundValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp); + + template + static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix const& A, std::vector const& x, std::vector const& b, std::vector& tmp); template - std::enable_if::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + typename std::enable_if::value && !NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + template + typename std::enable_if::value && NumberTraits::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; template - std::enable_if::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + typename std::enable_if::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 05a9de733..c5c8d7b82 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -232,6 +232,11 @@ namespace storm { return std::log10(number); } + template + ValueType trunc(ValueType const& number) { + return std::trunc(number); + } + template ValueType mod(ValueType const& first, ValueType const& second) { return std::fmod(first, second); @@ -368,6 +373,11 @@ namespace storm { return carl::log10(number); } + template<> + ClnRationalNumber trunc(ClnRationalNumber const& number) { + cln::truncate2(number); + } + template<> ClnRationalNumber mod(ClnRationalNumber const& first, ClnRationalNumber const& second) { STORM_LOG_ASSERT(isInteger(first) && isInteger(second), "Expecting integers in modulo computation."); @@ -523,6 +533,12 @@ namespace storm { GmpRationalNumber log10(GmpRationalNumber const& number) { return carl::log10(number); } + + template<> + GmpRationalNumber trunc(GmpRationalNumber const& number) { + // FIXME: precision issue. + return carl::rationalize(std::trunc(number.get_d())); + } template<> GmpRationalNumber mod(GmpRationalNumber const& first, GmpRationalNumber const& second) { diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index 8c47cc490..fbb251855 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -108,6 +108,9 @@ namespace storm { template ValueType log10(ValueType const& number); + template + ValueType trunc(ValueType const& number); + template ValueType mod(ValueType const& first, ValueType const& second);