Storm 1.10.0.1
A Modern Probabilistic Model Checker
Loading...
Searching...
No Matches
SparsePcaaQuery.cpp
Go to the documentation of this file.
2
5#include "storm/io/export.h"
16
18
19namespace storm {
20namespace modelchecker {
21namespace multiobjective {
22
23template<class SparseModelType, typename GeometryValueType>
25 if (storm::settings::getModule<storm::settings::modules::CoreSettings>().isShowStatisticsSet()) {
26 STORM_PRINT_AND_LOG("Pareto Curve Approximation Algorithm terminated after " << this->refinementSteps.size() << " refinement steps.\n");
27 }
28}
29
30template<class SparseModelType, typename GeometryValueType>
39
40template<class SparseModelType, typename GeometryValueType>
42 Point const& pointToBeSeparated) {
43 STORM_LOG_DEBUG("Searching a weight vector to seperate the point given by "
44 << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(pointToBeSeparated)) << ".");
45
46 if (underApproximation->isEmpty()) {
47 // In this case, every weight vector is separating
48 uint_fast64_t objIndex = diracWeightVectorsToBeChecked.getNextSetIndex(0) % pointToBeSeparated.size();
49 WeightVector result(pointToBeSeparated.size(), storm::utility::zero<GeometryValueType>());
50 result[objIndex] = storm::utility::one<GeometryValueType>();
51 diracWeightVectorsToBeChecked.set(objIndex, false);
52 return result;
53 }
54
55 // Reaching this point means that the underApproximation contains halfspaces. The seperating vector has to be the normal vector of one of these halfspaces.
56 // We pick one with maximal distance to the given point. However, Dirac weight vectors that only assign a non-zero weight to a single objective take
57 // precedence.
58 STORM_LOG_ASSERT(!underApproximation->contains(pointToBeSeparated),
59 "Tried to find a separating point but the point is already contained in the underApproximation");
60 std::vector<storm::storage::geometry::Halfspace<GeometryValueType>> halfspaces = underApproximation->getHalfspaces();
61 uint_fast64_t farestHalfspaceIndex = halfspaces.size();
62 GeometryValueType farestDistance = -storm::utility::one<GeometryValueType>();
63 bool foundSeparatingDiracVector = false;
64 for (uint_fast64_t halfspaceIndex = 0; halfspaceIndex < halfspaces.size(); ++halfspaceIndex) {
65 GeometryValueType distance = halfspaces[halfspaceIndex].euclideanDistance(pointToBeSeparated);
66 if (!storm::utility::isZero(distance)) {
67 storm::storage::BitVector nonZeroVectorEntries = ~storm::utility::vector::filterZero<GeometryValueType>(halfspaces[halfspaceIndex].normalVector());
68 bool isSingleObjectiveVector =
69 nonZeroVectorEntries.getNumberOfSetBits() == 1 && diracWeightVectorsToBeChecked.get(nonZeroVectorEntries.getNextSetIndex(0));
70 // Check if this halfspace is a better candidate than the current one
71 if ((!foundSeparatingDiracVector && isSingleObjectiveVector) ||
72 (foundSeparatingDiracVector == isSingleObjectiveVector && distance > farestDistance)) {
73 foundSeparatingDiracVector = foundSeparatingDiracVector || isSingleObjectiveVector;
74 farestHalfspaceIndex = halfspaceIndex;
75 farestDistance = distance;
76 }
77 }
78 }
79 if (foundSeparatingDiracVector) {
80 diracWeightVectorsToBeChecked &= storm::utility::vector::filterZero<GeometryValueType>(halfspaces[farestHalfspaceIndex].normalVector());
81 }
82
83 STORM_LOG_THROW(farestHalfspaceIndex < halfspaces.size(), storm::exceptions::UnexpectedException, "There is no seperating vector.");
84 STORM_LOG_DEBUG("Found separating weight vector: "
85 << storm::utility::vector::toString(storm::utility::vector::convertNumericVector<double>(halfspaces[farestHalfspaceIndex].normalVector()))
86 << ".");
87 return halfspaces[farestHalfspaceIndex].normalVector();
88}
89
90template<class SparseModelType, typename GeometryValueType>
92 // Normalize the direction vector so that the entries sum up to one
94 direction, storm::utility::one<GeometryValueType>() / std::accumulate(direction.begin(), direction.end(), storm::utility::zero<GeometryValueType>()));
95 weightVectorChecker->check(env, storm::utility::vector::convertNumericVector<typename SparseModelType::ValueType>(direction));
96 STORM_LOG_DEBUG("weighted objectives checker result (under approximation) is " << storm::utility::vector::toString(
97 storm::utility::vector::convertNumericVector<double>(weightVectorChecker->getUnderApproximationOfInitialStateResults())));
98 RefinementStep step;
99 step.weightVector = direction;
100 step.lowerBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getUnderApproximationOfInitialStateResults());
101 step.upperBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getOverApproximationOfInitialStateResults());
102 // For the minimizing objectives, we need to scale the corresponding entries with -1 as we want to consider the downward closure
103 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
104 if (storm::solver::minimize(this->objectives[objIndex].formula->getOptimalityType())) {
105 step.lowerBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
106 step.upperBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
107 }
108 }
109 if (produceScheduler) {
110 step.scheduler = weightVectorChecker->computeScheduler();
111 }
112 refinementSteps.push_back(std::move(step));
113
114 updateOverApproximation();
115 updateUnderApproximation();
116}
117
118template<class SparseModelType, typename GeometryValueType>
121 refinementSteps.back().weightVector, storm::utility::vector::dotProduct(refinementSteps.back().weightVector, refinementSteps.back().upperBoundPoint));
122
123 // Due to numerical issues, it might be the case that the updated overapproximation does not contain the underapproximation,
124 // e.g., when the new point is strictly contained in the underapproximation. Check if this is the case.
125 GeometryValueType maximumOffset = h.offset();
126 for (auto const& step : refinementSteps) {
127 maximumOffset = std::max(maximumOffset, storm::utility::vector::dotProduct(h.normalVector(), step.lowerBoundPoint));
128 }
129 if (maximumOffset > h.offset()) {
130 // We correct the issue by shifting the halfspace such that it contains the underapproximation
131 h.offset() = maximumOffset;
132 STORM_LOG_WARN("Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by "
133 << storm::utility::convertNumber<double>(h.invert().euclideanDistance(refinementSteps.back().upperBoundPoint)) << ".");
134 }
135 overApproximation = overApproximation->intersection(h);
136 STORM_LOG_DEBUG("Updated OverApproximation to " << overApproximation->toString(true));
137}
138
139template<class SparseModelType, typename GeometryValueType>
141 std::vector<Point> paretoPoints;
142 paretoPoints.reserve(refinementSteps.size());
143 for (auto const& step : refinementSteps) {
144 paretoPoints.push_back(step.lowerBoundPoint);
145 }
147 STORM_LOG_DEBUG("Updated UnderApproximation to " << underApproximation->toString(true));
148}
149
150template<class SparseModelType, typename GeometryValueType>
152 return env.modelchecker().multi().isMaxStepsSet() && this->refinementSteps.size() >= env.modelchecker().multi().getMaxSteps();
153}
154
155template<typename SparseModelType, typename GeometryValueType>
157 STORM_LOG_ERROR_COND(objectives.size() == 2, "Exporting plot requested but this is only implemented for the two-dimensional case.");
158
159 auto transformedUnderApprox = transformObjectivePolytopeToOriginal(this->objectives, underApproximation);
160 auto transformedOverApprox = transformObjectivePolytopeToOriginal(this->objectives, overApproximation);
161
162 // Get pareto points as well as a hyperrectangle that is used to guarantee that the resulting polytopes are bounded.
164 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()),
165 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()));
166 std::vector<std::vector<GeometryValueType>> paretoPoints;
167 paretoPoints.reserve(refinementSteps.size());
168 for (auto const& step : refinementSteps) {
169 paretoPoints.push_back(transformObjectiveValuesToOriginal(this->objectives, step.lowerBoundPoint));
170 boundaries.enlarge(paretoPoints.back());
171 }
172 auto underApproxVertices = transformedUnderApprox->getVertices();
173 for (auto const& v : underApproxVertices) {
174 boundaries.enlarge(v);
175 }
176 auto overApproxVertices = transformedOverApprox->getVertices();
177 for (auto const& v : overApproxVertices) {
178 boundaries.enlarge(v);
179 }
180
181 // Further enlarge the boundaries a little
182 storm::utility::vector::scaleVectorInPlace(boundaries.lowerBounds(), GeometryValueType(15) / GeometryValueType(10));
183 storm::utility::vector::scaleVectorInPlace(boundaries.upperBounds(), GeometryValueType(15) / GeometryValueType(10));
184
185 auto boundariesAsPolytope = boundaries.asPolytope();
186 std::vector<std::string> columnHeaders = {"x", "y"};
187
188 std::vector<std::vector<double>> pointsForPlotting;
190 underApproxVertices = transformedUnderApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
191 pointsForPlotting.reserve(underApproxVertices.size());
192 for (auto const& v : underApproxVertices) {
193 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
194 }
195 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathUnderApproximation().get(), pointsForPlotting, columnHeaders);
196 }
197
199 pointsForPlotting.clear();
200 overApproxVertices = transformedOverApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
201 pointsForPlotting.reserve(overApproxVertices.size());
202 for (auto const& v : overApproxVertices) {
203 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
204 }
205 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathOverApproximation().get(), pointsForPlotting, columnHeaders);
206 }
207
209 pointsForPlotting.clear();
210 pointsForPlotting.reserve(paretoPoints.size());
211 for (auto const& v : paretoPoints) {
212 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
213 }
214 storm::io::exportDataToCSVFile<double, std::string>(env.modelchecker().multi().getPlotPathParetoPoints().get(), pointsForPlotting, columnHeaders);
215 }
216}
217
218#ifdef STORM_HAVE_CARL
219template class SparsePcaaQuery<storm::models::sparse::Mdp<double>, storm::RationalNumber>;
220template class SparsePcaaQuery<storm::models::sparse::MarkovAutomaton<double>, storm::RationalNumber>;
221
222template class SparsePcaaQuery<storm::models::sparse::Mdp<storm::RationalNumber>, storm::RationalNumber>;
224#endif
225} // namespace multiobjective
226} // namespace modelchecker
227} // namespace storm
ModelCheckerEnvironment & modelchecker()
MultiObjectiveModelCheckerEnvironment & multi()
SparsePcaaQuery(preprocessing::SparseMultiObjectivePreprocessorResult< SparseModelType > &preprocessorResult)
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > overApproximation
WeightVector findSeparatingVector(Point const &pointToBeSeparated)
std::unique_ptr< PcaaWeightVectorChecker< SparseModelType > > weightVectorChecker
bool maxStepsPerformed(Environment const &env) const
void performRefinementStep(Environment const &env, WeightVector &&direction, bool produceScheduler)
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > underApproximation
void exportPlotOfCurrentApproximation(Environment const &env) const
std::vector< Objective< typename SparseModelType::ValueType > > objectives
static std::unique_ptr< PcaaWeightVectorChecker< ModelType > > create(preprocessing::SparseMultiObjectivePreprocessorResult< ModelType > const &preprocessorResult)
A bit vector that is internally represented as a vector of 64-bit values.
Definition BitVector.h:16
uint64_t getNextSetIndex(uint64_t startingIndex) const
Retrieves the index of the bit that is the next bit set to true in the bit vector.
uint64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
Halfspace< ValueType > invert() const
Definition Halfspace.h:67
std::vector< ValueType > const & normalVector() const
Definition Halfspace.h:112
ValueType const & offset() const
Definition Halfspace.h:120
void enlarge(std::vector< ValueType > const &point)
std::vector< ValueType > const & lowerBounds() const
std::shared_ptr< Polytope< ValueType > > asPolytope() const
std::vector< ValueType > const & upperBounds() const
static std::shared_ptr< Polytope< ValueType > > createDownwardClosure(std::vector< Point > const &points)
Creates the downward closure of the given points (i.e., the set { x | ex.
Definition Polytope.cpp:43
static std::shared_ptr< Polytope< ValueType > > createUniversalPolytope()
Creates the universal polytope (i.e., the set R^n)
Definition Polytope.cpp:27
static std::shared_ptr< Polytope< ValueType > > createEmptyPolytope()
Creates the empty polytope (i.e., emptyset)
Definition Polytope.cpp:32
#define STORM_LOG_WARN(message)
Definition logging.h:25
#define STORM_LOG_DEBUG(message)
Definition logging.h:18
#define STORM_LOG_ASSERT(cond, message)
Definition macros.h:11
#define STORM_LOG_ERROR_COND(cond, message)
Definition macros.h:52
#define STORM_LOG_THROW(cond, exception, message)
Definition macros.h:30
#define STORM_PRINT_AND_LOG(message)
Definition macros.h:68
std::vector< GeometryValueType > transformObjectiveValuesToOriginal(std::vector< Objective< ValueType > > objectives, std::vector< GeometryValueType > const &point)
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > transformObjectivePolytopeToOriginal(std::vector< Objective< ValueType > > objectives, std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > const &polytope)
bool constexpr minimize(OptimizationDirection d)
T dotProduct(std::vector< T > const &firstOperand, std::vector< T > const &secondOperand)
Computes the dot product (aka scalar product) and returns the result.
Definition vector.h:477
std::string toString(std::vector< ValueType > const &vector)
Output vector as string.
Definition vector.h:1149
void scaleVectorInPlace(std::vector< ValueType1 > &target, ValueType2 const &factor)
Multiplies each element of the given vector with the given factor and writes the result into the vect...
Definition vector.h:451
bool isZero(ValueType const &a)
Definition constants.cpp:41
LabParser.cpp.
std::optional< storm::storage::Scheduler< typename SparseModelType::ValueType > > scheduler