From 53db0b1f228cd0724052cc4fbb86fbacf72a629a Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Mon, 9 Mar 2020 13:10:21 +0100 Subject: [PATCH] Added AcyclicMinMaxLinearEquationSolver and AcyclicLinearEquationSolver which are optimized for many calls on an acyclic model. --- src/storm/settings/modules/CoreSettings.cpp | 4 +- .../modules/MinMaxEquationSolverSettings.cpp | 4 +- .../solver/AcyclicLinearEquationSolver.cpp | 117 ++++++++++++++++++ .../solver/AcyclicLinearEquationSolver.h | 40 ++++-- .../AcyclicMinMaxLinearEquationSolver.cpp | 100 +++++++++++++++ .../AcyclicMinMaxLinearEquationSolver.h | 50 ++++++++ src/storm/solver/LinearEquationSolver.cpp | 10 +- .../LinearEquationSolverRequirements.cpp | 30 ++++- .../solver/LinearEquationSolverRequirements.h | 7 ++ .../solver/MinMaxLinearEquationSolver.cpp | 5 + ...MinMaxLinearEquationSolverRequirements.cpp | 25 +++- .../MinMaxLinearEquationSolverRequirements.h | 9 +- src/storm/solver/SolverSelectionOptions.cpp | 4 + src/storm/solver/SolverSelectionOptions.h | 4 +- .../solver/helper/AcyclicSolverHelper.cpp | 117 ++++++++++++++++++ 15 files changed, 498 insertions(+), 28 deletions(-) create mode 100644 src/storm/solver/AcyclicLinearEquationSolver.cpp create mode 100644 src/storm/solver/AcyclicMinMaxLinearEquationSolver.cpp create mode 100644 src/storm/solver/AcyclicMinMaxLinearEquationSolver.h create mode 100644 src/storm/solver/helper/AcyclicSolverHelper.cpp diff --git a/src/storm/settings/modules/CoreSettings.cpp b/src/storm/settings/modules/CoreSettings.cpp index 19df21f3e..0f05c1aa6 100644 --- a/src/storm/settings/modules/CoreSettings.cpp +++ b/src/storm/settings/modules/CoreSettings.cpp @@ -40,7 +40,7 @@ namespace storm { this->addOption(storm::settings::OptionBuilder(moduleName, engineOptionName, false, "Sets which engine is used for model building and model checking.").setShortName(engineOptionShortName) .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the engine to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(engines)).setDefaultValueString("sparse").build()).build()); - std::vector linearEquationSolver = {"gmm++", "native", "eigen", "elimination", "topological"}; + std::vector linearEquationSolver = {"gmm++", "native", "eigen", "elimination", "topological", "acyclic"}; this->addOption(storm::settings::OptionBuilder(moduleName, eqSolverOptionName, false, "Sets which solver is preferred for solving systems of linear equations.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the solver to prefer.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(linearEquationSolver)).setDefaultValueString("topological").build()).build()); @@ -73,6 +73,8 @@ namespace storm { return storm::solver::EquationSolverType::Elimination; } else if (equationSolverName == "topological") { return storm::solver::EquationSolverType::Topological; + } else if (equationSolverName == "acyclic") { + return storm::solver::EquationSolverType::Acyclic; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown equation solver '" << equationSolverName << "'."); } diff --git a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp index 0014edb2d..875e9144d 100644 --- a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp +++ b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp @@ -21,7 +21,7 @@ namespace storm { const std::string MinMaxEquationSolverSettings::intervalIterationSymmetricUpdatesOptionName = "symmetricupdates"; MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) { - std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "lp", "linear-programming", "rs", "ratsearch", "ii", "interval-iteration", "svi", "sound-value-iteration", "ovi", "optimistic-value-iteration", "topological", "vi-to-pi"}; + std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "lp", "linear-programming", "rs", "ratsearch", "ii", "interval-iteration", "svi", "sound-value-iteration", "ovi", "optimistic-value-iteration", "topological", "vi-to-pi", "acyclic"}; this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, false, "Sets which min/max linear equation solving technique is preferred.").setIsAdvanced() .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(minMaxSolvingTechniques)).setDefaultValueString("topological").build()).build()); @@ -59,6 +59,8 @@ namespace storm { return storm::solver::MinMaxMethod::Topological; } else if (minMaxEquationSolvingTechnique == "vi-to-pi") { return storm::solver::MinMaxMethod::ViToPi; + } else if (minMaxEquationSolvingTechnique == "acyclic") { + return storm::solver::MinMaxMethod::Acyclic; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown min/max equation solving technique '" << minMaxEquationSolvingTechnique << "'."); diff --git a/src/storm/solver/AcyclicLinearEquationSolver.cpp b/src/storm/solver/AcyclicLinearEquationSolver.cpp new file mode 100644 index 000000000..d5a2940dd --- /dev/null +++ b/src/storm/solver/AcyclicLinearEquationSolver.cpp @@ -0,0 +1,117 @@ +#include "storm/solver/AcyclicLinearEquationSolver.h" + +#include "storm/solver/helper/AcyclicSolverHelper.cpp" + +#include "storm/utility/vector.h" + +namespace storm { + namespace solver { + + template + AcyclicLinearEquationSolver::AcyclicLinearEquationSolver() { + // Intentionally left empty. + } + + template + AcyclicLinearEquationSolver::AcyclicLinearEquationSolver(storm::storage::SparseMatrix const& A) { + this->setMatrix(A); + } + + template + AcyclicLinearEquationSolver::AcyclicLinearEquationSolver(storm::storage::SparseMatrix&& A) { + this->setMatrix(std::move(A)); + } + + template + void AcyclicLinearEquationSolver::setMatrix(storm::storage::SparseMatrix const& A) { + localA.reset(); + this->A = &A; + clearCache(); + } + + template + void AcyclicLinearEquationSolver::setMatrix(storm::storage::SparseMatrix&& A) { + localA = std::make_unique>(std::move(A)); + this->A = localA.get(); + clearCache(); + } + + template + uint64_t AcyclicLinearEquationSolver::getMatrixRowCount() const { + return this->A->getRowCount(); + } + + template + uint64_t AcyclicLinearEquationSolver::getMatrixColumnCount() const { + return this->A->getColumnCount(); + } + + template + bool AcyclicLinearEquationSolver::internalSolveEquations(Environment const& env, std::vector& x, std::vector const& b) const { + STORM_LOG_ASSERT(x.size() == this->A->getRowGroupCount(), "Provided x-vector has invalid size."); + STORM_LOG_ASSERT(b.size() == this->A->getRowCount(), "Provided b-vector has invalid size."); + + if (!multiplier) { + // We have not allocated cache memory, yet + rowOrdering = helper::computeTopologicalGroupOrdering(*this->A); + if (!rowOrdering) { + // It is not required to reorder the elements. + this->multiplier = storm::solver::MultiplierFactory().create(env, *this->A); + } else { + bFactors.clear(); + orderedMatrix = helper::createReorderedMatrix(*this->A, *rowOrdering, bFactors); + this->multiplier = storm::solver::MultiplierFactory().create(env, *orderedMatrix); + } + auxiliaryRowVector = std::vector(); + } + + std::vector const* bPtr = &b; + if (rowOrdering) { + STORM_LOG_ASSERT(rowOrdering->size() == b.size(), "b-vector has unexpected size."); + auxiliaryRowVector->resize(b.size()); + storm::utility::vector::selectVectorValues(*auxiliaryRowVector, *rowOrdering, b); + for (auto const& bFactor : bFactors) { + (*auxiliaryRowVector)[bFactor.first] *= bFactor.second; + } + bPtr = &auxiliaryRowVector.get(); + } + + this->multiplier->multiplyGaussSeidel(env, x, bPtr, true); + + if (!this->isCachingEnabled()) { + this->clearCache(); + } + return true; + } + + template + LinearEquationSolverProblemFormat AcyclicLinearEquationSolver::getEquationProblemFormat(storm::Environment const& env) const { + return LinearEquationSolverProblemFormat::FixedPointSystem; + } + + template + LinearEquationSolverRequirements AcyclicLinearEquationSolver::getRequirements(Environment const& env) const { + // Return the requirements of the underlying solver + LinearEquationSolverRequirements requirements; + requirements.requireAcyclic(); + return requirements; + } + + template + void AcyclicLinearEquationSolver::clearCache() const { + multiplier.reset(); + orderedMatrix = boost::none; + rowOrdering = boost::none; + auxiliaryRowVector = boost::none; + bFactors.clear(); + } + + // Explicitly instantiate the min max linear equation solver. + template class AcyclicLinearEquationSolver; + +#ifdef STORM_HAVE_CARL + template class AcyclicLinearEquationSolver; + template class AcyclicLinearEquationSolver; +#endif + } +} diff --git a/src/storm/solver/AcyclicLinearEquationSolver.h b/src/storm/solver/AcyclicLinearEquationSolver.h index e6e5fe3a3..e16931d09 100644 --- a/src/storm/solver/AcyclicLinearEquationSolver.h +++ b/src/storm/solver/AcyclicLinearEquationSolver.h @@ -1,6 +1,8 @@ #pragma once -#include "storm/solver/StandardMinMaxLinearEquationSolver.h" +#include +#include "storm/solver/Multiplier.h" +#include "storm/solver/LinearEquationSolver.h" namespace storm { @@ -13,30 +15,44 @@ namespace storm { * It is optimized for solving many instances of the equation system with the same underlying matrix. */ template - class AcyclicMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver { + class AcyclicLinearEquationSolver : public LinearEquationSolver { public: - AcyclicMinMaxLinearEquationSolver(); - AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A); - AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A); - - virtual ~AcyclicMinMaxLinearEquationSolver() { + AcyclicLinearEquationSolver(); + AcyclicLinearEquationSolver(storm::storage::SparseMatrix const& A); + AcyclicLinearEquationSolver(storm::storage::SparseMatrix&& A); + + virtual void setMatrix(storm::storage::SparseMatrix const& A) override; + virtual void setMatrix(storm::storage::SparseMatrix&& A) override; + + virtual ~AcyclicLinearEquationSolver() { } - + virtual void clearCache() const override; - virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const& env, boost::optional const& direction = boost::none, bool const& hasInitialScheduler = false) const override ; + virtual LinearEquationSolverProblemFormat getEquationProblemFormat(storm::Environment const& env) const override; + virtual LinearEquationSolverRequirements getRequirements(Environment const& env) const override; protected: - - virtual bool internalSolveEquations(storm::Environment const& env, OptimizationDirection d, std::vector& x, std::vector const& b) const override; + virtual bool internalSolveEquations(storm::Environment const& env, std::vector& x, std::vector const& b) const override; private: + + virtual uint64_t getMatrixRowCount() const override; + virtual uint64_t getMatrixColumnCount() const override; + + // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted + // when the solver is destructed. + std::unique_ptr> localA; + // A pointer to the original sparse matrix given to this solver. If the solver takes posession of the matrix + // the pointer refers to orderedMatrix. + storm::storage::SparseMatrix const* A; + // cached multiplier either with original matrix or ordered matrix mutable std::unique_ptr> multiplier; // cached matrix for the multiplier (only if different from original matrix) mutable boost::optional> orderedMatrix; // cached row group ordering (only if not identity) - mutable boost::optional> rowGroupOrdering; // A.rowGroupCount() entries + mutable boost::optional> rowOrdering; // A.rowGroupCount() entries // can be used if the entries in 'b' need to be reordered mutable boost::optional> auxiliaryRowVector; // A.rowCount() entries // contains factors applied to scale the entries of the 'b' vector diff --git a/src/storm/solver/AcyclicMinMaxLinearEquationSolver.cpp b/src/storm/solver/AcyclicMinMaxLinearEquationSolver.cpp new file mode 100644 index 000000000..3bee1aa3a --- /dev/null +++ b/src/storm/solver/AcyclicMinMaxLinearEquationSolver.cpp @@ -0,0 +1,100 @@ +#include "storm/solver/AcyclicMinMaxLinearEquationSolver.h" + +#include "storm/solver/helper/AcyclicSolverHelper.cpp" + +#include "storm/utility/vector.h" + +namespace storm { + namespace solver { + + template + AcyclicMinMaxLinearEquationSolver::AcyclicMinMaxLinearEquationSolver() { + // Intentionally left empty. + } + + template + AcyclicMinMaxLinearEquationSolver::AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A) : StandardMinMaxLinearEquationSolver(A) { + // Intentionally left empty. + } + + template + AcyclicMinMaxLinearEquationSolver::AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A) : StandardMinMaxLinearEquationSolver(std::move(A)) { + // Intentionally left empty. + } + + + template + bool AcyclicMinMaxLinearEquationSolver::internalSolveEquations(Environment const& env, OptimizationDirection dir, std::vector& x, std::vector const& b) const { + STORM_LOG_ASSERT(x.size() == this->A->getRowGroupCount(), "Provided x-vector has invalid size."); + STORM_LOG_ASSERT(b.size() == this->A->getRowCount(), "Provided b-vector has invalid size."); + // Allocate memory for the scheduler (if required) + if (this->isTrackSchedulerSet()) { + if (this->schedulerChoices) { + this->schedulerChoices->resize(this->A->getRowGroupCount()); + } else { + this->schedulerChoices = std::vector(this->A->getRowGroupCount()); + } + } + + if (!multiplier) { + // We have not allocated cache memory, yet + rowGroupOrdering = helper::computeTopologicalGroupOrdering(*this->A); + if (!rowGroupOrdering) { + // It is not required to reorder the elements. + this->multiplier = storm::solver::MultiplierFactory().create(env, *this->A); + } else { + bFactors.clear(); + orderedMatrix = helper::createReorderedMatrix(*this->A, *rowGroupOrdering, bFactors); + this->multiplier = storm::solver::MultiplierFactory().create(env, *orderedMatrix); + } + auxiliaryRowVector = std::vector(); + } + + std::vector const* bPtr = &b; + if (rowGroupOrdering) { + STORM_LOG_ASSERT(rowGroupOrdering->size() == b.size(), "b-vector has unexpected size."); + auxiliaryRowVector->resize(b.size()); + storm::utility::vector::selectVectorValues(*auxiliaryRowVector, *rowGroupOrdering, b); + for (auto const& bFactor : bFactors) { + (*auxiliaryRowVector)[bFactor.first] *= bFactor.second; + } + bPtr = &auxiliaryRowVector.get(); + } + + if (this->isTrackSchedulerSet()) { + this->multiplier->multiplyAndReduceGaussSeidel(env, dir, x, bPtr, &this->schedulerChoices.get(), true); + } else { + this->multiplier->multiplyAndReduceGaussSeidel(env, dir, x, bPtr, nullptr, true); + } + + if (!this->isCachingEnabled()) { + this->clearCache(); + } + return true; + } + + template + MinMaxLinearEquationSolverRequirements AcyclicMinMaxLinearEquationSolver::getRequirements(Environment const& env, boost::optional const& direction, bool const& hasInitialScheduler) const { + // Return the requirements of the underlying solver + MinMaxLinearEquationSolverRequirements requirements; + requirements.requireAcyclic(); + return requirements; + } + + template + void AcyclicMinMaxLinearEquationSolver::clearCache() const { + multiplier.reset(); + orderedMatrix = boost::none; + rowGroupOrdering = boost::none; + auxiliaryRowVector = boost::none; + bFactors.clear(); + } + + // Explicitly instantiate the min max linear equation solver. + template class AcyclicMinMaxLinearEquationSolver; + +#ifdef STORM_HAVE_CARL + template class AcyclicMinMaxLinearEquationSolver; +#endif + } +} diff --git a/src/storm/solver/AcyclicMinMaxLinearEquationSolver.h b/src/storm/solver/AcyclicMinMaxLinearEquationSolver.h new file mode 100644 index 000000000..f199c98ab --- /dev/null +++ b/src/storm/solver/AcyclicMinMaxLinearEquationSolver.h @@ -0,0 +1,50 @@ +#pragma once + +#include + +#include "storm/solver/Multiplier.h" +#include "storm/solver/StandardMinMaxLinearEquationSolver.h" + +namespace storm { + + class Environment; + + namespace solver { + + /*! + * This solver can be used on equation systems that are known to be acyclic. + * It is optimized for solving many instances of the equation system with the same underlying matrix. + */ + template + class AcyclicMinMaxLinearEquationSolver : public StandardMinMaxLinearEquationSolver { + public: + AcyclicMinMaxLinearEquationSolver(); + AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix const& A); + AcyclicMinMaxLinearEquationSolver(storm::storage::SparseMatrix&& A); + + virtual ~AcyclicMinMaxLinearEquationSolver() { + } + + virtual void clearCache() const override; + + virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const& env, boost::optional const& direction = boost::none, bool const& hasInitialScheduler = false) const override ; + + protected: + + virtual bool internalSolveEquations(storm::Environment const& env, OptimizationDirection d, std::vector& x, std::vector const& b) const override; + + private: + // cached multiplier either with original matrix or ordered matrix + mutable std::unique_ptr> multiplier; + // cached matrix for the multiplier (only if different from original matrix) + mutable boost::optional> orderedMatrix; + // cached row group ordering (only if not identity) + mutable boost::optional> rowGroupOrdering; // A.rowGroupCount() entries + // can be used if the entries in 'b' need to be reordered + mutable boost::optional> auxiliaryRowVector; // A.rowCount() entries + // contains factors applied to scale the entries of the 'b' vector + mutable std::vector> bFactors; + + }; + } +} diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index 5d2884d86..39aee7d89 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -8,6 +8,7 @@ #include "storm/solver/EigenLinearEquationSolver.h" #include "storm/solver/EliminationLinearEquationSolver.h" #include "storm/solver/TopologicalLinearEquationSolver.h" +#include "storm/solver/AcyclicLinearEquationSolver.h" #include "storm/utility/vector.h" @@ -88,7 +89,7 @@ namespace storm { EquationSolverType type = env.solver().getLinearEquationSolverType(); // Adjust the solver type if it is not supported by this value type - if (type != EquationSolverType::Eigen && type != EquationSolverType::Topological && (env.solver().isLinearEquationSolverTypeSetFromDefaultValue() || type == EquationSolverType::Gmmxx)) { + if (type != EquationSolverType::Eigen && type != EquationSolverType::Topological && type != EquationSolverType::Acyclic && (env.solver().isLinearEquationSolverTypeSetFromDefaultValue() || type == EquationSolverType::Gmmxx)) { STORM_LOG_INFO("Selecting '" + toString(EquationSolverType::Eigen) + "' as the linear equation solver since the previously selected one (" << toString(type) << ") does not support exact computations."); type = EquationSolverType::Eigen; } @@ -98,6 +99,7 @@ namespace storm { case EquationSolverType::Eigen: return std::make_unique>(); case EquationSolverType::Elimination: return std::make_unique>(); case EquationSolverType::Topological: return std::make_unique>(); + case EquationSolverType::Acyclic: return std::make_unique>(); default: STORM_LOG_THROW(false, storm::exceptions::InvalidEnvironmentException, "Unknown solver type."); return nullptr; @@ -118,6 +120,7 @@ namespace storm { case EquationSolverType::Eigen: return std::make_unique>(); case EquationSolverType::Elimination: return std::make_unique>(); case EquationSolverType::Topological: return std::make_unique>(); + case EquationSolverType::Acyclic: return std::make_unique>(); default: STORM_LOG_THROW(false, storm::exceptions::InvalidEnvironmentException, "Unknown solver type."); return nullptr; @@ -129,14 +132,14 @@ namespace storm { EquationSolverType type = env.solver().getLinearEquationSolverType(); // Adjust the solver type if none was specified and we want sound/exact computations - if (env.solver().isForceExact() && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological) { + if (env.solver().isForceExact() && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological && type != EquationSolverType::Acyclic) { if (env.solver().isLinearEquationSolverTypeSetFromDefaultValue()) { type = EquationSolverType::Eigen; STORM_LOG_INFO("Selecting '" + toString(type) + "' as the linear equation solver to guarantee exact results. If you want to override this, please explicitly specify a different solver."); } else { STORM_LOG_WARN("The selected solver does not yield exact results."); } - } else if (env.solver().isForceSoundness() && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological) { + } else if (env.solver().isForceSoundness() && type != EquationSolverType::Native && type != EquationSolverType::Eigen && type != EquationSolverType::Elimination && type != EquationSolverType::Topological && type != EquationSolverType::Acyclic) { if (env.solver().isLinearEquationSolverTypeSetFromDefaultValue()) { type = EquationSolverType::Native; STORM_LOG_INFO("Selecting '" + toString(type) + "' as the linear equation solver to guarantee sound results. If you want to override this, please explicitly specify a different solver."); @@ -151,6 +154,7 @@ namespace storm { case EquationSolverType::Eigen: return std::make_unique>(); case EquationSolverType::Elimination: return std::make_unique>(); case EquationSolverType::Topological: return std::make_unique>(); + case EquationSolverType::Acyclic: return std::make_unique>(); default: STORM_LOG_THROW(false, storm::exceptions::InvalidEnvironmentException, "Unknown solver type."); return nullptr; diff --git a/src/storm/solver/LinearEquationSolverRequirements.cpp b/src/storm/solver/LinearEquationSolverRequirements.cpp index e4075b1e1..64d333ab1 100644 --- a/src/storm/solver/LinearEquationSolverRequirements.cpp +++ b/src/storm/solver/LinearEquationSolverRequirements.cpp @@ -10,6 +10,11 @@ namespace storm { // Intentionally left empty. } + LinearEquationSolverRequirements& LinearEquationSolverRequirements::requireAcyclic(bool critical) { + acyclicRequirement.enable(critical); + return *this; + } + LinearEquationSolverRequirements& LinearEquationSolverRequirements::requireLowerBounds(bool critical) { lowerBoundsRequirement.enable(critical); return *this; @@ -26,6 +31,10 @@ namespace storm { return *this; } + SolverRequirement const& LinearEquationSolverRequirements::acyclic() const { + return acyclicRequirement; + } + SolverRequirement const& LinearEquationSolverRequirements::lowerBounds() const { return lowerBoundsRequirement; } @@ -36,12 +45,17 @@ namespace storm { SolverRequirement const& LinearEquationSolverRequirements::get(Element const& element) const { switch (element) { - case Element::LowerBounds: return lowerBounds(); break; - case Element::UpperBounds: return upperBounds(); break; + case Element::Acyclic: return acyclic(); + case Element::LowerBounds: return lowerBounds(); + case Element::UpperBounds: return upperBounds(); } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Unknown ElementType"); } + void LinearEquationSolverRequirements::clearAcyclic() { + acyclicRequirement.clear(); + } + void LinearEquationSolverRequirements::clearLowerBounds() { lowerBoundsRequirement.clear(); } @@ -51,17 +65,23 @@ namespace storm { } bool LinearEquationSolverRequirements::hasEnabledRequirement() const { - return lowerBoundsRequirement || upperBoundsRequirement; + return acyclicRequirement || lowerBoundsRequirement || upperBoundsRequirement; } bool LinearEquationSolverRequirements::hasEnabledCriticalRequirement() const { - return lowerBoundsRequirement.isCritical() || upperBoundsRequirement.isCritical(); + return acyclicRequirement.isCritical() || lowerBoundsRequirement.isCritical() || upperBoundsRequirement.isCritical(); } - std::string LinearEquationSolverRequirements::getEnabledRequirementsAsString() const { std::string res = "["; bool first = true; + if (acyclic()) { + if (!first) { res += ", "; } else {first = false;} + res += "acyclic"; + if (acyclic().isCritical()) { + res += "(mandatory)"; + } + } if (lowerBounds()) { if (!first) { res += ", "; } else {first = false;} res += "lowerBounds"; diff --git a/src/storm/solver/LinearEquationSolverRequirements.h b/src/storm/solver/LinearEquationSolverRequirements.h index a7e23d248..2903e8398 100644 --- a/src/storm/solver/LinearEquationSolverRequirements.h +++ b/src/storm/solver/LinearEquationSolverRequirements.h @@ -11,6 +11,9 @@ namespace storm { public: // The different requirements a solver can have. enum class Element { + // Requirements that are related to the graph structure of the model. + Acyclic, + // Requirements that are related to bounds for the actual solution. LowerBounds, UpperBounds @@ -18,14 +21,17 @@ namespace storm { LinearEquationSolverRequirements(); + LinearEquationSolverRequirements& requireAcyclic(bool critical = true); LinearEquationSolverRequirements& requireLowerBounds(bool critical = true); LinearEquationSolverRequirements& requireUpperBounds(bool critical = true); LinearEquationSolverRequirements& requireBounds(bool critical = true); + SolverRequirement const& acyclic() const; SolverRequirement const& lowerBounds() const; SolverRequirement const& upperBounds() const; SolverRequirement const& get(Element const& element) const; + void clearAcyclic(); void clearLowerBounds(); void clearUpperBounds(); @@ -39,6 +45,7 @@ namespace storm { std::string getEnabledRequirementsAsString() const; private: + SolverRequirement acyclicRequirement; SolverRequirement lowerBoundsRequirement; SolverRequirement upperBoundsRequirement; }; diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index 0378e0cce..c36e291c2 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -7,6 +7,7 @@ #include "storm/solver/TopologicalMinMaxLinearEquationSolver.h" #include "storm/solver/TopologicalCudaMinMaxLinearEquationSolver.h" #include "storm/solver/LpMinMaxLinearEquationSolver.h" +#include "storm/solver/AcyclicMinMaxLinearEquationSolver.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" @@ -211,6 +212,8 @@ namespace storm { result = std::make_unique>(); } else if (method == MinMaxMethod::LinearProgramming) { result = std::make_unique>(std::make_unique>()); + } else if (method == MinMaxMethod::Acyclic) { + result = std::make_unique>(); } else { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } @@ -226,6 +229,8 @@ namespace storm { result = std::make_unique>(std::make_unique>()); } else if (method == MinMaxMethod::LinearProgramming) { result = std::make_unique>(std::make_unique>()); + } else if (method == MinMaxMethod::Acyclic) { + result = std::make_unique>(); } else if (method == MinMaxMethod::Topological) { result = std::make_unique>(); } else { diff --git a/src/storm/solver/MinMaxLinearEquationSolverRequirements.cpp b/src/storm/solver/MinMaxLinearEquationSolverRequirements.cpp index 3d7cdc607..1fdfc172e 100644 --- a/src/storm/solver/MinMaxLinearEquationSolverRequirements.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolverRequirements.cpp @@ -10,6 +10,11 @@ namespace storm { // Intentionally left empty. } + MinMaxLinearEquationSolverRequirements& MinMaxLinearEquationSolverRequirements::requireAcyclic(bool critical) { + acyclicRequirement.enable(critical); + return *this; + } + MinMaxLinearEquationSolverRequirements& MinMaxLinearEquationSolverRequirements::requireUniqueSolution(bool critical) { uniqueSolutionRequirement.enable(critical); return *this; @@ -36,6 +41,10 @@ namespace storm { return *this; } + SolverRequirement const& MinMaxLinearEquationSolverRequirements::acyclic() const { + return acyclicRequirement; + } + SolverRequirement const& MinMaxLinearEquationSolverRequirements::uniqueSolution() const { return uniqueSolutionRequirement; } @@ -54,6 +63,7 @@ namespace storm { SolverRequirement const& MinMaxLinearEquationSolverRequirements::get(Element const& element) const { switch (element) { + case Element::Acyclic: return acyclic(); break; case Element::UniqueSolution: return uniqueSolution(); break; case Element::ValidInitialScheduler: return validInitialScheduler(); break; case Element::LowerBounds: return lowerBounds(); break; @@ -62,6 +72,10 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentException, "Unknown ElementType"); } + void MinMaxLinearEquationSolverRequirements::clearAcyclic() { + acyclicRequirement.clear(); + } + void MinMaxLinearEquationSolverRequirements::clearUniqueSolution() { uniqueSolutionRequirement.clear(); } @@ -84,16 +98,23 @@ namespace storm { } bool MinMaxLinearEquationSolverRequirements::hasEnabledRequirement() const { - return uniqueSolutionRequirement || validInitialSchedulerRequirement || lowerBoundsRequirement || upperBoundsRequirement; + return acyclicRequirement || uniqueSolutionRequirement || validInitialSchedulerRequirement || lowerBoundsRequirement || upperBoundsRequirement; } bool MinMaxLinearEquationSolverRequirements::hasEnabledCriticalRequirement() const { - return uniqueSolutionRequirement.isCritical() || validInitialSchedulerRequirement.isCritical() || lowerBoundsRequirement.isCritical() || upperBoundsRequirement.isCritical(); + return acyclicRequirement.isCritical() || uniqueSolutionRequirement.isCritical() || validInitialSchedulerRequirement.isCritical() || lowerBoundsRequirement.isCritical() || upperBoundsRequirement.isCritical(); } std::string MinMaxLinearEquationSolverRequirements::getEnabledRequirementsAsString() const { std::string res = "["; bool first = true; + if (acyclic()) { + if (!first) { res += ", "; } else {first = false;} + res += "Acyclic"; + if (acyclic().isCritical()) { + res += "(mandatory)"; + } + } if (uniqueSolution()) { if (!first) { res += ", "; } else {first = false;} res += "UniqueSolution"; diff --git a/src/storm/solver/MinMaxLinearEquationSolverRequirements.h b/src/storm/solver/MinMaxLinearEquationSolverRequirements.h index 38f2f1849..2232e8a6e 100644 --- a/src/storm/solver/MinMaxLinearEquationSolverRequirements.h +++ b/src/storm/solver/MinMaxLinearEquationSolverRequirements.h @@ -13,8 +13,9 @@ namespace storm { // The different requirements a solver can have. enum class Element { // Requirements that are related to the graph structure of the system. Note that the requirements in this - // category are to be interpreted incrementally in the following sense: whenever the system has no end - // components then automatically both requirements are fulfilled. + // category are to be interpreted incrementally in the following sense: whenever the system has a unique + // solution then a valid initial scheduler is no longer required. + Acyclic, UniqueSolution, ValidInitialScheduler, @@ -27,18 +28,21 @@ namespace storm { MinMaxLinearEquationSolverRequirements(LinearEquationSolverRequirements const& linearEquationSolverRequirements = LinearEquationSolverRequirements()); + MinMaxLinearEquationSolverRequirements& requireAcyclic(bool critical = true); MinMaxLinearEquationSolverRequirements& requireUniqueSolution(bool critical = true); MinMaxLinearEquationSolverRequirements& requireValidInitialScheduler(bool critical = true); MinMaxLinearEquationSolverRequirements& requireLowerBounds(bool critical = true); MinMaxLinearEquationSolverRequirements& requireUpperBounds(bool critical = true); MinMaxLinearEquationSolverRequirements& requireBounds(bool critical = true); + SolverRequirement const& acyclic() const; SolverRequirement const& uniqueSolution() const; SolverRequirement const& validInitialScheduler() const; SolverRequirement const& lowerBounds() const; SolverRequirement const& upperBounds() const; SolverRequirement const& get(Element const& element) const; + void clearAcyclic(); void clearUniqueSolution(); void clearValidInitialScheduler(); void clearLowerBounds(); @@ -54,6 +58,7 @@ namespace storm { std::string getEnabledRequirementsAsString() const; private: + SolverRequirement acyclicRequirement; SolverRequirement uniqueSolutionRequirement; SolverRequirement validInitialSchedulerRequirement; SolverRequirement lowerBoundsRequirement; diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp index 12e1cac02..74b181fd3 100644 --- a/src/storm/solver/SolverSelectionOptions.cpp +++ b/src/storm/solver/SolverSelectionOptions.cpp @@ -24,6 +24,8 @@ namespace storm { return "topologicalcuda"; case MinMaxMethod::ViToPi: return "vi-to-pi"; + case MinMaxMethod::Acyclic: + return "vi-to-pi"; } return "invalid"; } @@ -96,6 +98,8 @@ namespace storm { return "Elimination"; case EquationSolverType::Topological: return "Topological"; + case EquationSolverType::Acyclic: + return "Acyclic"; } return "invalid"; } diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h index 7c2c90693..3b0c2cd36 100644 --- a/src/storm/solver/SolverSelectionOptions.h +++ b/src/storm/solver/SolverSelectionOptions.h @@ -6,14 +6,14 @@ namespace storm { namespace solver { - ExtendEnumsWithSelectionField(MinMaxMethod, ValueIteration, PolicyIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, OptimisticValueIteration, TopologicalCuda, ViToPi) + ExtendEnumsWithSelectionField(MinMaxMethod, ValueIteration, PolicyIteration, LinearProgramming, Topological, RationalSearch, IntervalIteration, SoundValueIteration, OptimisticValueIteration, TopologicalCuda, ViToPi, Acyclic) ExtendEnumsWithSelectionField(MultiplierType, Native, Gmmxx) ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration) ExtendEnumsWithSelectionField(LraMethod, LinearProgramming, ValueIteration, GainBiasEquations, LraDistributionEquations) ExtendEnumsWithSelectionField(MaBoundedReachabilityMethod, Imca, UnifPlus) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk, Z3) - ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination, Topological) + ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination, Topological, Acyclic) ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat) ExtendEnumsWithSelectionField(NativeLinearEquationSolverMethod, Jacobi, GaussSeidel, SOR, WalkerChae, Power, SoundValueIteration, OptimisticValueIteration, IntervalIteration, RationalSearch) diff --git a/src/storm/solver/helper/AcyclicSolverHelper.cpp b/src/storm/solver/helper/AcyclicSolverHelper.cpp new file mode 100644 index 000000000..9bbdf6fbb --- /dev/null +++ b/src/storm/solver/helper/AcyclicSolverHelper.cpp @@ -0,0 +1,117 @@ +#include "storm/storage/SparseMatrix.h" +#include "storm/utility/constants.h" +#include "storm/utility/vector.h" +#include "storm/exceptions/InvalidStateException.h" +#include "storm/exceptions/InvalidEnvironmentException.h" +#include "storm/exceptions/UnexpectedException.h" +#include "storm/exceptions/UnmetRequirementException.h" + +namespace storm { + namespace solver { + + namespace helper { + + + /*! + * Returns a reordering of the matrix row(groups) and columns such that we can solve the (minmax or linear) equation system in one go. + * More precisely, let x be the result and i an arbitrary rowgroup index. Solving for rowgroup x[i] only requires knowledge of the result at rowgroups x[i+1], x[i+2], ... + */ + template + boost::optional> computeTopologicalGroupOrdering(storm::storage::SparseMatrix const& matrix) { + uint64_t numGroups = matrix.getRowGroupCount(); + bool orderedMatrixRequired = false; + std::vector result; + result.reserve(numGroups); + storm::storage::BitVector processed(numGroups, false); + storm::storage::BitVector visited(numGroups, false); + std::vector stack; + + // It's more likely that states without a successor are located at the end (due to the way we build the model). + // We therefore process the matrix in reverse order. + uint64_t startState = numGroups; + while (startState > 0) { + --startState; + // skip already processed states + if (processed.get(startState)) continue; + + // Now do a dfs from start state. + stack.push_back(startState); + while (!stack.empty()) { + uint64_t current = stack.back(); + if (visited.get(current)) { + // We are backtracking, so add this state now + result.push_back(current); + processed.set(current); + stack.pop_back(); + } else { + visited.set(current); + for (auto const& entry : matrix.getRowGroup(current)) { + if (!processed.get(entry.getColumn()) && !storm::utility::isZero(entry.getValue())) { + orderedMatrixRequired = true; + STORM_LOG_THROW(!visited.get(entry.getColumn()), storm::exceptions::UnmetRequirementException, "The model is not acyclic."); + stack.push_back(entry.getColumn()); + } + } + // If there are no successors to process, we will add the current state to the result in the next iteration. + } + } + } + // we will do backwards iterations, so the order has to be reversed + if (orderedMatrixRequired) { + std::reverse(result.begin(), result.end()); + return result; + } else { + return boost::none; + } + } + + /// reorders the row group such that the i'th row of the new matrix corresponds to the order[i]'th row of the source matrix. + /// Also eliminates selfloops p>0 and inserts 1/p into the bFactors + template + storm::storage::SparseMatrix createReorderedMatrix(storm::storage::SparseMatrix const& matrix, std::vector const& newToOrigIndexMap, std::vector>& bFactors) { + std::vector origToNewMap(newToOrigIndexMap.size(), std::numeric_limits::max()); + for (uint64_t i = 0; i < newToOrigIndexMap.size(); ++i) { + origToNewMap[newToOrigIndexMap[i]] = i; + } + + bool hasRowGrouping = !matrix.hasTrivialRowGrouping(); + storm::storage::SparseMatrixBuilder builder(matrix.getRowCount(), matrix.getColumnCount(), matrix.getNonzeroEntryCount(), false, hasRowGrouping, hasRowGrouping ? matrix.getRowGroupCount() : static_cast(0)); + uint64_t newRow = 0; + for (uint64_t newRowGroup = 0; newRowGroup < newToOrigIndexMap.size(); ++newRowGroup) { + auto const& origRowGroup = newToOrigIndexMap[newRowGroup]; + if (hasRowGrouping) { + builder.newRowGroup(newRowGroup); + } + for (uint64_t origRow = matrix.getRowGroupIndices()[origRowGroup]; origRow < matrix.getRowGroupIndices()[origRowGroup + 1]; ++origRow) { + for (auto const& entry : matrix.getRow(origRow)) { + if (storm::utility::isZero(entry.getValue())) { + continue; + } + if (entry.getColumn() == origRowGroup) { + if (storm::utility::isOne(entry.getValue())) { + // a one selfloop can only mean that there is never a non-zero value at the b vector for the current row. + // Let's "assert" this by considering infinity. (This is necessary to avoid division by zero) + bFactors.emplace_back(newRow, storm::utility::infinity()); + } + ValueType factor = storm::utility::one() / (storm::utility::one() - entry.getValue()); + bFactors.emplace_back(newRow, factor); + } + builder.addNextValue(newRow, origToNewMap[entry.getColumn()], entry.getValue()); + } + ++newRow; + } + } + auto result = builder.build(matrix.getRowCount(), matrix.getColumnCount(), matrix.getRowGroupCount()); + // apply the bFactors to the relevant rows + for (auto const& bFactor : bFactors) { + STORM_LOG_ASSERT(!storm::utility::isInfinity(bFactor.second) || storm::utility::isZero(result.getRowSum(bFactor.first)), "The input matrix does not seem to be probabilistic."); + for (auto& entry : result.getRow(bFactor.first)) { + entry.setValue(entry.getValue() * bFactor.second); + } + } + STORM_LOG_DEBUG("Reordered " << matrix.getDimensionsAsString() << " with " << bFactors.size() << " selfloop entries for acyclic solving."); + return result; + } + } + } +}