diff --git a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp index c13bb492d..c140e7520 100644 --- a/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp +++ b/src/storm/settings/modules/MinMaxEquationSolverSettings.cpp @@ -19,7 +19,7 @@ namespace storm { const std::string MinMaxEquationSolverSettings::absoluteOptionName = "absolute"; MinMaxEquationSolverSettings::MinMaxEquationSolverSettings() : ModuleSettings(moduleName) { - std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration"}; + std::vector minMaxSolvingTechniques = {"vi", "value-iteration", "pi", "policy-iteration", "acyclic"}; this->addOption(storm::settings::OptionBuilder(moduleName, solvingMethodOptionName, false, "Sets which min/max linear equation solving technique is preferred.") .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a min/max linear equation solving technique.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(minMaxSolvingTechniques)).setDefaultValueString("vi").build()).build()); @@ -36,6 +36,8 @@ namespace storm { return storm::solver::MinMaxMethod::ValueIteration; } else if (minMaxEquationSolvingTechnique == "policy-iteration" || minMaxEquationSolvingTechnique == "pi") { return storm::solver::MinMaxMethod::PolicyIteration; + } 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/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index 5f1e0c176..4bdb08f3e 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -170,7 +170,7 @@ namespace storm { std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { std::unique_ptr> result; auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); - if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration) { + if (method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic) { result = std::make_unique>(std::forward(matrix), std::make_unique>()); } else if (method == MinMaxMethod::Topological) { result = std::make_unique>(std::forward(matrix)); @@ -187,7 +187,7 @@ namespace storm { std::unique_ptr> GeneralMinMaxLinearEquationSolverFactory::selectSolver(MatrixType&& matrix) const { std::unique_ptr> result; auto method = storm::settings::getModule().getMinMaxEquationSolvingMethod(); - STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration, storm::exceptions::InvalidSettingsException, "For this data type only value iteration and policy iteration are available."); + STORM_LOG_THROW(method == MinMaxMethod::ValueIteration || method == MinMaxMethod::PolicyIteration || method == MinMaxMethod::Acyclic, storm::exceptions::InvalidSettingsException, "For this data type only value iteration, policy iteration, and acyclic value iteration are available."); return std::make_unique>(std::forward(matrix), std::make_unique>()); } #endif diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp index fe42f299a..2a39108ea 100644 --- a/src/storm/solver/SolverSelectionOptions.cpp +++ b/src/storm/solver/SolverSelectionOptions.cpp @@ -10,7 +10,8 @@ namespace storm { return "value"; case MinMaxMethod::Topological: return "topological"; - + case MinMaxMethod::Acyclic: + return "acyclic"; } return "invalid"; } diff --git a/src/storm/solver/SolverSelectionOptions.h b/src/storm/solver/SolverSelectionOptions.h index 32af0a4eb..94584c7c3 100644 --- a/src/storm/solver/SolverSelectionOptions.h +++ b/src/storm/solver/SolverSelectionOptions.h @@ -6,7 +6,7 @@ namespace storm { namespace solver { - ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, Topological) + ExtendEnumsWithSelectionField(MinMaxMethod, PolicyIteration, ValueIteration, Topological, Acyclic) ExtendEnumsWithSelectionField(GameMethod, PolicyIteration, ValueIteration) ExtendEnumsWithSelectionField(LpSolverType, Gurobi, Glpk, Z3) diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp index 0a93a70f4..627dcd131 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp @@ -29,6 +29,7 @@ namespace storm { switch (method) { case MinMaxMethod::ValueIteration: this->solutionMethod = SolutionMethod::ValueIteration; break; case MinMaxMethod::PolicyIteration: this->solutionMethod = SolutionMethod::PolicyIteration; break; + case MinMaxMethod::Acyclic: this->solutionMethod = SolutionMethod::Acyclic; break; default: STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unsupported technique."); } @@ -91,6 +92,8 @@ namespace storm { return solveEquationsValueIteration(dir, x, b); case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::PolicyIteration: return solveEquationsPolicyIteration(dir, x, b); + case StandardMinMaxLinearEquationSolverSettings::SolutionMethod::Acyclic: + return solveEquationsAcyclic(dir, x, b); default: STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "This solver does not implement the selected solution method"); } @@ -293,6 +296,113 @@ namespace storm { return status == Status::Converged || status == Status::TerminatedEarly; } + template + bool StandardMinMaxLinearEquationSolver::solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const { + uint64_t numGroups = this->A.getRowGroupCount(); + + // Allocate memory for the scheduler (if required) + if (this->isTrackSchedulerSet()) { + if (this->schedulerChoices) { + this->schedulerChoices->resize(numGroups); + } else { + this->schedulerChoices = std::vector(numGroups); + } + } + + // We now compute a topological sort of the row groups. + // If caching is enabled, it might be possible to obtain this sort from the cache. + if (this->isCachingEnabled()) { + if (rowGroupOrdering) { + for (auto const& group : *rowGroupOrdering) { + computeOptimalValueForRowGroup(group, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[group] : nullptr); + } + return true; + } else { + rowGroupOrdering = std::make_unique>(); + rowGroupOrdering->reserve(numGroups); + } + } + + auto transposedMatrix = this->A.transpose(true); + + // We store the groups that have already been processed, i.e., the groups for which x[group] was already set to the correct value. + storm::storage::BitVector processedGroups(numGroups, false); + // Furthermore, we keep track of all candidate groups for which we still need to check whether this group can be processed now. + // A group can be processed if all successors have already been processed. + // Notice that the BitVector candidates is considered in a reversed way, i.e., group i is a candidate iff candidates.get(numGroups - i - 1) is true. + // This is due to the observation that groups with higher indices usually need to be processed earlier. + storm::storage::BitVector candidates(numGroups, true); + uint64_t candidate = numGroups - 1; + for (uint64_t numCandidates = candidates.size(); numCandidates > 0; --numCandidates) { + candidates.set(numGroups - candidate - 1, false); + + // Check if the candidate row group has an unprocessed successor + bool hasUnprocessedSuccessor = false; + for (auto const& entry : this->A.getRowGroup(candidate)) { + if (!processedGroups.get(entry.getColumn())) { + hasUnprocessedSuccessor = true; + break; + } + } + + uint64_t nextCandidate = numGroups - candidates.getNextSetIndex(numGroups - candidate - 1 + 1) - 1; + + if (!hasUnprocessedSuccessor) { + // This candidate can be processed. + processedGroups.set(candidate); + computeOptimalValueForRowGroup(candidate, dir, x, b, this->isTrackSchedulerSet() ? &this->schedulerChoices.get()[candidate] : nullptr); + if (this->isCachingEnabled()) { + rowGroupOrdering->push_back(candidate); + } + + // Add new candidates + for (auto const& predecessorEntry : transposedMatrix.getRow(candidate)) { + uint64_t predecessor = predecessorEntry.getColumn(); + if (!candidates.get(numGroups - predecessor - 1)) { + candidates.set(numGroups - predecessor - 1, true); + nextCandidate = std::max(nextCandidate, predecessor); + ++numCandidates; + } + } + } + candidate = nextCandidate; + } + + assert(candidates.empty()); + STORM_LOG_THROW(processedGroups.full(), storm::exceptions::InvalidOperationException, "The MinMax equation system is not acyclic."); + return true; + } + + template + void StandardMinMaxLinearEquationSolver::computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice) const { + uint64_t row = this->A.getRowGroupIndices()[group]; + uint64_t groupEnd = this->A.getRowGroupIndices()[group + 1]; + assert(row != groupEnd); + + auto bIt = b.begin() + row; + ValueType& xi = x[group]; + xi = this->A.multiplyRowWithVector(row, x) + *bIt; + uint64_t optimalRow = row; + + for (++row, ++bIt; row < groupEnd; ++row, ++bIt) { + ValueType choiceVal = this->A.multiplyRowWithVector(row, x) + *bIt; + if (minimize(dir)) { + if (choiceVal < xi) { + xi = choiceVal; + optimalRow = row; + } + } else { + if (choiceVal > xi) { + xi = choiceVal; + optimalRow = row; + } + } + } + if (choice != nullptr) { + *choice = optimalRow - this->A.getRowGroupIndices()[group]; + } + } + template void StandardMinMaxLinearEquationSolver::repeatedMultiply(OptimizationDirection dir, std::vector& x, std::vector const* b, uint_fast64_t n) const { if(!linEqSolverA) { @@ -356,6 +466,7 @@ namespace storm { linEqSolverA.reset(); auxiliaryRowVector.reset(); auxiliaryRowGroupVector.reset(); + rowGroupOrdering.reset(); MinMaxLinearEquationSolver::clearCache(); } diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.h b/src/storm/solver/StandardMinMaxLinearEquationSolver.h index f7e1e40c6..050ce0bde 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.h +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.h @@ -12,7 +12,7 @@ namespace storm { StandardMinMaxLinearEquationSolverSettings(); enum class SolutionMethod { - ValueIteration, PolicyIteration + ValueIteration, PolicyIteration, Acyclic }; void setSolutionMethod(SolutionMethod const& solutionMethod); @@ -51,9 +51,12 @@ namespace storm { private: bool solveEquationsPolicyIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool solveEquationsValueIteration(OptimizationDirection dir, std::vector& x, std::vector const& b) const; + bool solveEquationsAcyclic(OptimizationDirection dir, std::vector& x, std::vector const& b) const; bool valueImproved(OptimizationDirection dir, ValueType const& value1, ValueType const& value2) const; + void computeOptimalValueForRowGroup(uint_fast64_t group, OptimizationDirection dir, std::vector& x, std::vector const& b, uint_fast64_t* choice = nullptr) const; + enum class Status { Converged, TerminatedEarly, MaximalIterationsExceeded, InProgress }; @@ -62,7 +65,8 @@ namespace storm { mutable std::unique_ptr> linEqSolverA; mutable std::unique_ptr> auxiliaryRowVector; // A.rowCount() entries mutable std::unique_ptr> auxiliaryRowGroupVector; // A.rowGroupCount() entries - + mutable std::unique_ptr> rowGroupOrdering; // A.rowGroupCount() entries + Status updateStatusIfNotConverged(Status status, std::vector const& x, uint64_t iterations) const; void reportStatus(Status status, uint64_t iterations) const;