Added LRA settings.

src/storm-pars/settings/ParsSettings.cpp

@@ -18,6 +18,7 @@
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 #include "storm/settings/modules/TopologicalEquationSolverSettings.h"
 #include "storm/settings/modules/EliminationSettings.h"
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
 #include "storm/settings/modules/GameSolverSettings.h"
 #include "storm/settings/modules/BisimulationSettings.h"
@@ -48,6 +49,7 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::NativeEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::TopologicalEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::EliminationSettings>();
+            storm::settings::addModule<storm::settings::modules::LongRunAverageSolverSettings>();
             storm::settings::addModule<storm::settings::modules::MinMaxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::GameSolverSettings>();
             storm::settings::addModule<storm::settings::modules::BisimulationSettings>();

src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp

@@ -14,6 +14,7 @@
 #include "storm/storage/StronglyConnectedComponentDecomposition.h"
 
 #include "storm/adapters/RationalFunctionAdapter.h"
+#include "storm/environment/solver/LongRunAverageSolverEnvironment.h"
 
 #include "storm/utility/macros.h"
 #include "storm/utility/vector.h"
@@ -559,17 +560,30 @@ namespace storm {
                 auto subm = rateMatrix.getSubmatrix(false,bsccStates,bsccStates);
                 std::cout << std::endl << subm << std::endl;
                 
-                std::cout << std::endl << "========== EQSYS =======" << std::endl;
 
-                // Catch the trivial case for bsccs of size 1
-                if (bscc.size() == 1) {
-                    std::cout << std::endl << "========== Result = " << valueGetter(*bscc.begin()) << " =======" << std::endl;
-
-                    return valueGetter(*bscc.begin());
+                // Catch the case where all values are the same (this includes the special case where the bscc is of size 1)
+                auto it = bscc.begin();
+                ValueType val = valueGetter(*it);
+                for (++it; it != bscc.end(); ++it) {
+                    if (valueGetter(*it) != val) {
+                        break;
+                    }
+                }
+                if (it == bscc.end()) {
+                    // All entries have the same LRA
+                    std::cout << std::endl << "==== TRIVIAL: RESULT=" << val << "=======" << std::endl;
+                    return val;
                 }
                 
-                return computeLongRunAveragesForBsccVi(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector);
-                return computeLongRunAveragesForBsccEqSys(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector);
+                std::cout << std::endl << "========== SOLVING =======" << std::endl;
+                storm::solver::LraMethod method = env.solver().lra().getDetLraMethod();
+                if (method == storm::solver::LraMethod::ValueIteration) {
+                    return computeLongRunAveragesForBsccVi<ValueType>(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector);
+                } else if (method == storm::solver::LraMethod::LraDistributionEquations) {
+                    //return computeLongRunAveragesForBsccLraDistrEqSys<ValueType>(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector);
+                }
+                STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to gain-bias-equations.");
+                return computeLongRunAveragesForBsccEqSys<ValueType>(env, bscc, rateMatrix, backwardTransitions, valueGetter, exitRateVector);
             }
             
             template <>
@@ -593,7 +607,7 @@ namespace storm {
                 }
                 // To ensure that the model is aperiodic, we need to make sure that every Markovian state gets a self loop.
                 // Hence, we increase the uniformization rate a little.
-                uniformizationRate += storm::utility::one<ValueType>(); // Todo: try other values such as *=1.01
+                uniformizationRate *= (storm::utility::one<ValueType>() + storm::utility::convertNumber<ValueType>(env.solver().lra().getAperiodicFactor()));
 
                 // Get the transitions of the submodel
                 typename storm::storage::SparseMatrix<ValueType> bsccMatrix = rateMatrix.getSubmatrix(true, bsccStates, bsccStates, true);
@@ -619,48 +633,69 @@ namespace storm {
                 std::vector<ValueType> markovianRewards;
                 markovianRewards.reserve(bsccMatrix.getRowCount());
                 for (auto const& state : bsccStates) {
-                    ValueType stateRewardScalingFactor = storm::utility::one<ValueType>() / uniformizationRate;
-                    // ValueType actionRewardScalingFactor = exitRateVector[state] / uniformizationRate; // action rewards are currently not supported
-                    markovianRewards.push_back(valueGetter(state) * stateRewardScalingFactor);
+                    markovianRewards.push_back(valueGetter(state) / uniformizationRate);
                 }
                 
                 // start the iterations
-                ValueType precision = storm::utility::convertNumber<ValueType>(1e-6) / uniformizationRate; // TODO env.solver().minMax().getPrecision()) / uniformizationRate;
-                bool relative = true; // TODO env.solver().minMax().getRelativeTerminationCriterion();
-                std::vector<ValueType> v(bsccMatrix.getRowCount(), storm::utility::zero<ValueType>());
-                std::vector<ValueType> w = v;
-                std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> solver;
-                while (true) {
+                ValueType precision = storm::utility::convertNumber<ValueType>(env.solver().lra().getPrecision()) / uniformizationRate;
+                bool relative = env.solver().lra().getRelativeTerminationCriterion();
+                if (!relative) {
+                    precision /= uniformizationRate;
+                }
+                std::vector<ValueType> x(bsccMatrix.getRowCount(), storm::utility::zero<ValueType>());
+                std::vector<ValueType> xPrime(x.size());
+                
+                auto multiplier = storm::solver::MultiplierFactory<ValueType>().create(env, bsccMatrix);
+                ValueType maxDiff, minDiff;
+                uint64_t iter = 0;
+                boost::optional<uint64_t> maxIter;
+                if (env.solver().lra().isMaximalIterationCountSet()) {
+                    maxIter = env.solver().lra().getMaximalIterationCount();
+                }
+                while (!maxIter.is_initialized() || iter < maxIter.get()) {
+                    ++iter;
                     // Compute the values for the markovian states. We also keep track of the maximal and minimal difference between two values (for convergence checking)
-                    auto vIt = v.begin();
-                    uint64_t row = 0;
-                    ValueType newValue = markovianRewards[row] + bsccMatrix.multiplyRowWithVector(row, w);
-                    ValueType maxDiff = newValue - *vIt;
-                    ValueType minDiff = maxDiff;
-                    *vIt = newValue;
-                    for (++vIt, ++row; row < bsccMatrix.getRowCount(); ++vIt, ++row) {
-                        newValue = markovianRewards[row] + bsccMatrix.multiplyRowWithVector(row, w);
-                        ValueType diff = newValue - *vIt;
+                    multiplier->multiply(env, x, &markovianRewards, xPrime);
+                    
+                    // update xPrime and check for convergence
+                    // to avoid large (and numerically unstable) x-values, we substract a reference value.
+                    auto xIt = x.begin();
+                    auto xPrimeIt = xPrime.begin();
+                    ValueType refVal = *xPrimeIt;
+                    maxDiff = *xPrimeIt - *xIt;
+                    minDiff = maxDiff;
+                    *xPrimeIt -= refVal;
+                    *xIt = *xPrimeIt;
+                    for (++xIt, ++xPrimeIt; xIt != x.end(); ++xIt, ++xPrimeIt) {
+                        ValueType diff = *xPrimeIt - *xIt;
                         maxDiff = std::max(maxDiff, diff);
                         minDiff = std::min(minDiff, diff);
-                        *vIt = newValue;
+                        *xPrimeIt -= refVal;
+                        *xIt = *xPrimeIt;
                     }
-
-                    // Check for convergence
-                    if ((maxDiff - minDiff) <= (relative ? (precision * (v.front() + minDiff)) : precision)) {
+                    
+                    // Check for convergence. The uniformization rate is already incorporated into the precision parameter
+                    if ((maxDiff - minDiff) <= (relative ? (precision * minDiff) : precision)) {
                         break;
                     }
                     
-                    // update the rhs of the equation system
-                    ValueType referenceValue = v.front();
-                    storm::utility::vector::applyPointwise<ValueType, ValueType>(v, w, [&referenceValue] (ValueType const& v_i) -> ValueType { return v_i - referenceValue; });
                 }
-                std::cout << std::endl << "========== VI Result = " << v.front() * uniformizationRate << " =======" << std::endl;
-                return v.front() * uniformizationRate;
+                return (maxDiff + minDiff) * uniformizationRate / (storm::utility::convertNumber<ValueType>(2.0));
+            }
+            
+            template <>
+            storm::RationalFunction SparseCtmcCslHelper::computeLongRunAveragesForBsccEqSys<storm::RationalFunction>(Environment const&, storm::storage::StronglyConnectedComponent const&, storm::storage::SparseMatrix<storm::RationalFunction> const&, storm::storage::SparseMatrix<storm::RationalFunction> const&, std::function<storm::RationalFunction (storm::storage::sparse::state_type const& state)> const&, std::vector<storm::RationalFunction> const*) {
+                STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The requested Method for LRA computation is not supported for parametric models.");
+                // This is not possible due to uniformization
             }
             
             template <typename ValueType>
             ValueType SparseCtmcCslHelper::computeLongRunAveragesForBsccEqSys(Environment const& env, storm::storage::StronglyConnectedComponent const& bscc, storm::storage::SparseMatrix<ValueType> const& rateMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::function<ValueType (storm::storage::sparse::state_type const& state)> const& valueGetter, std::vector<ValueType> const* exitRateVector) {
+                // We build the equation system as in Line 3 of Algorithm 3 from
+                // Kretinsky, Meggendorfer: Efficient Strategy Iteration for Mean Payoff in Markov Decision Processes (ATVA 2017)
+                // The first variable corresponds to the gain of the bscc whereas the subsequent variables yield the bias for each state s_1, s_2, ....
+                // No bias variable for s_0 is needed since it is always set to zero, yielding an nxn equation system matrix
+                // To properly deal with CMTC, we also need to uniformize the transitions.
                 
                 // Get a mapping from global state indices to local ones.
                 std::unordered_map<uint64_t, uint64_t> toLocalIndexMap;
@@ -670,104 +705,91 @@ namespace storm {
                     ++localIndex;
                 }
                 
-                // Build the equation system matrix
-                
+                // Get the uniformization rate
+                ValueType uniformizationRate = storm::utility::one<ValueType>();
+                if (exitRateVector) {
+                    uniformizationRate = *std::max_element(exitRateVector->begin(), exitRateVector->end());
+                }
                 
-                // Build the equation system matrix for A[s,s] =  R(s,s) - r(s) & A[s,s'] = R(s,s') (where s != s')
-                // x_0+...+x_n=1 & x*A=0  <=>  x_0+...+x_n=1 & A^t*x=0 [ <=> 1-x_1+...+x_n=x_0 & (1-A^t)*x = x ]
+                // Build the equation system matrix and vector.
                 storm::solver::GeneralLinearEquationSolverFactory<ValueType> linearEquationSolverFactory;
-                storm::storage::SparseMatrixBuilder<ValueType> builder(bscc.size(), bscc.size());
                 bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem;
-                // The first row asserts that the values sum up to one
+                storm::storage::SparseMatrixBuilder<ValueType> builder(bscc.size(), bscc.size());
+                std::vector<ValueType> eqsysVector;
+                eqsysVector.reserve(bscc.size());
+                // The first row asserts that the weighted bias variables and the reward at s_0 sum up to the gain
                 uint64_t row = 0;
-                if (isEquationSystemFormat) {
-                    for (uint64_t state = 0; state < bscc.size(); ++state) {
-                            builder.addNextValue(row, state, storm::utility::one<ValueType>());
-                    }
-                } else {
-                    for (uint64_t state = 1; state < bscc.size(); ++state) {
-                            builder.addNextValue(row, state, -storm::utility::one<ValueType>());
+                ValueType entryValue;
+                for (auto const& globalState : bscc) {
+                    // Coefficient for the gain variable
+                    if (isEquationSystemFormat) {
+                        // '1-0' in row 0 and -(-1) in other rows
+                        builder.addNextValue(row, 0, storm::utility::one<ValueType>());
+                    } else if (row > 0) {
+                        // No coeficient in row 0, othwerise substract the gain
+                        builder.addNextValue(row, 0, -storm::utility::one<ValueType>());
                     }
-                }
-                ++row;
-                // Build the equation system matrix
-                // We can skip the first state to make the equation system matrix square
-                auto stateIt = bscc.begin();
-                for (++stateIt; stateIt != bscc.end(); ++stateIt) {
-                    ValueType diagonalValue = (exitRateVector == nullptr) ? -storm::utility::one<ValueType>() : -(*exitRateVector)[*stateIt];
-                    if (!isEquationSystemFormat) {
-                        diagonalValue = storm::utility::one<ValueType>() - diagonalValue;
+                    // Compute weighted sum over successor state. As this is a BSCC, each successor state will again be in the BSCC.
+                    bool needDiagonalEntry = ((exitRateVector && (*exitRateVector)[globalState] != uniformizationRate) || isEquationSystemFormat) && (row > 0);
+                    ValueType diagonalValue;
+                    if (needDiagonalEntry) {
+                        if (isEquationSystemFormat) {
+                            diagonalValue = (*exitRateVector)[globalState] / uniformizationRate; // 1 - (1 - r(s)/unifrate)
+                        } else {
+                            diagonalValue = storm::utility::one<ValueType>() - (*exitRateVector)[globalState] / uniformizationRate;
+                        }
                     }
-                    bool insertedDiagonal = storm::utility::isZero(diagonalValue);
-                    for (auto const& backwardsEntry : backwardTransitions.getRow(*stateIt)) {
-                        auto localIndexIt = toLocalIndexMap.find(backwardsEntry.getColumn());
-                        if (localIndexIt != toLocalIndexMap.end()) {
-                            ValueType val = backwardsEntry.getValue();
-                            if (!isEquationSystemFormat) {
-                                val = -val;
-                            }
-                            uint64_t localIndex = localIndexIt->second;
-                            if (!insertedDiagonal && localIndex == row) {
-                                builder.addNextValue(row, localIndex, val + diagonalValue);
-                                insertedDiagonal = true;
-                            } else {
-                                if (!insertedDiagonal && localIndex > row) {
-                                    builder.addNextValue(row, row, diagonalValue);
-                                    insertedDiagonal = true;
-                                }
-                                builder.addNextValue(row, localIndex, val);
+                    for (auto const& entry : rateMatrix.getRow(globalState)) {
+                        uint64_t col = toLocalIndexMap[entry.getColumn()];
+                        if (col == 0) {
+                            //Skip transition to state_0. This corresponds to setting the bias of state_0 to zero
+                            continue;
+                        }
+                        entryValue = entry.getValue();
+                        if (exitRateVector) {
+                     //       entryValue /= uniformizationRate;
+                        }
+                        if (isEquationSystemFormat) {
+                            entryValue = -entryValue;
+                        }
+                        if (needDiagonalEntry && col >= row) {
+                            if (col == row) {
+                                entryValue += diagonalValue;
+                            } else { // col > row
+                                builder.addNextValue(row, row, diagonalValue);
                             }
+                            needDiagonalEntry = false;
                         }
+                        builder.addNextValue(row, col, entryValue);
                     }
-                    if (!insertedDiagonal) {
+                    if (needDiagonalEntry) {
                         builder.addNextValue(row, row, diagonalValue);
                     }
+                    if (exitRateVector) {
+                        eqsysVector.push_back(valueGetter(globalState) / uniformizationRate);
+                    } else {
+                        eqsysVector.push_back(valueGetter(globalState));
+                    }
                     ++row;
                 }
 
                 // Create a linear equation solver
                 auto matrix = builder.build();
                 std::cout << matrix << std::endl;
+                std::cout << "RHS: " << storm::utility::vector::toString(eqsysVector) << std::endl;
                 auto solver = linearEquationSolverFactory.create(env, std::move(matrix));
-                solver->setBounds(storm::utility::zero<ValueType>(), storm::utility::one<ValueType>());
                 // Check solver requirements.
                 auto requirements = solver->getRequirements(env);
                 STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked.");
+                // Todo: Find bounds on the bias variable. Just inserting the maximal value from the vector probably does not work.
                 
-                std::vector<ValueType> bsccEquationSystemRightSide(bscc.size(), storm::utility::zero<ValueType>());
-                bsccEquationSystemRightSide.front() = storm::utility::one<ValueType>();
-                std::vector<ValueType> bsccEquationSystemSolution(bscc.size(), storm::utility::one<ValueType>() / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size()));
-                std::cout << storm::utility::vector::toString(bsccEquationSystemRightSide) << std::endl;
-                solver->solveEquations(env, bsccEquationSystemSolution, bsccEquationSystemRightSide);
-                std::cout << storm::utility::vector::toString(bsccEquationSystemSolution) << std::endl;
-                // If exit rates were given, we need to 'fix' the results to also account for the timing behaviour.
-                if (false && exitRateVector != nullptr) {
-                    ValueType totalValue = storm::utility::zero<ValueType>();
-                    auto solIt = bsccEquationSystemSolution.begin();
-                    for (auto const& globalState : bscc) {
-                        totalValue += (*solIt) * (storm::utility::one<ValueType>() / (*exitRateVector)[globalState]);
-                        ++solIt;
-                    }
-                    assert(solIt == bsccEquationSystemSolution.end());
-                    solIt = bsccEquationSystemSolution.begin();
-                    for (auto const& globalState : bscc) {
-                        *solIt = ((*solIt) * (storm::utility::one<ValueType>() / (*exitRateVector)[globalState])) / totalValue;
-                        ++solIt;
-                    }
-                    assert(solIt == bsccEquationSystemSolution.end());
-                }
+                std::vector<ValueType> eqSysSol(bscc.size(), storm::utility::zero<ValueType>());
+                // Take the mean of the rewards as an initial guess for the gain
+                eqSysSol.front() = std::accumulate(eqsysVector.begin(), eqsysVector.end(), storm::utility::zero<ValueType>()) / storm::utility::convertNumber<ValueType, uint64_t>(bscc.size());
+                solver->solveEquations(env, eqSysSol, eqsysVector);
                 
-                // Calculate final LRA Value
-                ValueType result = storm::utility::zero<ValueType>();
-                auto solIt = bsccEquationSystemSolution.begin();
-                for (auto const& globalState : bscc) {
-                    result += valueGetter(globalState) * (*solIt);
-                    ++solIt;
-                }
-                assert(solIt == bsccEquationSystemSolution.end());
-                std::cout << std::endl << "========== Result = " << result << " =======" << std::endl;
-
-                return result;
+                std::cout << std::endl << "========== Result = " << eqSysSol.front() << "*" << uniformizationRate << " =======" << std::endl;
+                return eqSysSol.front();// * uniformizationRate;
             }
             
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>

src/storm/settings/SettingsManager.cpp

@@ -24,6 +24,7 @@
 #include "storm/settings/modules/GmmxxEquationSolverSettings.h"
 #include "storm/settings/modules/NativeEquationSolverSettings.h"
 #include "storm/settings/modules/EliminationSettings.h"
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
 #include "storm/settings/modules/MinMaxEquationSolverSettings.h"
 #include "storm/settings/modules/GameSolverSettings.h"
 #include "storm/settings/modules/BisimulationSettings.h"
@@ -654,6 +655,7 @@ namespace storm {
             storm::settings::addModule<storm::settings::modules::EigenEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::NativeEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::EliminationSettings>();
+            storm::settings::addModule<storm::settings::modules::LongRunAverageSolverSettings>();
             storm::settings::addModule<storm::settings::modules::MinMaxEquationSolverSettings>();
             storm::settings::addModule<storm::settings::modules::GameSolverSettings>();
             storm::settings::addModule<storm::settings::modules::BisimulationSettings>();

src/storm/settings/modules/LongRunAverageSolverSettings.cpp

@@ -0,0 +1,101 @@
+#include "storm/settings/modules/LongRunAverageSolverSettings.h"
+
+#include "storm/settings/Option.h"
+#include "storm/settings/ArgumentBuilder.h"
+#include "storm/settings/OptionBuilder.h"
+
+#include "storm/utility/macros.h"
+#include "storm/exceptions/IllegalArgumentValueException.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            const std::string LongRunAverageSolverSettings::moduleName = "lra";
+            const std::string LongRunAverageSolverSettings::detLraMethodOptionName = "detmethod";
+            const std::string LongRunAverageSolverSettings::nondetLraMethodOptionName = "nondetmethod";
+            const std::string LongRunAverageSolverSettings::maximalIterationsOptionName = "maxiter";
+            const std::string LongRunAverageSolverSettings::maximalIterationsOptionShortName = "i";
+            const std::string LongRunAverageSolverSettings::precisionOptionName = "precision";
+            const std::string LongRunAverageSolverSettings::absoluteOptionName = "absolute";
+            const std::string LongRunAverageSolverSettings::aperiodicFactorOptionName = "aperiodicfactor";
+
+            LongRunAverageSolverSettings::LongRunAverageSolverSettings() : ModuleSettings(moduleName) {
+                
+                std::vector<std::string> detLraMethods = {"gb", "gain-bias-equations", "distr", "lra-distribution-equations", "vi", "value-iteration"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, detLraMethodOptionName, true, "Sets which method is preferred for computing long run averages on deterministic models.").setIsAdvanced()
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(detLraMethods)).setDefaultValueString("gb").build()).build());
+                
+                std::vector<std::string> nondetLraMethods = {"vi", "value-iteration", "linear-programming", "lp"};
+                this->addOption(storm::settings::OptionBuilder(moduleName, nondetLraMethodOptionName, true, "Sets which method is preferred for computing long run averages on models with nondeterminism.").setIsAdvanced()
+                                .addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of a long run average computation method.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(nondetLraMethods)).setDefaultValueString("vi").build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").build()).build());
+                
+                this->addOption(storm::settings::OptionBuilder(moduleName, precisionOptionName, false, "The precision used for detecting convergence of iterative methods.").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The precision to achieve.").setDefaultValueDouble(1e-06).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+
+                this->addOption(storm::settings::OptionBuilder(moduleName, absoluteOptionName, false, "Sets whether the relative or the absolute error is considered for detecting convergence.").setIsAdvanced().build());
+
+                this->addOption(storm::settings::OptionBuilder(moduleName, aperiodicFactorOptionName, true, "If required by the selected method (e.g. vi), this factor controls how the system is made aperiodic").setIsAdvanced().addArgument(storm::settings::ArgumentBuilder::createDoubleArgument("value", "The factor.").setDefaultValueDouble(0.125).addValidatorDouble(ArgumentValidatorFactory::createDoubleRangeValidatorExcluding(0.0, 1.0)).build()).build());
+                
+            }
+            
+            storm::solver::LraMethod LongRunAverageSolverSettings::getDetLraMethod() const {
+                std::string lraMethodString = this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getValueAsString();
+                if (lraMethodString == "gain-bias-equations" || lraMethodString == "gb") {
+                    return storm::solver::LraMethod::GainBiasEquations;
+                }
+                if (lraMethodString == "lra-distribution-equations" || lraMethodString == "distr") {
+                    return storm::solver::LraMethod::LraDistributionEquations;
+                }
+                if (lraMethodString == "value-iteration" || lraMethodString == "vi") {
+                    return storm::solver::LraMethod::ValueIteration;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique for deterministic models:'" << lraMethodString << "'.");
+            }
+            
+            bool LongRunAverageSolverSettings::isDetLraMethodSetFromDefaultValue() const {
+                return !this->getOption(detLraMethodOptionName).getArgumentByName("name").getHasBeenSet() || this->getOption(detLraMethodOptionName).getArgumentByName("name").wasSetFromDefaultValue();
+            }
+            
+            storm::solver::LraMethod LongRunAverageSolverSettings::getNondetLraMethod() const {
+                std::string lraMethodString = this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getValueAsString();
+                if (lraMethodString == "value-iteration" || lraMethodString == "vi") {
+                    return storm::solver::LraMethod::ValueIteration;
+                } else if (lraMethodString == "linear-programming" || lraMethodString == "lp") {
+                    return storm::solver::LraMethod::LinearProgramming;
+                }
+                STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown lra solving technique for nondeterministic models:'" << lraMethodString << "'.");
+            }
+            
+            bool LongRunAverageSolverSettings::isNondetLraMethodSetFromDefaultValue() const {
+                return !this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getHasBeenSet() || this->getOption(nondetLraMethodOptionName).getArgumentByName("name").wasSetFromDefaultValue();
+            }
+            
+            bool LongRunAverageSolverSettings::isMaximalIterationCountSet() const {
+                return this->getOption(maximalIterationsOptionName).getHasOptionBeenSet();
+            }
+            
+            uint_fast64_t LongRunAverageSolverSettings::getMaximalIterationCount() const {
+                return this->getOption(maximalIterationsOptionName).getArgumentByName("count").getValueAsUnsignedInteger();
+            }
+            
+            bool LongRunAverageSolverSettings::isPrecisionSet() const {
+                return this->getOption(precisionOptionName).getHasOptionBeenSet();
+            }
+            
+            double LongRunAverageSolverSettings::getPrecision() const {
+                return this->getOption(precisionOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+            bool LongRunAverageSolverSettings::isRelativePrecision() const {
+                return this->getOption(absoluteOptionName).getHasOptionBeenSet();
+            }
+            
+            double LongRunAverageSolverSettings::getAperiodicFactor() const {
+                return this->getOption(aperiodicFactorOptionName).getArgumentByName("value").getValueAsDouble();
+            }
+            
+        }
+    }
+}

src/storm/settings/modules/LongRunAverageSolverSettings.h

@@ -0,0 +1,92 @@
+#pragma once
+
+#include "storm-config.h"
+#include "storm/settings/modules/ModuleSettings.h"
+
+#include "storm/solver/SolverSelectionOptions.h"
+#include "storm/solver/MultiplicationStyle.h"
+
+namespace storm {
+    namespace settings {
+        namespace modules {
+            
+            /*!
+             * This class represents the LRA solver settings.
+             */
+            class LongRunAverageSolverSettings : public ModuleSettings {
+            public:
+                
+                LongRunAverageSolverSettings();
+              
+                /*!
+                 * Retrieves the selected long run average method for deterministic models.
+                 */
+                storm::solver::LraMethod getDetLraMethod() const;
+                
+                /*!
+                 * Retrieves whether the LraMethod for deterministic models was set from a default value.
+                 */
+                bool isDetLraMethodSetFromDefaultValue() const;
+                
+                /*!
+                 * Retrieves the selected long run average method for nondeterministic models.
+                 */
+                storm::solver::LraMethod getNondetLraMethod() const;
+                
+                /*!
+                 * Retrieves whether the LraMethod for nondeterministic models was set from a default value.
+                 */
+                bool isNondetLraMethodSetFromDefaultValue() const;
+                
+                /*!
+                 * Retrieves whether a maximal iteration count for iterative methods was set.
+                 */
+                bool isMaximalIterationCountSet() const;
+                
+                /*!
+                 * Retrieves the maximal number of iterations to perform until giving up on converging.
+                 *
+                 * @return The maximal iteration count.
+                 */
+                uint_fast64_t getMaximalIterationCount() const;
+                
+                /*!
+                 * Retrieves whether the precision has been set.
+                 *
+                 * @return True iff the precision has been set.
+                 */
+                bool isPrecisionSet() const;
+                
+                /*!
+                 * Retrieves the precision that is used for detecting convergence.
+                 *
+                 * @return The precision to use for detecting convergence.
+                 */
+                double getPrecision() const;
+                
+                /*!
+                 * Retrieves whether the convergence criterion has been set to relative.
+                 */
+                bool isRelativePrecision() const;
+                
+                /*!
+                 * Retrieves a factor that describes how the system is made aperiodic (if necessary by the method)
+                 */
+                double getAperiodicFactor() const;
+                
+                // The name of the module.
+                static const std::string moduleName;
+                
+            private:
+                static const std::string detLraMethodOptionName;
+                static const std::string nondetLraMethodOptionName;
+                static const std::string maximalIterationsOptionName;
+                static const std::string maximalIterationsOptionShortName;
+                static const std::string precisionOptionName;
+                static const std::string absoluteOptionName;
+                static const std::string aperiodicFactorOptionName;
+            };
+            
+        }
+    }
+}

src/test/storm/modelchecker/prctl/dtmc/LraDtmcPrctlModelCheckerTest.cpp

@@ -52,7 +52,7 @@ namespace {
     
     class EigenRationalLUEnvironment {
     public:
-        typedef double ValueType;
+        typedef storm::RationalNumber ValueType;
         static const bool isExact = true;
         static storm::Environment createEnvironment() {
             storm::Environment env;
@@ -70,7 +70,7 @@ namespace {
             storm::Environment env;
             env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native);
             env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::SOR);
-            env.solver().native().setSorOmega(storm::utility::convertNumber<storm::RationalNumber>(0.9));
+            env.solver().native().setSorOmega(storm::utility::convertNumber<storm::RationalNumber>(0.8)); // A test fails if this is set to 0.9...
             env.solver().native().setPrecision(storm::utility::convertNumber<storm::RationalNumber>(1e-8));
             return env;
         }