diff --git a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp index 922adc4de..36e6eb9a2 100644 --- a/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp +++ b/src/storm/modelchecker/prctl/SparseDtmcPrctlModelChecker.cpp @@ -36,21 +36,32 @@ namespace storm { template bool SparseDtmcPrctlModelChecker::canHandle(CheckTask const& checkTask) const { storm::logic::Formula const& formula = checkTask.getFormula(); - return formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true)); + return formula.isInFragment(storm::logic::prctl().setLongRunAverageRewardFormulasAllowed(true).setLongRunAverageProbabilitiesAllowed(true).setConditionalProbabilityFormulasAllowed(true).setConditionalRewardFormulasAllowed(true).setOnlyEventuallyFormuluasInConditionalFormulasAllowed(true).setRewardBoundedUntilFormulasAllowed(true).setRewardBoundedCumulativeRewardFormulasAllowed(true).setMultiDimensionalBoundedUntilFormulasAllowed(true).setMultiDimensionalCumulativeRewardFormulasAllowed(true)); } template std::unique_ptr SparseDtmcPrctlModelChecker::computeBoundedUntilProbabilities(Environment const& env, CheckTask const& checkTask) { storm::logic::BoundedUntilFormula const& pathFormula = checkTask.getFormula(); - STORM_LOG_THROW(!pathFormula.hasLowerBound() && pathFormula.hasUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have single upper time bound."); - STORM_LOG_THROW(pathFormula.hasIntegerUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have discrete upper time bound."); - std::unique_ptr leftResultPointer = this->check(env, pathFormula.getLeftSubformula()); - std::unique_ptr rightResultPointer = this->check(env, pathFormula.getRightSubformula()); - ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); - ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); - std::vector numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeBoundedUntilProbabilities(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getNonStrictUpperBound(), *linearEquationSolverFactory, checkTask.getHint()); - std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); - return result; + if (pathFormula.isMultiDimensional() || pathFormula.getTimeBoundReference().isRewardBound()) { + STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Checking non-trivial bounded until probabilities can only be computed for the initial states of the model."); + storm::logic::OperatorInformation opInfo; + if (checkTask.isBoundSet()) { + opInfo.bound = checkTask.getBound(); + } + auto formula = std::make_shared(checkTask.getFormula().asSharedPointer(), opInfo); + auto numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeRewardBoundedValues(env, this->getModel(), formula, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + } else { + STORM_LOG_THROW(!pathFormula.hasLowerBound() && pathFormula.hasUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have single upper time bound."); + STORM_LOG_THROW(pathFormula.hasIntegerUpperBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have discrete upper time bound."); + std::unique_ptr leftResultPointer = this->check(env, pathFormula.getLeftSubformula()); + std::unique_ptr rightResultPointer = this->check(env, pathFormula.getRightSubformula()); + ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult(); + ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult(); + std::vector numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeStepBoundedUntilProbabilities(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), this->getModel().getBackwardTransitions(), leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), pathFormula.getNonStrictUpperBound(), *linearEquationSolverFactory, checkTask.getHint()); + std::unique_ptr result = std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + return result; + } } template @@ -85,9 +96,20 @@ namespace storm { template std::unique_ptr SparseDtmcPrctlModelChecker::computeCumulativeRewards(Environment const& env, storm::logic::RewardMeasureType, CheckTask const& checkTask) { storm::logic::CumulativeRewardFormula const& rewardPathFormula = checkTask.getFormula(); - STORM_LOG_THROW(rewardPathFormula.hasIntegerBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have a discrete time bound."); - std::vector numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeCumulativeRewards(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getRewardModel(""), rewardPathFormula.getNonStrictBound(), *linearEquationSolverFactory); - return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + if (rewardPathFormula.isMultiDimensional() || rewardPathFormula.getTimeBoundReference().isRewardBound()) { + STORM_LOG_THROW(checkTask.isOnlyInitialStatesRelevantSet(), storm::exceptions::InvalidOperationException, "Checking non-trivial bounded until probabilities can only be computed for the initial states of the model."); + storm::logic::OperatorInformation opInfo; + if (checkTask.isBoundSet()) { + opInfo.bound = checkTask.getBound(); + } + auto formula = std::make_shared(checkTask.getFormula().asSharedPointer(), checkTask.getRewardModel(), opInfo); + auto numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeRewardBoundedValues(env, this->getModel(), formula, *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + } else { + STORM_LOG_THROW(rewardPathFormula.hasIntegerBound(), storm::exceptions::InvalidPropertyException, "Formula needs to have a discrete time bound."); + std::vector numericResult = storm::modelchecker::helper::SparseDtmcPrctlHelper::computeCumulativeRewards(env, storm::solver::SolveGoal(this->getModel(), checkTask), this->getModel().getTransitionMatrix(), checkTask.isRewardModelSet() ? this->getModel().getRewardModel(checkTask.getRewardModel()) : this->getModel().getRewardModel(""), rewardPathFormula.getNonStrictBound(), *linearEquationSolverFactory); + return std::unique_ptr(new ExplicitQuantitativeCheckResult(std::move(numericResult))); + } } template diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp index 8d3f0f68d..fd52e4437 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.cpp @@ -15,6 +15,18 @@ #include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h" #include "storm/modelchecker/hints/ExplicitModelCheckerHint.h" #include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" +#include "storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h" + +#include "storm/environment/solver/SolverEnvironment.h" + +#include "storm/settings/SettingsManager.h" +#include "storm/settings/modules/GeneralSettings.h" +#include "storm/settings/modules/CoreSettings.h" +#include "storm/settings/modules/IOSettings.h" + +#include "storm/utility/Stopwatch.h" +#include "storm/utility/ProgressMeasurement.h" +#include "storm/utility/export.h" #include "storm/utility/macros.h" #include "storm/utility/ConstantsComparator.h" @@ -29,12 +41,12 @@ namespace storm { namespace modelchecker { namespace helper { template - std::vector SparseDtmcPrctlHelper::computeBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, ModelCheckerHint const& hint) { + std::vector SparseDtmcPrctlHelper::computeStepBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, ModelCheckerHint const& hint) { std::vector result(transitionMatrix.getRowCount(), storm::utility::zero()); // If we identify the states that have probability 0 of reaching the target states, we can exclude them in the further analysis. storm::storage::BitVector maybeStates; - + if (hint.isExplicitModelCheckerHint() && hint.template asExplicitModelCheckerHint().getComputeOnlyMaybeStates()) { maybeStates = hint.template asExplicitModelCheckerHint().getMaybeStates(); } else { @@ -68,8 +80,168 @@ namespace storm { return result; } + template + std::vector analyzeTrivialDtmcEpochModel(typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { + + std::vector epochResult; + epochResult.reserve(epochModel.epochInStates.getNumberOfSetBits()); + auto stepSolutionIt = epochModel.stepSolutions.begin(); + auto stepChoiceIt = epochModel.stepChoices.begin(); + for (auto const& state : epochModel.epochInStates) { + while (*stepChoiceIt < state) { + ++stepChoiceIt; + ++stepSolutionIt; + } + if (epochModel.objectiveRewardFilter.front().get(state)) { + if (*stepChoiceIt == state) { + epochResult.push_back(epochModel.objectiveRewards.front()[state] + *stepSolutionIt); + } else { + epochResult.push_back(epochModel.objectiveRewards.front()[state]); + } + } else { + if (*stepChoiceIt == state) { + epochResult.push_back(*stepSolutionIt); + } else { + epochResult.push_back(storm::utility::zero()); + } + } + } + return epochResult; + } + + template + std::vector analyzeNonTrivialDtmcEpochModel(Environment const& env, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& linEqSolver, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, boost::optional const& lowerBound, boost::optional const& upperBound) { + + // Update some data for the case that the Matrix has changed + if (epochModel.epochMatrixChanged) { + x.assign(epochModel.epochMatrix.getRowGroupCount(), storm::utility::zero()); + linEqSolver = linearEquationSolverFactory.create(env, epochModel.epochMatrix, storm::solver::LinearEquationSolverTask::SolveEquations); + linEqSolver->setCachingEnabled(true); + auto req = linEqSolver->getRequirements(env, storm::solver::LinearEquationSolverTask::SolveEquations); + if (lowerBound) { + linEqSolver->setLowerBound(lowerBound.get()); + req.clearLowerBounds(); + } + if (upperBound) { + linEqSolver->setUpperBound(upperBound.get()); + req.clearUpperBounds(); + } + STORM_LOG_THROW(req.empty(), storm::exceptions::UncheckedRequirementException, "At least one requirement was not checked."); + } + + // Prepare the right hand side of the equation system + b.assign(epochModel.epochMatrix.getRowCount(), storm::utility::zero()); + std::vector const& objectiveValues = epochModel.objectiveRewards.front(); + for (auto const& choice : epochModel.objectiveRewardFilter.front()) { + b[choice] = objectiveValues[choice]; + } + auto stepSolutionIt = epochModel.stepSolutions.begin(); + for (auto const& choice : epochModel.stepChoices) { + b[choice] += *stepSolutionIt; + ++stepSolutionIt; + } + assert(stepSolutionIt == epochModel.stepSolutions.end()); + + // Solve the minMax equation system + linEqSolver->solveEquations(env, x, b); + + return storm::utility::vector::filterVector(x, epochModel.epochInStates); + } + + template<> + std::map SparseDtmcPrctlHelper::computeRewardBoundedValues(Environment const& env, storm::models::sparse::Dtmc const& model, std::shared_ptr rewardBoundedFormula, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "The specified property is not supported by this value type."); + return std::map(); + } + template + std::map SparseDtmcPrctlHelper::computeRewardBoundedValues(Environment const& env, storm::models::sparse::Dtmc const& model, std::shared_ptr rewardBoundedFormula, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory) { + storm::utility::Stopwatch swAll(true), swBuild, swCheck; + + storm::modelchecker::helper::rewardbounded::MultiDimensionalRewardUnfolding rewardUnfolding(model, rewardBoundedFormula); + + // Get lower and upper bounds for the solution. + auto lowerBound = rewardUnfolding.getLowerObjectiveBound(); + auto upperBound = rewardUnfolding.getUpperObjectiveBound(); + + // Initialize epoch models + auto initEpoch = rewardUnfolding.getStartEpoch(); + auto epochOrder = rewardUnfolding.getEpochComputationOrder(initEpoch); + + // initialize data that will be needed for each epoch + std::vector x, b; + std::unique_ptr> linEqSolver; + Environment preciseEnv = env; + ValueType precision = rewardUnfolding.getRequiredEpochModelPrecision(initEpoch, storm::utility::convertNumber(storm::settings::getModule().getPrecision())); + preciseEnv.solver().setLinearEquationSolverPrecision(storm::utility::convertNumber(precision)); + + // In case of cdf export we store the necessary data. + std::vector> cdfData; + + // Set the correct equation problem format + rewardUnfolding.setEquationSystemFormatForEpochModel(linearEquationSolverFactory.getEquationProblemFormat(preciseEnv)); + bool convertToEquationSystem = linearEquationSolverFactory.getEquationProblemFormat(preciseEnv) == solver::LinearEquationSolverProblemFormat::EquationSystem; + + storm::utility::ProgressMeasurement progress("epochs"); + progress.setMaxCount(epochOrder.size()); + progress.startNewMeasurement(0); + uint64_t numCheckedEpochs = 0; + for (auto const& epoch : epochOrder) { + swBuild.start(); + auto& epochModel = rewardUnfolding.setCurrentEpoch(epoch); + swBuild.stop(); swCheck.start(); + // If the epoch matrix is empty we do not need to solve a linear equation system + if ((convertToEquationSystem && epochModel.epochMatrix.isIdentityMatrix()) || (!convertToEquationSystem && epochModel.epochMatrix.getEntryCount() == 0)) { + rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialDtmcEpochModel(epochModel)); + } else { + rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialDtmcEpochModel(preciseEnv, epochModel, x, b, linEqSolver, linearEquationSolverFactory, lowerBound, upperBound)); + } + swCheck.stop(); + if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { + std::vector cdfEntry; + for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) { + uint64_t offset = rewardUnfolding.getDimension(i).isUpperBounded ? 0 : 1; + cdfEntry.push_back(storm::utility::convertNumber(rewardUnfolding.getEpochManager().getDimensionOfEpoch(epoch, i) + offset) * rewardUnfolding.getDimension(i).scalingFactor); + } + cdfEntry.push_back(rewardUnfolding.getInitialStateResult(epoch)); + cdfData.push_back(std::move(cdfEntry)); + } + ++numCheckedEpochs; + progress.updateProgress(numCheckedEpochs); + } + + std::map result; + for (auto const& initState : model.getInitialStates()) { + result[initState] = rewardUnfolding.getInitialStateResult(initEpoch, initState); + } + + swAll.stop(); + + if (storm::settings::getModule().isExportCdfSet()) { + std::vector headers; + for (uint64_t i = 0; i < rewardUnfolding.getEpochManager().getDimensionCount(); ++i) { + headers.push_back(rewardUnfolding.getDimension(i).formula->toString()); + } + headers.push_back("Result"); + storm::utility::exportDataToCSVFile(storm::settings::getModule().getExportCdfDirectory() + "cdf.csv", cdfData, headers); + } + + if (storm::settings::getModule().isShowStatisticsSet()) { + STORM_PRINT_AND_LOG("---------------------------------" << std::endl); + STORM_PRINT_AND_LOG("Statistics:" << std::endl); + STORM_PRINT_AND_LOG("---------------------------------" << std::endl); + STORM_PRINT_AND_LOG(" #checked epochs: " << epochOrder.size() << "." << std::endl); + STORM_PRINT_AND_LOG(" overall Time: " << swAll << "." << std::endl); + STORM_PRINT_AND_LOG("Epoch Model building Time: " << swBuild << "." << std::endl); + STORM_PRINT_AND_LOG("Epoch Model checking Time: " << swCheck << "." << std::endl); + STORM_PRINT_AND_LOG("---------------------------------" << std::endl); + } + + return result; + } + + template std::vector SparseDtmcPrctlHelper::computeUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, bool qualitative, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, ModelCheckerHint const& hint) { std::vector result(transitionMatrix.getRowCount(), storm::utility::zero()); diff --git a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h index 0c0712723..ead15de74 100644 --- a/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h +++ b/src/storm/modelchecker/prctl/helper/SparseDtmcPrctlHelper.h @@ -5,6 +5,7 @@ #include +#include "storm/models/sparse/Dtmc.h" #include "storm/models/sparse/StandardRewardModel.h" #include "storm/modelchecker/hints/ModelCheckerHint.h" @@ -26,7 +27,10 @@ namespace storm { template > class SparseDtmcPrctlHelper { public: - static std::vector computeBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint()); + + static std::vector computeStepBoundedUntilProbabilities(Environment const& env, storm::solver::SolveGoal&& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, storm::storage::BitVector const& phiStates, storm::storage::BitVector const& psiStates, uint_fast64_t stepBound, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory, ModelCheckerHint const& hint = ModelCheckerHint()); + + static std::map computeRewardBoundedValues(Environment const& env, storm::models::sparse::Dtmc const& model, std::shared_ptr rewardBoundedFormula, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); static std::vector computeNextProbabilities(Environment const& env, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::BitVector const& nextStates, storm::solver::LinearEquationSolverFactory const& linearEquationSolverFactory); diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 8c7b0dc7e..10311c2a2 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -87,7 +87,7 @@ namespace storm { } template - std::vector analyzeTrivialEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { + std::vector analyzeTrivialMdpEpochModel(OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel) { // Assert that the epoch model is indeed trivial assert(epochModel.epochMatrix.getEntryCount() == 0); @@ -137,7 +137,7 @@ namespace storm { } template - std::vector analyzeNonTrivialEpochModel(Environment const& env, OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, boost::optional const& lowerBound, boost::optional const& upperBound) { + std::vector analyzeNonTrivialMdpEpochModel(Environment const& env, OptimizationDirection dir, typename rewardbounded::MultiDimensionalRewardUnfolding::EpochModel& epochModel, std::vector& x, std::vector& b, std::unique_ptr>& minMaxSolver, storm::solver::MinMaxLinearEquationSolverFactory const& minMaxLinearEquationSolverFactory, boost::optional const& lowerBound, boost::optional const& upperBound) { // Update some data for the case that the Matrix has changed if (epochModel.epochMatrixChanged) { @@ -215,9 +215,9 @@ namespace storm { swBuild.stop(); swCheck.start(); // If the epoch matrix is empty we do not need to solve a linear equation system if (epochModel.epochMatrix.getEntryCount() == 0) { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialEpochModel(dir, epochModel)); + rewardUnfolding.setSolutionForCurrentEpoch(analyzeTrivialMdpEpochModel(dir, epochModel)); } else { - rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialEpochModel(preciseEnv, dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, lowerBound, upperBound)); + rewardUnfolding.setSolutionForCurrentEpoch(analyzeNonTrivialMdpEpochModel(preciseEnv, dir, epochModel, x, b, minMaxSolver, minMaxLinearEquationSolverFactory, lowerBound, upperBound)); } swCheck.stop(); if (storm::settings::getModule().isExportCdfSet() && !rewardUnfolding.getEpochManager().hasBottomDimension(epoch)) { diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp index ec599fd37..efc3c81f7 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.cpp @@ -10,6 +10,8 @@ #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" #include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/Dtmc.h" #include "storm/transformer/EndComponentEliminator.h" @@ -24,14 +26,13 @@ namespace storm { namespace rewardbounded { template - MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::vector> const& objectives) : model(model), objectives(objectives) { + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::vector> const& objectives) : model(model), objectives(objectives) { initialize(); } template - MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::shared_ptr objectiveFormula) : model(model) { + MultiDimensionalRewardUnfolding::MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula) : model(model) { - STORM_LOG_THROW(objectiveFormula->hasOptimalityType(), storm::exceptions::InvalidPropertyException, "Formula needs to specify whether minimal or maximal values are to be computed on nondeterministic model."); if (objectiveFormula->isProbabilityOperatorFormula()) { if (objectiveFormula->getSubformula().isMultiObjectiveFormula()) { for (auto const& subFormula : objectiveFormula->getSubformula().asMultiObjectiveFormula().getSubformulas()) { @@ -91,7 +92,7 @@ namespace storm { // lower bounded until formulas with non-trivial left hand side are excluded as this would require some additional effort (in particular the ProductModel::transformMemoryState method). STORM_LOG_THROW(dimension.isUpperBounded || subformula.getLeftSubformula(dim).isTrueFormula(), storm::exceptions::NotSupportedException, "Lower bounded until formulas are only supported by this method if the left subformula is 'true'. Got " << subformula << " instead."); if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { - dimensionWiseEpochSteps.push_back(std::vector(model.getNumberOfChoices(), 1)); + dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); dimension.scalingFactor = storm::utility::one(); } else { STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); @@ -114,7 +115,7 @@ namespace storm { dimension.objectiveIndex = objIndex; dimension.isUpperBounded = true; if (subformula.getTimeBoundReference(dim).isTimeBound() || subformula.getTimeBoundReference(dim).isStepBound()) { - dimensionWiseEpochSteps.push_back(std::vector(model.getNumberOfChoices(), 1)); + dimensionWiseEpochSteps.push_back(std::vector(model.getTransitionMatrix().getRowCount(), 1)); dimension.scalingFactor = storm::utility::one(); } else { STORM_LOG_ASSERT(subformula.getTimeBoundReference(dim).isRewardBound(), "Unexpected type of time bound."); @@ -166,8 +167,8 @@ namespace storm { epochManager = EpochManager(dimensions.size()); // Convert the epoch steps to a choice-wise representation - epochSteps.reserve(model.getNumberOfChoices()); - for (uint64_t choice = 0; choice < model.getNumberOfChoices(); ++choice) { + epochSteps.reserve(model.getTransitionMatrix().getRowCount()); + for (uint64_t choice = 0; choice < model.getTransitionMatrix().getRowCount(); ++choice) { Epoch step; uint64_t dim = 0; for (auto const& dimensionSteps : dimensionWiseEpochSteps) { @@ -396,7 +397,7 @@ namespace storm { // std::cout << "Setting epoch class for epoch " << epochManager.toString(epoch) << std::endl; auto productObjectiveRewards = productModel->computeObjectiveRewards(epochClass, objectives); - storm::storage::BitVector stepChoices(productModel->getProduct().getNumberOfChoices(), false); + storm::storage::BitVector stepChoices(productModel->getProduct().getTransitionMatrix().getRowCount(), false); uint64_t choice = 0; for (auto const& step : productModel->getSteps()) { if (!epochManager.isZeroEpoch(step) && epochManager.getSuccessorEpoch(epoch, step) != epoch) { @@ -421,7 +422,7 @@ namespace storm { } } - storm::storage::BitVector zeroObjRewardChoices(productModel->getProduct().getNumberOfChoices(), true); + storm::storage::BitVector zeroObjRewardChoices(productModel->getProduct().getTransitionMatrix().getRowCount(), true); for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { if (violatedLowerBoundedDimensions.isDisjointFrom(objectiveDimensions[objIndex])) { zeroObjRewardChoices &= storm::utility::vector::filterZero(productObjectiveRewards[objIndex]); @@ -433,15 +434,68 @@ namespace storm { storm::storage::BitVector productInStates = productModel->getInStates(epochClass); // The epoch model only needs to consider the states that are reachable from a relevant state storm::storage::BitVector consideredStates = storm::utility::graph::getReachableStates(epochModel.epochMatrix, productInStates, allProductStates, ~allProductStates); - // std::cout << "numInStates = " << productInStates.getNumberOfSetBits() << std::endl; - // std::cout << "numConsideredStates = " << consideredStates.getNumberOfSetBits() << std::endl; // We assume that there is no end component in which objective reward is earned STORM_LOG_ASSERT(!storm::utility::graph::checkIfECWithChoiceExists(epochModel.epochMatrix, epochModel.epochMatrix.transpose(true), allProductStates, ~zeroObjRewardChoices & ~stepChoices), "There is a scheduler that yields infinite reward for one objective. This case should be excluded"); - auto ecElimResult = storm::transformer::EndComponentEliminator::transform(epochModel.epochMatrix, consideredStates, zeroObjRewardChoices & ~stepChoices, consideredStates); - epochModel.epochMatrix = std::move(ecElimResult.matrix); - epochModelToProductChoiceMap = std::move(ecElimResult.newToOldRowMapping); + // Create the epoch model matrix + std::vector productToEpochModelStateMapping; + if (model.isOfType(storm::models::ModelType::Dtmc)) { + assert(zeroObjRewardChoices.size() == productModel->getProduct().getNumberOfStates()); + assert(stepChoices.size() == productModel->getProduct().getNumberOfStates()); + STORM_LOG_ASSERT(equationSolverProblemFormatForEpochModel.is_initialized(), "Linear equation problem format was not set."); + bool convertToEquationSystem = equationSolverProblemFormatForEpochModel.get() == storm::solver::LinearEquationSolverProblemFormat::EquationSystem; + // For DTMCs we consider the subsystem induced by the considered states. + // The transitions for states with zero reward are filtered out to guarantee a unique solution of the eq-system. + auto backwardTransitions = epochModel.epochMatrix.transpose(true); + storm::storage::BitVector nonZeroRewardStates = storm::utility::graph::performProbGreater0(backwardTransitions, consideredStates, consideredStates & (~zeroObjRewardChoices | stepChoices)); + // If there is at least one considered state with reward zero, we have to add a 'zero-reward-state' to the epoch model. + bool requiresZeroRewardState = nonZeroRewardStates != consideredStates; + uint64_t numEpochModelStates = nonZeroRewardStates.getNumberOfSetBits(); + uint64_t zeroRewardInState = numEpochModelStates; + if (requiresZeroRewardState) { + ++numEpochModelStates; + } + storm::storage::SparseMatrixBuilder builder; + if (!nonZeroRewardStates.empty()) { + builder = storm::storage::SparseMatrixBuilder(epochModel.epochMatrix.getSubmatrix(true, nonZeroRewardStates, nonZeroRewardStates, convertToEquationSystem)); + } + if (requiresZeroRewardState) { + if (convertToEquationSystem) { + // add a diagonal entry + builder.addNextValue(zeroRewardInState, zeroRewardInState, storm::utility::zero()); + } + epochModel.epochMatrix = builder.build(numEpochModelStates, numEpochModelStates); + } else { + assert (!nonZeroRewardStates.empty()); + epochModel.epochMatrix = builder.build(); + } + if (convertToEquationSystem) { + epochModel.epochMatrix.convertToEquationSystem(); + } + + epochModelToProductChoiceMap.clear(); + epochModelToProductChoiceMap.reserve(numEpochModelStates); + productToEpochModelStateMapping.assign(nonZeroRewardStates.size(), zeroRewardInState); + for (auto const& productState : nonZeroRewardStates) { + productToEpochModelStateMapping[productState] = epochModelToProductChoiceMap.size(); + epochModelToProductChoiceMap.push_back(productState); + } + if (requiresZeroRewardState) { + uint64_t zeroRewardProductState = (consideredStates & ~nonZeroRewardStates).getNextSetIndex(0); + assert(zeroRewardProductState < consideredStates.size()); + epochModelToProductChoiceMap.push_back(zeroRewardProductState); + } + } else if (model.isOfType(storm::models::ModelType::Mdp)) { + // Eliminate zero-reward end components + auto ecElimResult = storm::transformer::EndComponentEliminator::transform(epochModel.epochMatrix, consideredStates, zeroObjRewardChoices & ~stepChoices, consideredStates); + epochModel.epochMatrix = std::move(ecElimResult.matrix); + epochModelToProductChoiceMap = std::move(ecElimResult.newToOldRowMapping); + productToEpochModelStateMapping = std::move(ecElimResult.oldToNewStateMapping); + } else { + STORM_LOG_THROW(false, storm::exceptions::UnexpectedException, "Unsupported model type."); + } + epochModel.stepChoices = storm::storage::BitVector(epochModel.epochMatrix.getRowCount(), false); for (uint64_t choice = 0; choice < epochModel.epochMatrix.getRowCount(); ++choice) { if (stepChoices.get(epochModelToProductChoiceMap[choice])) { @@ -466,14 +520,14 @@ namespace storm { epochModel.epochInStates = storm::storage::BitVector(epochModel.epochMatrix.getRowGroupCount(), false); for (auto const& productState : productInStates) { - STORM_LOG_ASSERT(ecElimResult.oldToNewStateMapping[productState] < epochModel.epochMatrix.getRowGroupCount(), "Selected product state does not exist in the epoch model."); - epochModel.epochInStates.set(ecElimResult.oldToNewStateMapping[productState], true); + STORM_LOG_ASSERT(productToEpochModelStateMapping[productState] < epochModel.epochMatrix.getRowGroupCount(), "Selected product state does not exist in the epoch model."); + epochModel.epochInStates.set(productToEpochModelStateMapping[productState], true); } std::vector toEpochModelInStatesMap(productModel->getProduct().getNumberOfStates(), std::numeric_limits::max()); std::vector epochModelStateToInStateMap = epochModel.epochInStates.getNumberOfSetBitsBeforeIndices(); for (auto const& productState : productInStates) { - toEpochModelInStatesMap[productState] = epochModelStateToInStateMap[ecElimResult.oldToNewStateMapping[productState]]; + toEpochModelInStatesMap[productState] = epochModelStateToInStateMap[productToEpochModelStateMapping[productState]]; } productStateToEpochModelInStateMap = std::make_shared const>(std::move(toEpochModelInStatesMap)); @@ -487,6 +541,13 @@ namespace storm { } + template + void MultiDimensionalRewardUnfolding::setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat) { + STORM_LOG_ASSERT(model.isOfType(storm::models::ModelType::Dtmc), "Trying to set the equation problem format although the model is not deterministic."); + equationSolverProblemFormatForEpochModel = eqSysFormat; + } + + template template::type> typename MultiDimensionalRewardUnfolding::SolutionType MultiDimensionalRewardUnfolding::getScaledSolution(SolutionType const& solution, ValueType const& scalingFactor) const { diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h index 11a15a82a..e1b96c64c 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/MultiDimensionalRewardUnfolding.h @@ -8,7 +8,8 @@ #include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h" #include "storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h" #include "storm/modelchecker/prctl/helper/rewardbounded/Dimension.h" -#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/Model.h" +#include "storm/solver/LinearEquationSolverProblemFormat.h" #include "storm/utility/vector.h" #include "storm/storage/memorystructure/MemoryStructure.h" #include "storm/utility/Stopwatch.h" @@ -43,8 +44,8 @@ namespace storm { * @param objectives The (preprocessed) objectives * */ - MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::vector> const& objectives); - MultiDimensionalRewardUnfolding(storm::models::sparse::Mdp const& model, std::shared_ptr objectiveFormula); + MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::vector> const& objectives); + MultiDimensionalRewardUnfolding(storm::models::sparse::Model const& model, std::shared_ptr objectiveFormula); ~MultiDimensionalRewardUnfolding() = default; @@ -53,6 +54,8 @@ namespace storm { EpochModel& setCurrentEpoch(Epoch const& epoch); + void setEquationSystemFormatForEpochModel(storm::solver::LinearEquationSolverProblemFormat eqSysFormat); + /*! * Returns the precision required for the analyzis of each epoch model in order to achieve the given overall precision */ @@ -106,7 +109,7 @@ namespace storm { EpochSolution const& getEpochSolution(std::map const& solutions, Epoch const& epoch); SolutionType const& getStateSolution(EpochSolution const& epochSolution, uint64_t const& productState); - storm::models::sparse::Mdp const& model; + storm::models::sparse::Model const& model; std::vector> objectives; std::unique_ptr> productModel; @@ -118,6 +121,9 @@ namespace storm { EpochModel epochModel; boost::optional currentEpoch; + // In case of DTMCs we have different options for the equation problem format the epoch model will have. + boost::optional equationSolverProblemFormatForEpochModel; + EpochManager epochManager; std::vector> dimensions; diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp index fa327d34e..af41c8a4f 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.cpp @@ -9,6 +9,8 @@ #include "storm/modelchecker/propositional/SparsePropositionalModelChecker.h" #include "storm/modelchecker/results/ExplicitQualitativeCheckResult.h" +#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/Dtmc.h" #include "storm/exceptions/UnexpectedException.h" #include "storm/exceptions/NotSupportedException.h" @@ -19,7 +21,7 @@ namespace storm { namespace rewardbounded { template - ProductModel::ProductModel(storm::models::sparse::Mdp const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()) { + ProductModel::ProductModel(storm::models::sparse::Model const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps) : dimensions(dimensions), objectiveDimensions(objectiveDimensions), epochManager(epochManager), memoryStateManager(dimensions.size()) { for (uint64_t dim = 0; dim < dimensions.size(); ++dim) { if (!dimensions[dim].memoryLabel) { @@ -34,7 +36,7 @@ namespace storm { storm::storage::SparseModelMemoryProduct productBuilder(memory.product(model)); setReachableProductStates(productBuilder, originalModelSteps, memoryStateMap); - product = productBuilder.build()->template as>(); + product = productBuilder.build(); uint64_t numModelStates = productBuilder.getOriginalModel().getNumberOfStates(); MemoryState upperMemStateBound = memoryStateManager.getUpperMemoryStateBound(); @@ -57,7 +59,7 @@ namespace storm { } // Map choice indices of the product to the state where it origins - choiceToStateMap.reserve(getProduct().getNumberOfChoices()); + choiceToStateMap.reserve(getProduct().getTransitionMatrix().getRowCount()); for (uint64_t productState = 0; productState < numProductStates; ++productState) { uint64_t groupSize = getProduct().getTransitionMatrix().getRowGroupSize(productState); for (uint64_t i = 0; i < groupSize; ++i) { @@ -66,7 +68,7 @@ namespace storm { } // Compute the epoch steps for the product - steps.resize(getProduct().getNumberOfChoices(), 0); + steps.resize(getProduct().getTransitionMatrix().getRowCount(), 0); for (uint64_t modelState = 0; modelState < numModelStates; ++modelState) { uint64_t numChoices = productBuilder.getOriginalModel().getTransitionMatrix().getRowGroupSize(modelState); uint64_t firstChoice = productBuilder.getOriginalModel().getTransitionMatrix().getRowGroupIndices()[modelState]; @@ -92,9 +94,9 @@ namespace storm { template - storm::storage::MemoryStructure ProductModel::computeMemoryStructure(storm::models::sparse::Mdp const& model, std::vector> const& objectives) const { + storm::storage::MemoryStructure ProductModel::computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives) const { - storm::modelchecker::SparsePropositionalModelChecker> mc(model); + storm::modelchecker::SparsePropositionalModelChecker> mc(model); // Create a memory structure that remembers whether (sub)objectives are satisfied storm::storage::MemoryStructure memory = storm::storage::MemoryStructureBuilder::buildTrivialMemoryStructure(model); @@ -285,7 +287,7 @@ namespace storm { } template - storm::models::sparse::Mdp const& ProductModel::getProduct() const { + storm::models::sparse::Model const& ProductModel::getProduct() const { return *product; } @@ -341,7 +343,7 @@ namespace storm { for (uint64_t objIndex = 0; objIndex < objectives.size(); ++objIndex) { auto const& formula = *objectives[objIndex].formula; if (formula.isProbabilityOperatorFormula()) { - storm::modelchecker::SparsePropositionalModelChecker> mc(getProduct()); + storm::modelchecker::SparsePropositionalModelChecker> mc(getProduct()); std::vector dimensionIndexMap; for (auto const& globalDimensionIndex : objectiveDimensions[objIndex]) { dimensionIndexMap.push_back(globalDimensionIndex); diff --git a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h index 4aaf90796..3f5b2a9a7 100644 --- a/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h +++ b/src/storm/modelchecker/prctl/helper/rewardbounded/ProductModel.h @@ -7,7 +7,7 @@ #include "storm/modelchecker/prctl/helper/rewardbounded/EpochManager.h" #include "storm/modelchecker/prctl/helper/rewardbounded/MemoryStateManager.h" #include "storm/modelchecker/prctl/helper/rewardbounded/Dimension.h" -#include "storm/models/sparse/Mdp.h" +#include "storm/models/sparse/Model.h" #include "storm/utility/vector.h" #include "storm/storage/memorystructure/MemoryStructure.h" #include "storm/storage/memorystructure/SparseModelMemoryProduct.h" @@ -26,9 +26,9 @@ namespace storm { typedef typename MemoryStateManager::MemoryState MemoryState; - ProductModel(storm::models::sparse::Mdp const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps); + ProductModel(storm::models::sparse::Model const& model, std::vector> const& objectives, std::vector> const& dimensions, std::vector const& objectiveDimensions, EpochManager const& epochManager, std::vector const& originalModelSteps); - storm::models::sparse::Mdp const& getProduct() const; + storm::models::sparse::Model const& getProduct() const; std::vector const& getSteps() const; bool productStateExists(uint64_t const& modelState, uint64_t const& memoryState) const; @@ -50,7 +50,7 @@ namespace storm { private: - storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Mdp const& model, std::vector> const& objectives) const; + storm::storage::MemoryStructure computeMemoryStructure(storm::models::sparse::Model const& model, std::vector> const& objectives) const; std::vector computeMemoryStateMap(storm::storage::MemoryStructure const& memory) const; @@ -66,7 +66,7 @@ namespace storm { EpochManager const& epochManager; MemoryStateManager memoryStateManager; - std::shared_ptr> product; + std::shared_ptr> product; std::vector steps; std::map reachableStates; std::map inStates;