fixed issue in rational search preventing convergence in many cases

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -776,9 +776,9 @@ namespace storm {
         template<typename ValueType>
         template<typename RationalType, typename ImpreciseType>
         bool NativeLinearEquationSolver<ValueType>::sharpen(uint64_t precision, storm::storage::SparseMatrix<RationalType> const& A, std::vector<ImpreciseType> const& x, std::vector<RationalType> const& b, std::vector<RationalType>& tmp) {
-            for (uint64_t p = 0; p <= precision; ++p) {
+            for (uint64_t p = 1; p <= precision; ++p) {
                 storm::utility::kwek_mehlhorn::sharpen(p, x, tmp);
-                
+
                 if (NativeLinearEquationSolver<RationalType>::isSolution(A, tmp, b)) {
                     return true;
                 }

src/storm/solver/SymbolicMinMaxLinearEquationSolver.cpp

@@ -142,9 +142,10 @@ namespace storm {
         template<typename RationalType, typename ImpreciseType>
         storm::dd::Add<DdType, RationalType> SymbolicMinMaxLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, OptimizationDirection dir, SymbolicMinMaxLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
             
-            // Storage for the rational sharpened vector.
+            // Storage for the rational sharpened vector and the power iteration intermediate vector.
             storm::dd::Add<DdType, RationalType> sharpenedX;
-            
+            storm::dd::Add<DdType, ImpreciseType> currentX = x;
+
             // The actual rational search.
             uint64_t overallIterations = 0;
             uint64_t valueIterationInvocations = 0;
@@ -153,7 +154,7 @@ namespace storm {
             bool relative = env.solver().minMax().getRelativeTerminationCriterion();
             SolverStatus status = SolverStatus::InProgress;
             while (status == SolverStatus::InProgress && overallIterations < maxIter) {
-                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
+                typename SymbolicMinMaxLinearEquationSolver<DdType, ImpreciseType>::ValueIterationResult viResult = impreciseSolver.performValueIteration(dir, currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter);
                 
                 ++valueIterationInvocations;
                 STORM_LOG_TRACE("Completed " << valueIterationInvocations << " value iteration invocations, the last one with precision " << precision << " completed in " << viResult.iterations << " iterations.");
@@ -170,6 +171,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
+                    currentX = viResult.values;
                     precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
                 }
             }

src/storm/solver/SymbolicNativeLinearEquationSolver.cpp

@@ -174,10 +174,10 @@ namespace storm {
             
             storm::dd::Add<DdType, RationalType> sharpenedX;
             
-            for (uint64_t p = 1; p < precision; ++p) {
+            for (uint64_t p = 1; p <= precision; ++p) {
                 sharpenedX = x.sharpenKwekMehlhorn(p);
+
                 isSolution = rationalSolver.isSolutionFixedPoint(sharpenedX, rationalB);
-                
                 if (isSolution) {
                     break;
                 }
@@ -190,8 +190,9 @@ namespace storm {
         template<typename RationalType, typename ImpreciseType>
         storm::dd::Add<DdType, RationalType> SymbolicNativeLinearEquationSolver<DdType, ValueType>::solveEquationsRationalSearchHelper(Environment const& env, SymbolicNativeLinearEquationSolver<DdType, RationalType> const& rationalSolver, SymbolicNativeLinearEquationSolver<DdType, ImpreciseType> const& impreciseSolver, storm::dd::Add<DdType, RationalType> const& rationalB, storm::dd::Add<DdType, ImpreciseType> const& x, storm::dd::Add<DdType, ImpreciseType> const& b) const {
             
-            // Storage for the rational sharpened vector.
+            // Storage for the rational sharpened vector and the power iteration intermediate vector.
             storm::dd::Add<DdType, RationalType> sharpenedX;
+            storm::dd::Add<DdType, ImpreciseType> currentX = x;
             
             // The actual rational search.
             uint64_t overallIterations = 0;
@@ -201,7 +202,7 @@ namespace storm {
             bool relative = env.solver().native().getRelativeTerminationCriterion();
             SolverStatus status = SolverStatus::InProgress;
             while (status == SolverStatus::InProgress && overallIterations < maxIter) {
-                typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(x, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
+                typename SymbolicNativeLinearEquationSolver<DdType, ImpreciseType>::PowerIterationResult result = impreciseSolver.performPowerIteration(currentX, b, storm::utility::convertNumber<ImpreciseType, ValueType>(precision), relative, maxIter - overallIterations);
                 
                 ++powerIterationInvocations;
                 STORM_LOG_TRACE("Completed " << powerIterationInvocations << " power iteration invocations, the last one with precision " << precision << " completed in " << result.iterations << " iterations.");
@@ -218,6 +219,7 @@ namespace storm {
                 if (isSolution) {
                     status = SolverStatus::Converged;
                 } else {
+                    currentX = result.values;
                     precision /= storm::utility::convertNumber<ValueType, uint64_t>(10);
                 }
             }

src/storm/storage/dd/Add.cpp

@@ -941,7 +941,7 @@ namespace storm {
         
         template<DdType LibraryType, typename ValueType>
         std::ostream& operator<<(std::ostream& out, Add<LibraryType, ValueType> const& add) {
-            out << "ADD with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
+            out << "ADD [" << add.getInternalAdd().getStringId() << "] with " << add.getNonZeroCount() << " nnz, " << add.getNodeCount() << " nodes, " << add.getLeafCount() << " leaves" << std::endl;
             std::vector<std::string> variableNames;
             for (auto const& variable : add.getContainedMetaVariables()) {
                 variableNames.push_back(variable.getName());
@@ -993,15 +993,19 @@ namespace storm {
         }
 		
         template class Add<storm::dd::DdType::CUDD, double>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, double> const& add);
         template class Add<storm::dd::DdType::CUDD, uint_fast64_t>;
 
         template class Add<storm::dd::DdType::Sylvan, double>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, double> const& add);
         template class Add<storm::dd::DdType::Sylvan, uint_fast64_t>;
 
 #ifdef STORM_HAVE_CARL
         template class Add<storm::dd::DdType::CUDD, storm::RationalNumber>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::CUDD, storm::RationalNumber> const& add);
 
         template class Add<storm::dd::DdType::Sylvan, storm::RationalNumber>;
+        template std::ostream& operator<<(std::ostream& out, Add<storm::dd::DdType::Sylvan, storm::RationalNumber> const& add);
 		template class Add<storm::dd::DdType::Sylvan, storm::RationalFunction>;
 
         template Add<storm::dd::DdType::CUDD, storm::RationalNumber> Add<storm::dd::DdType::CUDD, storm::RationalNumber>::toValueType<storm::RationalNumber>() const;

src/storm/storage/dd/cudd/InternalCuddAdd.cpp

@@ -420,6 +420,13 @@ namespace storm {
         DdNode* InternalAdd<DdType::CUDD, ValueType>::getCuddDdNode() const {
             return this->getCuddAdd().getNode();
         }
+        
+        template<typename ValueType>
+        std::string InternalAdd<DdType::CUDD, ValueType>::getStringId() const {
+            std::stringstream ss;
+            ss << this->getCuddDdNode();
+            return ss.str();
+        }
 
         template<typename ValueType>
         Odd InternalAdd<DdType::CUDD, ValueType>::createOdd(std::vector<uint_fast64_t> const& ddVariableIndices) const {

src/storm/storage/dd/cudd/InternalCuddAdd.h

@@ -629,6 +629,11 @@ namespace storm {
              */
             DdNode* getCuddDdNode() const;
             
+            /*!
+             * Retrieves a string representation of an ID for thid ADD.
+             */
+            std::string getStringId() const;
+            
         private:
             /*!
              * Performs a recursive step to perform the given function between the given DD-based vector and the given

src/storm/storage/dd/sylvan/InternalSylvanAdd.cpp

@@ -1151,6 +1151,13 @@ namespace storm {
             return sylvanMtbdd;
         }
         
+        template<typename ValueType>
+        std::string InternalAdd<DdType::Sylvan, ValueType>::getStringId() const {
+            std::stringstream ss;
+            ss << this->getSylvanMtbdd().GetMTBDD();
+            return ss.str();
+        }
+        
         template class InternalAdd<DdType::Sylvan, double>;
         template class InternalAdd<DdType::Sylvan, uint_fast64_t>;
 

src/storm/storage/dd/sylvan/InternalSylvanAdd.h

@@ -620,6 +620,8 @@ namespace storm {
              */
             static ValueType getValue(MTBDD const& node);
             
+            std::string getStringId() const;
+            
         private:
             /*!
              * Recursively builds the ODD from an ADD.