diff --git a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp index 528182ea5..6d82a6d79 100644 --- a/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp +++ b/src/storm/modelchecker/csl/helper/SparseCtmcCslHelper.cpp @@ -15,6 +15,7 @@ #include "storm/adapters/RationalFunctionAdapter.h" #include "storm/environment/solver/LongRunAverageSolverEnvironment.h" +#include "storm/environment/solver/TopologicalSolverEnvironment.h" #include "storm/utility/macros.h" #include "storm/utility/vector.h" @@ -575,13 +576,21 @@ namespace storm { } storm::solver::LraMethod method = env.solver().lra().getDetLraMethod(); + if (storm::NumberTraits::IsExact && env.solver().lra().isDetLraMethodSetFromDefault() && method == storm::solver::LraMethod::ValueIteration) { + method = storm::solver::LraMethod::GainBiasEquations; + STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee exact results. If you want to override this, please explicitly specify a different LRA method."); + } else if (env.solver().isForceSoundness() && env.solver().lra().isDetLraMethodSetFromDefault() && method != storm::solver::LraMethod::ValueIteration) { + method = storm::solver::LraMethod::ValueIteration; + STORM_LOG_INFO("Selecting " << storm::solver::toString(method) << " as the solution technique for long-run properties to guarantee sound results. If you want to override this, please explicitly specify a different LRA method."); + } + STORM_LOG_TRACE("Computing LRA for BSCC of size " << bscc.size() << " using '" << storm::solver::toString(method) << "'."); if (method == storm::solver::LraMethod::ValueIteration) { return computeLongRunAveragesForBsccVi(env, bscc, rateMatrix, valueGetter, exitRateVector); } else if (method == storm::solver::LraMethod::LraDistributionEquations) { // We only need the first element of the pair as the lra distribution is not relevant at this point. return computeLongRunAveragesForBsccLraDistr(env, bscc, rateMatrix, valueGetter, exitRateVector).first; } - STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to gain-bias-equations."); + STORM_LOG_WARN_COND(method == storm::solver::LraMethod::GainBiasEquations, "Unsupported lra method selected. Defaulting to " << storm::solver::toString(storm::solver::LraMethod::GainBiasEquations) << "."); // We don't need the bias values return computeLongRunAveragesForBsccGainBias(env, bscc, rateMatrix, valueGetter, exitRateVector).first; } @@ -705,12 +714,20 @@ namespace storm { ++localIndex; } + // Prepare an environment for the underlying equation solver + auto subEnv = env; + if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) { + // Topological solver does not make any sense since the BSCC is connected. + subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType()); + } + subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); + // Build the equation system matrix and vector. storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; - bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + bool isEquationSystemFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; storm::storage::SparseMatrixBuilder builder(bscc.size(), bscc.size()); - std::vector eqsysVector; - eqsysVector.reserve(bscc.size()); + std::vector 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; ValueType entryValue; @@ -724,11 +741,15 @@ namespace storm { builder.addNextValue(row, 0, -storm::utility::one()); } // Compute weighted sum over successor state. As this is a BSCC, each successor state will again be in the BSCC. - bool needDiagonalEntry = isEquationSystemFormat && (row > 0); - ValueType diagonalValue; - if (needDiagonalEntry) { - diagonalValue = exitRateVector ? (*exitRateVector)[globalState] : storm::utility::one(); + auto diagonalValue = storm::utility::zero(); + if (row > 0) { + if (isEquationSystemFormat) { + diagonalValue = exitRateVector ? (*exitRateVector)[globalState] : storm::utility::one(); + } else { + diagonalValue = storm::utility::one() - (exitRateVector ? (*exitRateVector)[globalState] : storm::utility::one()); + } } + bool needDiagonalEntry = !storm::utility::isZero(diagonalValue); for (auto const& entry : rateMatrix.getRow(globalState)) { uint64_t col = toLocalIndexMap[entry.getColumn()]; if (col == 0) { @@ -752,23 +773,21 @@ namespace storm { if (needDiagonalEntry) { builder.addNextValue(row, row, diagonalValue); } - eqsysVector.push_back(valueGetter(globalState)); + eqSysVector.push_back(valueGetter(globalState)); ++row; } // Create a linear equation solver - auto subEnv = env; - subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); auto solver = linearEquationSolverFactory.create(subEnv, builder.build()); // Check solver requirements. auto requirements = solver->getRequirements(subEnv); 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. + // Todo: Find bounds on the bias variables. Just inserting the maximal value from the vector probably does not work. std::vector eqSysSol(bscc.size(), storm::utility::zero()); // Take the mean of the rewards as an initial guess for the gain - eqSysSol.front() = std::accumulate(eqsysVector.begin(), eqsysVector.end(), storm::utility::zero()) / storm::utility::convertNumber(bscc.size()); - solver->solveEquations(subEnv, eqSysSol, eqsysVector); + //eqSysSol.front() = std::accumulate(eqSysVector.begin(), eqSysVector.end(), storm::utility::zero()) / storm::utility::convertNumber(bscc.size()); + solver->solveEquations(subEnv, eqSysSol, eqSysVector); ValueType gain = eqSysSol.front(); // insert bias value for state 0 @@ -787,9 +806,17 @@ namespace storm { // This method assumes that this BSCC consist of more than one state if (bscc.size() == 1) { - return { storm::utility::one(), {storm::utility::one()} }; + return { valueGetter(*bscc.begin()), {storm::utility::one()} }; } + // Prepare an environment for the underlying linear equation solver + auto subEnv = env; + if (subEnv.solver().getLinearEquationSolverType() == storm::solver::EquationSolverType::Topological) { + // Topological solver does not make any sense since the BSCC is connected. + subEnv.solver().setLinearEquationSolverType(subEnv.solver().topological().getUnderlyingEquationSolverType()); + } + subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); + // Get a mapping from global state indices to local ones as well as a bitvector containing states within the BSCC. std::unordered_map toLocalIndexMap; storm::storage::BitVector bsccStates(rateMatrix.getRowCount(), false); @@ -823,12 +850,11 @@ namespace storm { // Check whether we need the fixpoint characterization storm::solver::GeneralLinearEquationSolverFactory linearEquationSolverFactory; - bool isFixpointFormat = linearEquationSolverFactory.getEquationProblemFormat(env) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem; - + bool isFixpointFormat = linearEquationSolverFactory.getEquationProblemFormat(subEnv) == storm::solver::LinearEquationSolverProblemFormat::FixedPointSystem; if (isFixpointFormat) { // Add a 1 on the diagonal for (row = 0; row < auxMatrix.getRowCount(); ++row) { - for (auto& entry : auxMatrix) { + for (auto& entry : auxMatrix.getRow(row)) { if (entry.getColumn() == row) { entry.setValue(storm::utility::one() + entry.getValue()); } @@ -869,32 +895,16 @@ namespace storm { bsccEquationSystemRightSide.back() = storm::utility::one(); // Create a linear equation solver - auto subEnv = env; - subEnv.solver().setLinearEquationSolverPrecision(env.solver().lra().getPrecision(), env.solver().lra().getRelativeTerminationCriterion()); auto solver = linearEquationSolverFactory.create(subEnv, builder.build()); solver->setBounds(storm::utility::zero(), storm::utility::one()); // Check solver requirements. auto requirements = solver->getRequirements(subEnv); + requirements.clearLowerBounds(); + requirements.clearUpperBounds(); STORM_LOG_THROW(!requirements.hasEnabledCriticalRequirement(), storm::exceptions::UncheckedRequirementException, "Solver requirements " + requirements.getEnabledRequirementsAsString() + " not checked."); std::vector lraDistr(bscc.size(), storm::utility::one() / storm::utility::convertNumber(bscc.size())); solver->solveEquations(subEnv, lraDistr, bsccEquationSystemRightSide); - // 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(); - auto solIt = lraDistr.begin(); - for (auto const& globalState : bscc) { - totalValue += (*solIt) * (storm::utility::one() / (*exitRateVector)[globalState]); - ++solIt; - } - assert(solIt == lraDistr.end()); - solIt = lraDistr.begin(); - for (auto const& globalState : bscc) { - *solIt = ((*solIt) * (storm::utility::one() / (*exitRateVector)[globalState])) / totalValue; - ++solIt; - } - assert(solIt == lraDistr.end()); - } // Calculate final LRA Value ValueType result = storm::utility::zero(); diff --git a/src/storm/settings/modules/LongRunAverageSolverSettings.cpp b/src/storm/settings/modules/LongRunAverageSolverSettings.cpp index 9826bdca8..8337bb779 100644 --- a/src/storm/settings/modules/LongRunAverageSolverSettings.cpp +++ b/src/storm/settings/modules/LongRunAverageSolverSettings.cpp @@ -41,7 +41,7 @@ namespace storm { } storm::solver::LraMethod LongRunAverageSolverSettings::getDetLraMethod() const { - std::string lraMethodString = this->getOption(nondetLraMethodOptionName).getArgumentByName("name").getValueAsString(); + std::string lraMethodString = this->getOption(detLraMethodOptionName).getArgumentByName("name").getValueAsString(); if (lraMethodString == "gain-bias-equations" || lraMethodString == "gb") { return storm::solver::LraMethod::GainBiasEquations; } diff --git a/src/storm/solver/SolverSelectionOptions.cpp b/src/storm/solver/SolverSelectionOptions.cpp index 26001226f..f635b5a38 100644 --- a/src/storm/solver/SolverSelectionOptions.cpp +++ b/src/storm/solver/SolverSelectionOptions.cpp @@ -49,9 +49,13 @@ namespace storm { std::string toString(LraMethod m) { switch(m) { case LraMethod::LinearProgramming: - return "linearprogramming"; + return "linear-programming"; case LraMethod::ValueIteration: - return "valueiteration"; + return "value-iteration"; + case LraMethod::LraDistributionEquations: + return "lra-distribution-equations"; + case LraMethod::GainBiasEquations: + return "gain-bias-equations"; } return "invalid"; } diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 857a8d478..d37177b6c 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -2301,6 +2301,7 @@ namespace storm { template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + template double SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; @@ -2313,6 +2314,7 @@ namespace storm { template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + template float SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; @@ -2341,6 +2343,7 @@ namespace storm { template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + template storm::ClnRationalNumber SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif @@ -2351,6 +2354,7 @@ namespace storm { template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + template storm::GmpRationalNumber SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template bool SparseMatrix::isSubmatrixOf(SparseMatrix const& matrix) const; #endif @@ -2361,6 +2365,10 @@ namespace storm { template class SparseMatrixBuilder; template class SparseMatrix; template std::ostream& operator<<(std::ostream& out, SparseMatrix const& matrix); + template storm::RationalFunction SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; + template storm::RationalFunction SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; + template storm::RationalFunction SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; + template storm::RationalFunction SparseMatrix::getPointwiseProductRowSum(storm::storage::SparseMatrix const& otherMatrix, typename SparseMatrix::index_type const& row) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; template std::vector SparseMatrix::getPointwiseProductRowSumVector(storm::storage::SparseMatrix const& otherMatrix) const; diff --git a/src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp b/src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp index 629684560..b8d3b1fc7 100755 --- a/src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp +++ b/src/test/storm/modelchecker/csl/LraCtmcCslModelCheckerTest.cpp @@ -144,9 +144,8 @@ namespace { env.solver().setLinearEquationSolverType(storm::solver::EquationSolverType::Native); env.solver().native().setMethod(storm::solver::NativeLinearEquationSolverMethod::SOR); env.solver().native().setSorOmega(storm::utility::convertNumber(0.7)); // LRA computation fails for 0.9 - env.solver().native().setPrecision(storm::utility::convertNumber(1e-9)); - env.solver().native().setRelativeTerminationCriterion(false); env.solver().lra().setDetLraMethod(storm::solver::LraMethod::GainBiasEquations); + env.solver().lra().setPrecision(storm::utility::convertNumber(1e-9)); return env; } };