From afd2acd06b87bb6ac13efbaa67646a33fbfa7874 Mon Sep 17 00:00:00 2001 From: dehnert Date: Thu, 28 Sep 2017 15:52:13 +0200 Subject: [PATCH] more work on Kwek-Mehlhorn approach --- .../exceptions/PrecisionExceededException.h | 13 + .../IterativeMinMaxLinearEquationSolver.cpp | 304 +++++++++--------- .../IterativeMinMaxLinearEquationSolver.h | 21 +- .../solver/MinMaxLinearEquationSolver.cpp | 2 +- src/storm/utility/KwekMehlhorn.h | 68 ++++ src/storm/utility/NumberTraits.h | 6 + src/storm/utility/constants.cpp | 101 +++--- src/storm/utility/constants.h | 19 +- 8 files changed, 321 insertions(+), 213 deletions(-) create mode 100644 src/storm/exceptions/PrecisionExceededException.h create mode 100644 src/storm/utility/KwekMehlhorn.h diff --git a/src/storm/exceptions/PrecisionExceededException.h b/src/storm/exceptions/PrecisionExceededException.h new file mode 100644 index 000000000..e520d73bc --- /dev/null +++ b/src/storm/exceptions/PrecisionExceededException.h @@ -0,0 +1,13 @@ +#pragma once + +#include "storm/exceptions/BaseException.h" +#include "storm/exceptions/ExceptionMacros.h" + +namespace storm { + namespace exceptions { + + STORM_NEW_EXCEPTION(PrecisionExceededException) + + } // namespace exceptions +} // namespace storm + diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp index f181cb1d0..45b7479b8 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp @@ -6,6 +6,8 @@ #include "storm/settings/modules/GeneralSettings.h" #include "storm/settings/modules/MinMaxEquationSolverSettings.h" +#include "storm/utility/KwekMehlhorn.h" + #include "storm/utility/vector.h" #include "storm/utility/macros.h" #include "storm/exceptions/InvalidSettingsException.h" @@ -632,7 +634,10 @@ namespace storm { // If the value does not match the one in the values vector, the given vector is not a solution. if (!comparator.isEqual(groupValue, *valueIt)) { + std::cout << "group value for " << group << " is " << groupValue << " is not " << *valueIt << ", diff is " << (groupValue - *valueIt) << std::endl; return false; + } else { + std::cout << "group value for " << group << " is " << groupValue << " is actually " << *valueIt << std::endl; } } @@ -641,63 +646,13 @@ namespace storm { } template - std::pair findRational(ValueType const& alpha, ValueType const& beta, ValueType const& gamma, ValueType const& delta) { - ValueType alphaDivBeta = alpha / beta; - ValueType alphaDivBetaFloor = storm::utility::floor(alphaDivBeta); - ValueType gammaDivDeltaFloor = storm::utility::floor(gamma / delta); - - if (alphaDivBetaFloor == gammaDivDeltaFloor && !storm::utility::isInteger(alphaDivBeta)) { - std::pair subresult = findRational(delta, storm::utility::mod(gamma, delta), beta, storm::utility::mod(alpha, beta)); - auto result = std::make_pair(alphaDivBetaFloor * subresult.first + subresult.second, subresult.first); - - return result; - } else { - auto result = std::make_pair(storm::utility::ceil(alphaDivBeta), storm::utility::one()); - return result; - } - } - - 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); - } - - 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; + 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) { - std::pair result = findRational(numerator, powerOfTen, numerator + storm::utility::one(), powerOfTen); - return result.first / result.second; - } - - template - 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) { - 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)); - } - if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { + storm::utility::kwek_mehlhorn::sharpen(precision, x, tmp); + + if (IterativeMinMaxLinearEquationSolver::isSolution(dir, A, tmp, b)) { return true; } } @@ -710,8 +665,8 @@ namespace storm { } template - template - typename std::enable_if::value && !NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + template + typename std::enable_if::value && !NumberTraits::IsExact, 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(); @@ -727,52 +682,25 @@ namespace storm { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - std::vector* newX = auxiliaryRowGroupVector.get(); - std::vector* currentX = &x; - - Status status = Status::InProgress; - uint64_t overallIterations = 0; - ValueType precision = this->getSettings().getPrecision(); - this->startMeasureProgress(); - while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { - // Perform value iteration with the current precision. - 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() / 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; - } - } + // Forward the call to the core rational search routine. + bool converged = solveEquationsRationalSearchHelper(dir, *this, rationalA, rationalX, rationalB, *this->A, x, b, *auxiliaryRowGroupVector); - if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { - status = Status::MaximalIterationsExceeded; + // Translate back rational result to imprecise result. + auto targetIt = x.begin(); + for (auto it = rationalX.begin(), ite = rationalX.end(); it != ite; ++it, ++targetIt) { + *targetIt = storm::utility::convertNumber(*it); } - - reportStatus(status, overallIterations); - + if (!this->isCachingEnabled()) { - clearCache(); + this->clearCache(); } - return status == Status::Converged || status == Status::TerminatedEarly; + return converged; } template - template - typename std::enable_if::value && NumberTraits::IsExact, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + 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); @@ -783,95 +711,168 @@ namespace storm { auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); } - std::vector* newX = auxiliaryRowGroupVector.get(); - std::vector* currentX = &x; - - Status status = Status::InProgress; - uint64_t overallIterations = 0; - ValueType precision = this->getSettings().getPrecision(); - this->startMeasureProgress(); - while (status == Status::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) { - // Perform value iteration with the current precision. - 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, *newX); - - if (foundSolution) { - status = Status::Converged; - - // Swap the result in place. - if (&x == currentX) { - std::swap(*currentX, *newX); - } - } else { - precision = precision / 10; - } - } + // Forward the call to the core rational search routine. + bool converged = solveEquationsRationalSearchHelper(dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x); - if (status == Status::InProgress && overallIterations == this->getSettings().getMaximalNumberOfIterations()) { - status = Status::MaximalIterationsExceeded; - } - - reportStatus(status, overallIterations); - if (!this->isCachingEnabled()) { - clearCache(); + this->clearCache(); } - return status == Status::Converged || status == Status::TerminatedEarly; + return converged; } 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); + template + typename std::enable_if::value, bool>::type IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + // Translate A to its imprecise version. + storm::storage::SparseMatrix impreciseA = this->A->template toValueType(); - std::vector impreciseX(x.size()); + // Translate x to its imprecise version. + 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); + *targetIt = storm::utility::convertNumber(*sourceIt); } } - std::vector impreciseTmpX(x.size()); - std::vector impreciseB(b.size()); + + // Create temporary storage for an imprecise x. + std::vector impreciseTmpX(x.size()); + + // Translate b to its imprecise version. + std::vector impreciseB(b.size()); auto targetIt = impreciseB.begin(); for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) { - *targetIt = storm::utility::convertNumber(*sourceIt); + *targetIt = storm::utility::convertNumber(*sourceIt); } - std::vector* currentX = &impreciseX; - std::vector* newX = &impreciseTmpX; + // Create imprecise solver from the imprecise data. + IterativeMinMaxLinearEquationSolver impreciseSolver(std::make_unique>()); + impreciseSolver.setMatrix(impreciseA); + impreciseSolver.createLinearEquationSolver(); + impreciseSolver.setCachingEnabled(true); + + bool converged = false; + try { + // Forward the call to the core rational search routine. + converged = solveEquationsRationalSearchHelper(dir, impreciseSolver, *this->A, x, b, impreciseA, impreciseX, impreciseB, impreciseTmpX); + } catch (storm::exceptions::PrecisionExceededException const& e) { + STORM_LOG_WARN("Precision of double was exceeded, trying to recover by switching to rational arithmetic."); + + if (!auxiliaryRowGroupVector) { + auxiliaryRowGroupVector = std::make_unique>(this->A->getRowGroupCount()); + } + + // Translate the imprecise value iteration result to the one we are going to use from now on. + auto targetIt = auxiliaryRowGroupVector->begin(); + for (auto it = impreciseX.begin(), ite = impreciseX.end(); it != ite; ++it, ++targetIt) { + *targetIt = storm::utility::convertNumber(*it); + + std::cout << "converting vector[" << std::distance(impreciseX.begin(), it) << "] from " << *it << " to " << *targetIt << std::endl; + } + + + // Get rid of the superfluous data structures. + impreciseX = std::vector(); + impreciseTmpX = std::vector(); + impreciseB = std::vector(); + impreciseA = storm::storage::SparseMatrix(); + + if (!this->linEqSolverA) { + this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A); + this->linEqSolverA->setCachingEnabled(true); + } + + // Forward the call to the core rational search routine, but now with our value type as the imprecise value type. + converged = solveEquationsRationalSearchHelper(dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x); + } + if (!this->isCachingEnabled()) { + this->clearCache(); + } + + return converged; + } + + template + struct TemporaryHelper { + static std::vector* getFreeRational(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { + return &rationalX; + } + + static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { + // Nothing to do. + } + }; + + template + struct TemporaryHelper { + static std::vector* getFreeRational(std::vector& rationalX, std::vector*& currentX, std::vector*& newX) { + if (&rationalX == currentX) { + return newX; + } else { + return currentX; + } + } + + static void swapSolutions(std::vector& rationalX, std::vector*& rationalSolution, std::vector& x, std::vector*& currentX, std::vector*& newX) { + if (&rationalX == rationalSolution) { + // In this case, the rational solution is in place. + + // However, since the rational solution is no alias to current x, the imprecise solution is stored + // in current x and and rational x is not an alias to x, we can swap the contents of currentX to x. + std::swap(x, *currentX); + } else { + // Still, we may assume that the rational solution is not current x and is therefore new x. + std::swap(rationalX, *rationalSolution); + std::swap(x, *currentX); + } + } + }; + + template + template + bool IterativeMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const { + + std::vector* currentX = &x; + std::vector* newX = &tmpX; + Status status = Status::InProgress; uint64_t overallIterations = 0; + uint64_t valueIterationInvocations = 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); + typename IterativeMinMaxLinearEquationSolver::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, b, storm::utility::convertNumber(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations); + + // At this point, the result of the imprecise value iteration is stored in the (imprecise) current x. + + ++valueIterationInvocations; + STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << "."); // 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); + + // Make sure that currentX and rationalX are not aliased. + std::vector* freeRational = TemporaryHelper::getFreeRational(rationalX, currentX, newX); + + // Sharpen solution and place it in rationalX. + bool foundSolution = sharpen(dir, p, rationalA, *currentX, rationalB, *freeRational); + + // After sharpen, if a solution was found, it is contained in the current rational x. if (foundSolution) { status = Status::Converged; + + TemporaryHelper::swapSolutions(rationalX, freeRational, x, currentX, newX); } else { + // Increase the precision. precision = precision / 10; } } @@ -879,14 +880,9 @@ namespace storm { 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; } @@ -1019,7 +1015,7 @@ namespace storm { } template - void IterativeMinMaxLinearEquationSolver::reportStatus(Status status, uint64_t iterations) const { + void IterativeMinMaxLinearEquationSolver::reportStatus(Status status, uint64_t iterations) { switch (status) { case Status::Converged: STORM_LOG_INFO("Iterative solver converged after " << iterations << " iterations."); break; case Status::TerminatedEarly: STORM_LOG_INFO("Iterative solver terminated early after " << iterations << " iterations."); break; diff --git a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h index 974faa018..f7d062b4f 100644 --- a/src/storm/solver/IterativeMinMaxLinearEquationSolver.h +++ b/src/storm/solver/IterativeMinMaxLinearEquationSolver.h @@ -73,15 +73,18 @@ namespace storm { bool solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsRationalSearch(OptimizationDirection dir, std::vector& x, std::vector const& b) const; - 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 + 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 - 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 - typename std::enable_if::value, 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 + 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, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + + template + bool solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver const& impreciseSolver, storm::storage::SparseMatrix const& rationalA, std::vector& rationalX, std::vector const& rationalB, storm::storage::SparseMatrix const& A, std::vector& x, std::vector const& b, std::vector& tmpX) const; void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; @@ -111,7 +114,7 @@ namespace storm { mutable std::unique_ptr> rowGroupOrdering; // A.rowGroupCount() entries Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations, SolverGuarantee const& guarantee) const; - void reportStatus(Status status, uint64_t iterations) const; + static void reportStatus(Status status, uint64_t iterations); /// The settings of this solver. IterativeMinMaxLinearEquationSolverSettings settings; diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index 9a7a3f357..3cf349c69 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -189,7 +189,7 @@ namespace storm { template void MinMaxLinearEquationSolverFactory::setMinMaxMethod(MinMaxMethod const& newMethod) { - STORM_LOG_WARN_COND(!(std::is_same::value) || newMethod == MinMaxMethod::PolicyIteration, "The selected solution method does not guarantee exact result. Consider using policy iteration."); + STORM_LOG_WARN_COND(!(std::is_same::value) || newMethod == MinMaxMethod::PolicyIteration, "The selected solution method does not guarantee exact results. Consider using policy iteration."); method = newMethod; } diff --git a/src/storm/utility/KwekMehlhorn.h b/src/storm/utility/KwekMehlhorn.h new file mode 100644 index 000000000..30c48cba3 --- /dev/null +++ b/src/storm/utility/KwekMehlhorn.h @@ -0,0 +1,68 @@ +#pragma once + +#include "storm/utility/constants.h" +#include "storm/utility/macros.h" + +#include "storm/exceptions/PrecisionExceededException.h" + +namespace storm { + namespace utility{ + namespace kwek_mehlhorn { + + template + std::pair findRational(IntegerType const& alpha, IntegerType const& beta, IntegerType const& gamma, IntegerType const& delta) { + IntegerType alphaDivBetaFloor = alpha / beta; + IntegerType gammaDivDeltaFloor = gamma / delta; + IntegerType alphaModBeta = storm::utility::mod(alpha, beta); + + if (alphaDivBetaFloor == gammaDivDeltaFloor && !storm::utility::isZero(alphaModBeta)) { + std::pair subresult = findRational(delta, storm::utility::mod(gamma, delta), beta, alphaModBeta); + auto result = std::make_pair(alphaDivBetaFloor * subresult.first + subresult.second, subresult.first); + + return result; + } else { + auto result = std::make_pair(storm::utility::isZero(alphaModBeta) ? alphaDivBetaFloor : alphaDivBetaFloor + storm::utility::one(), storm::utility::one()); + return result; + } + } + + template + std::pair::IntegerType, typename NumberTraits::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) { + typedef typename NumberTraits::IntegerType IntegerType; + + IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber(10ull), precision); + IntegerType truncated = storm::utility::trunc(value * powerOfTen); + return std::make_pair(truncated, powerOfTen); + } + + template + std::pair::IntegerType, typename NumberTraits::IntegerType> truncateToRational(double const& value, uint64_t precision) { + STORM_LOG_THROW(precision < 17, storm::exceptions::PrecisionExceededException, "Exceeded precision of double, consider switching to rational numbers."); + double powerOfTen = std::pow(10, precision); + double truncated = storm::utility::trunc(value * powerOfTen); + return std::make_pair(truncated, powerOfTen); + } + + template + RationalType findRational(uint64_t precision, ImpreciseType const& value) { + typedef typename NumberTraits::IntegerType IntegerType; + + std::pair truncatedFraction = truncateToRational(value, precision); + std::pair result = findRational(truncatedFraction.first, truncatedFraction.second, truncatedFraction.first + storm::utility::one(), truncatedFraction.second); + + // Convert one of the arguments to a rational type to not get integer division. + return storm::utility::convertNumber(result.first) / result.second; + } + + template + void sharpen(uint64_t precision, std::vector const& input, std::vector& output) { + for (uint64_t index = 0; index < input.size(); ++index) { + ImpreciseType integer = storm::utility::floor(input[index]); + ImpreciseType fraction = input[index] - integer; + output[index] = storm::utility::convertNumber(integer) + findRational(precision, fraction); + } + } + + } + } +} diff --git a/src/storm/utility/NumberTraits.h b/src/storm/utility/NumberTraits.h index c52dad448..081ce54bc 100644 --- a/src/storm/utility/NumberTraits.h +++ b/src/storm/utility/NumberTraits.h @@ -13,6 +13,8 @@ namespace storm { struct NumberTraits { static const bool SupportsExponential = true; static const bool IsExact = false; + + typedef uint64_t IntegerType; }; #if defined(STORM_HAVE_CLN) @@ -20,6 +22,8 @@ namespace storm { struct NumberTraits { static const bool SupportsExponential = false; static const bool IsExact = true; + + typedef cln::cl_I IntegerType; }; #endif @@ -28,6 +32,8 @@ namespace storm { struct NumberTraits { static const bool SupportsExponential = false; static const bool IsExact = true; + + typedef mpz_class IntegerType; }; #endif diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index c5c8d7b82..7b0bb1baf 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -12,6 +12,8 @@ #include "storm/adapters/RationalFunctionAdapter.h" +#include "storm/utility/NumberTraits.h" + #include "storm/utility/macros.h" #include "storm/exceptions/NotSupportedException.h" @@ -233,25 +235,15 @@ namespace storm { } template - ValueType trunc(ValueType const& number) { - return std::trunc(number); + typename NumberTraits::IntegerType trunc(ValueType const& number) { + return static_cast::IntegerType>(std::trunc(number)); } - template - ValueType mod(ValueType const& first, ValueType const& second) { + template + IntegerType mod(IntegerType const& first, IntegerType const& second) { return std::fmod(first, second); } - template - ValueType numerator(ValueType const& number) { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported."); - } - - template - ValueType denominator(ValueType const& number) { - STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported."); - } - template std::string to_string(ValueType const& value) { std::stringstream ss; @@ -306,11 +298,6 @@ namespace storm { return carl::toInt(number); } - template<> - ClnRationalNumber convertNumber(ClnRationalNumber const& number) { - return number; - } - template<> ClnRationalNumber convertNumber(double const& number) { return carl::rationalize(number); @@ -321,12 +308,23 @@ namespace storm { return carl::rationalize(number); } + template<> + ClnRationalNumber convertNumber(NumberTraits::IntegerType const& number) { + return ClnRationalNumber(number); + } + template<> ClnRationalNumber convertNumber(uint_fast64_t const& number) { STORM_LOG_ASSERT(static_cast(number) == number, "Rationalizing failed, because the number is too large."); return carl::rationalize(static_cast(number)); } + template<> + typename NumberTraits::IntegerType convertNumber(uint_fast64_t const& number){ + STORM_LOG_ASSERT(static_cast(number) == number, "Conversion failed, because the number is too large."); + return NumberTraits::IntegerType(static_cast(number)); + } + template<> ClnRationalNumber convertNumber(int_fast64_t const& number) { STORM_LOG_ASSERT(static_cast(number) == number, "Rationalizing failed, because the number is too large."); @@ -374,14 +372,18 @@ namespace storm { } template<> - ClnRationalNumber trunc(ClnRationalNumber const& number) { - cln::truncate2(number); + typename NumberTraits::IntegerType trunc(ClnRationalNumber const& number) { + return cln::truncate1(number); + } + + template<> + typename NumberTraits::IntegerType mod(NumberTraits::IntegerType const& first, NumberTraits::IntegerType const& second) { + return carl::mod(first, second); } template<> - ClnRationalNumber mod(ClnRationalNumber const& first, ClnRationalNumber const& second) { - STORM_LOG_ASSERT(isInteger(first) && isInteger(second), "Expecting integers in modulo computation."); - return carl::mod(carl::getNum(first), carl::getNum(second)); + typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { + return carl::pow(value, exponent); } template<> @@ -390,12 +392,12 @@ namespace storm { } template<> - ClnRationalNumber numerator(ClnRationalNumber const& number) { + NumberTraits::IntegerType numerator(ClnRationalNumber const& number) { return carl::getNum(number); } template<> - ClnRationalNumber denominator(ClnRationalNumber const& number) { + NumberTraits::IntegerType denominator(ClnRationalNumber const& number) { return carl::getDenom(number); } #endif @@ -467,11 +469,6 @@ namespace storm { return carl::toInt(number); } - template<> - GmpRationalNumber convertNumber(GmpRationalNumber const& number){ - return number; - } - template<> GmpRationalNumber convertNumber(double const& number){ return carl::rationalize(number); @@ -487,7 +484,18 @@ namespace storm { STORM_LOG_ASSERT(static_cast(number) == number, "Rationalizing failed, because the number is too large."); return carl::rationalize(static_cast(number)); } + + template<> + GmpRationalNumber convertNumber(NumberTraits::IntegerType const& number) { + return GmpRationalNumber(number); + } + template<> + typename NumberTraits::IntegerType convertNumber(uint_fast64_t const& number){ + STORM_LOG_ASSERT(static_cast(number) == number, "Conversion failed, because the number is too large."); + return NumberTraits::IntegerType(static_cast(number)); + } + template<> GmpRationalNumber convertNumber(int_fast64_t const& number){ STORM_LOG_ASSERT(static_cast(number) == number, "Rationalizing failed, because the number is too large."); @@ -535,15 +543,19 @@ namespace storm { } template<> - GmpRationalNumber trunc(GmpRationalNumber const& number) { + typename NumberTraits::IntegerType trunc(GmpRationalNumber const& number) { // FIXME: precision issue. - return carl::rationalize(std::trunc(number.get_d())); + return typename NumberTraits::IntegerType(static_cast(std::trunc(number.get_d()))); } - + template<> - GmpRationalNumber mod(GmpRationalNumber const& first, GmpRationalNumber const& second) { - STORM_LOG_ASSERT(isInteger(first) && isInteger(second), "Expecting integers in modulo computation."); - return carl::mod(carl::getNum(first), carl::getNum(second)); + typename NumberTraits::IntegerType mod(typename NumberTraits::IntegerType const& first, typename NumberTraits::IntegerType const& second) { + return carl::mod(first, second); + } + + template<> + typename NumberTraits::IntegerType pow(typename NumberTraits::IntegerType const& value, uint_fast64_t exponent) { + return carl::pow(value, exponent); } template<> @@ -552,12 +564,12 @@ namespace storm { } template<> - GmpRationalNumber numerator(GmpRationalNumber const& number) { + typename NumberTraits::IntegerType numerator(GmpRationalNumber const& number) { return carl::getNum(number); } template<> - GmpRationalNumber denominator(GmpRationalNumber const& number) { + typename NumberTraits::IntegerType denominator(GmpRationalNumber const& number) { return carl::getDenom(number); } #endif @@ -787,8 +799,7 @@ namespace storm { template double ceil(double const& number); template double log(double const& number); template double log10(double const& number); - template double denominator(double const& number); - template double numerator(double const& number); + template typename NumberTraits::IntegerType trunc(double const& number); template double mod(double const& first, double const& second); template std::string to_string(double const& value); @@ -860,14 +871,18 @@ namespace storm { #if defined(STORM_HAVE_CLN) // Instantiations for (CLN) rational number. template storm::ClnRationalNumber one(); + template NumberTraits::IntegerType one(); template storm::ClnRationalNumber zero(); template storm::ClnRationalNumber infinity(); template bool isOne(storm::ClnRationalNumber const& value); + template bool isZero(NumberTraits::IntegerType const& value); template bool isZero(storm::ClnRationalNumber const& value); template bool isConstant(storm::ClnRationalNumber const& value); template bool isInfinity(storm::ClnRationalNumber const& value); template bool isInteger(storm::ClnRationalNumber const& number); template double convertNumber(storm::ClnRationalNumber const& number); + template storm::NumberTraits::IntegerType convertNumber(storm::NumberTraits::IntegerType const& number); + template storm::NumberTraits::IntegerType convertNumber(uint64_t const& number); template uint_fast64_t convertNumber(storm::ClnRationalNumber const& number); template storm::ClnRationalNumber convertNumber(double const& number); template storm::ClnRationalNumber convertNumber(storm::ClnRationalNumber const& number); @@ -896,15 +911,19 @@ namespace storm { #if defined(STORM_HAVE_GMP) // Instantiations for (GMP) rational number. template storm::GmpRationalNumber one(); + template NumberTraits::IntegerType one(); template storm::GmpRationalNumber zero(); template storm::GmpRationalNumber infinity(); template bool isOne(storm::GmpRationalNumber const& value); template bool isZero(storm::GmpRationalNumber const& value); + template bool isZero(NumberTraits::IntegerType const& value); template bool isConstant(storm::GmpRationalNumber const& value); template bool isInfinity(storm::GmpRationalNumber const& value); template bool isInteger(storm::GmpRationalNumber const& number); template double convertNumber(storm::GmpRationalNumber const& number); template uint_fast64_t convertNumber(storm::GmpRationalNumber const& number); + template storm::NumberTraits::IntegerType convertNumber(storm::NumberTraits::IntegerType const& number); + template storm::NumberTraits::IntegerType convertNumber(uint64_t const& number); template storm::GmpRationalNumber convertNumber(double const& number); template storm::GmpRationalNumber convertNumber(storm::GmpRationalNumber const& number); template GmpRationalNumber convertNumber(std::string const& number); diff --git a/src/storm/utility/constants.h b/src/storm/utility/constants.h index fbb251855..f2ce01183 100644 --- a/src/storm/utility/constants.h +++ b/src/storm/utility/constants.h @@ -22,6 +22,9 @@ namespace storm { namespace storage { template class MatrixEntry; } + + template + struct NumberTraits; namespace utility { @@ -109,17 +112,17 @@ namespace storm { ValueType log10(ValueType const& number); template - ValueType trunc(ValueType const& number); + typename NumberTraits::IntegerType trunc(ValueType const& number); - template - ValueType mod(ValueType const& first, ValueType const& second); + template + typename NumberTraits::IntegerType numerator(RationalType const& number); - template - ValueType numerator(ValueType const& number); - - template - ValueType denominator(ValueType const& number); + template + typename NumberTraits::IntegerType denominator(RationalType const& number); + template + IntegerType mod(IntegerType const& first, IntegerType const& second); + template std::string to_string(ValueType const& value); }