initial implementation for quick and sound vi for DTMCs

src/storm/settings/modules/NativeEquationSolverSettings.cpp

@@ -25,7 +25,7 @@ namespace storm {
             const std::string NativeEquationSolverSettings::powerMethodMultiplicationStyleOptionName = "powmult";
 
             NativeEquationSolverSettings::NativeEquationSolverSettings() : ModuleSettings(moduleName) {
-                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae", "power", "ratsearch" };
+                std::vector<std::string> methods = { "jacobi", "gaussseidel", "sor", "walkerchae", "power", "ratsearch", "qpower" };
                 this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(methods)).setDefaultValueString("jacobi").build()).build());
                 
                 this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build());
@@ -63,6 +63,8 @@ namespace storm {
                     return storm::solver::NativeLinearEquationSolverMethod::Power;
                 } else if (linearEquationSystemTechniqueAsString == "ratsearch") {
                     return storm::solver::NativeLinearEquationSolverMethod::RationalSearch;
+                } else if (linearEquationSystemTechniqueAsString == "qpower") {
+                    return storm::solver::NativeLinearEquationSolverMethod::QuickPower;
                 }
                 STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected.");
             }

src/storm/solver/NativeLinearEquationSolver.cpp

@@ -555,6 +555,99 @@ namespace storm {
             return converged;
         }
         
