Browse Source

rational search for symbolic linear equation solvers

tempestpy_adaptions
dehnert 7 years ago
parent
commit
f2e581b3df
  1. 4
      resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h
  2. 2
      resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h
  3. 18
      resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c
  4. 14
      resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h
  5. 10
      src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp
  6. 4
      src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp
  7. 6
      src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp
  8. 115
      src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp
  9. 15
      src/storm/solver/SymbolicEliminationLinearEquationSolver.h
  10. 94
      src/storm/solver/SymbolicEquationSolver.cpp
  11. 54
      src/storm/solver/SymbolicEquationSolver.h
  12. 81
      src/storm/solver/SymbolicLinearEquationSolver.cpp
  13. 48
      src/storm/solver/SymbolicLinearEquationSolver.h
  14. 33
      src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp
  15. 12
      src/storm/solver/SymbolicMinMaxLinearEquationSolver.h
  16. 236
      src/storm/solver/SymbolicNativeLinearEquationSolver.cpp
  17. 62
      src/storm/solver/SymbolicNativeLinearEquationSolver.h
  18. 6
      src/storm/storage/dd/bisimulation/QuotientExtractor.cpp
  19. 5
      src/storm/utility/constants.cpp

4
resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h

@ -164,8 +164,10 @@ TASK_DECL_3(BDD, mtbdd_min_abstract_representative, MTBDD, MTBDD, uint32_t);
TASK_DECL_3(BDD, mtbdd_max_abstract_representative, MTBDD, MTBDD, uint32_t);
#define mtbdd_max_abstract_representative(a, vars) (CALL(mtbdd_max_abstract_representative, a, vars, 0))
// A version of unary apply that performs no caching. This is needed of the argument is actually used by the unary operation,
// but may not be used to identify cache entries, for example if the argument is a pointer.
TASK_DECL_3(MTBDD, mtbdd_uapply_nocache, MTBDD, mtbdd_uapply_op, size_t);
#define mtbdd_uapply_nocache(dd, op, param) CALL(mtbdd_uapply_nocache, dd, op, param)
#define mtbdd_uapply_nocache(dd, op, param) (CALL(mtbdd_uapply_nocache, dd, op, param))
#ifdef __cplusplus
}

2
resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h

@ -74,7 +74,7 @@ TASK_DECL_3(MTBDD, sylvan_storm_rational_function_and_exists, MTBDD, MTBDD, MTBD
#define sylvan_storm_rational_function_abstract_max(dd, v) mtbdd_abstract(dd, v, TASK(sylvan_storm_rational_function_abstract_op_max))
TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_to_double, MTBDD, size_t)
#define sylvan_storm_rational_function_to_double(a) mtbdd_uapply_nocache(a, TASK(sylvan_storm_rational_function_op_to_double), 0)
#define sylvan_storm_rational_function_to_double(a) mtbdd_uapply(a, TASK(sylvan_storm_rational_function_op_to_double), 0)
TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_threshold, MTBDD, size_t*)
TASK_DECL_2(MTBDD, sylvan_storm_rational_function_op_strict_threshold, MTBDD, size_t*)

18
resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c

@ -583,8 +583,8 @@ TASK_IMPL_3(MTBDD, sylvan_storm_rational_number_and_exists, MTBDD, a, MTBDD, b,
return result;
}
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t*, svalue) {
storm_rational_number_ptr value = (storm_rational_number_ptr)svalue;
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t, svalue) {
storm_rational_number_ptr value = (storm_rational_number_ptr)(void*)svalue;
if (mtbdd_isleaf(a)) {
storm_rational_number_ptr ma = mtbdd_getstorm_rational_number_ptr(a);
@ -594,8 +594,8 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, a, size_t*,
return mtbdd_invalid;
}
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, size_t*, svalue) {
storm_rational_number_ptr value = (storm_rational_number_ptr)svalue;
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, size_t, svalue) {
storm_rational_number_ptr value = (storm_rational_number_ptr)(void*)svalue;
if (mtbdd_isleaf(a)) {
storm_rational_number_ptr ma = mtbdd_getstorm_rational_number_ptr(a);
@ -606,16 +606,6 @@ TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, a, s
return mtbdd_invalid;
}
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_threshold, MTBDD, dd, storm_rational_number_ptr, value)
{
return mtbdd_uapply(dd, TASK(sylvan_storm_rational_number_op_threshold), (size_t*)value);
}
TASK_IMPL_2(MTBDD, sylvan_storm_rational_number_strict_threshold, MTBDD, dd, storm_rational_number_ptr, value)
{
return mtbdd_uapply(dd, TASK(sylvan_storm_rational_number_op_strict_threshold), (size_t*)value);
}
TASK_IMPL_1(MTBDD, sylvan_storm_rational_number_minimum, MTBDD, a) {
/* Check terminal case */
if (a == mtbdd_false) return mtbdd_false;

14
resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h

@ -74,16 +74,12 @@ TASK_DECL_3(MTBDD, sylvan_storm_rational_number_and_exists, MTBDD, MTBDD, MTBDD)
#define sylvan_storm_rational_number_abstract_max(dd, v) mtbdd_abstract(dd, v, TASK(sylvan_storm_rational_number_abstract_op_max))
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_to_double, MTBDD, size_t)
#define sylvan_storm_rational_number_to_double(a) mtbdd_uapply_nocache(a, TASK(sylvan_storm_rational_number_op_to_double), 0)
#define sylvan_storm_rational_number_to_double(a) mtbdd_uapply(a, TASK(sylvan_storm_rational_number_op_to_double), 0)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, size_t*)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, size_t*)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_threshold, MTBDD, storm_rational_number_ptr);
#define sylvan_storm_rational_number_threshold(dd, value) CALL(sylvan_storm_rational_number_strict_threshold, dd, value)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_strict_threshold, MTBDD, storm_rational_number_ptr);
#define sylvan_storm_rational_number_strict_threshold(dd, value) CALL(sylvan_storm_rational_number_strict_threshold, dd, value)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_threshold, MTBDD, size_t)
TASK_DECL_2(MTBDD, sylvan_storm_rational_number_op_strict_threshold, MTBDD, size_t)
#define sylvan_storm_rational_number_threshold(dd, value) mtbdd_uapply_nocache(dd, TASK(sylvan_storm_rational_number_op_threshold), (size_t)(void*)value)
#define sylvan_storm_rational_number_strict_threshold(dd, value) mtbdd_uapply_nocache(dd, TASK(sylvan_storm_rational_number_op_strict_threshold), (size_t)(void*)value)
TASK_DECL_3(MTBDD, sylvan_storm_rational_number_equal_norm_d, MTBDD, MTBDD, storm_rational_number_ptr);
#define sylvan_storm_rational_number_equal_norm_d(a, b, epsilon) CALL(sylvan_storm_rational_number_equal_norm_d, a, b, epsilon)

