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<typename ValueType>
- std::pair<ValueType, ValueType> 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<ValueType>(gamma / delta);
-
- if (alphaDivBetaFloor == gammaDivDeltaFloor && !storm::utility::isInteger(alphaDivBeta)) {
- std::pair<ValueType, ValueType> 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<ValueType>());
- return result;
- }
- }
-
- template<typename PreciseValueType, typename ImpreciseValueType>
- 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<PreciseValueType>(10), precision);
- PreciseValueType truncated = storm::utility::trunc<PreciseValueType>(value * powerOfTen);
- return truncated / powerOfTen;
- }
-
- template<typename PreciseValueType>
- 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<PreciseValueType>(truncated) / storm::utility::convertNumber<PreciseValueType>(powerOfTen);
- }
-
- template<typename PreciseValueType, typename ImpreciseValueType>
- PreciseValueType findRational(uint64_t precision, ImpreciseValueType const& value) {
- PreciseValueType truncatedFraction = truncateToRational<PreciseValueType>(value, precision);
-
- PreciseValueType ten = storm::utility::convertNumber<PreciseValueType>(10);
- PreciseValueType powerOfTen = storm::utility::pow<PreciseValueType>(ten, precision);
- PreciseValueType multiplier = powerOfTen / storm::utility::denominator(truncatedFraction);
- PreciseValueType numerator = storm::utility::numerator(truncatedFraction) * multiplier;
+ template<typename RationalType, typename ImpreciseType>
+ bool IterativeMinMaxLinearEquationSolver<ValueType>::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp) {
- std::pair<PreciseValueType, PreciseValueType> result = findRational<PreciseValueType>(numerator, powerOfTen, numerator + storm::utility::one<PreciseValueType>(), powerOfTen);
- return result.first / result.second;
- }
-
- template<typename ValueType>
- template<typename PreciseValueType, typename ImpreciseValueType>
- bool IterativeMinMaxLinearEquationSolver<ValueType>::sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<PreciseValueType> const& A, std::vector<ImpreciseValueType> const& x, std::vector<PreciseValueType> const& b, std::vector<PreciseValueType>& tmp) {
for (uint64_t p = 1; p < precision; ++p) {
- for (uint64_t state = 0; state < x.size(); ++state) {
- PreciseValueType rationalX = storm::utility::convertNumber<PreciseValueType, ImpreciseValueType>(x[state]);
- PreciseValueType integer = storm::utility::floor(rationalX);
- PreciseValueType fraction = rationalX - integer;
-
- tmp[state] = integer + findRational<PreciseValueType>(p, storm::utility::convertNumber<ImpreciseValueType, PreciseValueType>(fraction));
- }
- if (IterativeMinMaxLinearEquationSolver<PreciseValueType>::isSolution(dir, A, tmp, b)) {
+ storm::utility::kwek_mehlhorn::sharpen(precision, x, tmp);
+
+ if (IterativeMinMaxLinearEquationSolver<RationalType>::isSolution(dir, A, tmp, b)) {
return true;
}
}
@@ -710,8 +665,8 @@ namespace storm {
}
template<typename ValueType>
- template<typename ImpreciseValueType>
- typename std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value && !NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+ template<typename ImpreciseType>
+ typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
// Create a rational representation of the input so we can check for a proper solution later.
storm::storage::SparseMatrix<storm::RationalNumber> rationalA = this->A->template toValueType<storm::RationalNumber>();
@@ -727,52 +682,25 @@ namespace storm {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
- std::vector<ValueType>* newX = auxiliaryRowGroupVector.get();
- std::vector<ValueType>* 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<uint64_t>(storm::utility::ceil(storm::utility::log10<storm::RationalNumber>(storm::utility::one<storm::RationalNumber>() / storm::utility::convertNumber<storm::RationalNumber>(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<ValueType>(*rationalIt);
- }
- } else {
- precision = precision / 10;
- }
- }
+ // Forward the call to the core rational search routine.
+ bool converged = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(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<ValueType>(*it);
}
-
- reportStatus(status, overallIterations);
-
+
if (!this->isCachingEnabled()) {
- clearCache();
+ this->clearCache();
}
- return status == Status::Converged || status == Status::TerminatedEarly;
+ return converged;
}
template<typename ValueType>
- template<typename ImpreciseValueType>
- typename std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value && NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+ template<typename ImpreciseType>
+ typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && NumberTraits<ValueType>::IsExact, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
if (!this->linEqSolverA) {
this->linEqSolverA = this->linearEquationSolverFactory->create(*this->A);
@@ -783,95 +711,168 @@ namespace storm {
auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
}
- std::vector<ValueType>* newX = auxiliaryRowGroupVector.get();
- std::vector<ValueType>* 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<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / 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<ValueType, ImpreciseType>(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<typename ValueType>
- template<typename ImpreciseValueType>
- typename std::enable_if<!std::is_same<ValueType, ImpreciseValueType>::value, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
- // Create an imprecise solver with caching enabled.
- IterativeMinMaxLinearEquationSolver<ImpreciseValueType> impreciseSolver(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ImpreciseValueType>>());
- impreciseSolver.setMatrix(this->A->template toValueType<ImpreciseValueType>());
- impreciseSolver.createLinearEquationSolver();
- impreciseSolver.setCachingEnabled(true);
+ template<typename ImpreciseType>
+ typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, bool>::type IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+ // Translate A to its imprecise version.
+ storm::storage::SparseMatrix<ImpreciseType> impreciseA = this->A->template toValueType<ImpreciseType>();
- std::vector<ImpreciseValueType> impreciseX(x.size());
+ // Translate x to its imprecise version.
+ std::vector<ImpreciseType> impreciseX(x.size());
{
std::vector<ValueType> tmp(x.size());
this->createLowerBoundsVector(tmp);
auto targetIt = impreciseX.begin();
for (auto sourceIt = tmp.begin(); targetIt != impreciseX.end(); ++targetIt, ++sourceIt) {
- *targetIt = storm::utility::convertNumber<ImpreciseValueType, ValueType>(*sourceIt);
+ *targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
}
- std::vector<ImpreciseValueType> impreciseTmpX(x.size());
- std::vector<ImpreciseValueType> impreciseB(b.size());
+
+ // Create temporary storage for an imprecise x.
+ std::vector<ImpreciseType> impreciseTmpX(x.size());
+
+ // Translate b to its imprecise version.
+ std::vector<ImpreciseType> impreciseB(b.size());
auto targetIt = impreciseB.begin();
for (auto sourceIt = b.begin(); targetIt != impreciseB.end(); ++targetIt, ++sourceIt) {
- *targetIt = storm::utility::convertNumber<ImpreciseValueType, ValueType>(*sourceIt);
+ *targetIt = storm::utility::convertNumber<ImpreciseType, ValueType>(*sourceIt);
}
- std::vector<ImpreciseValueType>* currentX = &impreciseX;
- std::vector<ImpreciseValueType>* newX = &impreciseTmpX;
+ // Create imprecise solver from the imprecise data.
+ IterativeMinMaxLinearEquationSolver<ImpreciseType> impreciseSolver(std::make_unique<storm::solver::GeneralLinearEquationSolverFactory<ImpreciseType>>());
+ impreciseSolver.setMatrix(impreciseA);
+ impreciseSolver.createLinearEquationSolver();
+ impreciseSolver.setCachingEnabled(true);
+
+ bool converged = false;
+ try {
+ // Forward the call to the core rational search routine.
+ converged = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(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<std::vector<ValueType>>(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<ValueType>(*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<ImpreciseType>();
+ impreciseTmpX = std::vector<ImpreciseType>();
+ impreciseB = std::vector<ImpreciseType>();
+ impreciseA = storm::storage::SparseMatrix<ImpreciseType>();
+
+ 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<ValueType, ValueType>(dir, *this, *this->A, x, b, *this->A, *auxiliaryRowGroupVector, b, x);
+ }
+ if (!this->isCachingEnabled()) {
+ this->clearCache();
+ }
+
+ return converged;
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ struct TemporaryHelper {
+ static std::vector<RationalType>* getFreeRational(std::vector<RationalType>& rationalX, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
+ return &rationalX;
+ }
+
+ static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<ImpreciseType>& x, std::vector<ImpreciseType>*& currentX, std::vector<ImpreciseType>*& newX) {
+ // Nothing to do.
+ }
+ };
+
+ template<typename RationalType>
+ struct TemporaryHelper<RationalType, RationalType> {
+ static std::vector<RationalType>* getFreeRational(std::vector<RationalType>& rationalX, std::vector<RationalType>*& currentX, std::vector<RationalType>*& newX) {
+ if (&rationalX == currentX) {
+ return newX;
+ } else {
+ return currentX;
+ }
+ }
+
+ static void swapSolutions(std::vector<RationalType>& rationalX, std::vector<RationalType>*& rationalSolution, std::vector<RationalType>& x, std::vector<RationalType>*& currentX, std::vector<RationalType>*& 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<typename ValueType>
+ template<typename RationalType, typename ImpreciseType>
+ bool IterativeMinMaxLinearEquationSolver<ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver<ImpreciseType> const& impreciseSolver, storm::storage::SparseMatrix<RationalType> const& rationalA, std::vector<RationalType>& rationalX, std::vector<RationalType> const& rationalB, storm::storage::SparseMatrix<ImpreciseType> const& A, std::vector<ImpreciseType>& x, std::vector<ImpreciseType> const& b, std::vector<ImpreciseType>& tmpX) const {
+
+ std::vector<ImpreciseType>* currentX = &x;
+ std::vector<ImpreciseType>* 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<ImpreciseValueType>::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, impreciseB, storm::utility::convertNumber<ImpreciseValueType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), SolverGuarantee::LessOrEqual, overallIterations);
+ typename IterativeMinMaxLinearEquationSolver<ImpreciseType>::ValueIterationResult result = impreciseSolver.performValueIteration(dir, currentX, newX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(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<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
- bool foundSolution = sharpen(dir, p, *this->A, *currentX, b, x);
+
+ // Make sure that currentX and rationalX are not aliased.
+ std::vector<RationalType>* freeRational = TemporaryHelper<RationalType, ImpreciseType>::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<RationalType, ImpreciseType>::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<typename ValueType>
- void IterativeMinMaxLinearEquationSolver<ValueType>::reportStatus(Status status, uint64_t iterations) const {
+ void IterativeMinMaxLinearEquationSolver<ValueType>::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<ValueType>& x, std::vector<ValueType> const& b) const;
bool solveEquationsRationalSearch(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
- template<typename PreciseValueType, typename ImpreciseValueType>
- static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<PreciseValueType> const& A, std::vector<ImpreciseValueType> const& x, std::vector<PreciseValueType> const& b, std::vector<PreciseValueType>& tmp);
+ template<typename RationalType, typename ImpreciseType>
+ static bool sharpen(storm::OptimizationDirection dir, uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp);
- template<typename ImpreciseValueType>
- typename std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value && !NumberTraits<ValueType>::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
- template<typename ImpreciseValueType>
- typename std::enable_if<std::is_same<ValueType, ImpreciseValueType>::value && NumberTraits<ValueType>::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
- template<typename ImpreciseValueType>
- typename std::enable_if<!std::is_same<ValueType, ImpreciseValueType>::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+ template<typename ImpreciseType>
+ typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !NumberTraits<ValueType>::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+ template<typename ImpreciseType>
+ typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && NumberTraits<ValueType>::IsExact, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+ template<typename ImpreciseType>
+ typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, bool>::type solveEquationsRationalSearchHelper(OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+
+ template<typename RationalType, typename ImpreciseType>
+ bool solveEquationsRationalSearchHelper(OptimizationDirection dir, IterativeMinMaxLinearEquationSolver<ImpreciseType> const& impreciseSolver, storm::storage::SparseMatrix<RationalType> const& rationalA, std::vector<RationalType>& rationalX, std::vector<RationalType> const& rationalB, storm::storage::SparseMatrix<ImpreciseType> const& A, std::vector<ImpreciseType>& x, std::vector<ImpreciseType> const& b, std::vector<ImpreciseType>& tmpX) const;
void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector<ValueType>& x, std::vector<ValueType> const& b, uint_fast64_t* choice = nullptr) const;
@@ -111,7 +114,7 @@ namespace storm {
mutable std::unique_ptr<std::vector<uint64_t>> rowGroupOrdering; // A.rowGroupCount() entries
Status updateStatusIfNotConverged(Status status, std::vector<ValueType> 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<ValueType> 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<typename ValueType>
void MinMaxLinearEquationSolverFactory<ValueType>::setMinMaxMethod(MinMaxMethod const& newMethod) {
- STORM_LOG_WARN_COND(!(std::is_same<ValueType, storm::RationalNumber>::value) || newMethod == MinMaxMethod::PolicyIteration, "The selected solution method does not guarantee exact result. Consider using policy iteration.");
+ STORM_LOG_WARN_COND(!(std::is_same<ValueType, storm::RationalNumber>::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<typename IntegerType>
+ std::pair<IntegerType, IntegerType> 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<IntegerType, IntegerType> 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<IntegerType>(), storm::utility::one<IntegerType>());
+ return result;
+ }
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::IntegerType> truncateToRational(ImpreciseType const& value, uint64_t precision) {
+ typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
+
+ IntegerType powerOfTen = storm::utility::pow(storm::utility::convertNumber<IntegerType>(10ull), precision);
+ IntegerType truncated = storm::utility::trunc<RationalType>(value * powerOfTen);
+ return std::make_pair(truncated, powerOfTen);
+ }
+
+ template<typename RationalType>
+ std::pair<typename NumberTraits<RationalType>::IntegerType, typename NumberTraits<RationalType>::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<double>(value * powerOfTen);
+ return std::make_pair(truncated, powerOfTen);
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ RationalType findRational(uint64_t precision, ImpreciseType const& value) {
+ typedef typename NumberTraits<RationalType>::IntegerType IntegerType;
+
+ std::pair<IntegerType, IntegerType> truncatedFraction = truncateToRational<RationalType>(value, precision);
+ std::pair<IntegerType, IntegerType> result = findRational<IntegerType>(truncatedFraction.first, truncatedFraction.second, truncatedFraction.first + storm::utility::one<IntegerType>(), truncatedFraction.second);
+
+ // Convert one of the arguments to a rational type to not get integer division.
+ return storm::utility::convertNumber<RationalType>(result.first) / result.second;
+ }
+
+ template<typename RationalType, typename ImpreciseType>
+ void sharpen(uint64_t precision, std::vector<ImpreciseType> const& input, std::vector<RationalType>& 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<RationalType>(integer) + findRational<RationalType>(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<double> {
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<storm::ClnRationalNumber> {
static const bool SupportsExponential = false;
static const bool IsExact = true;
+
+ typedef cln::cl_I IntegerType;
};
#endif
@@ -28,6 +32,8 @@ namespace storm {
struct NumberTraits<storm::GmpRationalNumber> {
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<typename ValueType>
- ValueType trunc(ValueType const& number) {
- return std::trunc(number);
+ typename NumberTraits<ValueType>::IntegerType trunc(ValueType const& number) {
+ return static_cast<typename NumberTraits<ValueType>::IntegerType>(std::trunc(number));
}
- template<typename ValueType>
- ValueType mod(ValueType const& first, ValueType const& second) {
+ template<typename IntegerType>
+ IntegerType mod(IntegerType const& first, IntegerType const& second) {
return std::fmod(first, second);
}
- template<typename ValueType>
- ValueType numerator(ValueType const& number) {
- STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported.");
- }
-
- template<typename ValueType>
- ValueType denominator(ValueType const& number) {
- STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "Not supported.");
- }
-
template<typename ValueType>
std::string to_string(ValueType const& value) {
std::stringstream ss;
@@ -306,11 +298,6 @@ namespace storm {
return carl::toInt<carl::uint>(number);
}
- template<>
- ClnRationalNumber convertNumber(ClnRationalNumber const& number) {
- return number;
- }
-
template<>
ClnRationalNumber convertNumber(double const& number) {
return carl::rationalize<ClnRationalNumber>(number);
@@ -321,12 +308,23 @@ namespace storm {
return carl::rationalize<ClnRationalNumber>(number);
}
+ template<>
+ ClnRationalNumber convertNumber(NumberTraits<ClnRationalNumber>::IntegerType const& number) {
+ return ClnRationalNumber(number);
+ }
+
template<>
ClnRationalNumber convertNumber(uint_fast64_t const& number) {
STORM_LOG_ASSERT(static_cast<carl::uint>(number) == number, "Rationalizing failed, because the number is too large.");
return carl::rationalize<ClnRationalNumber>(static_cast<carl::uint>(number));
}
+ template<>
+ typename NumberTraits<ClnRationalNumber>::IntegerType convertNumber(uint_fast64_t const& number){
+ STORM_LOG_ASSERT(static_cast<unsigned long int>(number) == number, "Conversion failed, because the number is too large.");
+ return NumberTraits<ClnRationalNumber>::IntegerType(static_cast<unsigned long int>(number));
+ }
+
template<>
ClnRationalNumber convertNumber(int_fast64_t const& number) {
STORM_LOG_ASSERT(static_cast<carl::sint>(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<ClnRationalNumber>::IntegerType trunc(ClnRationalNumber const& number) {
+ return cln::truncate1(number);
+ }
+
+ template<>
+ typename NumberTraits<ClnRationalNumber>::IntegerType mod(NumberTraits<ClnRationalNumber>::IntegerType const& first, NumberTraits<ClnRationalNumber>::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<ClnRationalNumber>::IntegerType pow(typename NumberTraits<ClnRationalNumber>::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<ClnRationalNumber>::IntegerType numerator(ClnRationalNumber const& number) {
return carl::getNum(number);
}
template<>
- ClnRationalNumber denominator(ClnRationalNumber const& number) {
+ NumberTraits<ClnRationalNumber>::IntegerType denominator(ClnRationalNumber const& number) {
return carl::getDenom(number);
}
#endif
@@ -467,11 +469,6 @@ namespace storm {
return carl::toInt<carl::uint>(number);
}
- template<>
- GmpRationalNumber convertNumber(GmpRationalNumber const& number){
- return number;
- }
-
template<>
GmpRationalNumber convertNumber(double const& number){
return carl::rationalize<GmpRationalNumber>(number);
@@ -487,7 +484,18 @@ namespace storm {
STORM_LOG_ASSERT(static_cast<carl::uint>(number) == number, "Rationalizing failed, because the number is too large.");
return carl::rationalize<GmpRationalNumber>(static_cast<carl::uint>(number));
}
+
+ template<>
+ GmpRationalNumber convertNumber(NumberTraits<GmpRationalNumber>::IntegerType const& number) {
+ return GmpRationalNumber(number);
+ }
+ template<>
+ typename NumberTraits<GmpRationalNumber>::IntegerType convertNumber(uint_fast64_t const& number){
+ STORM_LOG_ASSERT(static_cast<unsigned long int>(number) == number, "Conversion failed, because the number is too large.");
+ return NumberTraits<GmpRationalNumber>::IntegerType(static_cast<unsigned long int>(number));
+ }
+
template<>
GmpRationalNumber convertNumber(int_fast64_t const& number){
STORM_LOG_ASSERT(static_cast<carl::sint>(number) == number, "Rationalizing failed, because the number is too large.");
@@ -535,15 +543,19 @@ namespace storm {
}
template<>
- GmpRationalNumber trunc(GmpRationalNumber const& number) {
+ typename NumberTraits<GmpRationalNumber>::IntegerType trunc(GmpRationalNumber const& number) {
// FIXME: precision issue.
- return carl::rationalize<GmpRationalNumber>(std::trunc(number.get_d()));
+ return typename NumberTraits<GmpRationalNumber>::IntegerType(static_cast<unsigned long int>(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<GmpRationalNumber>::IntegerType mod(typename NumberTraits<GmpRationalNumber>::IntegerType const& first, typename NumberTraits<GmpRationalNumber>::IntegerType const& second) {
+ return carl::mod(first, second);
+ }
+
+ template<>
+ typename NumberTraits<GmpRationalNumber>::IntegerType pow(typename NumberTraits<GmpRationalNumber>::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<GmpRationalNumber>::IntegerType numerator(GmpRationalNumber const& number) {
return carl::getNum(number);
}
template<>
- GmpRationalNumber denominator(GmpRationalNumber const& number) {
+ typename NumberTraits<GmpRationalNumber>::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<double>::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<storm::ClnRationalNumber>::IntegerType one();
template storm::ClnRationalNumber zero();
template storm::ClnRationalNumber infinity();
template bool isOne(storm::ClnRationalNumber const& value);
+ template bool isZero(NumberTraits<storm::ClnRationalNumber>::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<ClnRationalNumber>::IntegerType convertNumber(storm::NumberTraits<ClnRationalNumber>::IntegerType const& number);
+ template storm::NumberTraits<ClnRationalNumber>::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<storm::GmpRationalNumber>::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<storm::GmpRationalNumber>::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<GmpRationalNumber>::IntegerType convertNumber(storm::NumberTraits<GmpRationalNumber>::IntegerType const& number);
+ template storm::NumberTraits<GmpRationalNumber>::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<typename IndexType, typename ValueType> class MatrixEntry;
}
+
+ template<typename RationalType>
+ struct NumberTraits;
namespace utility {
@@ -109,17 +112,17 @@ namespace storm {
ValueType log10(ValueType const& number);
template<typename ValueType>
- ValueType trunc(ValueType const& number);
+ typename NumberTraits<ValueType>::IntegerType trunc(ValueType const& number);
- template<typename ValueType>
- ValueType mod(ValueType const& first, ValueType const& second);
+ template<typename RationalType>
+ typename NumberTraits<RationalType>::IntegerType numerator(RationalType const& number);
- template<typename ValueType>
- ValueType numerator(ValueType const& number);
-
- template<typename ValueType>
- ValueType denominator(ValueType const& number);
+ template<typename RationalType>
+ typename NumberTraits<RationalType>::IntegerType denominator(RationalType const& number);
+ template<typename IntegerType>
+ IntegerType mod(IntegerType const& first, IntegerType const& second);
+
template<typename ValueType>
std::string to_string(ValueType const& value);
}