sp
/
tempest
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

more work on Kwek-Mehlhorn approach

tempestpy_adaptions
dehnert 7 years ago
parent
254bc05e94
commit
afd2acd06b
8 changed files with 321 additions and 213 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 13
      src/storm/exceptions/PrecisionExceededException.h
  2. 304
      src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
  3. 21
      src/storm/solver/IterativeMinMaxLinearEquationSolver.h
  4. 2
      src/storm/solver/MinMaxLinearEquationSolver.cpp
  5. 68
      src/storm/utility/KwekMehlhorn.h
  6. 6
      src/storm/utility/NumberTraits.h
  7. 101
      src/storm/utility/constants.cpp
  8. 19
      src/storm/utility/constants.h

13
src/storm/exceptions/PrecisionExceededException.h
View File

@ -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

304
src/storm/solver/IterativeMinMaxLinearEquationSolver.cpp
View File

@ -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;

21
src/storm/solver/IterativeMinMaxLinearEquationSolver.h
View File

@ -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;

2
src/storm/solver/MinMaxLinearEquationSolver.cpp
View File

@ -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;
}

68
src/storm/utility/KwekMehlhorn.h
View File

@ -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);
}
}
}
}
}

6
src/storm/utility/NumberTraits.h
View File

@ -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

101
src/storm/utility/constants.cpp
View File

@ -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);

19
src/storm/utility/constants.h
View File

@ -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);
}

Write Preview
Loading…
Cancel
Save