10
src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp

@ -52,10 +52,13 @@ namespace storm {
// Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
// for solving the equation system (i.e. compute (I-A)).
submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
if (linearEquationSolverFactory.getEquationProblemFormat() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem) {
submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
}
// Solve the equation system.
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().template getAddZero<ValueType>(), subvector);
return statesWithProbability01.second.template toAdd<ValueType>() + result;
@ -168,10 +171,13 @@ namespace storm {
// Finally cut away all columns targeting non-maybe states and convert the matrix into the matrix needed
// for solving the equation system (i.e. compute (I-A)).
submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs());
submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
if (linearEquationSolverFactory.getEquationProblemFormat() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem) {
submatrix = (model.getRowColumnIdentity() * maybeStatesAdd) - submatrix;
}
// Solve the equation system.
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs());
solver->setLowerBound(storm::utility::zero<ValueType>());
storm::dd::Add<DdType, ValueType> result = solver->solveEquations(model.getManager().template getAddZero<ValueType>(), subvector);
return infinityStates.ite(model.getManager().getConstant(storm::utility::infinity<ValueType>()), result);

4
src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp

@ -87,11 +87,13 @@ namespace storm {
initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::UntilProbabilities, model, transitionMatrix, maybeStates, statesWithProbability01.second);
requirements.clearValidInitialScheduler();
}
requirements.clearBounds();
STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
}
if (initialScheduler) {
solver->setInitialScheduler(initialScheduler.get());
}
solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
solver->setRequirementsChecked();
storm::dd::Add<DdType, ValueType> result = solver->solveEquations(dir, model.getManager().template getAddZero<ValueType>(), subvector);
@ -244,11 +246,13 @@ namespace storm {
initialScheduler = computeValidSchedulerHint(storm::solver::EquationSystemType::ReachabilityRewards, model, transitionMatrix, maybeStates, targetStates);
requirements.clearValidInitialScheduler();
}
requirements.clearLowerBounds();
STORM_LOG_THROW(requirements.empty(), storm::exceptions::UncheckedRequirementException, "Could not establish requirements of solver.");
}
if (initialScheduler) {
solver->setInitialScheduler(initialScheduler.get());
}
solver->setLowerBound(storm::utility::zero<ValueType>());
solver->setRequirementsChecked();
storm::dd::Add<DdType, ValueType> result = solver->solveEquations(dir, model.getManager().template getAddZero<ValueType>(), subvector);

6
src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp

