diff --git a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h b/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h index 449a98c3e..d249c5564 100644 --- a/resources/3rdparty/sylvan/src/sylvan_mtbdd_storm.h +++ b/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 } diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h b/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h index 888a80d51..b6fdfeed2 100644 --- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_function.h +++ b/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*) diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c index 22f9ef053..261156ba7 100644 --- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.c +++ b/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; diff --git a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h b/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h index 227a9831d..070f03e1a 100644 --- a/resources/3rdparty/sylvan/src/sylvan_storm_rational_number.h +++ b/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) diff --git a/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp index 86f624459..ffc602aaa 100644 --- a/src/storm/modelchecker/prctl/helper/SymbolicDtmcPrctlHelper.cpp +++ b/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> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + solver->setBounds(storm::utility::zero(), storm::utility::one()); storm::dd::Add result = solver->solveEquations(model.getManager().template getAddZero(), subvector); return statesWithProbability01.second.template toAdd() + 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> solver = linearEquationSolverFactory.create(submatrix, maybeStates, model.getRowVariables(), model.getColumnVariables(), model.getRowColumnMetaVariablePairs()); + solver->setLowerBound(storm::utility::zero()); storm::dd::Add result = solver->solveEquations(model.getManager().template getAddZero(), subvector); return infinityStates.ite(model.getManager().getConstant(storm::utility::infinity()), result); diff --git a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp index bca6bec0f..518f50acc 100644 --- a/src/storm/modelchecker/prctl/helper/SymbolicMdpPrctlHelper.cpp +++ b/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(), storm::utility::one()); solver->setRequirementsChecked(); storm::dd::Add result = solver->solveEquations(dir, model.getManager().template getAddZero(), 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()); solver->setRequirementsChecked(); storm::dd::Add result = solver->solveEquations(dir, model.getManager().template getAddZero(), subvector); diff --git a/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp b/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp index 6f5597562..b73190bab 100644 --- a/src/storm/modelchecker/results/SymbolicQuantitativeCheckResult.cpp +++ b/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()) && reachableStates).getNonZeroCount() << std::endl; - std::cout << "states one: " << (values.equals(values.getDdManager().template getAddOne()) && 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()) && 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) { diff --git a/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp b/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp index 1cdc67743..7443ecdcf 100644 --- a/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp +++ b/src/storm/solver/SymbolicEliminationLinearEquationSolver.cpp @@ -10,6 +10,11 @@ namespace storm { namespace solver { + template + SymbolicEliminationLinearEquationSolver::SymbolicEliminationLinearEquationSolver() : SymbolicLinearEquationSolver() { + // Intentionally left empty. + } + template SymbolicEliminationLinearEquationSolver::SymbolicEliminationLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : SymbolicEliminationLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) { this->setMatrix(A); @@ -17,6 +22,64 @@ namespace storm { template SymbolicEliminationLinearEquationSolver::SymbolicEliminationLinearEquationSolver(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : SymbolicLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) { + this->createInternalData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + } + + template + storm::dd::Add SymbolicEliminationLinearEquationSolver::solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const { + storm::dd::DdManager& ddManager = x.getDdManager(); + + // Build diagonal BDD over new meta variables. + storm::dd::Bdd diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs); + diagonal &= this->getAllRows(); + diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs); + + storm::dd::Add rowsAdd = this->getAllRows().swapVariables(rowRowMetaVariablePairs).template toAdd(); + storm::dd::Add diagonalAdd = diagonal.template toAdd(); + + // Move the matrix to the new meta variables. + storm::dd::Add matrix = this->A.swapVariables(oldNewMetaVariablePairs); + + // Initialize solution over the new meta variables. + storm::dd::Add 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 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(), 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 + LinearEquationSolverProblemFormat SymbolicEliminationLinearEquationSolver::getEquationProblemFormat() const { + return LinearEquationSolverProblemFormat::FixedPointSystem; + } + + template + void SymbolicEliminationLinearEquationSolver::createInternalData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) { storm::dd::DdManager& 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::Add SymbolicEliminationLinearEquationSolver::solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const { - storm::dd::DdManager& ddManager = x.getDdManager(); - - // Build diagonal BDD over new meta variables. - storm::dd::Bdd diagonal = storm::utility::dd::getRowColumnDiagonal(ddManager, this->rowColumnMetaVariablePairs); - diagonal &= this->allRows; - diagonal = diagonal.swapVariables(this->oldNewMetaVariablePairs); - - storm::dd::Add rowsAdd = this->allRows.swapVariables(rowRowMetaVariablePairs).template toAdd(); - storm::dd::Add diagonalAdd = diagonal.template toAdd(); - - // Revert the conversion to an equation system and move it to the new meta variables. - storm::dd::Add matrix = diagonalAdd - this->A.swapVariables(oldNewMetaVariablePairs); + void SymbolicEliminationLinearEquationSolver::setData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) { - // Initialize solution over the new meta variables. - storm::dd::Add solution = b.swapVariables(oldNewMetaVariablePairs); + // Call superclass function. + SymbolicLinearEquationSolver::setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); - // As long as there are transitions, we eliminate them. - uint64_t iterations = 0; - while (!matrix.isZero()) { - // Determine inverse loop probabilies. - storm::dd::Add 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(), 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); + // Now create new variables as needed. + this->createInternalData(this->getAllRows(), this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs); } - + template - std::unique_ptr> SymbolicEliminationLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { - return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + std::unique_ptr> SymbolicEliminationLinearEquationSolverFactory::create() const { + return std::make_unique>(); } template diff --git a/src/storm/solver/SymbolicEliminationLinearEquationSolver.h b/src/storm/solver/SymbolicEliminationLinearEquationSolver.h index 45b0f0b12..2e1ce5a4b 100644 --- a/src/storm/solver/SymbolicEliminationLinearEquationSolver.h +++ b/src/storm/solver/SymbolicEliminationLinearEquationSolver.h @@ -14,13 +14,26 @@ namespace storm { template class SymbolicEliminationLinearEquationSolver : public SymbolicLinearEquationSolver { public: + /*! + * Constructs a symbolic linear equation solver. + * + * @param settings The settings to use. + */ + SymbolicEliminationLinearEquationSolver(); + SymbolicEliminationLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs); SymbolicEliminationLinearEquationSolver(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs); virtual storm::dd::Add solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const override; + virtual void setData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) override; + + virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override; + private: + void createInternalData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs); + std::vector> oldToNewMapping; std::set newRowVariables; std::set newColumnVariables; @@ -38,7 +51,7 @@ namespace storm { public: using SymbolicLinearEquationSolverFactory::create; - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const override; + virtual std::unique_ptr> create() const override; SymbolicEliminationLinearEquationSolverSettings& getSettings(); SymbolicEliminationLinearEquationSolverSettings const& getSettings() const; diff --git a/src/storm/solver/SymbolicEquationSolver.cpp b/src/storm/solver/SymbolicEquationSolver.cpp new file mode 100644 index 000000000..a76f94873 --- /dev/null +++ b/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 + SymbolicEquationSolver::SymbolicEquationSolver(storm::dd::Bdd const& allRows) : allRows(allRows) { + // Intentionally left empty. + } + + template + storm::dd::DdManager& SymbolicEquationSolver::getDdManager() const { + return this->allRows.getDdManager(); + } + + template + void SymbolicEquationSolver::setAllRows(storm::dd::Bdd const& allRows) { + this->allRows = allRows; + } + + template + storm::dd::Bdd const& SymbolicEquationSolver::getAllRows() const { + return this->allRows; + } + + template + void SymbolicEquationSolver::setLowerBounds(storm::dd::Add const& lowerBounds) { + this->lowerBounds = lowerBounds; + } + + template + void SymbolicEquationSolver::setLowerBound(ValueType const& lowerBound) { + this->lowerBound = lowerBound; + } + + template + void SymbolicEquationSolver::setUpperBounds(storm::dd::Add const& upperBounds) { + this->upperBounds = upperBounds; + } + + template + void SymbolicEquationSolver::setUpperBound(ValueType const& upperBound) { + this->upperBound = upperBound; + } + + template + void SymbolicEquationSolver::setBounds(ValueType const& lowerBound, ValueType const& upperBound) { + setLowerBound(lowerBound); + setUpperBound(upperBound); + } + + template + void SymbolicEquationSolver::setBounds(storm::dd::Add const& lowerBounds, storm::dd::Add const& upperBounds) { + setLowerBounds(lowerBounds); + setUpperBounds(upperBounds); + } + + template + storm::dd::Add SymbolicEquationSolver::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()); + } + } + + template + storm::dd::Add SymbolicEquationSolver::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()); + } + } + + template class SymbolicEquationSolver; + template class SymbolicEquationSolver; + +#ifdef STORM_HAVE_CARL + template class SymbolicEquationSolver; + template class SymbolicEquationSolver; + + template class SymbolicEquationSolver; +#endif + } +} diff --git a/src/storm/solver/SymbolicEquationSolver.h b/src/storm/solver/SymbolicEquationSolver.h new file mode 100644 index 000000000..2e6322cfc --- /dev/null +++ b/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 + class SymbolicEquationSolver { + public: + SymbolicEquationSolver() = default; + SymbolicEquationSolver(storm::dd::Bdd const& allRows); + + void setLowerBounds(storm::dd::Add const& lowerBounds); + void setLowerBound(ValueType const& lowerBound); + void setUpperBounds(storm::dd::Add const& upperBounds); + void setUpperBound(ValueType const& lowerBound); + void setBounds(ValueType const& lowerBound, ValueType const& upperBound); + void setBounds(storm::dd::Add const& lowerBounds, storm::dd::Add const& upperBounds); + + /*! + * Retrieves a vector of lower bounds for all values (if any lower bounds are known). + */ + storm::dd::Add getLowerBounds() const; + + /*! + * Retrieves a vector of upper bounds for all values (if any lower bounds are known). + */ + storm::dd::Add getUpperBounds() const; + + protected: + storm::dd::DdManager& getDdManager() const; + + void setAllRows(storm::dd::Bdd const& allRows); + storm::dd::Bdd const& getAllRows() const; + + // The relevant rows to this equation solver. + storm::dd::Bdd allRows; + + private: + // Lower bounds (if given). + boost::optional> lowerBounds; + boost::optional lowerBound; + + // Upper bounds (if given). + boost::optional> upperBounds; + boost::optional upperBound; + }; + + } +} diff --git a/src/storm/solver/SymbolicLinearEquationSolver.cpp b/src/storm/solver/SymbolicLinearEquationSolver.cpp index 41620bbf8..fb6e6059e 100644 --- a/src/storm/solver/SymbolicLinearEquationSolver.cpp +++ b/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 + SymbolicLinearEquationSolver::SymbolicLinearEquationSolver() { + // Intentionally left empty. + } + template SymbolicLinearEquationSolver::SymbolicLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : SymbolicLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs) { this->setMatrix(A); } template - SymbolicLinearEquationSolver::SymbolicLinearEquationSolver(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : allRows(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) { + SymbolicLinearEquationSolver::SymbolicLinearEquationSolver(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) : SymbolicEquationSolver(allRows), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs) { // Intentionally left empty. } @@ -45,16 +51,30 @@ namespace storm { return xCopy; } + template + LinearEquationSolverProblemFormat SymbolicLinearEquationSolver::getEquationProblemFormat() const { + return LinearEquationSolverProblemFormat::EquationSystem; + } + + template + LinearEquationSolverRequirements SymbolicLinearEquationSolver::getRequirements() const { + // Return empty requirements by default. + return LinearEquationSolverRequirements(); + } + template void SymbolicLinearEquationSolver::setMatrix(storm::dd::Add const& newA) { this->A = newA; } template - storm::dd::DdManager& SymbolicLinearEquationSolver::getDdManager() const { - return this->allRows.getDdManager(); + void SymbolicLinearEquationSolver::setData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) { + this->setAllRows(allRows); + this->rowMetaVariables = rowMetaVariables; + this->columnMetaVariables = columnMetaVariables; + this->rowColumnMetaVariablePairs = rowColumnMetaVariablePairs; } - + template std::unique_ptr> SymbolicLinearEquationSolverFactory::create(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { std::unique_ptr> solver = this->create(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); @@ -63,55 +83,68 @@ namespace storm { } template - std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { + std::unique_ptr> SymbolicLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { + std::unique_ptr> solver = this->create(); + solver->setData(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + return solver; + } + + template + LinearEquationSolverProblemFormat SymbolicLinearEquationSolverFactory::getEquationProblemFormat() const { + return this->create()->getEquationProblemFormat(); + } + + template + LinearEquationSolverRequirements SymbolicLinearEquationSolverFactory::getRequirements() const { + return this->create()->getRequirements(); + } + + template + std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create() const { storm::solver::EquationSolverType equationSolver = storm::settings::getModule().getEquationSolver(); switch (equationSolver) { - case storm::solver::EquationSolverType::Elimination: return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + case storm::solver::EquationSolverType::Elimination: return std::make_unique>(); break; - case storm::solver::EquationSolverType::Native: return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + case storm::solver::EquationSolverType::Native: return std::make_unique>(); 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>(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>(); } } template - std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { + std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create() const { auto const& coreSettings = storm::settings::getModule(); 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 (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; } - 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."); - } - switch (equationSolver) { - case storm::solver::EquationSolverType::Elimination: return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + case storm::solver::EquationSolverType::Elimination: return std::make_unique>(); break; case storm::solver::EquationSolverType::Native: - return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + return std::make_unique>(); 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>(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>(); } } template - std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { + std::unique_ptr> GeneralSymbolicLinearEquationSolverFactory::create() const { storm::solver::EquationSolverType equationSolver = storm::settings::getModule().getEquationSolver(); switch (equationSolver) { - case storm::solver::EquationSolverType::Elimination: return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs); + case storm::solver::EquationSolverType::Elimination: return std::make_unique>(); 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>(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>(); } } diff --git a/src/storm/solver/SymbolicLinearEquationSolver.h b/src/storm/solver/SymbolicLinearEquationSolver.h index 9b3c3b490..d8d2e00ca 100644 --- a/src/storm/solver/SymbolicLinearEquationSolver.h +++ b/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 class Bdd; - } namespace solver { - template - 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 - class SymbolicLinearEquationSolver { + class SymbolicLinearEquationSolver : public SymbolicEquationSolver { 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 multiply(storm::dd::Add const& x, storm::dd::Add 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 const& newA); + virtual void setData(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs); protected: - storm::dd::DdManager& getDdManager() const; - // The matrix defining the coefficients of the linear equation system. storm::dd::Add A; - // A BDD characterizing all rows of the equation system. - storm::dd::Bdd allRows; - // The row variables. std::set rowMetaVariables; @@ -103,15 +110,20 @@ namespace storm { std::set columnMetaVariables; // The pairs of meta variables used for renaming. - std::vector> const& rowColumnMetaVariablePairs; + std::vector> rowColumnMetaVariablePairs; }; template class SymbolicLinearEquationSolverFactory { public: - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const = 0; + std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const; std::unique_ptr> create(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const; + + LinearEquationSolverProblemFormat getEquationProblemFormat() const; + LinearEquationSolverRequirements getRequirements() const; + + virtual std::unique_ptr> create() const = 0; }; template @@ -119,7 +131,7 @@ namespace storm { public: using SymbolicLinearEquationSolverFactory::create; - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const; + virtual std::unique_ptr> create() const override; }; template @@ -127,7 +139,7 @@ namespace storm { public: using SymbolicLinearEquationSolverFactory::create; - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const; + virtual std::unique_ptr> create() const override; }; template @@ -135,7 +147,7 @@ namespace storm { public: using SymbolicLinearEquationSolverFactory::create; - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const; + virtual std::unique_ptr> create() const override; }; } // namespace solver diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp index 490f341e3..4cb05ee90 100644 --- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp @@ -78,12 +78,12 @@ namespace storm { } template - SymbolicMinMaxLinearEquationSolver::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings const& settings) : settings(settings), requirementsChecked(false) { + SymbolicMinMaxLinearEquationSolver::SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings const& settings) : SymbolicEquationSolver(), settings(settings), requirementsChecked(false) { // Intentionally left empty. } template - SymbolicMinMaxLinearEquationSolver::SymbolicMinMaxLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, storm::dd::Bdd const& illegalMask, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& choiceVariables, std::vector> const& rowColumnMetaVariablePairs, std::unique_ptr>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings const& settings) : A(A), allRows(allRows), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity()), A.getDdManager().template getAddZero())), rowMetaVariables(rowMetaVariables), columnMetaVariables(columnMetaVariables), choiceVariables(choiceVariables), rowColumnMetaVariablePairs(rowColumnMetaVariablePairs), linearEquationSolverFactory(std::move(linearEquationSolverFactory)), settings(settings), requirementsChecked(false) { + SymbolicMinMaxLinearEquationSolver::SymbolicMinMaxLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, storm::dd::Bdd const& illegalMask, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& choiceVariables, std::vector> const& rowColumnMetaVariablePairs, std::unique_ptr>&& linearEquationSolverFactory, SymbolicMinMaxLinearEquationSolverSettings const& settings) : SymbolicEquationSolver(allRows), A(A), illegalMask(illegalMask), illegalMaskAdd(illegalMask.ite(A.getDdManager().getConstant(storm::utility::infinity()), A.getDdManager().template getAddZero())), 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 - typename SymbolicMinMaxLinearEquationSolver::ValueIterationResult SymbolicMinMaxLinearEquationSolver::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b, ValueType const& precision, bool relativeTerminationCriterion) const { + typename SymbolicMinMaxLinearEquationSolver::ValueIterationResult SymbolicMinMaxLinearEquationSolver::performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const { // Set up local variables. storm::dd::Add 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 localXAsColumn = localX.swapVariables(this->rowColumnMetaVariablePairs); storm::dd::Add tmp = this->A.multiplyMatrix(localXAsColumn, this->columnMetaVariables); @@ -178,7 +178,7 @@ namespace storm { template template - storm::dd::Add SymbolicMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver const& rationalSolver, SymbolicMinMaxLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalX, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const { + storm::dd::Add SymbolicMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver const& rationalSolver, SymbolicMinMaxLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const { // Storage for the rational sharpened vector. storm::dd::Add 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::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber(precision), this->getSettings().getRelativeTerminationCriterion()); + typename SymbolicMinMaxLinearEquationSolver::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber(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 template typename std::enable_if::value && storm::NumberTraits::IsExact, storm::dd::Add>::type SymbolicMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b) const { - - return solveEquationsRationalSearchHelper(dir, *this, *this, x, b, x, b); + return solveEquationsRationalSearchHelper(dir, *this, *this, b, this->getLowerBounds(), b); } template template typename std::enable_if::value && !storm::NumberTraits::IsExact, storm::dd::Add>::type SymbolicMinMaxLinearEquationSolver::solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b) const { - storm::dd::Add rationalX = x.template toValueType(); storm::dd::Add rationalB = b.template toValueType(); SymbolicMinMaxLinearEquationSolver rationalSolver(this->A.template toValueType(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique>()); - storm::dd::Add rationalResult = solveEquationsRationalSearchHelper(dir, rationalSolver, *this, rationalX, rationalB, x, b); + storm::dd::Add rationalResult = solveEquationsRationalSearchHelper(dir, rationalSolver, *this, rationalB, this->getLowerBounds(), b); return rationalResult.template toValueType(); } @@ -248,17 +246,18 @@ namespace storm { // First try to find a solution using the imprecise value type. storm::dd::Add rationalResult; + storm::dd::Add impreciseX; try { - storm::dd::Add impreciseX = x.template toValueType(); + impreciseX = this->getLowerBounds().template toValueType(); storm::dd::Add impreciseB = b.template toValueType(); SymbolicMinMaxLinearEquationSolver impreciseSolver(this->A.template toValueType(), this->allRows, this->illegalMask, this->rowMetaVariables, this->columnMetaVariables, this->choiceVariables, this->rowColumnMetaVariablePairs, std::make_unique>()); - rationalResult = solveEquationsRationalSearchHelper(dir, *this, impreciseSolver, x, b, impreciseX, impreciseB); + rationalResult = solveEquationsRationalSearchHelper(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(dir, *this, *this, x, b, x, b); + rationalResult = solveEquationsRationalSearchHelper(dir, *this, *this, b, impreciseX.template toValueType(), b); } return rationalResult.template toValueType(); } @@ -271,14 +270,16 @@ namespace storm { template storm::dd::Add SymbolicMinMaxLinearEquationSolver::solveEquationsValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b) const { // Set up the environment. - storm::dd::Add localX = x; + storm::dd::Add 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::SolutionMethod::ValueIteration) { + requirements.requireLowerBounds(); if (equationSystemType == EquationSystemType::ReachabilityRewards) { if (!direction || direction.get() == OptimizationDirection::Minimize) { requirements.requireValidInitialScheduler(); } } } else if (this->getSettings().getSolutionMethod() == SymbolicMinMaxLinearEquationSolverSettings::SolutionMethod::RationalSearch) { + requirements.requireLowerBounds(); if (equationSystemType == EquationSystemType::ReachabilityRewards) { if (!direction || direction.get() == OptimizationDirection::Minimize) { requirements.requireNoEndComponents(); diff --git a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h index 408faf21e..4e5fadbfd 100644 --- a/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h +++ b/src/storm/solver/SymbolicMinMaxLinearEquationSolver.h @@ -6,6 +6,7 @@ #include #include +#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 - class SymbolicMinMaxLinearEquationSolver { + class SymbolicMinMaxLinearEquationSolver : public SymbolicEquationSolver { public: SymbolicMinMaxLinearEquationSolver(SymbolicMinMaxLinearEquationSolverSettings const& settings = SymbolicMinMaxLinearEquationSolverSettings()); @@ -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 const& x, storm::dd::Add const& b) const; @@ -156,7 +157,7 @@ namespace storm { static storm::dd::Add sharpen(OptimizationDirection dir, uint64_t precision, SymbolicMinMaxLinearEquationSolver const& rationalSolver, storm::dd::Add const& x, storm::dd::Add const& rationalB, bool& isSolution); template - storm::dd::Add solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver const& rationalSolver, SymbolicMinMaxLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalX, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const; + storm::dd::Add solveEquationsRationalSearchHelper(OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver const& rationalSolver, SymbolicMinMaxLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const; template typename std::enable_if::value && storm::NumberTraits::IsExact, storm::dd::Add>::type solveEquationsRationalSearchHelper(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b) const; template @@ -177,15 +178,12 @@ namespace storm { storm::dd::Add values; }; - ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add const& b, ValueType const& precision, bool relativeTerminationCriterion) const; + ValueIterationResult performValueIteration(storm::solver::OptimizationDirection const& dir, storm::dd::Add const& x, storm::dd::Add 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 A; - // A BDD characterizing all rows of the equation system. - storm::dd::Bdd allRows; - // A BDD characterizing the illegal choices. storm::dd::Bdd illegalMask; diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp b/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp index 89f75bd41..5f3613100 100644 --- a/src/storm/solver/SymbolicNativeLinearEquationSolver.cpp +++ b/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(); // 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::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(settings.getPrecision()); relative = settings.getConvergenceCriterion() == storm::settings::modules::NativeEquationSolverSettings::ConvergenceCriterion::Relative; @@ -38,6 +64,11 @@ namespace storm { this->relative = value; } + template + void SymbolicNativeLinearEquationSolverSettings::setSolutionMethod(SolutionMethod const& method) { + this->method = method; + } + template ValueType SymbolicNativeLinearEquationSolverSettings::getPrecision() const { return precision; @@ -52,6 +83,16 @@ namespace storm { bool SymbolicNativeLinearEquationSolverSettings::getRelativeTerminationCriterion() const { return relative; } + + template + typename SymbolicNativeLinearEquationSolverSettings::SolutionMethod SymbolicNativeLinearEquationSolverSettings::getSolutionMethod() const { + return this->method; + } + + template + SymbolicNativeLinearEquationSolver::SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings const& settings) : SymbolicLinearEquationSolver(), settings(settings) { + // Intentionally left empty. + } template SymbolicNativeLinearEquationSolver::SymbolicNativeLinearEquationSolver(storm::dd::Add const& A, storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs, SymbolicNativeLinearEquationSolverSettings const& settings) : SymbolicNativeLinearEquationSolver(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings) { @@ -65,6 +106,18 @@ namespace storm { template storm::dd::Add SymbolicNativeLinearEquationSolver::solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const { + if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::Jacobi) { + return solveEquationsJacobi(x, b); + } else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::Power) { + return solveEquationsPower(x, b); + } else if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::RationalSearch) { + return solveEquationsRationalSearch(x, b); + } + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The selected solution technique is not supported."); + } + + template + storm::dd::Add SymbolicNativeLinearEquationSolver::solveEquationsJacobi(storm::dd::Add const& x, storm::dd::Add const& b) const { storm::dd::DdManager& manager = this->getDdManager(); // Start by computing the Jacobi decomposition of the matrix A. @@ -74,7 +127,7 @@ namespace storm { storm::dd::Add lu = diagonal.ite(manager.template getAddZero(), this->A); storm::dd::Add diagonalAdd = diagonal.template toAdd(); storm::dd::Add diag = diagonalAdd.multiplyMatrix(this->A, this->columnMetaVariables); - + storm::dd::Add scaledLu = lu / diag; storm::dd::Add scaledB = b / diag; @@ -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 + typename SymbolicNativeLinearEquationSolver::PowerIterationResult SymbolicNativeLinearEquationSolver::performPowerIteration(storm::dd::Add const& x, storm::dd::Add const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const { + + // Set up additional environment variables. + storm::dd::Add currentX = x; + uint_fast64_t iterations = 0; + SolverStatus status = SolverStatus::InProgress; + + while (status == SolverStatus::InProgress && iterations < maximalIterations) { + storm::dd::Add currentXAsColumn = currentX.swapVariables(this->rowColumnMetaVariablePairs); + storm::dd::Add 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::Add SymbolicNativeLinearEquationSolver::solveEquationsPower(storm::dd::Add const& x, storm::dd::Add 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 + bool SymbolicNativeLinearEquationSolver::isSolutionFixedPoint(storm::dd::Add const& x, storm::dd::Add const& b) const { + storm::dd::Add xAsColumn = x.swapVariables(this->rowColumnMetaVariablePairs); + storm::dd::Add tmp = this->A.multiplyMatrix(xAsColumn, this->columnMetaVariables); + tmp += b; + + return x == tmp; + } + + template + template + storm::dd::Add SymbolicNativeLinearEquationSolver::sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver const& rationalSolver, storm::dd::Add const& x, storm::dd::Add const& rationalB, bool& isSolution) { + + storm::dd::Add sharpenedX; + + for (uint64_t p = 1; p < precision; ++p) { + sharpenedX = x.sharpenKwekMehlhorn(precision); + isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB); + + if (isSolution) { + break; + } + } + + return sharpenedX; + } + + template + template + storm::dd::Add SymbolicNativeLinearEquationSolver::solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver const& rationalSolver, SymbolicNativeLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const { + + // Storage for the rational sharpened vector. + storm::dd::Add 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::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber(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(storm::utility::ceil(storm::utility::log10(storm::utility::one() / precision))); + + bool isSolution = false; + sharpenedX = sharpen(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 + template + typename std::enable_if::value && storm::NumberTraits::IsExact, storm::dd::Add>::type SymbolicNativeLinearEquationSolver::solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const { + return solveEquationsRationalSearchHelper(*this, *this, b, this->getLowerBounds(), b); + } + + template + template + typename std::enable_if::value && !storm::NumberTraits::IsExact, storm::dd::Add>::type SymbolicNativeLinearEquationSolver::solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const { + + storm::dd::Add rationalB = b.template toValueType(); + SymbolicNativeLinearEquationSolver rationalSolver(this->A.template toValueType(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs); + + storm::dd::Add rationalResult = solveEquationsRationalSearchHelper(rationalSolver, *this, rationalB, this->getLowerBounds(), b); + return rationalResult.template toValueType(); + } + + template + template + typename std::enable_if::value, storm::dd::Add>::type SymbolicNativeLinearEquationSolver::solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const { + + // First try to find a solution using the imprecise value type. + storm::dd::Add rationalResult; + storm::dd::Add impreciseX; + try { + impreciseX = this->getLowerBounds().template toValueType(); + storm::dd::Add impreciseB = b.template toValueType(); + SymbolicNativeLinearEquationSolver impreciseSolver(this->A.template toValueType(), this->allRows, this->rowMetaVariables, this->columnMetaVariables, this->rowColumnMetaVariablePairs); + + rationalResult = solveEquationsRationalSearchHelper(*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(*this, *this, b, impreciseX.template toValueType(), b); + } + return rationalResult.template toValueType(); + } + + template + storm::dd::Add SymbolicNativeLinearEquationSolver::solveEquationsRationalSearch(storm::dd::Add const& x, storm::dd::Add const& b) const { + return solveEquationsRationalSearchHelper(x, b); + } + + template + LinearEquationSolverProblemFormat SymbolicNativeLinearEquationSolver::getEquationProblemFormat() const { + if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::Jacobi) { + return LinearEquationSolverProblemFormat::EquationSystem; + } + return LinearEquationSolverProblemFormat::FixedPointSystem; + } + + template + LinearEquationSolverRequirements SymbolicNativeLinearEquationSolver::getRequirements() const { + LinearEquationSolverRequirements requirements; + + if (this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::Power || this->getSettings().getSolutionMethod() == SymbolicNativeLinearEquationSolverSettings::SolutionMethod::RationalSearch) { + requirements.requireLowerBounds(); + } + + return requirements; + } + template SymbolicNativeLinearEquationSolverSettings const& SymbolicNativeLinearEquationSolver::getSettings() const { return settings; } template - std::unique_ptr> SymbolicNativeLinearEquationSolverFactory::create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const { - return std::make_unique>(allRows, rowMetaVariables, columnMetaVariables, rowColumnMetaVariablePairs, settings); + std::unique_ptr> SymbolicNativeLinearEquationSolverFactory::create() const { + return std::make_unique>(settings); } template diff --git a/src/storm/solver/SymbolicNativeLinearEquationSolver.h b/src/storm/solver/SymbolicNativeLinearEquationSolver.h index 23b166a7c..b04b22fc5 100644 --- a/src/storm/solver/SymbolicNativeLinearEquationSolver.h +++ b/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 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 class SymbolicNativeLinearEquationSolver : public SymbolicLinearEquationSolver { public: + /*! + * Constructs a symbolic linear equation solver. + * + * @param settings The settings to use. + */ + SymbolicNativeLinearEquationSolver(SymbolicNativeLinearEquationSolverSettings const& settings = SymbolicNativeLinearEquationSolverSettings()); + /*! * 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 solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const; - + virtual storm::dd::Add solveEquations(storm::dd::Add const& x, storm::dd::Add const& b) const override; + + virtual LinearEquationSolverProblemFormat getEquationProblemFormat() const override; + virtual LinearEquationSolverRequirements getRequirements() const override; + SymbolicNativeLinearEquationSolverSettings const& getSettings() const; private: + storm::dd::Add solveEquationsJacobi(storm::dd::Add const& x, storm::dd::Add const& b) const; + storm::dd::Add solveEquationsPower(storm::dd::Add const& x, storm::dd::Add const& b) const; + storm::dd::Add solveEquationsRationalSearch(storm::dd::Add const& x, storm::dd::Add const& b) const; + + /*! + * Determines whether the given vector x satisfies x = Ax + b. + */ + bool isSolutionFixedPoint(storm::dd::Add const& x, storm::dd::Add const& b) const; + + template + static storm::dd::Add sharpen(uint64_t precision, SymbolicNativeLinearEquationSolver const& rationalSolver, storm::dd::Add const& x, storm::dd::Add const& rationalB, bool& isSolution); + + template + storm::dd::Add solveEquationsRationalSearchHelper(SymbolicNativeLinearEquationSolver const& rationalSolver, SymbolicNativeLinearEquationSolver const& impreciseSolver, storm::dd::Add const& rationalB, storm::dd::Add const& x, storm::dd::Add const& b) const; + template + typename std::enable_if::value && storm::NumberTraits::IsExact, storm::dd::Add>::type solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const; + template + typename std::enable_if::value && !storm::NumberTraits::IsExact, storm::dd::Add>::type solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const; + template + typename std::enable_if::value, storm::dd::Add>::type solveEquationsRationalSearchHelper(storm::dd::Add const& x, storm::dd::Add const& b) const; + + template + friend class SymbolicNativeLinearEquationSolver; + + struct PowerIterationResult { + PowerIterationResult(SolverStatus status, uint64_t iterations, storm::dd::Add const& values) : status(status), iterations(iterations), values(values) { + // Intentionally left empty. + } + + SolverStatus status; + uint64_t iterations; + storm::dd::Add values; + }; + + PowerIterationResult performPowerIteration(storm::dd::Add const& x, storm::dd::Add const& b, ValueType const& precision, bool relativeTerminationCriterion, uint64_t maximalIterations) const; + // The settings to use. SymbolicNativeLinearEquationSolverSettings settings; }; @@ -87,7 +145,7 @@ namespace storm { public: using SymbolicLinearEquationSolverFactory::create; - virtual std::unique_ptr> create(storm::dd::Bdd const& allRows, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::vector> const& rowColumnMetaVariablePairs) const override; + virtual std::unique_ptr> create() const override; SymbolicNativeLinearEquationSolverSettings& getSettings(); SymbolicNativeLinearEquationSolverSettings const& getSettings() const; diff --git a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp b/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp index 6385c1d93..273ecb8fc 100644 --- a/src/storm/storage/dd/bisimulation/QuotientExtractor.cpp +++ b/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()).isZero(), "Illegal entries in quotient matrix."); + if (std::is_same::value) { + STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one()).isZero(), "Illegal entries in quotient matrix."); + } else { + STORM_LOG_ASSERT(quotientTransitionMatrix.greater(storm::utility::one() + storm::utility::convertNumber(1e-6)).isZero(), "Illegal entries in quotient matrix."); + } STORM_LOG_ASSERT(quotientTransitionMatrix.sumAbstract(blockPrimeVariableSet).equalModuloPrecision(quotientTransitionMatrix.notZero().existsAbstract(blockPrimeVariableSet).template toAdd(), ValueType(1e-6)), "Illegal non-probabilistic matrix."); STORM_LOG_TRACE("Quotient transition matrix extracted in " << std::chrono::duration_cast(end - start).count() << "ms."); diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 7fb973653..77665f6e5 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -708,6 +708,11 @@ namespace storm { return number; } + template<> + RationalFunction convertNumber(std::string const& number){ + return RationalFunction(convertNumber(number)); + } + template<> RationalFunction convertNumber(storm::storage::sparse::state_type const& number) { return RationalFunction(convertNumber(number));