+        
+        
+        template<typename ValueType>
+        bool NativeLinearEquationSolver<ValueType>::solveEquationsQuickSoundPower(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
+            STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (QuickPower)");
+            // Prepare the solution vectors.
+            if (!this->cachedRowVector) {
+                this->cachedRowVector = std::make_unique<std::vector<ValueType>>(getMatrixRowCount(), storm::utility::zero<ValueType>());
+            } else {
+                this->cachedRowVector->assign(getMatrixRowCount(), storm::utility::zero<ValueType>());
+            }
+            if (!this->cachedRowVector2) {
+                this->cachedRowVector2 = std::make_unique<std::vector<ValueType>>(getMatrixRowCount(), storm::utility::one<ValueType>());
+            } else {
+                this->cachedRowVector2->assign(getMatrixRowCount(), storm::utility::one<ValueType>());
+            }
+
+            std::vector<ValueType>* stepBoundedValues = this->cachedRowVector.get();
+            std::vector<ValueType>* stepBoundedStayProbabilities = this->cachedRowVector2.get();
+
+            std::vector<ValueType>* tmp = &x;
+            
+            bool converged = false;
+            bool terminate = false;
+            uint64_t iterations = 0;
+            bool doConvergenceCheck = true;
+            bool useDiffs = this->hasRelevantValues();
+            ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().native().getPrecision());
+            ValueType lowerValueBound, upperValueBound;
+            bool relative = env.solver().native().getRelativeTerminationCriterion();
+            if (!relative) {
+                precision *= storm::utility::convertNumber<ValueType>(2.0);
+            }
+            uint64_t maxIter = env.solver().native().getMaximalNumberOfIterations();
+            this->startMeasureProgress();
+            auto firstProb1EntryIt = stepBoundedStayProbabilities->begin();
+            while (!converged && !terminate && iterations < maxIter) {
+                this->multiplier.multAdd(*this->A, *stepBoundedValues, &b, *tmp);
+                std::swap(tmp, stepBoundedValues);
+                this->multiplier.multAdd(*this->A, *stepBoundedStayProbabilities, nullptr, *tmp);
+                std::swap(tmp, stepBoundedStayProbabilities);
+                for (; firstProb1EntryIt != stepBoundedStayProbabilities->end(); ++firstProb1EntryIt) {
+                    static_assert(NumberTraits<ValueType>::IsExact || std::is_same<ValueType, double>::value, "Considered ValueType not handled.");
+                    if (NumberTraits<ValueType>::IsExact) {
+                        if (storm::utility::isOne(*firstProb1EntryIt)) {
+                            break;
+                        }
+                    } else {
+                        if (storm::utility::isAlmostOne(storm::utility::convertNumber<double>(*firstProb1EntryIt))) {
+                            break;
+                        }
+                    }
+                }
+                if (firstProb1EntryIt == stepBoundedStayProbabilities->end()) {
+                    auto valIt = stepBoundedValues->begin();
+                    auto valIte = stepBoundedValues->end();
+                    auto probIt = stepBoundedStayProbabilities->begin();
+                    STORM_LOG_ASSERT(!storm::utility::isOne(*probIt), "Did not expect staying-probability 1 at this point.");
+                    lowerValueBound = *valIt / (storm::utility::one<ValueType>() - *probIt);
+                    upperValueBound = lowerValueBound;
+                    ValueType largestStayProb = *probIt;
+                    for (; valIt != valIte; ++valIt, ++probIt) {
+                        ValueType currentBound = *valIt / (storm::utility::one<ValueType>() - *probIt);
+                        lowerValueBound = std::min(lowerValueBound, currentBound);
+                        upperValueBound = std::max(upperValueBound, currentBound);
+                        largestStayProb = std::max(largestStayProb, *probIt);
+                    }
+                    STORM_LOG_ASSERT(!relative, "Relative termination criterion not implemented currently.");
+                    converged = largestStayProb * (upperValueBound - lowerValueBound) < precision;
+                }
+                
+                // Potentially show progress.
+                this->showProgressIterative(iterations);
+                
+                // Set up next iteration.
+                ++iterations;
+
+            }
+            
+            // Finally set up the solution vector
+            ValueType meanBound = (upperValueBound - lowerValueBound) / storm::utility::convertNumber<ValueType>(2.0);
+            storm::utility::vector::applyPointwise(*stepBoundedValues, *stepBoundedStayProbabilities, x, [&meanBound] (ValueType const& v, ValueType const& p) { return v + p * meanBound; });
+            
+            if (!this->isCachingEnabled()) {
+                clearCache();
+            }
+            
+            this->logIterations(converged, terminate, iterations);
+
+            return converged;
+
+        }
+        
         template<typename ValueType>
         bool NativeLinearEquationSolver<ValueType>::solveEquationsRationalSearch(Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const {
             return solveEquationsRationalSearchHelper<double>(env, x, b);
@@ -828,7 +921,7 @@ namespace storm {
                 } else {
                     STORM_LOG_WARN("The selected solution method does not guarantee exact results.");
                 }
-            } else if (env.solver().isForceSoundness() && method != NativeLinearEquationSolverMethod::Power && method != NativeLinearEquationSolverMethod::RationalSearch) {
+            } else if (env.solver().isForceSoundness() && method != NativeLinearEquationSolverMethod::Power && method != NativeLinearEquationSolverMethod::RationalSearch && method != NativeLinearEquationSolverMethod::QuickPower) {
                 if (env.solver().native().isMethodSetFromDefault()) {
                     method = NativeLinearEquationSolverMethod::Power;
                     STORM_LOG_INFO("Selecting '" + toString(method) + "' as the solution technique to guarantee sound results. If you want to override this, please explicitly specify a different method.");
@@ -857,6 +950,8 @@ namespace storm {
                     } else {
                         return this->solveEquationsPower(env, x, b);
                     }
+                case NativeLinearEquationSolverMethod::QuickPower:
+                        return this->solveEquationsQuickSoundPower(env, x, b);
                 case NativeLinearEquationSolverMethod::RationalSearch:
                     return this->solveEquationsRationalSearch(env, x, b);
             }
@@ -921,7 +1016,7 @@ namespace storm {
         template<typename ValueType>
         LinearEquationSolverProblemFormat NativeLinearEquationSolver<ValueType>::getEquationProblemFormat(Environment const& env) const {
             auto method = getMethod(env, storm::NumberTraits<ValueType>::IsExact);
-            if (method == NativeLinearEquationSolverMethod::Power || method == NativeLinearEquationSolverMethod::RationalSearch) {
+            if (method == NativeLinearEquationSolverMethod::Power || method == NativeLinearEquationSolverMethod::RationalSearch || method == NativeLinearEquationSolverMethod::QuickPower) {
                 return LinearEquationSolverProblemFormat::FixedPointSystem;
             } else {
                 return LinearEquationSolverProblemFormat::EquationSystem;

src/storm/solver/NativeLinearEquationSolver.h

@@ -71,6 +71,7 @@ namespace storm {
             virtual bool solveEquationsWalkerChae(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             virtual bool solveEquationsPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             virtual bool solveEquationsSoundPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
+            virtual bool solveEquationsQuickSoundPower(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
             virtual bool solveEquationsRationalSearch(storm::Environment const& env, std::vector<ValueType>& x, std::vector<ValueType> const& b) const;
 
             template<typename RationalType, typename ImpreciseType>

src/storm/solver/SolverSelectionOptions.cpp

@@ -78,6 +78,8 @@ namespace storm {
                     return "Power";
                 case NativeLinearEquationSolverMethod::RationalSearch:
                     return "RationalSearch";
+                case NativeLinearEquationSolverMethod::QuickPower:
+                    return "QuickPower";
             }
             return "invalid";
         }

src/storm/solver/SolverSelectionOptions.h

@@ -14,7 +14,7 @@ namespace storm {
         ExtendEnumsWithSelectionField(EquationSolverType, Native, Gmmxx, Eigen, Elimination)
         ExtendEnumsWithSelectionField(SmtSolverType, Z3, Mathsat)
         
-        ExtendEnumsWithSelectionField(NativeLinearEquationSolverMethod, Jacobi, GaussSeidel, SOR, WalkerChae, Power, RationalSearch)
+        ExtendEnumsWithSelectionField(NativeLinearEquationSolverMethod, Jacobi, GaussSeidel, SOR, WalkerChae, Power, RationalSearch, QuickPower)
         ExtendEnumsWithSelectionField(GmmxxLinearEquationSolverMethod, Bicgstab, Qmr, Gmres)
         ExtendEnumsWithSelectionField(GmmxxLinearEquationSolverPreconditioner, Ilu, Diagonal, None)
         ExtendEnumsWithSelectionField(EigenLinearEquationSolverMethod, SparseLU, Bicgstab, DGmres, Gmres)

src/storm/utility/constants.cpp

@@ -48,7 +48,11 @@ namespace storm {
         bool isAlmostZero(double const& a) {
             return a < 1e-12 && a > -1e-12;
         }
-
+        
+        bool isAlmostOne(double const& a) {
+            return a < (1.0 + 1e-12) && a > (1.0 - 1e-12);
+        }
+        
         template<typename ValueType>
         bool isConstant(ValueType const&) {
             return true;

src/storm/utility/constants.h

@@ -44,6 +44,8 @@ namespace storm {
         bool isZero(ValueType const& a);
 
         bool isAlmostZero(double const& a);
+
+        bool isAlmostOne(double const& a);
         
         template<typename ValueType>
         bool isConstant(ValueType const& a);

src/test/storm/modelchecker/DtmcPrctlModelCheckerTest.cpp

@@ -222,6 +222,24 @@ namespace {
         }
     };
 
+    class SparseNativeQuickSoundPowerEnvironment {
+    public:
+        static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
+        static const storm::settings::modules::CoreSettings::Engine engine = storm::settings::modules::CoreSettings::Engine::Sparse;
+        static const bool isExact = false;
+        typedef double ValueType;
+        typedef storm::models::sparse::Dtmc<ValueType> ModelType;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setForceSoundness(true);
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
+            env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::QuickPower);
+            env.solver().native().setRelativeTerminationCriterion(false);
+            env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-6));
+            return env;
+        }
+    };
+
     class SparseNativeRationalSearchEnvironment {
     public:
         static const storm::dd::DdType ddType = storm::dd::DdType::Sylvan; // unused for sparse models
@@ -447,6 +465,7 @@ namespace {
             SparseNativeSorEnvironment,
             SparseNativePowerEnvironment,
             SparseNativeSoundPowerEnvironment,
+            SparseNativeQuickSoundPowerEnvironment,
             SparseNativeRationalSearchEnvironment,
             HybridSylvanGmmxxGmresEnvironment,
             HybridCuddNativeJacobiEnvironment,

src/test/storm/solver/LinearEquationSolverTest.cpp

@@ -38,6 +38,21 @@ namespace {
         }
     };
     
+    class NativeDoubleQuickSoundPowerEnvironment {
+    public:
+        typedef double ValueType;
+        static const bool isExact = false;
+        static storm::Environment createEnvironment() {
+            storm::Environment env;
+            env.solver().setForceSoundness(true);
+            env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
+            env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::QuickPower);
+            env.solver().native().setRelativeTerminationCriterion(false);
+            env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber, std::string>("1e-6"));
+            return env;
+        }
+    };
+    
     class NativeDoubleJacobiEnvironment {
     public:
         typedef double ValueType;
@@ -265,6 +280,7 @@ namespace {
     typedef ::testing::Types<
             NativeDoublePowerEnvironment,
             NativeDoubleSoundPowerEnvironment,
+            NativeDoubleQuickSoundPowerEnvironment,
             NativeDoubleJacobiEnvironment,
             NativeDoubleGaussSeidelEnvironment,
             NativeDoubleSorEnvironment,