Storm
A Modern Probabilistic Model Checker
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TopologicalMinMaxLinearEquationSolver.cpp
Go to the documentation of this file.
2
5
15
16namespace storm {
17namespace solver {
18
19template<typename ValueType, typename SolutionType>
23
24template<typename ValueType, typename SolutionType>
29
30template<typename ValueType, typename SolutionType>
35
36template<typename ValueType, typename SolutionType>
38 bool adaptPrecision) const {
39 storm::Environment subEnv(env);
40 subEnv.solver().minMax().setMethod(env.solver().topological().getUnderlyingMinMaxMethod(),
42 if (adaptPrecision) {
43 STORM_LOG_ASSERT(this->longestSccChainSize, "Did not compute the longest SCC chain size although it is needed.");
44 storm::RationalNumber subEnvPrec =
45 subEnv.solver().minMax().getPrecision() / storm::utility::convertNumber<storm::RationalNumber>(this->longestSccChainSize.get());
46 subEnv.solver().minMax().setPrecision(subEnvPrec);
47 }
48 return subEnv;
49}
50
51template<typename ValueType, typename SolutionType>
53 std::vector<SolutionType>& x,
54 std::vector<ValueType> const& b) const {
55 STORM_LOG_ASSERT(x.size() == this->A->getRowGroupCount(), "Provided x-vector has invalid size.");
56 STORM_LOG_ASSERT(b.size() == this->A->getRowCount(), "Provided b-vector has invalid size.");
57
58 // For sound computations we need to increase the precision in each SCC
59 bool needAdaptPrecision = env.solver().isForceSoundness();
60
61 if (!this->sortedSccDecomposition || (needAdaptPrecision && !this->longestSccChainSize)) {
62 STORM_LOG_TRACE("Creating SCC decomposition.");
63 storm::utility::Stopwatch sccSw(true);
64 createSortedSccDecomposition(needAdaptPrecision);
65 sccSw.stop();
66 STORM_LOG_INFO("SCC decomposition computed in "
67 << sccSw << ". Found " << this->sortedSccDecomposition->size() << " SCC(s) containing a total of " << x.size()
68 << " states. Average SCC size is "
69 << static_cast<double>(this->A->getRowGroupCount()) / static_cast<double>(this->sortedSccDecomposition->size()) << ".");
70 }
71
72 // We do not need to adapt the precision if all SCCs are trivial (i.e., the system is acyclic)
73 needAdaptPrecision = needAdaptPrecision && (this->sortedSccDecomposition->size() != this->A->getRowGroupCount());
74
75 storm::Environment sccSolverEnvironment = getEnvironmentForUnderlyingSolver(env, needAdaptPrecision);
76
77 if (this->longestSccChainSize) {
78 STORM_LOG_INFO("Longest SCC chain size is " << this->longestSccChainSize.get());
79 }
80
81 bool returnValue = true;
82 if (this->sortedSccDecomposition->size() == 1 && (!this->choiceFixedForRowGroup || this->choiceFixedForRowGroup.get().empty())) {
83 // Handle the case where there is just one large SCC, as there are no fixed choices for states, we solve it like this
84 if (auto const& scc = *this->sortedSccDecomposition->begin(); scc.size() == 1) {
85 // Catch the trivial case where the whole system is just a single state.
86 if (this->isTrackSchedulerSet()) {
87 this->schedulerChoices = std::vector<uint64_t>(1);
88 }
89 returnValue = solveTrivialScc(*scc.begin(), dir, x, b);
90 } else {
91 returnValue = solveFullyConnectedEquationSystem(sccSolverEnvironment, dir, x, b);
92 }
93 } else {
94 // Solve each SCC individually
95 if (this->isTrackSchedulerSet()) {
96 if (this->schedulerChoices) {
97 this->schedulerChoices.get().resize(x.size());
98 } else {
99 this->schedulerChoices = std::vector<uint64_t>(x.size());
100 }
101 }
102 std::optional<storm::storage::BitVector> newRelevantValues;
103 if (env.solver().topological().isExtendRelevantValues() && this->hasRelevantValues() &&
104 this->sortedSccDecomposition->size() < this->A->getRowGroupCount()) {
105 newRelevantValues = this->getRelevantValues();
106 // Extend the relevant values towards those that have an incoming transition from another SCC
107 std::vector<uint64_t> rowGroupToScc = this->sortedSccDecomposition->computeStateToSccIndexMap(this->A->getRowGroupCount());
108 for (uint64_t rowGroup = 0; rowGroup < this->A->getRowGroupCount(); ++rowGroup) {
109 auto currScc = rowGroupToScc[rowGroup];
110 for (auto const& successor : this->A->getRowGroup(rowGroup)) {
111 if (rowGroupToScc[successor.getColumn()] != currScc) {
112 newRelevantValues->set(successor.getColumn(), true);
113 }
114 }
115 }
116 }
117 storm::storage::BitVector sccRowGroupsAsBitVector(x.size(), false);
118 storm::storage::BitVector sccRowsAsBitVector(b.size(), false);
119 uint64_t sccIndex = 0;
120 storm::utility::ProgressMeasurement progress("states");
121 progress.setMaxCount(x.size());
122 progress.startNewMeasurement(0);
123 for (auto const& scc : *this->sortedSccDecomposition) {
124 if (scc.size() == 1) {
125 returnValue = solveTrivialScc(*scc.begin(), dir, x, b) && returnValue;
126 } else {
127 STORM_LOG_TRACE("Solving SCC of size " << scc.size() << ".");
128 sccRowGroupsAsBitVector.clear();
129 sccRowsAsBitVector.clear();
130 for (auto const& group : scc) { // Group refers to state
131 sccRowGroupsAsBitVector.set(group, true);
132
133 if (!this->choiceFixedForRowGroup || !this->choiceFixedForRowGroup.get()[group]) {
134 for (uint64_t row = this->A->getRowGroupIndices()[group]; row < this->A->getRowGroupIndices()[group + 1]; ++row) {
135 sccRowsAsBitVector.set(row, true);
136 }
137 } else {
138 auto row = this->A->getRowGroupIndices()[group] + this->getInitialScheduler()[group];
139 sccRowsAsBitVector.set(row, true);
140 STORM_LOG_TRACE("Fixing state " << group << " to choice " << this->getInitialScheduler()[group] << ".");
141 }
142 }
143 returnValue = solveScc(sccSolverEnvironment, dir, sccRowGroupsAsBitVector, sccRowsAsBitVector, x, b, newRelevantValues) && returnValue;
144 }
145 ++sccIndex;
146 progress.updateProgress(sccIndex);
148 STORM_LOG_WARN("Topological solver aborted after analyzing " << sccIndex << "/" << this->sortedSccDecomposition->size() << " SCCs.");
149 break;
150 }
151 }
152
153 // If requested, we store the scheduler for retrieval.
154 if (this->isTrackSchedulerSet()) {
155 if (!auxiliaryRowGroupVector) {
156 auxiliaryRowGroupVector = std::make_unique<std::vector<ValueType>>(this->A->getRowGroupCount());
157 }
158 this->schedulerChoices = std::vector<uint_fast64_t>(this->A->getRowGroupCount());
159 this->A->multiplyAndReduce(dir, this->A->getRowGroupIndices(), x, &b, *auxiliaryRowGroupVector.get(), &this->schedulerChoices.get());
160 }
161 }
162
163 if (!this->isCachingEnabled()) {
164 clearCache();
165 }
166
167 return returnValue;
168}
169
170template<typename ValueType, typename SolutionType>
172 // Obtain the scc decomposition
173 this->sortedSccDecomposition = std::make_unique<storm::storage::StronglyConnectedComponentDecomposition<ValueType>>(
175 if (needLongestChainSize) {
176 this->longestSccChainSize = this->sortedSccDecomposition->getMaxSccDepth() + 1;
177 }
178}
179
180template<typename ValueType, typename SolutionType>
181bool TopologicalMinMaxLinearEquationSolver<ValueType, SolutionType>::solveTrivialScc(uint64_t const& sccState, OptimizationDirection dir,
182 std::vector<ValueType>& globalX,
183 std::vector<ValueType> const& globalB) const {
184 ValueType& xi = globalX[sccState];
185 if (this->choiceFixedForRowGroup && this->choiceFixedForRowGroup.get()[sccState]) {
186 // if the choice in the scheduler is fixed we only update for the fixed choice
187 uint_fast64_t row = this->A->getRowGroupIndices()[sccState] + this->getInitialScheduler()[sccState];
188 ValueType rowValue = globalB[row];
189 bool hasDiagonalEntry = false;
190 ValueType denominator;
191 for (auto const& entry : this->A->getRow(row)) {
192 if (entry.getColumn() == sccState) {
193 hasDiagonalEntry = true;
194 denominator = storm::utility::one<ValueType>() - entry.getValue();
195 } else {
196 rowValue += entry.getValue() * globalX[entry.getColumn()];
197 }
198 }
199 if (hasDiagonalEntry) {
202 "State " << sccState << " has a selfloop with probability '1-(" << denominator << ")'. This could be an indication for numerical issues.");
203 if (storm::utility::isZero(denominator)) {
204 // In this case we have a selfloop on this state. This can never an optimal choice:
205 // When minimizing, we are looking for the largest fixpoint (which will never be attained by this action)
206 // When maximizing, this choice reflects probability zero (non-optimal) or reward infinity (should already be handled during preprocessing).
207 } else {
208 rowValue /= denominator;
209 }
210 }
211 xi = std::move(rowValue);
212 } else {
213 bool firstRow = true;
214 uint64_t bestRow;
215 for (uint64_t row = this->A->getRowGroupIndices()[sccState]; row < this->A->getRowGroupIndices()[sccState + 1]; ++row) {
216 ValueType rowValue = globalB[row];
217 bool hasDiagonalEntry = false;
219 for (auto const& entry : this->A->getRow(row)) {
220 if (entry.getColumn() == sccState) {
221 hasDiagonalEntry = true;
222 denominator = storm::utility::one<ValueType>() - entry.getValue();
223 } else {
224 rowValue += entry.getValue() * globalX[entry.getColumn()];
225 }
226 }
227 if (hasDiagonalEntry) {
230 "State " << sccState << " has a selfloop with probability '1-(" << denominator << ")'. This could be an indication for numerical issues.");
231 if (storm::utility::isZero(denominator)) {
232 // In this case we have a selfloop on this state. This can never an optimal choice:
233 // When minimizing, we are looking for the largest fixpoint (which will never be attained by this action)
234 // When maximizing, this choice reflects probability zero (non-optimal) or reward infinity (should already be handled during preprocessing).
235 continue;
236 } else {
237 rowValue /= denominator;
238 }
239 }
240 if (firstRow) {
241 xi = std::move(rowValue);
242 bestRow = row;
243 firstRow = false;
244 } else {
245 if (minimize(dir)) {
246 if (rowValue < xi) {
247 xi = std::move(rowValue);
248 bestRow = row;
249 }
250 } else {
251 if (rowValue > xi) {
252 xi = std::move(rowValue);
253 bestRow = row;
254 }
255 }
256 }
257 }
258 if (this->isTrackSchedulerSet()) {
259 this->schedulerChoices.get()[sccState] = bestRow - this->A->getRowGroupIndices()[sccState];
260 }
261 STORM_LOG_THROW(!firstRow, storm::exceptions::UnexpectedException, "Empty row group in MinMax equation system.");
262 }
263 return true;
264}
265
266template<typename ValueType, typename SolutionType>
267bool TopologicalMinMaxLinearEquationSolver<ValueType, SolutionType>::solveFullyConnectedEquationSystem(storm::Environment const& sccSolverEnvironment,
268 OptimizationDirection dir, std::vector<SolutionType>& x,
269 std::vector<ValueType> const& b) const {
270 STORM_LOG_ASSERT(!this->choiceFixedForRowGroup || this->choiceFixedForRowGroup.get().empty(),
271 "Expecting no fixed choices for states when solving the fully connected equation system");
272 if (!this->sccSolver) {
273 this->sccSolver = GeneralMinMaxLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment);
274 this->sccSolver->setCachingEnabled(true);
275 }
276 this->sccSolver->setMatrix(*this->A);
277 this->sccSolver->setHasUniqueSolution(this->hasUniqueSolution());
278 this->sccSolver->setHasNoEndComponents(this->hasNoEndComponents());
279 this->sccSolver->setTrackScheduler(this->isTrackSchedulerSet());
280 if (this->hasInitialScheduler()) {
281 auto choices = this->getInitialScheduler();
282 this->sccSolver->setInitialScheduler(std::move(choices));
283 }
284 if (this->hasRelevantValues()) {
285 this->sccSolver->setRelevantValues(this->getRelevantValues());
286 }
287 auto req = this->sccSolver->getRequirements(sccSolverEnvironment, dir);
288 this->sccSolver->setBoundsFromOtherSolver(*this);
289 if (req.upperBounds() && this->hasUpperBound()) {
290 req.clearUpperBounds();
291 }
292 if (req.lowerBounds() && this->hasLowerBound()) {
293 req.clearLowerBounds();
294 }
295 if (req.validInitialScheduler() && this->hasInitialScheduler()) {
296 req.clearValidInitialScheduler();
297 }
298 if (req.uniqueSolution() && this->hasUniqueSolution()) {
299 req.clearUniqueSolution();
300 }
301 STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException,
302 "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked.");
303 this->sccSolver->setRequirementsChecked(true);
304
305 bool res = this->sccSolver->solveEquations(sccSolverEnvironment, dir, x, b);
306 if (this->isTrackSchedulerSet()) {
307 this->schedulerChoices = this->sccSolver->getSchedulerChoices();
308 }
309 return res;
310}
311
312template<typename ValueType, typename SolutionType>
313bool TopologicalMinMaxLinearEquationSolver<ValueType, SolutionType>::solveScc(storm::Environment const& sccSolverEnvironment, OptimizationDirection dir,
314 storm::storage::BitVector const& sccRowGroups,
315 storm::storage::BitVector const& sccRows, std::vector<ValueType>& globalX,
316 std::vector<ValueType> const& globalB,
317 std::optional<storm::storage::BitVector> const& globalRelevantValues) const {
318 // Set up the SCC solver
319 if (!this->sccSolver) {
320 this->sccSolver = GeneralMinMaxLinearEquationSolverFactory<ValueType>().create(sccSolverEnvironment);
321 this->sccSolver->setCachingEnabled(true);
322 }
323 this->sccSolver->setHasUniqueSolution(this->hasUniqueSolution());
324 this->sccSolver->setHasNoEndComponents(this->hasNoEndComponents());
325 this->sccSolver->setTrackScheduler(this->isTrackSchedulerSet());
326 if (globalRelevantValues) {
327 this->sccSolver->setRelevantValues((*globalRelevantValues) % sccRowGroups);
328 }
329
331 if (this->choiceFixedForRowGroup) {
332 // Obtain choiceFixedForState bitvector containing only the states of the scc.
333 storm::storage::BitVector choiceFixedForStateSCC = this->choiceFixedForRowGroup.get() % sccRowGroups;
334 sccA = this->A->getSubmatrix(false, sccRows, sccRowGroups);
335
336 // initial scheduler
337 if (this->hasInitialScheduler()) {
338 std::vector<uint_fast64_t> sccInitChoices = storm::utility::vector::filterVector(this->getInitialScheduler(), sccRowGroups);
339 // As we removed the entries where the choice was fixed, we need to change the scheduler.
340 // We set the scheduler to 0 for those states.
341 storm::utility::vector::setVectorValues<uint_fast64_t>(sccInitChoices, choiceFixedForStateSCC, 0);
342 this->sccSolver->setInitialScheduler(std::move(sccInitChoices));
343 }
344
345 } else {
346 sccA = this->A->getSubmatrix(true, sccRowGroups, sccRowGroups);
347
348 // initial scheduler
349 if (this->hasInitialScheduler()) {
350 auto sccInitChoices = storm::utility::vector::filterVector(this->getInitialScheduler(), sccRowGroups);
351 this->sccSolver->setInitialScheduler(std::move(sccInitChoices));
352 }
353 }
354
355 this->sccSolver->setMatrix(std::move(sccA));
356
357 // x Vector
358 auto sccX = storm::utility::vector::filterVector(globalX, sccRowGroups);
359
360 // b Vector
361 std::vector<ValueType> sccB;
362 sccB.reserve(sccRows.getNumberOfSetBits());
363 for (auto row : sccRows) {
364 ValueType bi = globalB[row];
365 for (auto const& entry : this->A->getRow(row)) {
366 if (!sccRowGroups.get(entry.getColumn())) {
367 bi += entry.getValue() * globalX[entry.getColumn()];
368 }
369 }
370 sccB.push_back(std::move(bi));
371 }
372
373 auto req = this->sccSolver->getRequirements(sccSolverEnvironment, dir);
374 this->sccSolver->clearBounds();
375 // lower/upper bounds
377 this->sccSolver->setLowerBound(this->getLowerBound());
378 req.clearLowerBounds();
380 this->sccSolver->setLowerBounds(storm::utility::vector::filterVector(this->getLowerBounds(), sccRowGroups));
381 req.clearLowerBounds();
382 }
384 this->sccSolver->setUpperBound(this->getUpperBound());
385 req.clearUpperBounds();
387 this->sccSolver->setUpperBounds(storm::utility::vector::filterVector(this->getUpperBounds(), sccRowGroups));
388 req.clearUpperBounds();
389 }
390
391 // Requirements
392 if (req.validInitialScheduler() && (this->hasInitialScheduler() || this->hasNoEndComponents())) {
393 req.clearValidInitialScheduler();
394 }
395 if (req.uniqueSolution() && this->hasUniqueSolution()) {
396 req.clearUniqueSolution();
397 }
398 STORM_LOG_THROW(!req.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException,
399 "Solver requirements " + req.getEnabledRequirementsAsString() + " not checked.");
400 this->sccSolver->setRequirementsChecked(true);
401
402 // Invoke scc solver
403 bool res = this->sccSolver->solveEquations(sccSolverEnvironment, dir, sccX, sccB);
404
405 // Set Scheduler choices
406 if (this->isTrackSchedulerSet()) {
407 storm::utility::vector::setVectorValues(this->schedulerChoices.get(), sccRowGroups, this->sccSolver->getSchedulerChoices());
408 }
409
410 // Set solution
411 storm::utility::vector::setVectorValues(globalX, sccRowGroups, sccX);
412
413 return res;
414}
415
416template<typename ValueType, typename SolutionType>
418 Environment const& env, boost::optional<storm::solver::OptimizationDirection> const& direction, bool const& hasInitialScheduler) const {
419 // Return the requirements of the underlying solver
420 return GeneralMinMaxLinearEquationSolverFactory<ValueType>().getRequirements(getEnvironmentForUnderlyingSolver(env), this->hasUniqueSolution(),
421 this->hasNoEndComponents(), direction, hasInitialScheduler,
422 this->isTrackSchedulerSet());
423}
424
425template<typename ValueType, typename SolutionType>
427 sortedSccDecomposition.reset();
428 longestSccChainSize = boost::none;
429 sccSolver.reset();
430 auxiliaryRowGroupVector.reset();
432}
433
434// Explicitly instantiate the min max linear equation solver.
437
438} // namespace solver
439} // namespace storm
SolverEnvironment & solver()
TopologicalSolverEnvironment & topological()
storm::solver::MinMaxMethod const & getUnderlyingMinMaxMethod() const
MinMaxLinearEquationSolverRequirements getRequirements(Environment const &env, bool hasUniqueSolution=false, bool hasNoEndComponents=false, boost::optional< storm::solver::OptimizationDirection > const &direction=boost::none, bool hasInitialScheduler=false, bool trackScheduler=false) const
Retrieves the requirements of the solver that would be created when calling create() right now.
virtual void clearCache() const
Clears the currently cached data that has been stored during previous calls of the solver.
virtual bool internalSolveEquations(storm::Environment const &env, OptimizationDirection d, std::vector< SolutionType > &x, std::vector< ValueType > const &b) const override
virtual void clearCache() const override
Clears the currently cached data that has been stored during previous calls of the solver.
virtual MinMaxLinearEquationSolverRequirements getRequirements(Environment const &env, boost::optional< storm::solver::OptimizationDirection > const &direction=boost::none, bool const &hasInitialScheduler=false) const override
Retrieves the requirements of this solver for solving equations with the current settings.
A bit vector that is internally represented as a vector of 64-bit values.
Definition BitVector.h:18
void clear()
Removes all set bits from the bit vector.
void set(uint_fast64_t index, bool value=true)
Sets the given truth value at the given index.
uint_fast64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
bool get(uint_fast64_t index) const
Retrieves the truth value of the bit at the given index and performs a bound check.
A class that holds a possibly non-square matrix in the compressed row storage format.
SparseMatrix getSubmatrix(bool useGroups, storm::storage::BitVector const &rowConstraint, storm::storage::BitVector const &columnConstraint, bool insertDiagonalEntries=false, storm::storage::BitVector const &makeZeroColumns=storm::storage::BitVector()) const
Creates a submatrix of the current matrix by dropping all rows and columns whose bits are not set to ...
A class that provides convenience operations to display run times.
bool updateProgress(uint64_t count)
Updates the progress to the current count and prints it if the delay passed.
void setMaxCount(uint64_t maxCount)
Sets the maximal possible count.
void startNewMeasurement(uint64_t startCount)
Starts a new measurement, dropping all progress information collected so far.
A class that provides convenience operations to display run times.
Definition Stopwatch.h:14
void stop()
Stop stopwatch and add measured time to total time.
Definition Stopwatch.cpp:42
#define STORM_LOG_INFO(message)
Definition logging.h:29
#define STORM_LOG_WARN(message)
Definition logging.h:30
#define STORM_LOG_TRACE(message)
Definition logging.h:17
#define STORM_LOG_ASSERT(cond, message)
Definition macros.h:11
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
#define STORM_LOG_WARN_COND_DEBUG(cond, message)
Definition macros.h:18
SFTBDDChecker::ValueType ValueType
bool constexpr minimize(OptimizationDirection d)
bool isTerminate()
Check whether the program should terminate (due to some abort signal).
void setVectorValues(std::vector< T > &vector, storm::storage::BitVector const &positions, std::vector< T > const &values)
Sets the provided values at the provided positions in the given vector.
Definition vector.h:83
std::vector< Type > filterVector(std::vector< Type > const &in, storm::storage::BitVector const &filter)
Definition vector.h:1213
NumberTraits< RationalType >::IntegerType denominator(RationalType const &number)
bool isAlmostZero(ValueType const &a)
Definition constants.cpp:56
bool isZero(ValueType const &a)
Definition constants.cpp:41
LabParser.cpp.
Definition cli.cpp:18
StronglyConnectedComponentDecompositionOptions & computeSccDepths(bool value=true)
Sets if scc depths can be retrieved.
StronglyConnectedComponentDecompositionOptions & forceTopologicalSort(bool value=true)
Enforces that the returned SCCs are sorted in a topological order.