20namespace modelchecker {
21namespace multiobjective {
23template<
class SparseModelType,
typename GeometryValueType>
26 STORM_PRINT_AND_LOG(
"Pareto Curve Approximation Algorithm terminated after " << this->refinementSteps.size() <<
" refinement steps.\n");
30template<
class SparseModelType,
typename GeometryValueType>
32 : originalModel(preprocessorResult.originalModel), originalFormula(preprocessorResult.originalFormula), objectives(preprocessorResult.objectives) {
40template<
class SparseModelType,
typename GeometryValueType>
42 Point const& pointToBeSeparated) {
43 STORM_LOG_DEBUG(
"Searching a weight vector to seperate the point given by "
46 if (underApproximation->isEmpty()) {
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);
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);
68 bool isSingleObjectiveVector =
71 if ((!foundSeparatingDiracVector && isSingleObjectiveVector) ||
72 (foundSeparatingDiracVector == isSingleObjectiveVector && distance > farestDistance)) {
73 foundSeparatingDiracVector = foundSeparatingDiracVector || isSingleObjectiveVector;
74 farestHalfspaceIndex = halfspaceIndex;
75 farestDistance = distance;
79 if (foundSeparatingDiracVector) {
80 diracWeightVectorsToBeChecked &= storm::utility::vector::filterZero<GeometryValueType>(halfspaces[farestHalfspaceIndex].normalVector());
83 STORM_LOG_THROW(farestHalfspaceIndex < halfspaces.size(), storm::exceptions::UnexpectedException,
"There is no seperating vector.");
87 return halfspaces[farestHalfspaceIndex].normalVector();
90template<
class SparseModelType,
typename GeometryValueType>
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));
97 storm::utility::vector::convertNumericVector<double>(weightVectorChecker->getUnderApproximationOfInitialStateResults())));
100 step.
lowerBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getUnderApproximationOfInitialStateResults());
101 step.
upperBoundPoint = storm::utility::vector::convertNumericVector<GeometryValueType>(weightVectorChecker->getOverApproximationOfInitialStateResults());
103 for (uint_fast64_t objIndex = 0; objIndex < this->objectives.size(); ++objIndex) {
105 step.
lowerBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
106 step.
upperBoundPoint[objIndex] *= -storm::utility::one<GeometryValueType>();
109 refinementSteps.push_back(std::move(step));
111 updateOverApproximation();
112 updateUnderApproximation();
115template<
class SparseModelType,
typename GeometryValueType>
122 GeometryValueType maximumOffset = h.
offset();
123 for (
auto const& step : refinementSteps) {
126 if (maximumOffset > h.
offset()) {
128 h.
offset() = maximumOffset;
129 STORM_LOG_WARN(
"Numerical issues: The overapproximation would not contain the underapproximation. Hence, a halfspace is shifted by "
130 << storm::utility::convertNumber<double>(h.
invert().euclideanDistance(refinementSteps.back().upperBoundPoint)) <<
".");
132 overApproximation = overApproximation->intersection(h);
133 STORM_LOG_DEBUG(
"Updated OverApproximation to " << overApproximation->toString(
true));
136template<
class SparseModelType,
typename GeometryValueType>
138 std::vector<Point> paretoPoints;
139 paretoPoints.reserve(refinementSteps.size());
140 for (
auto const& step : refinementSteps) {
141 paretoPoints.push_back(step.lowerBoundPoint);
144 STORM_LOG_DEBUG(
"Updated UnderApproximation to " << underApproximation->toString(
true));
147template<
class SparseModelType,
typename GeometryValueType>
152template<
typename SparseModelType,
typename GeometryValueType>
154 STORM_LOG_ERROR_COND(objectives.size() == 2,
"Exporting plot requested but this is only implemented for the two-dimensional case.");
161 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()),
162 std::vector<GeometryValueType>(objectives.size(), storm::utility::zero<GeometryValueType>()));
163 std::vector<std::vector<GeometryValueType>> paretoPoints;
164 paretoPoints.reserve(refinementSteps.size());
165 for (
auto const& step : refinementSteps) {
167 boundaries.
enlarge(paretoPoints.back());
169 auto underApproxVertices = transformedUnderApprox->getVertices();
170 for (
auto const& v : underApproxVertices) {
173 auto overApproxVertices = transformedOverApprox->getVertices();
174 for (
auto const& v : overApproxVertices) {
182 auto boundariesAsPolytope = boundaries.
asPolytope();
183 std::vector<std::string> columnHeaders = {
"x",
"y"};
185 std::vector<std::vector<double>> pointsForPlotting;
187 underApproxVertices = transformedUnderApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
188 pointsForPlotting.reserve(underApproxVertices.size());
189 for (
auto const& v : underApproxVertices) {
190 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
196 pointsForPlotting.clear();
197 overApproxVertices = transformedOverApprox->intersection(boundariesAsPolytope)->getVerticesInClockwiseOrder();
198 pointsForPlotting.reserve(overApproxVertices.size());
199 for (
auto const& v : overApproxVertices) {
200 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
206 pointsForPlotting.clear();
207 pointsForPlotting.reserve(paretoPoints.size());
208 for (
auto const& v : paretoPoints) {
209 pointsForPlotting.push_back(storm::utility::vector::convertNumericVector<double>(v));
215#ifdef STORM_HAVE_CARL
ModelCheckerEnvironment & modelchecker()
MultiObjectiveModelCheckerEnvironment & multi()
uint64_t const & getMaxSteps() const
boost::optional< std::string > getPlotPathOverApproximation() const
bool isMaxStepsSet() const
boost::optional< std::string > getPlotPathParetoPoints() const
boost::optional< std::string > getPlotPathUnderApproximation() const
void performRefinementStep(Environment const &env, WeightVector &&direction)
SparsePcaaQuery(preprocessing::SparseMultiObjectivePreprocessorResult< SparseModelType > &preprocessorResult)
void updateUnderApproximation()
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
virtual ~SparsePcaaQuery()
storm::storage::BitVector diracWeightVectorsToBeChecked
std::vector< GeometryValueType > WeightVector
std::shared_ptr< storm::storage::geometry::Polytope< GeometryValueType > > underApproximation
void updateOverApproximation()
void exportPlotOfCurrentApproximation(Environment const &env) const
std::vector< Objective< typename SparseModelType::ValueType > > objectives
std::vector< GeometryValueType > Point
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.
uint_fast64_t getNextSetIndex(uint_fast64_t startingIndex) const
Retrieves the index of the bit that is the next bit set to true in the bit vector.
uint_fast64_t getNumberOfSetBits() const
Returns the number of bits that are set to true in this bit vector.
Halfspace< ValueType > invert() const
std::vector< ValueType > const & normalVector() const
ValueType const & offset() const
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.
static std::shared_ptr< Polytope< ValueType > > createUniversalPolytope()
Creates the universal polytope (i.e., the set R^n)
static std::shared_ptr< Polytope< ValueType > > createEmptyPolytope()
Creates the empty polytope (i.e., emptyset)
#define STORM_LOG_WARN(message)
#define STORM_LOG_DEBUG(message)
#define STORM_LOG_ASSERT(cond, message)
#define STORM_LOG_ERROR_COND(cond, message)
#define STORM_LOG_THROW(cond, exception, message)
#define STORM_PRINT_AND_LOG(message)
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)
SettingsType const & getModule()
Get module.
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.
std::string toString(std::vector< ValueType > const &vector)
Output vector as string.
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...
bool isZero(ValueType const &a)
WeightVector weightVector