@ -28,13 +28,7 @@ namespace storm {
if (comparisonType == storm::logic::ComparisonType::Less) {
states &= values.less(bound);
} else if (comparisonType == storm::logic::ComparisonType::LessEqual) {
std::cout << "states zero: " << (values.equals(values.getDdManager().template getAddZero<ValueType>()) && reachableStates).getNonZeroCount() << std::endl;
std::cout << "states one: " << (values.equals(values.getDdManager().template getAddOne<ValueType>()) && reachableStates).getNonZeroCount() << std::endl;
std::cout << "states before: " << states.getNonZeroCount() << std::endl;
values.exportToDot("vals.dot");
std::cout << "total: " <<((values.equals(values.getDdManager().template getAddOne<ValueType>()) && states) || (values.lessOrEqual(bound) && states)).getNonZeroCount() << std::endl;
states &= values.lessOrEqual(bound);
std::cout << "states after: " << states.getNonZeroCount() << std::endl;
} else if (comparisonType == storm::logic::ComparisonType::Greater) {
states &= values.greater(bound);
} else if (comparisonType == storm::logic::ComparisonType::GreaterEqual) {

115
src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp

@ -10,6 +10,11 @@
namespace storm {
namespace solver {
template<storm::dd::DdType DdType, typename ValueType>
SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver() : SymbolicLinearEquationSolver<DdType, ValueType>() {
// Intentionally left empty.
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : SymbolicEliminationLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
this->setMatrix(A);
@ -17,6 +22,64 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType>
SymbolicEliminationLinearEquationSolver<DdType, ValueType>::SymbolicEliminationLinearEquationSolver(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : SymbolicLinearEquationSolver<DdType, ValueType>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
this->createInternalData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicEliminationLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::DdManager<DdType>& ddManager = x.getDdManager();
// Build diagonal BDD over new meta variables.
storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs);
diagonal &= this->getAllRows();
diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs);
storm::dd::Add<DdType, ValueType> rowsAdd = this->getAllRows().swapVariables(rowRowMetaVariablePairs).template toAdd<ValueType>();
storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
// Move the matrix to the new meta variables.
storm::dd::Add<DdType, ValueType> matrix = this->A.swapVariables(oldNewMetaVariablePairs);
// Initialize solution over the new meta variables.
storm::dd::Add<DdType, ValueType> solution = b.swapVariables(oldNewMetaVariablePairs);
// As long as there are transitions, we eliminate them.
uint64_t iterations = 0;
while (!matrix.isZero()) {
// Determine inverse loop probabilies.
storm::dd::Add<DdType, ValueType> inverseLoopProbabilities = rowsAdd / (rowsAdd - (diagonalAdd * matrix).sumAbstract(newColumnVariables));
// Scale all transitions with the inverse loop probabilities.
matrix *= inverseLoopProbabilities;
solution *= inverseLoopProbabilities;
// Delete diagonal elements, i.e. remove self-loops.
matrix = diagonal.ite(ddManager.template getAddZero<ValueType>(), matrix);
// Update the one-step probabilities.
solution += (matrix * solution.swapVariables(newRowColumnMetaVariablePairs)).sumAbstract(newColumnVariables);
// Shortcut all transitions by eliminating one intermediate step.
matrix = matrix.multiplyMatrix(matrix.permuteVariables(shiftMetaVariablePairs), newColumnVariables);
matrix = matrix.swapVariables(columnHelperMetaVariablePairs);
++iterations;
STORM_LOG_TRACE("Completed iteration " << iterations << " of elimination process.");
}
STORM_LOG_INFO("Elimination completed in " << iterations << " iterations.");
return solution.swapVariables(rowRowMetaVariablePairs);
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverProblemFormat SymbolicEliminationLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
return LinearEquationSolverProblemFormat::FixedPointSystem;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEliminationLinearEquationSolver<DdType, ValueType>::createInternalData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) {
storm::dd::DdManager<DdType>& ddManager = allRows.getDdManager();
// Create triple-layered meta variables for all original meta variables. We will use them later in the elimination process.
@ -61,56 +124,18 @@ namespace storm {
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicEliminationLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::DdManager<DdType>& ddManager = x.getDdManager();
// Build diagonal BDD over new meta variables.
storm::dd::Bdd<DdType> diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs);
diagonal &= this->allRows;
diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs);
storm::dd::Add<DdType, ValueType> rowsAdd = this->allRows.swapVariables(rowRowMetaVariablePairs).template toAdd<ValueType>();
storm::dd::Add<DdType, ValueType> diagonalAdd = diagonal.template toAdd<ValueType>();
// Revert the conversion to an equation system and move it to the new meta variables.
storm::dd::Add<DdType, ValueType> matrix = diagonalAdd - this->A.swapVariables(oldNewMetaVariablePairs);
// Initialize solution over the new meta variables.
storm::dd::Add<DdType, ValueType> solution = b.swapVariables(oldNewMetaVariablePairs);
// As long as there are transitions, we eliminate them.
uint64_t iterations = 0;
while (!matrix.isZero()) {
// Determine inverse loop probabilies.
storm::dd::Add<DdType, ValueType> inverseLoopProbabilities = rowsAdd / (rowsAdd - (diagonalAdd * matrix).sumAbstract(newColumnVariables));
// Scale all transitions with the inverse loop probabilities.
matrix *= inverseLoopProbabilities;
solution *= inverseLoopProbabilities;
// Delete diagonal elements, i.e. remove self-loops.
matrix = diagonal.ite(ddManager.template getAddZero<ValueType>(), matrix);
void SymbolicEliminationLinearEquationSolver<DdType, ValueType>::setData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) {
// Update the one-step probabilities.
solution += (matrix * solution.swapVariables(newRowColumnMetaVariablePairs)).sumAbstract(newColumnVariables);
// Call superclass function.
SymbolicLinearEquationSolver<DdType, ValueType>::setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
// Shortcut all transitions by eliminating one intermediate step.
matrix = matrix.multiplyMatrix(matrix.permuteVariables(shiftMetaVariablePairs), newColumnVariables);
matrix = matrix.swapVariables(columnHelperMetaVariablePairs);
++iterations;
STORM_LOG_TRACE("Completed iteration " << iterations << " of elimination process.");
}
STORM_LOG_INFO("Elimination completed in " << iterations << " iterations.");
return solution.swapVariables(rowRowMetaVariablePairs);
// Now create new variables as needed.
this->createInternalData(this->getAllRows(), this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
}
template<storm::dd::DdType DdType, typename ValueType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicEliminationLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
return std::make_unique<SymbolicEliminationLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicEliminationLinearEquationSolverFactory<DdType, ValueType>::create() const {
return std::make_unique<SymbolicEliminationLinearEquationSolver<DdType, ValueType>>();
}
template<storm::dd::DdType DdType, typename ValueType>

15
src/storm/solver/SymbolicEliminationLinearEquationSolver.h

@ -14,13 +14,26 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType = double>
class SymbolicEliminationLinearEquationSolver : public SymbolicLinearEquationSolver<DdType, ValueType> {
public:
/*!
* Constructs a symbolic linear equation solver.
*
* @param settings The settings to use.
*/
SymbolicEliminationLinearEquationSolver();
SymbolicEliminationLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
SymbolicEliminationLinearEquationSolver(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const override;
virtual void setData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) override;
virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
private:
void createInternalData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
std::vector<std::vector<storm::expressions::Variable>> oldToNewMapping;
std::set<storm::expressions::Variable> newRowVariables;
std::set<storm::expressions::Variable> newColumnVariables;
@ -38,7 +51,7 @@ namespace storm {
public:
using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const override;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
SymbolicEliminationLinearEquationSolverSettings<ValueType>& getSettings();
SymbolicEliminationLinearEquationSolverSettings<ValueType> const& getSettings() const;

94
src/storm/solver/SymbolicEquationSolver.cpp

@ -0,0 +1,94 @@
#include "storm/solver/SymbolicEquationSolver.h"
#include "storm/adapters/RationalNumberAdapter.h"
#include "storm/adapters/RationalFunctionAdapter.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/UnmetRequirementException.h"
namespace storm {
namespace solver {
template<storm::dd::DdType DdType, typename ValueType>
SymbolicEquationSolver<DdType, ValueType>::SymbolicEquationSolver(storm::dd::Bdd<DdType> const& allRows) : allRows(allRows) {
// Intentionally left empty.
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::DdManager<DdType>& SymbolicEquationSolver<DdType, ValueType>::getDdManager() const {
return this->allRows.getDdManager();
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setAllRows(storm::dd::Bdd<DdType> const& allRows) {
this->allRows = allRows;
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Bdd<DdType> const& SymbolicEquationSolver<DdType, ValueType>::getAllRows() const {
return this->allRows;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setLowerBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds) {
this->lowerBounds = lowerBounds;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setLowerBound(ValueType const& lowerBound) {
this->lowerBound = lowerBound;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setUpperBounds(storm::dd::Add<DdType, ValueType> const& upperBounds) {
this->upperBounds = upperBounds;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setUpperBound(ValueType const& upperBound) {
this->upperBound = upperBound;
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setBounds(ValueType const& lowerBound, ValueType const& upperBound) {
setLowerBound(lowerBound);
setUpperBound(upperBound);
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicEquationSolver<DdType, ValueType>::setBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds, storm::dd::Add<DdType, ValueType> const& upperBounds) {
setLowerBounds(lowerBounds);
setUpperBounds(upperBounds);
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicEquationSolver<DdType, ValueType>::getLowerBounds() const {
STORM_LOG_THROW(lowerBound || lowerBounds, storm::exceptions::UnmetRequirementException, "Requiring lower bounds, but did not get any.");
if (lowerBounds) {
return lowerBounds.get();
} else {
return this->allRows.ite(this->allRows.getDdManager().getConstant(lowerBound.get()), this->allRows.getDdManager().template getAddZero<ValueType>());
}
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicEquationSolver<DdType, ValueType>::getUpperBounds() const {
STORM_LOG_THROW(upperBound || upperBounds, storm::exceptions::UnmetRequirementException, "Requiring upper bounds, but did not get any.");
if (upperBounds) {
return upperBounds.get();
} else {
return this->allRows.ite(this->allRows.getDdManager().getConstant(upperBound.get()), this->allRows.getDdManager().template getAddZero<ValueType>());
}
}
template class SymbolicEquationSolver<storm::dd::DdType::CUDD, double>;
template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, double>;
#ifdef STORM_HAVE_CARL
template class SymbolicEquationSolver<storm::dd::DdType::CUDD, storm::RationalNumber>;
template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, storm::RationalNumber>;
template class SymbolicEquationSolver<storm::dd::DdType::Sylvan, storm::RationalFunction>;
#endif
}
}

54
src/storm/solver/SymbolicEquationSolver.h

@ -0,0 +1,54 @@
#pragma once
#include "storm/storage/dd/DdType.h"
#include "storm/storage/dd/DdManager.h"
#include "storm/storage/dd/Bdd.h"
namespace storm {
namespace solver {
template<storm::dd::DdType DdType, typename ValueType = double>
class SymbolicEquationSolver {
public:
SymbolicEquationSolver() = default;
SymbolicEquationSolver(storm::dd::Bdd<DdType> const& allRows);
void setLowerBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds);
void setLowerBound(ValueType const& lowerBound);
void setUpperBounds(storm::dd::Add<DdType, ValueType> const& upperBounds);
void setUpperBound(ValueType const& lowerBound);
void setBounds(ValueType const& lowerBound, ValueType const& upperBound);
void setBounds(storm::dd::Add<DdType, ValueType> const& lowerBounds, storm::dd::Add<DdType, ValueType> const& upperBounds);
/*!
* Retrieves a vector of lower bounds for all values (if any lower bounds are known).
*/
storm::dd::Add<DdType, ValueType> getLowerBounds() const;
/*!
* Retrieves a vector of upper bounds for all values (if any lower bounds are known).
*/
storm::dd::Add<DdType, ValueType> getUpperBounds() const;
protected:
storm::dd::DdManager<DdType>& getDdManager() const;
void setAllRows(storm::dd::Bdd<DdType> const& allRows);
storm::dd::Bdd<DdType> const& getAllRows() const;
// The relevant rows to this equation solver.
storm::dd::Bdd<DdType> allRows;
private:
// Lower bounds (if given).
boost::optional<storm::dd::Add<DdType, ValueType>> lowerBounds;
boost::optional<ValueType> lowerBound;
// Upper bounds (if given).
boost::optional<storm::dd::Add<DdType, ValueType>> upperBounds;
boost::optional<ValueType> upperBound;
};
}
}

81
src/storm/solver/SymbolicLinearEquationSolver.cpp

@ -13,19 +13,25 @@
#include "storm/settings/modules/CoreSettings.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/UnmetRequirementException.h"
#include "storm/adapters/RationalFunctionAdapter.h"
namespace storm {
namespace solver {
template<storm::dd::DdType DdType, typename ValueType>
SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver() {
// Intentionally left empty.
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : SymbolicLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) {
this->setMatrix(A);
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
SymbolicLinearEquationSolver<DdType, ValueType>::SymbolicLinearEquationSolver(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) : SymbolicEquationSolver<DdType, ValueType>(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) {
// Intentionally left empty.
}
@ -45,14 +51,28 @@ namespace storm {
return xCopy;
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverProblemFormat SymbolicLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
return LinearEquationSolverProblemFormat::EquationSystem;
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverRequirements SymbolicLinearEquationSolver<DdType, ValueType>::getRequirements() const {
// Return empty requirements by default.
return LinearEquationSolverRequirements();
}
template<storm::dd::DdType DdType, typename ValueType>
void SymbolicLinearEquationSolver<DdType, ValueType>::setMatrix(storm::dd::Add<DdType, ValueType> const& newA) {
this->A = newA;
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::DdManager<DdType>& SymbolicLinearEquationSolver<DdType, ValueType>::getDdManager() const {
return this->allRows.getDdManager();
void SymbolicLinearEquationSolver<DdType, ValueType>::setData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) {
this->setAllRows(allRows);
this->rowMetaVariables = rowMetaVariables;
this->columnMetaVariables = columnMetaVariables;
this->rowColumnMetaVariablePairs = rowColumnMetaVariablePairs;
}
template<storm::dd::DdType DdType, typename ValueType>
@ -63,55 +83,68 @@ namespace storm {
}
template<storm::dd::DdType DdType, typename ValueType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> GeneralSymbolicLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> solver = this->create();
solver->setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
return solver;
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverProblemFormat SymbolicLinearEquationSolverFactory<DdType, ValueType>::getEquationProblemFormat() const {
return this->create()->getEquationProblemFormat();
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverRequirements SymbolicLinearEquationSolverFactory<DdType, ValueType>::getRequirements() const {
return this->create()->getRequirements();
}
template<storm::dd::DdType DdType, typename ValueType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> GeneralSymbolicLinearEquationSolverFactory<DdType, ValueType>::create() const {
storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::CoreSettings>().getEquationSolver();
switch (equationSolver) {
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, ValueType>>();
break;
case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
case storm::solver::EquationSolverType::Native: return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>();
break;
default:
STORM_LOG_WARN("The selected equation solver is not available in the dd engine. Falling back to native solver.");
return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
STORM_LOG_INFO("The selected equation solver is not available in the dd engine. Falling back to native solver.");
return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, ValueType>>();
}
}
template<storm::dd::DdType DdType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create() const {
auto const& coreSettings = storm::settings::getModule<storm::settings::modules::CoreSettings>();
storm::solver::EquationSolverType equationSolver = coreSettings.getEquationSolver();
if (coreSettings.isEquationSolverSetFromDefaultValue() && equationSolver != storm::solver::EquationSolverType::Elimination) {
STORM_LOG_WARN("Selecting the elimination solver to guarantee exact results. If you want to override this, please explicitly specify a different equation solver.");
equationSolver = storm::solver::EquationSolverType::Elimination;
}
if (equationSolver != storm::solver::EquationSolverType::Elimination) {
STORM_LOG_WARN("The chosen equation solver does not guarantee precise results despite using exact arithmetic. Consider using the elimination solver instead.");
if (coreSettings.isEquationSolverSetFromDefaultValue() && equationSolver != storm::solver::EquationSolverType::Native && equationSolver != storm::solver::EquationSolverType::Elimination) {
STORM_LOG_INFO("Selecting the native solver to provide a method that guarantees exact results. If you want to override this, please explicitly specify a different equation solver.");
equationSolver = storm::solver::EquationSolverType::Native;
}
switch (equationSolver) {
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>();
break;
case storm::solver::EquationSolverType::Native:
return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
return std::make_unique<storm::solver::SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber>>();
break;
default:
STORM_LOG_WARN("The selected equation solver is not available in the dd engine. Falling back to elimination solver.");
return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
STORM_LOG_INFO("The selected equation solver is not available in the dd engine. Falling back to elimination solver.");
return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalNumber>>();
}
}
template<storm::dd::DdType DdType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create() const {
storm::solver::EquationSolverType equationSolver = storm::settings::getModule<storm::settings::modules::CoreSettings>().getEquationSolver();
switch (equationSolver) {
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
case storm::solver::EquationSolverType::Elimination: return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>();
break;
default:
STORM_LOG_WARN("The selected equation solver is not available in the DD setting. Falling back to elimination solver.");
return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs);
STORM_LOG_INFO("The selected equation solver is not available in the DD setting. Falling back to elimination solver.");
return std::make_unique<storm::solver::SymbolicEliminationLinearEquationSolver<DdType, storm::RationalFunction>>();
}
}

48
src/storm/solver/SymbolicLinearEquationSolver.h

@ -8,6 +8,10 @@
#include "storm/storage/dd/DdManager.h"
#include "storm/storage/dd/DdType.h"
#include "storm/solver/SymbolicEquationSolver.h"
#include "storm/solver/LinearEquationSolverProblemFormat.h"
#include "storm/solver/LinearEquationSolverRequirements.h"
#include "storm/adapters/RationalFunctionAdapter.h"
namespace storm {
@ -17,24 +21,19 @@ namespace storm {
template<storm::dd::DdType Type>
class Bdd;
}
namespace solver {
template<storm::dd::DdType DdType, typename ValueType = double>
class SymbolicLinearEquationSolverSettings {
public:
// Currently empty.
};
/*!
* An interface that represents an abstract symbolic linear equation solver. In addition to solving a system of
* linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
*/
template<storm::dd::DdType DdType, typename ValueType = double>
class SymbolicLinearEquationSolver {
class SymbolicLinearEquationSolver : public SymbolicEquationSolver<DdType, ValueType> {
public:
SymbolicLinearEquationSolver();
/*!
* Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
*
@ -85,17 +84,25 @@ namespace storm {
*/
virtual storm::dd::Add<DdType, ValueType> multiply(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const* b = nullptr, uint_fast64_t n = 1) const;
/*!
* Retrieves the format in which this solver expects to solve equations. If the solver expects the equation
* system format, it solves Ax = b. If it it expects a fixed point format, it solves Ax + b = x.
*/
virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const;
/*!
* Retrieves the requirements of the solver under the current settings. Note that these requirements only
* apply to solving linear equations and not to the matrix vector multiplications.
*/
virtual LinearEquationSolverRequirements getRequirements() const;
void setMatrix(storm::dd::Add<DdType, ValueType> const& newA);
virtual void setData(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs);
protected:
storm::dd::DdManager<DdType>& getDdManager() const;
// The matrix defining the coefficients of the linear equation system.
storm::dd::Add<DdType, ValueType> A;
// A BDD characterizing all rows of the equation system.
storm::dd::Bdd<DdType> allRows;
// The row variables.
std::set<storm::expressions::Variable> rowMetaVariables;
@ -103,15 +110,20 @@ namespace storm {
std::set<storm::expressions::Variable> columnMetaVariables;
// The pairs of meta variables used for renaming.
std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs;
std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> rowColumnMetaVariablePairs;
};
template<storm::dd::DdType DdType, typename ValueType>
class SymbolicLinearEquationSolverFactory {
public:
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const = 0;
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
LinearEquationSolverProblemFormat getEquationProblemFormat() const;
LinearEquationSolverRequirements getRequirements() const;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const = 0;
};
template<storm::dd::DdType DdType, typename ValueType>
@ -119,7 +131,7 @@ namespace storm {
public:
using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
};
template<storm::dd::DdType DdType>
@ -127,7 +139,7 @@ namespace storm {
public:
using SymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>::create;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalNumber>> create() const override;
};
template<storm::dd::DdType DdType>
@ -135,7 +147,7 @@ namespace storm {
public:
using SymbolicLinearEquationSolverFactory<DdType, storm::RationalFunction>::create;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, storm::RationalFunction>> create() const override;
};
} // namespace solver

33
src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp

@ -78,12 +78,12 @@ namespace storm {
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : settings(settings), requirementsChecked(false) {
SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : SymbolicEquationSolver<DdType, ValueType>(), settings(settings), requirementsChecked(false) {
// Intentionally left empty.
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::unique_ptr<SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : A(A), allRows(allRows), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity<ValueType>()), A.getDdManager().template getAddZero<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), settings(settings), requirementsChecked(false) {
SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::SymbolicMinMaxLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, storm::dd::Bdd<DdType> const& illegalMask, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::set<storm::expressions::Variable> const& choiceVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, std::unique_ptr<SymbolicLinearEquationSolverFactory<DdType, ValueType>>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings) : SymbolicEquationSolver<DdType, ValueType>(allRows), A(A), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity<ValueType>()), A.getDdManager().template getAddZero<ValueType>())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), settings(settings), requirementsChecked(false) {
// Intentionally left empty.
}
@ -104,7 +104,7 @@ namespace storm {
}
template<storm::dd::DdType DdType, typename ValueType>
typename SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::ValueIterationResult SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion) const {
typename SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::ValueIterationResult SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const {
// Set up local variables.
storm::dd::Add<DdType, ValueType> localX = x;
@ -112,7 +112,7 @@ namespace storm {
// Value iteration loop.
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && iterations < this->settings.getMaximalNumberOfIterations()) {
while (status == SolverStatus::InProgress && iterations < maximalIterations) {
// Compute tmp = A * x + b
storm::dd::Add<DdType, ValueType> localXAsColumn = localX.swapVariables(this->rowColumnMetaVariablePairs);
storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(localXAsColumn, this->columnMetaVariables);
@ -178,7 +178,7 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType>
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalX, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
// Storage for the rational sharpened vector.
storm::dd::Add<DdType, RationalType> sharpenedX;
@ -189,7 +189,7 @@ namespace storm {
ValueType precision = this->getSettings().getPrecision();
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion());
typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
++valueIterationInvocations;
STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << viResult.iterations << " iterations.");
@ -226,19 +226,17 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
return solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, x, b, x, b);
return solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, b, this->getLowerBounds(), b);
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::Add<DdType, storm::RationalNumber> rationalX = x.template toValueType<storm::RationalNumber>();
storm::dd::Add<DdType, storm::RationalNumber> rationalB = b.template toValueType<storm::RationalNumber>();
SymbolicMinMaxLinearEquationSolver<DdType, storm::RationalNumber> rationalSolver(this->A.template toValueType<storm::RationalNumber>(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique<GeneralSymbolicLinearEquationSolverFactory<DdType, storm::RationalNumber>>());
storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(dir, rationalSolver, *this, rationalX, rationalB, x, b);
storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(dir, rationalSolver, *this, rationalB, this->getLowerBounds(), b);
return rationalResult.template toValueType<ValueType>();
}
@ -248,17 +246,18 @@ namespace storm {
// First try to find a solution using the imprecise value type.
storm::dd::Add<DdType, ValueType> rationalResult;
storm::dd::Add<DdType, ImpreciseType> impreciseX;
try {
storm::dd::Add<DdType, ImpreciseType> impreciseX = x.template toValueType<ImpreciseType>();
impreciseX = this->getLowerBounds().template toValueType<ImpreciseType>();
storm::dd::Add<DdType, ImpreciseType> impreciseB = b.template toValueType<ImpreciseType>();
SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> impreciseSolver(this->A.template toValueType<ImpreciseType>(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique<GeneralSymbolicLinearEquationSolverFactory<DdType, ImpreciseType>>());
rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(dir, *this, impreciseSolver, x, b, impreciseX, impreciseB);
rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(dir, *this, impreciseSolver, b, impreciseX, impreciseB);
} catch (storm::exceptions::PrecisionExceededException const& e) {
STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
// Fall back to precise value type if the precision of the imprecise value type was exceeded.
rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, x, b, x, b);
rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(dir, *this, *this, b, impreciseX.template toValueType<ValueType>(), b);
}
return rationalResult.template toValueType<ValueType>();
}
@ -271,14 +270,16 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
// Set up the environment.
storm::dd::Add<DdType, ValueType> localX = x;
storm::dd::Add<DdType, ValueType> localX;
// If we were given an initial scheduler, we take its solution as the starting point.
if (this->hasInitialScheduler()) {
localX = solveEquationsWithScheduler(this->getInitialScheduler(), x, b);
} else {
localX = this->getLowerBounds();
}
ValueIterationResult viResult = performValueIteration(dir, localX, b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion());
ValueIterationResult viResult = performValueIteration(dir, localX, b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
if (viResult.status == SolverStatus::Converged) {
STORM_LOG_INFO("Iterative solver (value iteration) converged in " << viResult.iterations << " iterations.");
@ -419,12 +420,14 @@ namespace storm {
}
}
} else if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::ValueIteration) {
requirements.requireLowerBounds();
if (equationSystemType == EquationSystemType::ReachabilityRewards) {
if (!direction || direction.get() == OptimizationDirection::Minimize) {
requirements.requireValidInitialScheduler();
}
}
} else if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
requirements.requireLowerBounds();
if (equationSystemType == EquationSystemType::ReachabilityRewards) {
if (!direction || direction.get() == OptimizationDirection::Minimize) {
requirements.requireNoEndComponents();

12
src/storm/solver/SymbolicMinMaxLinearEquationSolver.h

@ -6,6 +6,7 @@
#include <vector>
#include <boost/optional.hpp>
#include "storm/solver/SymbolicEquationSolver.h"
#include "storm/solver/OptimizationDirection.h"
#include "storm/solver/SymbolicLinearEquationSolver.h"
#include "storm/solver/EquationSystemType.h"
@ -58,7 +59,7 @@ namespace storm {
* linear equations, the functionality to repeatedly multiply a matrix with a given vector is provided.
*/
template<storm::dd::DdType DdType, typename ValueType>
class SymbolicMinMaxLinearEquationSolver {
class SymbolicMinMaxLinearEquationSolver : public SymbolicEquationSolver<DdType, ValueType> {
public:
SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings<ValueType> const& settings = SymbolicMinMaxLinearEquationSolverSettings<ValueType>());
@ -140,7 +141,7 @@ namespace storm {
bool isRequirementsCheckedSet() const;
/*!
* Determines whether the given vector x satisfies x = Ax + b.
* Determines whether the given vector x satisfies x = min/max Ax + b.
*/
bool isSolution(OptimizationDirection dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
@ -156,7 +157,7 @@ namespace storm {
static storm::dd::Add<DdType, RationalType> sharpen(OptimizationDirection dir, uint64_t precision, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution);
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalX, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
template<typename ImpreciseType>
@ -177,15 +178,12 @@ namespace storm {
storm::dd::Add<DdType, ValueType> values;
};
ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion) const;
ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const;
protected:
// The matrix defining the coefficients of the linear equation system.
storm::dd::Add<DdType, ValueType> A;
// A BDD characterizing all rows of the equation system.
storm::dd::Bdd<DdType> allRows;
// A BDD characterizing the illegal choices.
storm::dd::Bdd<DdType> illegalMask;

236
src/storm/solver/SymbolicNativeLinearEquationSolver.cpp

@ -9,6 +9,11 @@
#include "storm/settings/SettingsManager.h"
#include "storm/settings/modules/NativeEquationSolverSettings.h"
#include "storm/utility/KwekMehlhorn.h"
#include "storm/utility/macros.h"
#include "storm/exceptions/NotSupportedException.h"
#include "storm/exceptions/PrecisionExceededException.h"
namespace storm {
namespace solver {
@ -18,6 +23,27 @@ namespace storm {
storm::settings::modules::NativeEquationSolverSettings const& settings = storm::settings::getModule<storm::settings::modules::NativeEquationSolverSettings>();
// Get appropriate settings.
switch (settings.getLinearEquationSystemMethod()) {
case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Power:
this->method = SolutionMethod::Power;
break;
case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::RationalSearch:
this->method = SolutionMethod::RationalSearch;
break;
case storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::Jacobi:
default:
this->method = SolutionMethod::Jacobi;
break;
}
// Adjust the method if none was specified and we are using rational numbers.
if (std::is_same<ValueType, storm::RationalNumber>::value) {
if (settings.isLinearEquationSystemTechniqueSetFromDefaultValue() && this->method != SolutionMethod::RationalSearch) {
STORM_LOG_INFO("Selecting the rational search as the solution technique to guarantee exact results. If you want to override this, please explicitly specify a different method.");
this->method = SolutionMethod::RationalSearch;
}
}
maximalNumberOfIterations = settings.getMaximalIterationCount();
precision = storm::utility::convertNumber<ValueType>(settings.getPrecision());
relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative;
@ -38,6 +64,11 @@ namespace storm {
this->relative = value;
}
template<typename ValueType>
void SymbolicNativeLinearEquationSolverSettings<ValueType>::setSolutionMethod(SolutionMethod const& method) {
this->method = method;
}
template<typename ValueType>
ValueType SymbolicNativeLinearEquationSolverSettings<ValueType>::getPrecision() const {
return precision;
@ -53,6 +84,16 @@ namespace storm {
return relative;
}
template<typename ValueType>
typename SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod SymbolicNativeLinearEquationSolverSettings<ValueType>::getSolutionMethod() const {
return this->method;
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicNativeLinearEquationSolver<DdType, ValueType>::SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings) : SymbolicLinearEquationSolver<DdType, ValueType>(), settings(settings) {
// Intentionally left empty.
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicNativeLinearEquationSolver<DdType, ValueType>::SymbolicNativeLinearEquationSolver(storm::dd::Add<DdType, ValueType> const& A, storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs, SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings) : SymbolicNativeLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings) {
this->setMatrix(A);
@ -65,6 +106,18 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
return solveEquationsJacobi(x, b);
} else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Power) {
return solveEquationsPower(x, b);
} else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
return solveEquationsRationalSearch(x, b);
}
STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The selected solution technique is not supported.");
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsJacobi(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::DdManager<DdType>& manager = this->getDdManager();
// Start by computing the Jacobi decomposition of the matrix A.
@ -97,22 +150,197 @@ namespace storm {
}
if (converged) {
STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations.");
STORM_LOG_INFO("Iterative solver (jacobi) converged in " << iterationCount << " iterations.");
} else {
STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations.");
STORM_LOG_WARN("Iterative solver (jacobi) did not converge in " << iterationCount << " iterations.");
}
return xCopy;
}
template<storm::dd::DdType DdType, typename ValueType>
typename SymbolicNativeLinearEquationSolver<DdType, ValueType>::PowerIterationResult SymbolicNativeLinearEquationSolver<DdType, ValueType>::performPowerIteration(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const {
// Set up additional environment variables.
storm::dd::Add<DdType, ValueType> currentX = x;
uint_fast64_t iterations = 0;
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && iterations < maximalIterations) {
storm::dd::Add<DdType, ValueType> currentXAsColumn = currentX.swapVariables(this->rowColumnMetaVariablePairs);
storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(currentXAsColumn, this->columnMetaVariables) + b;
// Now check if the process already converged within our precision.
if (tmp.equalModuloPrecision(currentX, precision, relativeTerminationCriterion)) {
status = SolverStatus::Converged;
}
// Set up next iteration.
++iterations;
currentX = tmp;
}
return PowerIterationResult(status, iterations, currentX);
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsPower(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
PowerIterationResult result = performPowerIteration(this->getLowerBounds(), b, this->getSettings().getPrecision(), this->getSettings().getRelativeTerminationCriterion(), this->getSettings().getMaximalNumberOfIterations());
if (result.status == SolverStatus::Converged) {
STORM_LOG_INFO("Iterative solver (power iteration) converged in " << result.iterations << " iterations.");
} else {
STORM_LOG_WARN("Iterative solver (power iteration) did not converge in " << result.iterations << " iterations.");
}
return result.values;
}
template<storm::dd::DdType DdType, typename ValueType>
bool SymbolicNativeLinearEquationSolver<DdType, ValueType>::isSolutionFixedPoint(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::Add<DdType, ValueType> xAsColumn = x.swapVariables(this->rowColumnMetaVariablePairs);
storm::dd::Add<DdType, ValueType> tmp = this->A.multiplyMatrix(xAsColumn, this->columnMetaVariables);
tmp += b;
return x == tmp;
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution) {
storm::dd::Add<DdType, RationalType> sharpenedX;
for (uint64_t p = 1; p < precision; ++p) {
sharpenedX = x.sharpenKwekMehlhorn(precision);
isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
if (isSolution) {
break;
}
}
return sharpenedX;
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
// Storage for the rational sharpened vector.
storm::dd::Add<DdType, RationalType> sharpenedX;
// The actual rational search.
uint64_t overallIterations = 0;
uint64_t powerIterationInvocations = 0;
ValueType precision = this->getSettings().getPrecision();
SolverStatus status = SolverStatus::InProgress;
while (status == SolverStatus::InProgress && overallIterations < this->getSettings().getMaximalNumberOfIterations()) {
typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), this->getSettings().getRelativeTerminationCriterion(), this->settings.getMaximalNumberOfIterations());
++powerIterationInvocations;
STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
// Count the iterations.
overallIterations += result.iterations;
// Compute maximal precision until which to sharpen.
uint64_t p = storm::utility::convertNumber<uint64_t>(storm::utility::ceil(storm::utility::log10<ValueType>(storm::utility::one<ValueType>() / precision)));
bool isSolution = false;
sharpenedX = sharpen<RationalType, ImpreciseType>(p, rationalSolver, result.values, rationalB, isSolution);
if (isSolution) {
status = SolverStatus::Converged;
} else {
precision = precision / 100;
}
}
if (status == SolverStatus::InProgress) {
status = SolverStatus::MaximalIterationsExceeded;
}
if (status == SolverStatus::Converged) {
STORM_LOG_INFO("Iterative solver (rational search) converged in " << overallIterations << " iterations.");
} else {
STORM_LOG_WARN("Iterative solver (rational search) did not converge in " << overallIterations << " iterations.");
}
return sharpenedX;
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
return solveEquationsRationalSearchHelper<ValueType, ValueType>(*this, *this, b, this->getLowerBounds(), b);
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
storm::dd::Add<DdType, storm::RationalNumber> rationalB = b.template toValueType<storm::RationalNumber>();
SymbolicNativeLinearEquationSolver<DdType, storm::RationalNumber> rationalSolver(this->A.template toValueType<storm::RationalNumber>(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
storm::dd::Add<DdType, storm::RationalNumber> rationalResult = solveEquationsRationalSearchHelper<storm::RationalNumber, ImpreciseType>(rationalSolver, *this, rationalB, this->getLowerBounds(), b);
return rationalResult.template toValueType<ValueType>();
}
template<storm::dd::DdType DdType, typename ValueType>
template<typename ImpreciseType>
typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, storm::dd::Add<DdType, ValueType>>::type SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
// First try to find a solution using the imprecise value type.
storm::dd::Add<DdType, ValueType> rationalResult;
storm::dd::Add<DdType, ImpreciseType> impreciseX;
try {
impreciseX = this->getLowerBounds().template toValueType<ImpreciseType>();
storm::dd::Add<DdType, ImpreciseType> impreciseB = b.template toValueType<ImpreciseType>();
SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> impreciseSolver(this->A.template toValueType<ImpreciseType>(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs);
rationalResult = solveEquationsRationalSearchHelper<ValueType, ImpreciseType>(*this, impreciseSolver, b, impreciseX, impreciseB);
} catch (storm::exceptions::PrecisionExceededException const& e) {
STORM_LOG_WARN("Precision of value type was exceeded, trying to recover by switching to rational arithmetic.");
// Fall back to precise value type if the precision of the imprecise value type was exceeded.
rationalResult = solveEquationsRationalSearchHelper<ValueType, ValueType>(*this, *this, b, impreciseX.template toValueType<ValueType>(), b);
}
return rationalResult.template toValueType<ValueType>();
}
template<storm::dd::DdType DdType, typename ValueType>
storm::dd::Add<DdType, ValueType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearch(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const {
return solveEquationsRationalSearchHelper<double>(x, b);
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverProblemFormat SymbolicNativeLinearEquationSolver<DdType, ValueType>::getEquationProblemFormat() const {
if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Jacobi) {
return LinearEquationSolverProblemFormat::EquationSystem;
}
return LinearEquationSolverProblemFormat::FixedPointSystem;
}
template<storm::dd::DdType DdType, typename ValueType>
LinearEquationSolverRequirements SymbolicNativeLinearEquationSolver<DdType, ValueType>::getRequirements() const {
LinearEquationSolverRequirements requirements;
if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::Power || this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings<ValueType>::SolutionMethod::RationalSearch) {
requirements.requireLowerBounds();
}
return requirements;
}
template<storm::dd::DdType DdType, typename ValueType>
SymbolicNativeLinearEquationSolverSettings<ValueType> const& SymbolicNativeLinearEquationSolver<DdType, ValueType>::getSettings() const {
return settings;
}
template<storm::dd::DdType DdType, typename ValueType>
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicNativeLinearEquationSolverFactory<DdType, ValueType>::create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const {
return std::make_unique<SymbolicNativeLinearEquationSolver<DdType, ValueType>>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings);
std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> SymbolicNativeLinearEquationSolverFactory<DdType, ValueType>::create() const {
return std::make_unique<SymbolicNativeLinearEquationSolver<DdType, ValueType>>(settings);
}
template<storm::dd::DdType DdType, typename ValueType>

62
src/storm/solver/SymbolicNativeLinearEquationSolver.h

@ -1,6 +1,9 @@
#pragma once
#include "storm/solver/SymbolicLinearEquationSolver.h"
#include "storm/solver/SolverStatus.h"
#include "storm/utility/NumberTraits.h"
namespace storm {
namespace solver {
@ -8,17 +11,26 @@ namespace storm {
template<typename ValueType>
class SymbolicNativeLinearEquationSolverSettings {
public:
enum class SolutionMethod {
Jacobi, Power, RationalSearch
};
SymbolicNativeLinearEquationSolverSettings();
void setPrecision(ValueType precision);
void setMaximalNumberOfIterations(uint64_t maximalNumberOfIterations);
void setRelativeTerminationCriterion(bool value);
void setSolutionMethod(SolutionMethod const& method);
ValueType getPrecision() const;
uint64_t getMaximalNumberOfIterations() const;
bool getRelativeTerminationCriterion() const;
SolutionMethod getSolutionMethod() const;
private:
// The selected solution method.
SolutionMethod method;
// The required precision for the iterative methods.
ValueType precision;
@ -37,6 +49,13 @@ namespace storm {
template<storm::dd::DdType DdType, typename ValueType = double>
class SymbolicNativeLinearEquationSolver : public SymbolicLinearEquationSolver<DdType, ValueType> {
public:
/*!
* Constructs a symbolic linear equation solver.
*
* @param settings The settings to use.
*/
SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings<ValueType> const& settings = SymbolicNativeLinearEquationSolverSettings<ValueType>());
/*!
* Constructs a symbolic linear equation solver with the given meta variable sets and pairs.
*
@ -73,11 +92,50 @@ namespace storm {
* @param b The right-hand side of the equation system. Its length must be equal to the number of rows of A.
* @return The solution of the equation system.
*/
virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
virtual storm::dd::Add<DdType, ValueType> solveEquations(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const override;
virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override;
virtual LinearEquationSolverRequirements getRequirements() const override;
SymbolicNativeLinearEquationSolverSettings<ValueType> const& getSettings() const;
private:
storm::dd::Add<DdType, ValueType> solveEquationsJacobi(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
storm::dd::Add<DdType, ValueType> solveEquationsPower(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
storm::dd::Add<DdType, ValueType> solveEquationsRationalSearch(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
/*!
* Determines whether the given vector x satisfies x = Ax + b.
*/
bool isSolutionFixedPoint(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
template<typename RationalType, typename ImpreciseType>
static storm::dd::Add<DdType, RationalType> sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, RationalType> const& rationalB, bool& isSolution);
template<typename RationalType, typename ImpreciseType>
storm::dd::Add<DdType, RationalType> solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const;
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
template<typename ImpreciseType>
typename std::enable_if<std::is_same<ValueType, ImpreciseType>::value && !storm::NumberTraits<ValueType>::IsExact, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
template<typename ImpreciseType>
typename std::enable_if<!std::is_same<ValueType, ImpreciseType>::value, storm::dd::Add<DdType, ValueType>>::type solveEquationsRationalSearchHelper(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b) const;
template<storm::dd::DdType DdTypePrime, typename ValueTypePrime>
friend class SymbolicNativeLinearEquationSolver;
struct PowerIterationResult {
PowerIterationResult(SolverStatus status, uint64_t iterations, storm::dd::Add<DdType, ValueType> const& values) : status(status), iterations(iterations), values(values) {
// Intentionally left empty.
}
SolverStatus status;
uint64_t iterations;
storm::dd::Add<DdType, ValueType> values;
};
PowerIterationResult performPowerIteration(storm::dd::Add<DdType, ValueType> const& x, storm::dd::Add<DdType, ValueType> const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const;
// The settings to use.
SymbolicNativeLinearEquationSolverSettings<ValueType> settings;
};
@ -87,7 +145,7 @@ namespace storm {
public:
using SymbolicLinearEquationSolverFactory<DdType, ValueType>::create;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create(storm::dd::Bdd<DdType> const& allRows, std::set<storm::expressions::Variable> const& rowMetaVariables, std::set<storm::expressions::Variable> const& columnMetaVariables, std::vector<std::pair<storm::expressions::Variable, storm::expressions::Variable>> const& rowColumnMetaVariablePairs) const override;
virtual std::unique_ptr<storm::solver::SymbolicLinearEquationSolver<DdType, ValueType>> create() const override;
SymbolicNativeLinearEquationSolverSettings<ValueType>& getSettings();
SymbolicNativeLinearEquationSolverSettings<ValueType> const& getSettings() const;

6
src/storm/storage/dd/bisimulation/QuotientExtractor.cpp

@ -984,7 +984,11 @@ namespace storm {
end = std::chrono::high_resolution_clock::now();
// Check quotient matrix for sanity.
STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>()).isZero(), "Illegal entries in quotient matrix.");
if (std::is_same<ValueType, storm::RationalNumber>::value) {
STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>()).isZero(), "Illegal entries in quotient matrix.");
} else {
STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(1e-6)).isZero(), "Illegal entries in quotient matrix.");
}
STORM_LOG_ASSERT(quotientTransitionMatrix.sumAbstract(blockPrimeVariableSet).equalModuloPrecision(quotientTransitionMatrix.notZero().existsAbstract(blockPrimeVariableSet).template toAdd<ValueType>(), ValueType(1e-6)), "Illegal non-probabilistic matrix.");
STORM_LOG_TRACE("Quotient transition matrix extracted in " << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() << "ms.");

5
src/storm/utility/constants.cpp

@ -708,6 +708,11 @@ namespace storm {
return number;
}
template<>
RationalFunction convertNumber(std::string const& number){
return RationalFunction(convertNumber<RationalFunctionCoefficient>(number));
}
template<>
RationalFunction convertNumber(storm::storage::sparse::state_type const& number) {
return RationalFunction(convertNumber<RationalFunctionCoefficient>(number));

Loading…
Cancel
Save