diff --git a/src/storm/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp index 8870456cd..af3bbc2cf 100644 --- a/src/storm/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/HybridDtmcPrctlHelper.cpp @@ -73,6 +73,7 @@ namespace storm { std::vector b = subvector.toVector(odd); std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitSubmatrix)); + solver->setBounds(storm::utility::zero(), storm::utility::one()); solver->solveEquations(x, b); // Return a hybrid check result that stores the numerical values explicitly. @@ -140,7 +141,6 @@ namespace storm { } } - template std::unique_ptr HybridDtmcPrctlHelper::computeInstantaneousRewards(storm::models::symbolic::Model const& model, storm::dd::Add const& transitionMatrix, RewardModelType const& rewardModel, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { // Only compute the result if the model has at least one reward this->getModel(). @@ -238,6 +238,7 @@ namespace storm { // Now solve the resulting equation system. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitSubmatrix)); + solver->setLowerBound(storm::utility::zero()); solver->solveEquations(x, b); // Return a hybrid check result that stores the numerical values explicitly. diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 5d84bd9fb..f768a2428 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -94,6 +94,7 @@ namespace storm { // Now solve the created system of linear equations. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(submatrix)); + solver->setBounds(storm::utility::zero(), storm::utility::one()); solver->solveEquations(x, b); // Set values of resulting vector according to result. @@ -217,6 +218,7 @@ namespace storm { // Now solve the resulting equation system. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(submatrix)); + solver->setLowerBound(storm::utility::zero()); solver->solveEquations(x, b); // Set values of resulting vector according to result. diff --git a/src/storm/modelchecker/results/ExplicitQuantitativeCheckResult.cpp b/src/storm/modelchecker/results/ExplicitQuantitativeCheckResult.cpp index a4dd5f96c..5ea662a57 100644 --- a/src/storm/modelchecker/results/ExplicitQuantitativeCheckResult.cpp +++ b/src/storm/modelchecker/results/ExplicitQuantitativeCheckResult.cpp @@ -122,10 +122,12 @@ namespace storm { ValueType sum = storm::utility::zero(); if (this->isResultForAllStates()) { for (auto& element : boost::get(values)) { + STORM_LOG_THROW(element != storm::utility::infinity(), storm::exceptions::InvalidOperationException, "Cannot compute the sum of values containing infinity."); sum += element; } } else { for (auto& element : boost::get(values)) { + STORM_LOG_THROW(element.second != storm::utility::infinity(), storm::exceptions::InvalidOperationException, "Cannot compute the sum of values containing infinity."); sum += element.second; } } @@ -139,11 +141,13 @@ namespace storm { ValueType sum = storm::utility::zero(); if (this->isResultForAllStates()) { for (auto& element : boost::get(values)) { + STORM_LOG_THROW(element != storm::utility::infinity(), storm::exceptions::InvalidOperationException, "Cannot compute the average of values containing infinity."); sum += element; } return sum / boost::get(values).size(); } else { for (auto& element : boost::get(values)) { + STORM_LOG_THROW(element.second != storm::utility::infinity(), storm::exceptions::InvalidOperationException, "Cannot compute the average of values containing infinity."); sum += element.second; } return sum / boost::get(values).size(); @@ -168,17 +172,47 @@ namespace storm { template void print(std::ostream& out, ValueType const& value) { - out << value; - if (std::is_same::value) { - out << " (approx. " << storm::utility::convertNumber(value) << ")"; + if (value == storm::utility::infinity()) { + out << "inf"; + } else { + out << value; + if (std::is_same::value) { + out << " (approx. " << storm::utility::convertNumber(value) << ")"; + } } } template - void printApproxRange(std::ostream& out, ValueType const& min, ValueType const& max) { + void printRange(std::ostream& out, ValueType const& min, ValueType const& max) { + out << "["; + if (min == storm::utility::infinity()) { + out << "inf"; + } else { + out << min; + } + out << ", "; + if (max == storm::utility::infinity()) { + out << "inf"; + } else { + out << max; + } + out << "]"; if (std::is_same::value) { - out << " (approx. [" << storm::utility::convertNumber(min) << ", " << storm::utility::convertNumber(max) << "])"; + out << " (approx. ["; + if (min == storm::utility::infinity()) { + out << "inf"; + } else { + out << storm::utility::convertNumber(min); + } + out << ", "; + if (max == storm::utility::infinity()) { + out << "inf"; + } else { + out << storm::utility::convertNumber(max); + } + out << "])"; } + out << " (range)"; } template @@ -228,9 +262,7 @@ namespace storm { if (printAsRange) { std::pair minmax = this->getMinMax(); - out << "[" << minmax.first << ", " << minmax.second << "]"; - printApproxRange(out, minmax.first, minmax.second); - out << " (range)"; + printRange(out, minmax.first, minmax.second); } return out; diff --git a/src/storm/solver/EigenLinearEquationSolver.cpp b/src/storm/solver/EigenLinearEquationSolver.cpp index bf80f4db9..b05fe6add 100644 --- a/src/storm/solver/EigenLinearEquationSolver.cpp +++ b/src/storm/solver/EigenLinearEquationSolver.cpp @@ -5,6 +5,7 @@ #include "storm/settings/SettingsManager.h" #include "storm/settings/modules/EigenEquationSolverSettings.h" +#include "storm/utility/vector.h" #include "storm/utility/macros.h" #include "storm/exceptions/InvalidSettingsException.h" @@ -230,6 +231,9 @@ namespace storm { } } + // Make sure that all results conform to the bounds. + storm::utility::vector::clip(x, this->lowerBound, this->upperBound); + // Check if the solver converged and issue a warning otherwise. if (converged) { STORM_LOG_DEBUG("Iterative solver converged after " << numberOfIterations << " iterations."); @@ -240,7 +244,7 @@ namespace storm { } } - return false; + return true; } template diff --git a/src/storm/solver/GmmxxLinearEquationSolver.cpp b/src/storm/solver/GmmxxLinearEquationSolver.cpp index c92817ac2..31c358a22 100644 --- a/src/storm/solver/GmmxxLinearEquationSolver.cpp +++ b/src/storm/solver/GmmxxLinearEquationSolver.cpp @@ -13,6 +13,7 @@ #include "storm/solver/NativeLinearEquationSolver.h" #include "storm/utility/gmm.h" +#include "storm/utility/vector.h" namespace storm { namespace solver { @@ -191,6 +192,9 @@ namespace storm { clearCache(); } + // Make sure that all results conform to the bounds. + storm::utility::vector::clip(x, this->lowerBound, this->upperBound); + // Check if the solver converged and issue a warning otherwise. if (iter.converged()) { STORM_LOG_DEBUG("Iterative solver converged after " << iter.get_iteration() << " iterations."); @@ -202,6 +206,9 @@ namespace storm { } else if (method == GmmxxLinearEquationSolverSettings::SolutionMethod::Jacobi) { uint_fast64_t iterations = solveLinearEquationSystemWithJacobi(x, b); + // Make sure that all results conform to the bounds. + storm::utility::vector::clip(x, this->lowerBound, this->upperBound); + // Check if the solver converged and issue a warning otherwise. if (iterations < this->getSettings().getMaximalNumberOfIterations()) { STORM_LOG_DEBUG("Iterative solver converged after " << iterations << " iterations."); diff --git a/src/storm/solver/LinearEquationSolver.cpp b/src/storm/solver/LinearEquationSolver.cpp index ca1eeda94..487b77f89 100644 --- a/src/storm/solver/LinearEquationSolver.cpp +++ b/src/storm/solver/LinearEquationSolver.cpp @@ -74,6 +74,22 @@ namespace storm { cachedRowVector.reset(); } + template + void LinearEquationSolver::setLowerBound(ValueType const& value) { + lowerBound = value; + } + + template + void LinearEquationSolver::setUpperBound(ValueType const& value) { + upperBound = value; + } + + template + void LinearEquationSolver::setBounds(ValueType const& lower, ValueType const& upper) { + setLowerBound(lower); + setUpperBound(upper); + } + template std::unique_ptr> LinearEquationSolverFactory::create(storm::storage::SparseMatrix&& matrix) const { return create(matrix); diff --git a/src/storm/solver/LinearEquationSolver.h b/src/storm/solver/LinearEquationSolver.h index 4912070d7..b94626cd2 100644 --- a/src/storm/solver/LinearEquationSolver.h +++ b/src/storm/solver/LinearEquationSolver.h @@ -70,7 +70,7 @@ namespace storm { /*! * Sets whether some of the generated data during solver calls should be cached. - * This possibly decreases the runtime of subsequent calls but also increases memory consumption. + * This possibly increases the runtime of subsequent calls but also increases memory consumption. */ void setCachingEnabled(bool value) const; @@ -83,11 +83,32 @@ namespace storm { * Clears the currently cached data that has been stored during previous calls of the solver. */ virtual void clearCache() const; + + /*! + * Sets a lower bound for the solution that can potentially used by the solver. + */ + void setLowerBound(ValueType const& value); + + /*! + * Sets an upper bound for the solution that can potentially used by the solver. + */ + void setUpperBound(ValueType const& value); + + /*! + * Sets bounds for the solution that can potentially used by the solver. + */ + void setBounds(ValueType const& lower, ValueType const& upper); protected: // auxiliary storage. If set, this vector has getMatrixRowCount() entries. mutable std::unique_ptr> cachedRowVector; - + + // A lower bound if one was set. + boost::optional lowerBound; + + // An upper bound if one was set. + boost::optional upperBound; + private: /*! * Retrieves the row count of the matrix associated with this solver. @@ -101,7 +122,6 @@ namespace storm { /// Whether some of the generated data during solver calls should be cached. mutable bool cachingEnabled; - }; template diff --git a/src/storm/solver/MinMaxLinearEquationSolver.cpp b/src/storm/solver/MinMaxLinearEquationSolver.cpp index e7a332e13..02269a248 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/MinMaxLinearEquationSolver.cpp @@ -98,6 +98,21 @@ namespace storm { // Intentionally left empty. } + template + void MinMaxLinearEquationSolver::setLowerBound(ValueType const& value) { + lowerBound = value; + } + + template + void MinMaxLinearEquationSolver::setUpperBound(ValueType const& value) { + upperBound = value; + } + + template + void MinMaxLinearEquationSolver::setBounds(ValueType const& lower, ValueType const& upper) { + setLowerBound(lower); + setUpperBound(upper); + } template MinMaxLinearEquationSolverFactory::MinMaxLinearEquationSolverFactory(bool trackScheduler) : trackScheduler(trackScheduler) { diff --git a/src/storm/solver/MinMaxLinearEquationSolver.h b/src/storm/solver/MinMaxLinearEquationSolver.h index 592f00067..473438a38 100644 --- a/src/storm/solver/MinMaxLinearEquationSolver.h +++ b/src/storm/solver/MinMaxLinearEquationSolver.h @@ -145,6 +145,20 @@ namespace storm { */ virtual void clearCache() const; + /*! + * Sets a lower bound for the solution that can potentially used by the solver. + */ + void setLowerBound(ValueType const& value); + + /*! + * Sets an upper bound for the solution that can potentially used by the solver. + */ + void setUpperBound(ValueType const& value); + + /*! + * Sets bounds for the solution that can potentially used by the solver. + */ + void setBounds(ValueType const& lower, ValueType const& upper); protected: /// The optimization direction to use for calls to functions that do not provide it explicitly. Can also be unset. @@ -156,6 +170,12 @@ namespace storm { /// The scheduler (if it could be successfully generated). mutable boost::optional> scheduler; + // A lower bound if one was set. + boost::optional lowerBound; + + // An upper bound if one was set. + boost::optional upperBound; + private: /// Whether some of the generated data during solver calls should be cached. bool cachingEnabled; diff --git a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp index 7a6e54df5..c9c405341 100644 --- a/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp +++ b/src/storm/solver/StandardMinMaxLinearEquationSolver.cpp @@ -112,6 +112,12 @@ namespace storm { // Create a solver that we will use throughout the procedure. We will modify the matrix in each iteration. auto solver = linearEquationSolverFactory->create(std::move(submatrix)); + if (this->lowerBound) { + solver->setLowerBound(this->lowerBound.get()); + } + if (this->upperBound) { + solver->setUpperBound(this->upperBound.get()); + } solver->setCachingEnabled(true); Status status = Status::InProgress; diff --git a/src/storm/utility/constants.cpp b/src/storm/utility/constants.cpp index 7eeae8106..a6326d5e8 100644 --- a/src/storm/utility/constants.cpp +++ b/src/storm/utility/constants.cpp @@ -226,6 +226,26 @@ namespace storm { return std::make_pair(min, max); } + template<> + std::pair minmax(std::vector const& values) { + assert(!values.empty()); + storm::RationalNumber min = values.front(); + storm::RationalNumber max = values.front(); + for (auto const& vt : values) { + if (vt == storm::utility::infinity()) { + max = vt; + } else { + if (vt < min) { + min = vt; + } + if (vt > max) { + max = vt; + } + } + } + return std::make_pair(min, max); + } + template<> storm::RationalFunction minimum(std::vector const&) { STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Minimum for rational functions is not defined."); @@ -251,7 +271,7 @@ namespace storm { STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Maximum/maximum for rational functions is not defined."); } - template< typename K, typename ValueType> + template std::pair minmax(std::map const& values) { assert(!values.empty()); ValueType min = values.begin()->second; @@ -267,6 +287,26 @@ namespace storm { return std::make_pair(min, max); } + template<> + std::pair minmax(std::map const& values) { + assert(!values.empty()); + storm::RationalNumber min = values.begin()->second; + storm::RationalNumber max = values.begin()->second; + for (auto const& vt : values) { + if (vt.second == storm::utility::infinity()) { + max = vt.second; + } else { + if (vt.second < min) { + min = vt.second; + } + if (vt.second > max) { + max = vt.second; + } + } + } + return std::make_pair(min, max); + } + template<> storm::RationalFunction minimum(std::map const&) { STORM_LOG_THROW(false, storm::exceptions::InvalidArgumentException, "Minimum for rational functions is not defined."); @@ -274,14 +314,7 @@ namespace storm { template< typename K, typename ValueType> ValueType minimum(std::map const& values) { - assert(!values.empty()); - ValueType min = values.begin()->second; - for (auto const& vt : values) { - if (vt.second < min) { - min = vt.second; - } - } - return min; + return minmax(values).first; } template<> @@ -291,14 +324,7 @@ namespace storm { template ValueType maximum(std::map const& values) { - assert(!values.empty()); - ValueType max = values.begin()->second; - for (auto const& vt : values) { - if (vt.second > max) { - max = vt.second; - } - } - return max; + return minmax(values).second; } #ifdef STORM_HAVE_CARL diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index dc22fa310..578a942a8 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -12,6 +12,8 @@ #include #include +#include + #include "storm/storage/BitVector.h" #include "storm/utility/constants.h" #include "storm/utility/macros.h" @@ -738,6 +740,20 @@ namespace storm { return true; } + /*! + * Takes the input vector and ensures that all entries conform to the bounds. + */ + template + void clip(std::vector& x, boost::optional const& lowerBound, boost::optional const& upperBound) { + for (auto& entry : x) { + if (lowerBound && entry < lowerBound.get()) { + entry = lowerBound.get(); + } else if (upperBound && entry > upperBound.get()) { + entry = upperBound.get(); + } + } + } + /*! * Takes the given offset vector and applies the given contraint. That is, it produces another offset vector that contains * the relative offsets of the entries given by the constraint.