convex reduction

src/storm-pomdp/generator/NondeterministicBeliefTracker.cpp

@@ -2,6 +2,7 @@
 #include "storm-pomdp/generator/NondeterministicBeliefTracker.h"
 #include "storm/utility/ConstantsComparator.h"
 #include "storm/storage/geometry/nativepolytopeconversion/QuickHull.h"
+#include "storm/storage/geometry/ReduceVertexCloud.h"
 #include "storm/utility/vector.h"
 
 namespace storm {
@@ -198,6 +199,12 @@ namespace storm {
             return manager->getNumberOfStates();
         }
 
+        template<typename ValueType>
+        std::map<uint64_t, ValueType> const& SparseBeliefState<ValueType>::getBeliefMap() const {
+            return belief;
+        }
+
+
         template<typename ValueType>
         void SparseBeliefState<ValueType>::setSupport(storm::storage::BitVector& support) const {
             assert(false);
@@ -486,45 +493,33 @@ namespace storm {
         }
 
 //
-//        template<typename ValueType, typename BeliefState>
-//        void NondeterministicBeliefTracker<ValueType, BeliefState>::reduce() {
-//            std::cout << "reduce" << std::endl;
-//            storm::storage::BitVector support(beliefs.begin()->getSupportSize());
-//            std::cout << "supp generated" << std::endl;
-//            for(auto const& belief : beliefs) {
-//                belief.setSupport(support);
-//            }
-//            std::cout << "Support:" << support << std::endl;
-//
-//            if (beliefs.size() <= support.getNumberOfSetBits()) {
-//                return;
-//            }
-//
-//            std::vector<typename storm::storage::geometry::QuickHull<ValueType>::EigenVector> points;
-//            for(auto const& bel : beliefs) {
-//                std::cout << bel.toEigenVector(support) << std::endl;
-//                points.push_back(bel.toEigenVector(support));
-//            }
-//
-//            storm::storage::geometry::QuickHull<ValueType> qh;
-//            std::cout << "Run QH (dim=" << points[0].size() << ")" << std::endl;
-//            qh.generateHalfspacesFromPoints(points, true);
-//            std::cout << "End QH" << std::endl;
-//            auto filteredPoints = qh.getRelevantVertices();
-//            //if (filteredPoints.size() < points.size()) {
-//            //    for(auto const& fp : filteredPoints) {
-//            //        beliefs.emplace(manager,filteredPoints);
-//            //    }
-//            //}
-//        }
+        template<typename ValueType, typename BeliefState>
+        uint64_t NondeterministicBeliefTracker<ValueType, BeliefState>::reduce() {
+            std::shared_ptr<storm::utility::solver::SmtSolverFactory> solverFactory = std::make_shared<storm::utility::solver::Z3SmtSolverFactory>();
+            storm::storage::geometry::ReduceVertexCloud<ValueType> rvc(solverFactory);
+            std::vector<std::map<uint64_t, ValueType>> points;
+            std::vector<typename std::unordered_set<BeliefState>::iterator> iterators;
+            for (auto it = beliefs.begin(); it != beliefs.end(); ++it) {
+                // TODO get rid of the getBeliefMap function.
+                points.push_back(it->getBeliefMap());
+                iterators.push_back(it);
+            }
+            storm::storage::BitVector eliminate = ~rvc.eliminate(points, pomdp.getNumberOfStates());
+
+            auto selectedIterators = storm::utility::vector::filterVector(iterators, eliminate);
+            for (auto iter : selectedIterators) {
+                beliefs.erase(iter);
+            }
+            return eliminate.getNumberOfSetBits();
+        }
 
 
         template class SparseBeliefState<double>;
         template bool operator==(SparseBeliefState<double> const&, SparseBeliefState<double> const&);
         template class NondeterministicBeliefTracker<double, SparseBeliefState<double>>;
-        template class ObservationDenseBeliefState<double>;
-        template bool operator==(ObservationDenseBeliefState<double> const&, ObservationDenseBeliefState<double> const&);
-        template class NondeterministicBeliefTracker<double, ObservationDenseBeliefState<double>>;
+        //template class ObservationDenseBeliefState<double>;
+        //template bool operator==(ObservationDenseBeliefState<double> const&, ObservationDenseBeliefState<double> const&);
+        //template class NondeterministicBeliefTracker<double, ObservationDenseBeliefState<double>>;
 
 
     }

src/storm-pomdp/generator/NondeterministicBeliefTracker.h

@@ -46,6 +46,7 @@ namespace storm {
             Eigen::Matrix<ValueType, Eigen::Dynamic, 1> toEigenVector(storm::storage::BitVector const& support) const;
             uint64_t getSupportSize() const;
             void setSupport(storm::storage::BitVector&) const;
+            std::map<uint64_t, ValueType> const& getBeliefMap() const;
 
             friend bool operator==<>(SparseBeliefState<ValueType> const& lhs, SparseBeliefState<ValueType> const& rhs);
         private:
@@ -102,7 +103,7 @@ namespace storm {
             uint32_t getCurrentObservation() const;
             ValueType getCurrentRisk(bool max=true);
             void setRisk(std::vector<ValueType> const& risk);
-            //void reduce();
+            uint64_t reduce();
 
         private:
 

src/storm-pomdp/transformer/ObservationTraceUnfolder.cpp

@@ -59,15 +59,15 @@ namespace storm {
 
                 for (auto const& unfoldedToOldEntry : unfoldedToOld) {
                     transitionMatrixBuilder.newRowGroup(newRowGroupStart);
-                    std::cout << "\tconsider new state " << unfoldedToOldEntry.first << std::endl;
+                    //std::cout << "\tconsider new state " << unfoldedToOldEntry.first << std::endl;
                     assert(step == 0 || newRowCount == transitionMatrixBuilder.getLastRow() + 1);
                     svbuilder.addState(unfoldedToOldEntry.first, {}, {static_cast<int64_t>(unfoldedToOldEntry.second)});
                     uint64_t oldRowIndexStart = model.getNondeterministicChoiceIndices()[unfoldedToOldEntry.second];
                     uint64_t oldRowIndexEnd = model.getNondeterministicChoiceIndices()[unfoldedToOldEntry.second+1];
 
                     for (uint64_t oldRowIndex = oldRowIndexStart; oldRowIndex != oldRowIndexEnd; oldRowIndex++) {
-                        std::cout << "\t\tconsider old action " << oldRowIndex << std::endl;
-                        std::cout << "\t\tconsider new row nr " << newRowCount << std::endl;
+                        //std::cout << "\t\tconsider old action " << oldRowIndex << std::endl;
+                        //std::cout << "\t\tconsider new row nr " << newRowCount << std::endl;
 
                         ValueType resetProb = storm::utility::zero<ValueType>();
                         // We first find the reset probability
@@ -76,14 +76,14 @@ namespace storm {
                                 resetProb += oldRowEntry.getValue();
                             }
                         }
-                        std::cout << "\t\t\t add reset" << std::endl;
+                        //std::cout << "\t\t\t add reset" << std::endl;
 
                         // Add the resets
                         if (resetProb != storm::utility::zero<ValueType>()) {
                             transitionMatrixBuilder.addNextValue(newRowCount, 0, resetProb);
                         }
 
-                        std::cout << "\t\t\t add other transitions..." << std::endl;
+                        //std::cout << "\t\t\t add other transitions..." << std::endl;
 
                         // Now, we build the outgoing transitions.
                         for (auto const &oldRowEntry : model.getTransitionMatrix().getRow(oldRowIndex)) {
@@ -101,7 +101,7 @@ namespace storm {
                             } else {
                                 column = entryIt->second;
                             }
-                            std::cout << "\t\t\t\t transition to " << column << std::endl;
+                            //std::cout << "\t\t\t\t transition to " << column << std::endl;
                             transitionMatrixBuilder.addNextValue(newRowCount, column,
                                                                    oldRowEntry.getValue());
                         }
@@ -145,7 +145,6 @@ namespace storm {
 
             storm::storage::sparse::ModelComponents<ValueType> components;
             components.transitionMatrix = transitionMatrixBuilder.build();
-            std::cout << components.transitionMatrix << std::endl;
 
             STORM_LOG_ASSERT(components.transitionMatrix.getRowGroupCount() == targetState + 1, "Expect row group count (" << components.transitionMatrix.getRowGroupCount() << ") one more as target state index " << targetState << ")");
 

src/storm/storage/geometry/ReduceVertexCloud.cpp

@@ -0,0 +1,119 @@
+#include "storm/storage/geometry/ReduceVertexCloud.h"
+#include "storm/utility/Stopwatch.h"
+#undef _DEBUG_REDUCE_VERTEX_CLOUD
+
+
+namespace storm {
+    namespace storage {
+        namespace geometry  {
+
+            template<typename ValueType>
+            std::string toString(std::map<uint64_t, ValueType> const& point) {
+                std::stringstream  sstr;
+                bool first = true;
+                for (auto const& entry : point) {
+                    if (first) {
+                        first = false;
+                    } else {
+                        sstr << ", ";
+                    }
+                    sstr << entry.first << " : " << entry.second;
+                }
+                return sstr.str();
+            }
+
+            template<typename ValueType>
+            storm::storage::BitVector ReduceVertexCloud<ValueType>::eliminate(std::vector<std::map<uint64_t, ValueType>> const& input, uint64_t maxdimension) {
+                std::shared_ptr<storm::expressions::ExpressionManager> expressionManager = std::make_shared<storm::expressions::ExpressionManager>();
+                std::vector<storm::storage::BitVector> supports;
+                std::vector<storm::expressions::Variable> weightVariables;
+                std::vector<storm::expressions::Expression> weightVariableExpressions;
+
+                for (uint64_t pointIndex = 0; pointIndex < input.size(); ++pointIndex) {
+
+                    // Compute the support vectors to quickly determine which input points could be relevant.
+                    supports.push_back(storm::storage::BitVector(maxdimension));
+                    for (auto const& entry : input[pointIndex]) {
+                        supports.back().set(entry.first, true);
+                    }
+                    // Add a weight variable for each input point
+                    weightVariables.push_back(expressionManager->declareRationalVariable("w_"+ std::to_string(pointIndex)));
+                    // For convenience and performance, obtain the expression.
+                    weightVariableExpressions.push_back(weightVariables.back().getExpression());
+                }
+
+                std::unique_ptr<storm::solver::SmtSolver> smtSolver = smtSolverFactory->create(*expressionManager);
+                for (auto const& weightVariableExpr : weightVariableExpressions) {
+                    //smtSolver->add((weightVariableExpr == expressionManager->rational(0.0)) || (weightVariableExpr > expressionManager->rational(0.00001)));
+                    smtSolver->add((weightVariableExpr >= expressionManager->rational(0.0)));
+                    smtSolver->add(weightVariableExpr < expressionManager->rational(1.0));
+                }
+                smtSolver->add(storm::expressions::sum(weightVariableExpressions) <= expressionManager->rational(1.0 + wiggle));
+                smtSolver->add(storm::expressions::sum(weightVariableExpressions) >= expressionManager->rational(1 - wiggle));
+
+                storm::utility::Stopwatch solverTime;
+                storm::utility::Stopwatch totalTime(true);
+                storm::storage::BitVector vertices(input.size());
+                for (uint64_t pointIndex = 0; pointIndex < input.size(); ++pointIndex) {
+                    std::cout << pointIndex << " out of " << input.size() << std::endl;
+                    smtSolver->push();
+                    std::map<uint64_t, std::vector<storm::expressions::Expression>> dimensionTerms;
+                    for (auto const& entry : input[pointIndex]) {
+                        dimensionTerms[entry.first] = {expressionManager->rational(-entry.second)};
+                    }
+                    for (uint64_t potentialSupport = 0; potentialSupport < input.size(); ++potentialSupport) {
+                        if (pointIndex == potentialSupport) {
+                            smtSolver->add(weightVariableExpressions[potentialSupport] == expressionManager->rational(0.0));
+                        } else if (potentialSupport < pointIndex && !vertices.get(potentialSupport)) {
+                            smtSolver->add(weightVariableExpressions[potentialSupport] == expressionManager->rational(0.0));
+                        } else if (supports[potentialSupport].isSubsetOf(supports[pointIndex])) {
+                            for (auto const& entry : input[potentialSupport]) {
+                                dimensionTerms[entry.first].push_back(weightVariableExpressions[potentialSupport] * expressionManager->rational(entry.second));
+                            }
+                        } else {
+                            smtSolver->add(weightVariableExpressions[potentialSupport] == expressionManager->rational(0.0));
+                        }
+                    }
+                    for (auto const& entry : dimensionTerms) {
+                        smtSolver->add(storm::expressions::sum(entry.second) == expressionManager->rational(0.0));
+                    }
+
+                    solverTime.start();
+                    auto result = smtSolver->check();
+                    solverTime.stop();
+                    if (result == storm::solver::SmtSolver::CheckResult::Unsat) {
+#ifdef _DEBUG_REDUCE_VERTEX_CLOUD
+                        if (input[pointIndex].size() == 2) {
+                            std::cout << "point " << toString(input[pointIndex]) << " is a vertex:";
+                            std::cout << smtSolver->getSmtLibString() << std::endl;
+                        }
+#endif
+                        vertices.set(pointIndex, true);
+                    }
+#ifdef _DEBUG_REDUCE_VERTEX_CLOUD
+                    else
+                    {
+                        std::cout << "point " << toString(input[pointIndex]) << " is a convex combination of ";
+                        auto val = smtSolver->getModelAsValuation();
+                        uint64_t varIndex = 0;
+                        for (auto const& wvar : weightVariables) {
+                            if (!storm::utility::isZero(val.getRationalValue(wvar))) {
+                                std::cout << toString(input[varIndex]) << " (weight: " <<  val.getRationalValue(wvar) << ")";
+                            }
+                            varIndex++;
+                        }
+                        std::cout << std::endl;
+                    }
+#endif
+                    smtSolver->pop();
+                    std::cout << "Solver time " << solverTime.getTimeInMilliseconds() << std::endl;
+                    std::cout << "Total time " << totalTime.getTimeInMilliseconds() << std::endl;
+                }
+                return vertices;
+
+            }
+
+            template class ReduceVertexCloud<double>;
+        }
+    }
+}

src/storm/storage/geometry/ReduceVertexCloud.h

@@ -0,0 +1,29 @@
+#pragma once
+
+#include "storm/storage/BitVector.h"
+#include "storm/solver/SmtSolver.h"
+#include "storm/utility/solver.h"
+
+namespace storm {
+    namespace storage {
+        namespace geometry {
+            template<typename ValueType>
+            class ReduceVertexCloud {
+            public:
+                ReduceVertexCloud(std::shared_ptr<storm::utility::solver::SmtSolverFactory>& smtSolverFactory, ValueType wiggle = storm::utility::convertNumber<ValueType>(0.001))
+                : smtSolverFactory(smtSolverFactory), wiggle(wiggle)
+                {
+
+                }
+
+                storm::storage::BitVector eliminate(std::vector<std::map<uint64_t, ValueType>> const& input, uint64_t maxdimension);
+
+            private:
+                std::shared_ptr<storm::utility::solver::SmtSolverFactory>& smtSolverFactory;
+                ValueType wiggle;
+
+            };
+
+        }
+    }
+}