diff --git a/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp b/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp index 1f6eaa8a6..9241e2e87 100644 --- a/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp +++ b/src/storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.cpp @@ -19,6 +19,7 @@ namespace storm { template ValueType BaierUpperRewardBoundsComputer::computeUpperBound() { + STORM_LOG_TRACE("Computing upper reward bounds using variant-2 of Baier et al."); std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now(); std::vector stateToScc(transitionMatrix.getRowGroupCount()); @@ -34,7 +35,7 @@ namespace storm { ++sccIndex; } } - + // The states that we still need to assign a value. storm::storage::BitVector remainingStates(transitionMatrix.getRowGroupCount(), true); diff --git a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp index f3789d31b..cdc0b9611 100644 --- a/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp +++ b/src/storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.cpp @@ -23,6 +23,7 @@ namespace storm { template std::vector DsMpiDtmcUpperRewardBoundsComputer::computeUpperBounds() { + STORM_LOG_TRACE("Computing upper reward bounds using DS-MPI."); std::chrono::high_resolution_clock::time_point start = std::chrono::high_resolution_clock::now(); sweep(); ValueType lambda = computeLambda(); diff --git a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp index 4d6f3426b..c43e3e66a 100644 --- a/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/HybridMdpPrctlHelper.cpp @@ -14,6 +14,8 @@ #include "storm/models/symbolic/StandardRewardModel.h" #include "storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h" +#include "storm/modelchecker/prctl/helper/DsMpiUpperRewardBoundsComputer.h" +#include "storm/modelchecker/prctl/helper/BaierUpperRewardBoundsComputer.h" #include "storm/modelchecker/results/SymbolicQualitativeCheckResult.h" #include "storm/modelchecker/results/SymbolicQuantitativeCheckResult.h" #include "storm/modelchecker/results/HybridQuantitativeCheckResult.h" @@ -32,6 +34,7 @@ namespace storm { boost::optional> ecInformation; boost::optional> initialScheduler; storm::storage::BitVector properMaybeStates; + boost::optional> oneStepTargetProbabilities; }; template @@ -66,15 +69,17 @@ namespace storm { } template - void eliminateExtendedStatesFromExplicitRepresentation(std::pair, std::vector>& explicitRepresentation, std::vector& scheduler, storm::storage::BitVector const& properMaybeStates) { - // Eliminate superfluous entries from the scheduler. - uint64_t position = 0; - for (auto state : properMaybeStates) { - scheduler[position] = scheduler[state]; - position++; + void eliminateExtendedStatesFromExplicitRepresentation(std::pair, std::vector>& explicitRepresentation, boost::optional>& scheduler, storm::storage::BitVector const& properMaybeStates) { + if (scheduler) { + // Eliminate superfluous entries from the scheduler. + uint64_t position = 0; + for (auto state : properMaybeStates) { + scheduler.get()[position] = scheduler.get()[state]; + position++; + } + scheduler.get().resize(properMaybeStates.getNumberOfSetBits()); + scheduler.get().shrink_to_fit(); } - scheduler.resize(properMaybeStates.getNumberOfSetBits()); - scheduler.shrink_to_fit(); // Treat the matrix. explicitRepresentation.first = explicitRepresentation.first.getSubmatrix(true, properMaybeStates, properMaybeStates); @@ -103,7 +108,7 @@ namespace storm { oneStepProbabilities = std::move(subvector); } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); - eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates); oneStepProbabilities = explicitRepresentation.first.getConstrainedRowGroupSumVector(solverRequirementsData.properMaybeStates, targetStates); } } @@ -192,14 +197,11 @@ namespace storm { storm::dd::Add subvector = submatrix * prob1StatesAsColumn; subvector = subvector.sumAbstract(model.getColumnVariables()); - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. - std::vector rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd().sumAbstract(model.getNondeterminismVariables()).toVector(odd); - // Finally cut away all columns targeting non-maybe states. submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); // Translate the symbolic matrix/vector to their explicit representations and solve the equation system. - explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); + explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd); if (requirements.requiresValidInitialScheduler()) { solverRequirementsData.initialScheduler = computeValidInitialSchedulerForUntilProbabilities(explicitRepresentation.first, explicitRepresentation.second); @@ -280,9 +282,6 @@ namespace storm { storm::dd::Add prob1StatesAsColumn = psiStates.template toAdd().swapVariables(model.getRowColumnMetaVariablePairs()); storm::dd::Add subvector = (submatrix * prob1StatesAsColumn).sumAbstract(model.getColumnVariables()); - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. - std::vector rowGroupSizes = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd().sumAbstract(model.getNondeterminismVariables()).toVector(odd); - // Finally cut away all columns targeting non-maybe states. submatrix *= maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); @@ -290,7 +289,7 @@ namespace storm { std::vector x(maybeStates.getNonZeroCount(), storm::utility::zero()); // Translate the symbolic matrix/vector to their explicit representations. - std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); + std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd); std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); solver->repeatedMultiply(dir, x, &explicitRepresentation.second, stepBound); @@ -338,12 +337,8 @@ namespace storm { // Create the solution vector. std::vector x(model.getNumberOfStates(), storm::utility::zero()); - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. - storm::dd::Add stateActionAdd = (transitionMatrix.notZero().existsAbstract(model.getColumnVariables()) || totalRewardVector.notZero()).template toAdd(); - std::vector rowGroupSizes = stateActionAdd.sumAbstract(model.getNondeterminismVariables()).toVector(odd); - // Translate the symbolic matrix/vector to their explicit representations. - std::pair, std::vector> explicitRepresentation = transitionMatrix.toMatrixVector(totalRewardVector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); + std::pair, std::vector> explicitRepresentation = transitionMatrix.toMatrixVector(totalRewardVector, model.getNondeterminismVariables(), odd, odd); // Perform the matrix-vector multiplication. std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); @@ -386,7 +381,12 @@ namespace storm { } template - void eliminateEndComponentsAndTargetStatesReachabilityRewards(std::pair, std::vector>& explicitRepresentation, SolverRequirementsData& solverRequirementsData) { + std::vector computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(storm::storage::SparseMatrix const& extendedMatrix, storm::storage::BitVector const& properMaybeStates, storm::storage::BitVector const& targetStates) { + return extendedMatrix.getConstrainedRowGroupSumVector(properMaybeStates, targetStates); + } + + template + void eliminateEndComponentsAndTargetStatesReachabilityRewards(std::pair, std::vector>& explicitRepresentation, SolverRequirementsData& solverRequirementsData, storm::storage::BitVector const& targetStates, bool computeOneStepTargetProbabilities) { // Get easier handles to the data. auto& transitionMatrix = explicitRepresentation.first; @@ -431,12 +431,31 @@ namespace storm { // Only do more work if there are actually end-components. if (doDecomposition && !endComponentDecomposition.empty()) { STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " EC(s)."); - std::vector subvector; - SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, rewardVector, solverRequirementsData.properMaybeStates, transitionMatrix, subvector); + if (computeOneStepTargetProbabilities) { + solverRequirementsData.oneStepTargetProbabilities = std::vector(solverRequirementsData.properMaybeStates.getNumberOfSetBits(), storm::utility::zero()); + } + + std::vector subvector(solverRequirementsData.properMaybeStates.getNumberOfSetBits(), storm::utility::zero()); + SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, solverRequirementsData.properMaybeStates, computeOneStepTargetProbabilities ? &targetStates : nullptr, nullptr, &rewardVector, transitionMatrix, computeOneStepTargetProbabilities ? &solverRequirementsData.oneStepTargetProbabilities.get() : nullptr, &subvector); rewardVector = std::move(subvector); } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); - eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); + if (computeOneStepTargetProbabilities) { + solverRequirementsData.oneStepTargetProbabilities = computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates); + } + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates); + } + } + + template + void setUpperRewardBounds(storm::solver::MinMaxLinearEquationSolver& solver, storm::OptimizationDirection const& direction, storm::storage::SparseMatrix const& submatrix, std::vector const& choiceRewards, std::vector const& oneStepTargetProbabilities) { + // For the min-case, we use DS-MPI, for the max-case variant 2 of the Baier et al. paper (CAV'17). + if (direction == storm::OptimizationDirection::Minimize) { + DsMpiMdpUpperRewardBoundsComputer dsmpi(submatrix, choiceRewards, oneStepTargetProbabilities); + solver.setUpperBounds(dsmpi.computeUpperBounds()); + } else { + BaierUpperRewardBoundsComputer baier(submatrix, choiceRewards, oneStepTargetProbabilities); + solver.setUpperBound(baier.computeUpperBound()); } } @@ -484,6 +503,11 @@ namespace storm { clearedRequirements.clearValidInitialScheduler(); } clearedRequirements.clearLowerBounds(); + if (clearedRequirements.requiresUpperBounds()) { + STORM_LOG_DEBUG("Computing upper bounds, because the solver requires it."); + extendMaybeStates = true; + clearedRequirements.clearUpperBounds(); + } STORM_LOG_THROW(clearedRequirements.empty(), storm::exceptions::UncheckedRequirementException, "Cannot establish requirements for solver."); } @@ -500,50 +524,54 @@ namespace storm { // (a) transitions from non-maybe states, and // (b) the choices in the transition matrix that lead to a state that is neither a maybe state // nor a target state ('infinity choices'). - storm::dd::Add choiceFilterAdd = (transitionMatrixBdd && maybeStatesWithTargetStates.renameVariables(model.getRowVariables(), model.getColumnVariables())).existsAbstract(model.getColumnVariables()).template toAdd(); - storm::dd::Add submatrix = transitionMatrix * maybeStatesAdd * choiceFilterAdd; + storm::dd::Add choiceFilterAdd = maybeStatesAdd * (transitionMatrixBdd && maybeStatesWithTargetStates.renameVariables(model.getRowVariables(), model.getColumnVariables())).existsAbstract(model.getColumnVariables()).template toAdd(); + storm::dd::Add submatrix = transitionMatrix * choiceFilterAdd; // Then compute the reward vector to use in the computation. - storm::dd::Add subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, submatrix, model.getColumnVariables()); - if (!rewardModel.hasStateActionRewards() && !rewardModel.hasTransitionRewards()) { - // If the reward model neither has state-action nor transition rewards, we need to multiply - // it with the legal nondetermism encodings in each state. - subvector *= choiceFilterAdd; - } - - // Before cutting the non-maybe columns, we need to compute the sizes of the row groups. - storm::dd::Add stateActionAdd = submatrix.notZero().existsAbstract(model.getColumnVariables()).template toAdd(); - std::vector rowGroupSizes = stateActionAdd.sumAbstract(model.getNondeterminismVariables()).toVector(odd); + storm::dd::Add subvector = rewardModel.getTotalRewardVector(maybeStatesAdd, choiceFilterAdd, submatrix, model.getColumnVariables()); // Finally cut away all columns targeting non-maybe states (or non-(maybe or target) states, respectively). submatrix *= extendMaybeStates ? maybeStatesWithTargetStates.swapVariables(model.getRowColumnMetaVariablePairs()).template toAdd() : maybeStatesAdd.swapVariables(model.getRowColumnMetaVariablePairs()); // Translate the symbolic matrix/vector to their explicit representations. - std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, std::move(rowGroupSizes), model.getNondeterminismVariables(), odd, odd); + std::pair, std::vector> explicitRepresentation = submatrix.toMatrixVector(subvector, model.getNondeterminismVariables(), odd, odd); // Fulfill the solver's requirements. SolverRequirementsData solverRequirementsData; - if (requirements.requiresNoEndComponents() || requirements.requiresValidInitialScheduler()) { + if (extendMaybeStates) { storm::storage::BitVector targetStates = computeTargetStatesForReachabilityRewardsFromExplicitRepresentation(explicitRepresentation.first); solverRequirementsData.properMaybeStates = ~targetStates; if (requirements.requiresNoEndComponents()) { - eliminateEndComponentsAndTargetStatesReachabilityRewards(explicitRepresentation, solverRequirementsData); + eliminateEndComponentsAndTargetStatesReachabilityRewards(explicitRepresentation, solverRequirementsData, targetStates, requirements.requiresUpperBounds()); } else { - // Compute a valid initial scheduler. - solverRequirementsData.initialScheduler = computeValidInitialSchedulerForReachabilityRewards(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates); + if (requirements.requiresValidInitialScheduler()) { + // Compute a valid initial scheduler. + solverRequirementsData.initialScheduler = computeValidInitialSchedulerForReachabilityRewards(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates); + } + + if (requirements.requiresUpperBounds()) { + solverRequirementsData.oneStepTargetProbabilities = computeOneStepTargetProbabilitiesFromExtendedExplicitRepresentation(explicitRepresentation.first, solverRequirementsData.properMaybeStates, targetStates); + } // Since we needed the transitions to target states to be translated as well for the computation // of the scheduler, we have to get rid of them now. - eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler.get(), solverRequirementsData.properMaybeStates); + eliminateExtendedStatesFromExplicitRepresentation(explicitRepresentation, solverRequirementsData.initialScheduler, solverRequirementsData.properMaybeStates); } } - + // Create the solution vector. std::vector x(explicitRepresentation.first.getRowGroupCount(), storm::utility::zero()); // Now solve the resulting equation system. - std::unique_ptr> solver = linearEquationSolverFactory.create(std::move(explicitRepresentation.first)); + std::unique_ptr> solver = linearEquationSolverFactory.create(); + + // If the solver requires upper bounds, compute them now. + if (requirements.requiresUpperBounds()) { + setUpperRewardBounds(*solver, dir, explicitRepresentation.first, explicitRepresentation.second, solverRequirementsData.oneStepTargetProbabilities.get()); + } + + solver->setMatrix(std::move(explicitRepresentation.first)); // Move the scheduler to the solver. if (solverRequirementsData.initialScheduler) { diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp index f9d46ae94..62dba3778 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.cpp @@ -30,6 +30,19 @@ namespace storm { // (3) Compute number of states not in MECs. numberOfMaybeStatesNotInEc = maybeStateToEc.size() - numberOfMaybeStatesInEc; + + // (4) Compute the number of maybe states that are not in ECs before every maybe state. + maybeStatesNotInEcBefore = std::vector(maybeStateToEc.size()); + + uint64_t count = 0; + auto resultIt = maybeStatesNotInEcBefore.begin(); + for (auto const& e : maybeStateToEc) { + *resultIt = count; + if (e == NOT_IN_EC) { + ++count; + } + ++resultIt; + } } template @@ -43,20 +56,8 @@ namespace storm { } template - std::vector SparseMdpEndComponentInformation::getNumberOfMaybeStatesNotInEcBeforeIndices() const { - std::vector result(maybeStateToEc.size()); - - uint64_t count = 0; - auto resultIt = result.begin(); - for (auto const& e : maybeStateToEc) { - *resultIt = count; - if (e != NOT_IN_EC) { - ++count; - } - ++resultIt; - } - - return result; + std::vector const& SparseMdpEndComponentInformation::getNumberOfMaybeStatesNotInEcBeforeIndices() const { + return maybeStatesNotInEcBefore; } template @@ -71,7 +72,11 @@ namespace storm { template uint64_t SparseMdpEndComponentInformation::getRowGroupAfterElimination(uint64_t state) const { - return numberOfMaybeStatesNotInEc + getEc(state); + if (this->isStateInEc(state)) { + return numberOfMaybeStatesNotInEc + getEc(state); + } else { + return maybeStatesNotInEcBefore[maybeStatesBefore[state]]; + } } template @@ -90,20 +95,22 @@ namespace storm { SparseMdpEndComponentInformation result(endComponentDecomposition, maybeStates); // (1) Compute the number of maybe states not in ECs before any other maybe state. - std::vector maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices(); + std::vector const& maybeStatesNotInEcBefore = result.getNumberOfMaybeStatesNotInEcBeforeIndices(); uint64_t numberOfStates = result.numberOfMaybeStatesNotInEc + result.numberOfEc; - STORM_LOG_TRACE("Found " << numberOfStates << " states, " << result.numberOfMaybeStatesNotInEc << " not in ECs, " << result.numberOfMaybeStatesInEc << " in ECs and " << result.numberOfEc << " ECs."); + STORM_LOG_TRACE("Creating " << numberOfStates << " states from " << result.numberOfMaybeStatesNotInEc << " states not in ECs and " << result.numberOfMaybeStatesInEc << " states in " << result.numberOfEc << " ECs."); // Prepare builder and vector storage. storm::storage::SparseMatrixBuilder builder(0, numberOfStates, 0, true, true, numberOfStates); STORM_LOG_ASSERT((sumColumns && columnSumVector) || (!sumColumns && !columnSumVector), "Expecting a bit vector for which columns to sum iff there is a column sum result vector."); if (columnSumVector) { - columnSumVector->resize(numberOfStates); + // Clearing column sum vector as we do not know the number of choices at this point. + columnSumVector->resize(0); } STORM_LOG_ASSERT((summand && summandResultVector) || (!summand && !summandResultVector), "Expecting summand iff there is a summand result vector."); if (summandResultVector) { - summandResultVector->resize(numberOfStates); + // Clearing summand result vector as we do not know the number of choices at this point. + summandResultVector->resize(0); } std::vector> ecValuePairs; @@ -121,14 +128,17 @@ namespace storm { ecValuePairs.clear(); if (summand) { - (*summandResultVector)[currentRow] += (*summand)[row]; + summandResultVector->emplace_back((*summand)[row]); + } + if (columnSumVector) { + columnSumVector->emplace_back(storm::utility::zero()); } for (auto const& e : transitionMatrix.getRow(row)) { if (sumColumns && sumColumns->get(e.getColumn())) { - (*columnSumVector)[currentRow] += e.getValue(); + columnSumVector->back() += e.getValue(); } else if (maybeStates.get(e.getColumn())) { // If the target state of the transition is not contained in an EC, we can just add the entry. - if (result.isStateInEc(e.getColumn())) { + if (!result.isStateInEc(e.getColumn())) { builder.addNextValue(currentRow, maybeStatesNotInEcBefore[result.maybeStatesBefore[e.getColumn()]], e.getValue()); } else { // Otherwise, we store the information that the state can go to a certain EC. @@ -170,14 +180,17 @@ namespace storm { ecValuePairs.clear(); if (summand) { - (*summandResultVector)[currentRow] += (*summand)[row]; + summandResultVector->emplace_back((*summand)[row]); + } + if (columnSumVector) { + columnSumVector->emplace_back(storm::utility::zero()); } for (auto const& e : transitionMatrix.getRow(row)) { if (sumColumns && sumColumns->get(e.getColumn())) { - (*columnSumVector)[currentRow] += e.getValue(); + columnSumVector->back() += e.getValue(); } else if (maybeStates.get(e.getColumn())) { // If the target state of the transition is not contained in an EC, we can just add the entry. - if (result.isStateInEc(e.getColumn())) { + if (!result.isStateInEc(e.getColumn())) { builder.addNextValue(currentRow, maybeStatesNotInEcBefore[result.maybeStatesBefore[e.getColumn()]], e.getValue()); } else { // Otherwise, we store the information that the state can go to a certain EC. diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h index a8a752305..513ffb49c 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h +++ b/src/storm/modelchecker/prctl/helper/SparseMdpEndComponentInformation.h @@ -29,7 +29,7 @@ namespace storm { * Retrieves for each maybe state the number of maybe states not contained in ECs with an index smaller * than the requested one. */ - std::vector getNumberOfMaybeStatesNotInEcBeforeIndices() const; + std::vector const& getNumberOfMaybeStatesNotInEcBeforeIndices() const; /*! * Retrieves the total number of maybe states on in ECs. @@ -65,6 +65,7 @@ namespace storm { // Data about end components. std::vector maybeStatesBefore; + std::vector maybeStatesNotInEcBefore; uint64_t numberOfMaybeStatesInEc; uint64_t numberOfMaybeStatesNotInEc; uint64_t numberOfEc; diff --git a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp index 37192633a..4a5688f5c 100644 --- a/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp +++ b/src/storm/modelchecker/prctl/helper/SparseMdpPrctlHelper.cpp @@ -449,7 +449,7 @@ namespace storm { } template - void computeFixedPointSystemUntilProbabilities(storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix& submatrix, std::vector& b) { + void computeFixedPointSystemUntilProbabilities(storm::solver::SolveGoal& goal, storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix& submatrix, std::vector& b) { // First, we can eliminate the rows and columns from the original transition probability matrix for states // whose probabilities are already known. submatrix = transitionMatrix.getSubmatrix(true, qualitativeStateSets.maybeStates, qualitativeStateSets.maybeStates, false); @@ -457,10 +457,13 @@ namespace storm { // Prepare the right-hand side of the equation system. For entry i this corresponds to // the accumulated probability of going from state i to some state that has probability 1. b = transitionMatrix.getConstrainedRowGroupSumVector(qualitativeStateSets.maybeStates, qualitativeStateSets.statesWithProbability1); + + // If the solve goal has relevant values, we need to adjust them. + goal.restrictRelevantValues(qualitativeStateSets.maybeStates); } template - boost::optional> computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix& submatrix, std::vector& b) { + boost::optional> computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(storm::solver::SolveGoal& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsUntilProbabilities const& qualitativeStateSets, storm::storage::SparseMatrix& submatrix, std::vector& b) { // Start by computing the states that are in MECs. storm::storage::MaximalEndComponentDecomposition endComponentDecomposition(transitionMatrix, backwardTransitions, qualitativeStateSets.maybeStates); @@ -468,10 +471,26 @@ namespace storm { // Only do more work if there are actually end-components. if (!endComponentDecomposition.empty()) { STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " EC(s)."); - return SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, &b, nullptr); + SparseMdpEndComponentInformation result = SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, &qualitativeStateSets.statesWithProbability1, nullptr, nullptr, submatrix, &b, nullptr); + + // If the solve goal has relevant values, we need to adjust them. + if (goal.hasRelevantValues()) { + storm::storage::BitVector newRelevantValues(submatrix.getRowGroupCount()); + for (auto state : goal.relevantValues()) { + if (qualitativeStateSets.maybeStates.get(state)) { + newRelevantValues.set(result.getRowGroupAfterElimination(state)); + } + } + if (!newRelevantValues.empty()) { + goal.setRelevantValues(std::move(newRelevantValues)); + } + } + + return result; } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); - computeFixedPointSystemUntilProbabilities(transitionMatrix, qualitativeStateSets, submatrix, b); + computeFixedPointSystemUntilProbabilities(goal, transitionMatrix, qualitativeStateSets, submatrix, b); + return boost::none; } } @@ -516,17 +535,16 @@ namespace storm { // If the hint information tells us that we have to eliminate MECs, we do so now. boost::optional> ecInformation; if (hintInformation.getEliminateEndComponents()) { - ecInformation = computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, submatrix, b); + ecInformation = computeFixedPointSystemUntilProbabilitiesEliminateEndComponents(goal, transitionMatrix, backwardTransitions, qualitativeStateSets, submatrix, b); // Make sure we are not supposed to produce a scheduler if we actually eliminate end components. STORM_LOG_THROW(!ecInformation || !ecInformation.get().getEliminatedEndComponents() || !produceScheduler, storm::exceptions::NotSupportedException, "Producing schedulers is not supported if end-components need to be eliminated for the solver."); } else { // Otherwise, we compute the standard equations. - computeFixedPointSystemUntilProbabilities(transitionMatrix, qualitativeStateSets, submatrix, b); + computeFixedPointSystemUntilProbabilities(goal, transitionMatrix, qualitativeStateSets, submatrix, b); } // Now compute the results for the maybe states. - goal.restrictRelevantValues(qualitativeStateSets.maybeStates); MaybeStateResult resultForMaybeStates = computeValuesForMaybeStates(std::move(goal), std::move(submatrix), b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation); // If we eliminated end components, we need to extract the result differently. @@ -725,7 +743,7 @@ namespace storm { } template - void computeFixedPointSystemReachabilityRewards(storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, std::vector* oneStepTargetProbabilities = nullptr) { + void computeFixedPointSystemReachabilityRewards(storm::solver::SolveGoal& goal, storm::storage::SparseMatrix const& transitionMatrix, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, std::vector* oneStepTargetProbabilities = nullptr) { // Remove rows and columns from the original transition probability matrix for states whose reward values are already known. // If there are infinity states, we additionally have to remove choices of maybeState that lead to infinity. if (qualitativeStateSets.infinityStates.empty()) { @@ -742,10 +760,13 @@ namespace storm { (*oneStepTargetProbabilities) = transitionMatrix.getConstrainedRowSumVector(*selectedChoices, targetStates); } } + + // If the solve goal has relevant values, we need to adjust them. + goal.restrictRelevantValues(qualitativeStateSets.maybeStates); } template - boost::optional> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, boost::optional>& oneStepTargetProbabilities) { + boost::optional> computeFixedPointSystemReachabilityRewardsEliminateEndComponents(storm::solver::SolveGoal& goal, storm::storage::SparseMatrix const& transitionMatrix, storm::storage::SparseMatrix const& backwardTransitions, QualitativeStateSetsReachabilityRewards const& qualitativeStateSets, storm::storage::BitVector const& targetStates, boost::optional const& selectedChoices, std::function(uint_fast64_t, storm::storage::SparseMatrix const&, storm::storage::BitVector const&)> const& totalStateRewardVectorGetter, storm::storage::SparseMatrix& submatrix, std::vector& b, boost::optional>& oneStepTargetProbabilities) { // Start by computing the choices with reward 0, as we only want ECs within this fragment. storm::storage::BitVector zeroRewardChoices(transitionMatrix.getRowCount()); @@ -789,10 +810,25 @@ namespace storm { // Only do more work if there are actually end-components. if (doDecomposition && !endComponentDecomposition.empty()) { STORM_LOG_DEBUG("Eliminating " << endComponentDecomposition.size() << " ECs."); - return SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, oneStepTargetProbabilities ? &targetStates : nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr, &b); + SparseMdpEndComponentInformation result = SparseMdpEndComponentInformation::eliminateEndComponents(endComponentDecomposition, transitionMatrix, qualitativeStateSets.maybeStates, oneStepTargetProbabilities ? &targetStates : nullptr, selectedChoices ? &selectedChoices.get() : nullptr, &rewardVector, submatrix, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr, &b); + + // If the solve goal has relevant values, we need to adjust them. + if (goal.hasRelevantValues()) { + storm::storage::BitVector newRelevantValues(submatrix.getRowGroupCount()); + for (auto state : goal.relevantValues()) { + if (qualitativeStateSets.maybeStates.get(state)) { + newRelevantValues.set(result.getRowGroupAfterElimination(state)); + } + } + if (!newRelevantValues.empty()) { + goal.setRelevantValues(std::move(newRelevantValues)); + } + } + + return result; } else { STORM_LOG_DEBUG("Not eliminating ECs as there are none."); - computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr); + computeFixedPointSystemReachabilityRewards(goal, transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr); return boost::none; } } @@ -862,13 +898,13 @@ namespace storm { // If the hint information tells us that we have to eliminate MECs, we do so now. boost::optional> ecInformation; if (hintInformation.getEliminateEndComponents()) { - ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(transitionMatrix, backwardTransitions, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities); + ecInformation = computeFixedPointSystemReachabilityRewardsEliminateEndComponents(goal, transitionMatrix, backwardTransitions, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities); // Make sure we are not supposed to produce a scheduler if we actually eliminate end components. STORM_LOG_THROW(!ecInformation || !ecInformation.get().getEliminatedEndComponents() || !produceScheduler, storm::exceptions::NotSupportedException, "Producing schedulers is not supported if end-components need to be eliminated for the solver."); } else { // Otherwise, we compute the standard equations. - computeFixedPointSystemReachabilityRewards(transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr); + computeFixedPointSystemReachabilityRewards(goal, transitionMatrix, qualitativeStateSets, targetStates, selectedChoices, totalStateRewardVectorGetter, submatrix, b, oneStepTargetProbabilities ? &oneStepTargetProbabilities.get() : nullptr); } // If we need to compute upper bounds, do so now. @@ -878,7 +914,6 @@ namespace storm { } // Now compute the results for the maybe states. - goal.restrictRelevantValues(qualitativeStateSets.maybeStates); MaybeStateResult resultForMaybeStates = computeValuesForMaybeStates(std::move(goal), std::move(submatrix), b, produceScheduler, minMaxLinearEquationSolverFactory, hintInformation); // If we eliminated end components, we need to extract the result differently. diff --git a/src/storm/models/symbolic/StandardRewardModel.cpp b/src/storm/models/symbolic/StandardRewardModel.cpp index 8fe296308..ffd94ad74 100644 --- a/src/storm/models/symbolic/StandardRewardModel.cpp +++ b/src/storm/models/symbolic/StandardRewardModel.cpp @@ -101,6 +101,21 @@ namespace storm { return result; } + template + storm::dd::Add StandardRewardModel::getTotalRewardVector(storm::dd::Add const& stateFilterAdd, storm::dd::Add const& choiceFilterAdd, storm::dd::Add const& transitionMatrix, std::set const& columnVariables) const { + storm::dd::Add result = transitionMatrix.getDdManager().template getAddZero(); + if (this->hasStateRewards()) { + result += stateFilterAdd * choiceFilterAdd * optionalStateRewardVector.get(); + } + if (this->hasStateActionRewards()) { + result += choiceFilterAdd * optionalStateActionRewardVector.get(); + } + if (this->hasTransitionRewards()) { + result += (transitionMatrix * this->getTransitionRewardMatrix()).sumAbstract(columnVariables); + } + return result; + } + template storm::dd::Add StandardRewardModel::getTotalRewardVector(storm::dd::Add const& transitionMatrix, std::set const& columnVariables) const { storm::dd::Add result = transitionMatrix.getDdManager().template getAddZero(); diff --git a/src/storm/models/symbolic/StandardRewardModel.h b/src/storm/models/symbolic/StandardRewardModel.h index a7ab321f4..c16f144e5 100644 --- a/src/storm/models/symbolic/StandardRewardModel.h +++ b/src/storm/models/symbolic/StandardRewardModel.h @@ -161,6 +161,19 @@ namespace storm { */ storm::dd::Add getTotalRewardVector(storm::dd::Add const& filterAdd, storm::dd::Add const& transitionMatrix, std::set const& columnVariables) const; + /*! + * Creates a vector representing the complete reward vector based on the state-, state-action- and + * transition-based rewards in the reward model. The state- and state-action rewards are filtered with + * the given filter DD. + * + * @param stateFilterAdd The DD used for filtering the state rewards. + * @param choiceFilterAdd The DD used for filtering the state action rewards. + * @param transitionMatrix The matrix that is used to weight the values of the transition reward matrix. + * @param columnVariables The column variables of the model. + * @return The full state-action reward vector. + */ + storm::dd::Add getTotalRewardVector(storm::dd::Add const& stateFilterAdd, storm::dd::Add const& choiceFilterAdd, storm::dd::Add const& transitionMatrix, std::set const& columnVariables) const; + /*! * Creates a vector representing the complete reward vector based on the state-, state-action- and * transition-based rewards in the reward model. diff --git a/src/storm/solver/GmmxxMultiplier.cpp b/src/storm/solver/GmmxxMultiplier.cpp index 706ad795e..b05ef8bcc 100644 --- a/src/storm/solver/GmmxxMultiplier.cpp +++ b/src/storm/solver/GmmxxMultiplier.cpp @@ -80,11 +80,19 @@ namespace storm { uint64_t choice; for (auto row_group_it = rowGroupIndices.end() - 2, row_group_ite = rowGroupIndices.begin() - 1; row_group_it != row_group_ite; --row_group_it, --choice_it, --target_it) { - T currentValue = b ? *add_it : storm::utility::zero(); - --add_it; + if (choices) { + *choice_it = 0; + } + + T currentValue = storm::utility::zero(); // Only multiply and reduce if the row group is not empty. if (*row_group_it != *(row_group_it + 1)) { + if (b) { + currentValue = *add_it; + --add_it; + } + currentValue += vect_sp(gmm::linalg_traits::row(itr), x); if (choices) { @@ -94,7 +102,7 @@ namespace storm { --itr; - for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr) { + for (uint64_t row = *row_group_it + 1, rowEnd = *(row_group_it + 1); row < rowEnd; ++row, --itr, --add_it) { T newValue = b ? *add_it : storm::utility::zero(); newValue += vect_sp(gmm::linalg_traits::row(itr), x); @@ -108,9 +116,6 @@ namespace storm { *choice_it = choice; } } - if (b) { - --add_it; - } } } else if (choices) { *choice_it = 0; @@ -165,14 +170,19 @@ namespace storm { auto resultIt = result.begin() + range.begin(); for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { - T currentValue = b ? *bIt : storm::utility::zero(); - ++bIt; if (choices) { *choiceIt = 0; } + T currentValue = storm::utility::zero(); + // Only multiply and reduce if the row group is not empty. if (*groupIt != *(groupIt + 1)) { + if (b) { + currentValue = *bIt; + ++bIt; + } + ++itr; for (auto itre = mat_row_const_begin(matrix) + *(groupIt + 1); itr != itre; ++itr) { diff --git a/src/storm/solver/SolveGoal.cpp b/src/storm/solver/SolveGoal.cpp index cd331b2fe..8913e1807 100644 --- a/src/storm/solver/SolveGoal.cpp +++ b/src/storm/solver/SolveGoal.cpp @@ -116,6 +116,11 @@ namespace storm { relevantValueVector = relevantValueVector.get() % filter; } } + + template + void SolveGoal::setRelevantValues(storm::storage::BitVector&& values) { + relevantValueVector = std::move(values); + } template class SolveGoal; diff --git a/src/storm/solver/SolveGoal.h b/src/storm/solver/SolveGoal.h index 3d0cab095..e8093b073 100644 --- a/src/storm/solver/SolveGoal.h +++ b/src/storm/solver/SolveGoal.h @@ -82,6 +82,7 @@ namespace storm { storm::storage::BitVector const& relevantValues() const; void restrictRelevantValues(storm::storage::BitVector const& filter); + void setRelevantValues(storm::storage::BitVector&& values); private: boost::optional optimizationDirection; diff --git a/src/storm/solver/StandardGameSolver.cpp b/src/storm/solver/StandardGameSolver.cpp index f522b3ffa..12f17f345 100644 --- a/src/storm/solver/StandardGameSolver.cpp +++ b/src/storm/solver/StandardGameSolver.cpp @@ -222,7 +222,7 @@ namespace storm { // Proceed with the iterations as long as the method did not converge or reach the maximum number of iterations. uint64_t iterations = 0; - + Status status = Status::InProgress; while (status == Status::InProgress) { multiplyAndReduce(player1Dir, player2Dir, *currentX, &b, *linEqSolverPlayer2Matrix, multiplyResult, reducedMultiplyResult, *newX); @@ -291,7 +291,7 @@ namespace storm { void StandardGameSolver::multiplyAndReduce(OptimizationDirection player1Dir, OptimizationDirection player2Dir, std::vector& x, std::vector const* b, storm::solver::LinearEquationSolver const& linEqSolver, std::vector& multiplyResult, std::vector& p2ReducedMultiplyResult, std::vector& p1ReducedMultiplyResult) const { linEqSolver.multiply(x, b, multiplyResult); - + storm::utility::vector::reduceVectorMinOrMax(player2Dir, multiplyResult, p2ReducedMultiplyResult, player2Matrix.getRowGroupIndices()); uint_fast64_t pl1State = 0; diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 928e190a4..9588b415c 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -1565,20 +1565,21 @@ namespace storm { } for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) { - ValueType currentValue = summand ? *summandIt : storm::utility::zero(); - ++summandIt; + ValueType currentValue = storm::utility::zero(); if (choices) { *choiceIt = 0; } // Only multiply and reduce if there is at least one row in the group. if (*rowGroupIt < *(rowGroupIt + 1)) { + if (summand) { + currentValue = *summandIt; + ++summandIt; + } + for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { currentValue += elementIt->getValue() * vector[elementIt->getColumn()]; } - if (choices) { - *choiceIt = 0; - } ++rowIt; @@ -1598,6 +1599,8 @@ namespace storm { ++summandIt; } } + } else if (choices) { + *choiceIt = 0; } // Finally write value to target vector. @@ -1627,21 +1630,27 @@ namespace storm { } for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt) { - ValueType currentValue = summand ? *summandIt : storm::utility::zero(); - --summandIt; - + if (choices) { + *choiceIt = 0; + } + ValueType currentValue = storm::utility::zero(); + // Only multiply and reduce if there is at least one row in the group. if (*rowGroupIt < *(rowGroupIt + 1)) { + if (summand) { + currentValue = *summandIt; + --summandIt; + } + for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) { currentValue += elementIt->getValue() * vector[elementIt->getColumn()]; } if (choices) { *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt; } - --rowIt; - for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, ++i) { + for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, ++i, --summandIt) { ValueType newValue = summand ? *summandIt : storm::utility::zero(); for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) { newValue += elementIt->getValue() * vector[elementIt->getColumn()]; @@ -1653,12 +1662,7 @@ namespace storm { *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt; } } - if (summand) { - --summandIt; - } } - } else if (choices) { - *choiceIt = 0; } // Finally write value to target vector. @@ -1703,14 +1707,19 @@ namespace storm { auto resultIt = result.begin() + range.begin(); for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) { - ValueType currentValue = summand ? *summandIt : storm::utility::zero(); - ++summandIt; if (choices) { *choiceIt = 0; } + ValueType currentValue = storm::utility::zero(); + // Only multiply and reduce if there is at least one row in the group. if (*groupIt < *(groupIt + 1)) { + if (summand) { + currentValue = *summandIt; + ++summandIt; + } + for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) { currentValue += elementIt->getValue() * x[elementIt->getColumn()]; } diff --git a/src/storm/storage/dd/Add.cpp b/src/storm/storage/dd/Add.cpp index f01680610..a7129ced8 100644 --- a/src/storm/storage/dd/Add.cpp +++ b/src/storm/storage/dd/Add.cpp @@ -472,6 +472,67 @@ namespace storm { internalAdd.composeWithExplicitVector(rowOdd, ddVariableIndices, result, std::plus()); return result; } + + template + std::vector Add::toVector(storm::dd::Add const& matrix, std::vector const& rowGroupIndices, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const { + std::vector ddRowVariableIndices; + std::vector ddColumnVariableIndices; + std::vector ddGroupVariableIndices; + std::set rowAndColumnMetaVariables; + + for (auto const& variable : rowMetaVariables) { + DdMetaVariable const& metaVariable = this->getDdManager().getMetaVariable(variable); + for (auto const& ddVariable : metaVariable.getDdVariables()) { + ddRowVariableIndices.push_back(ddVariable.getIndex()); + } + rowAndColumnMetaVariables.insert(variable); + } + std::sort(ddRowVariableIndices.begin(), ddRowVariableIndices.end()); + for (auto const& variable : columnMetaVariables) { + DdMetaVariable const& metaVariable = this->getDdManager().getMetaVariable(variable); + for (auto const& ddVariable : metaVariable.getDdVariables()) { + ddColumnVariableIndices.push_back(ddVariable.getIndex()); + } + rowAndColumnMetaVariables.insert(variable); + } + std::sort(ddColumnVariableIndices.begin(), ddColumnVariableIndices.end()); + for (auto const& variable : groupMetaVariables) { + DdMetaVariable const& metaVariable = this->getDdManager().getMetaVariable(variable); + for (auto const& ddVariable : metaVariable.getDdVariables()) { + ddGroupVariableIndices.push_back(ddVariable.getIndex()); + } + } + std::sort(ddGroupVariableIndices.begin(), ddGroupVariableIndices.end()); + + Bdd columnVariableCube = Bdd::getCube(this->getDdManager(), columnMetaVariables); + + // Copy the row group indices so we can modify them. + std::vector mutableRowGroupIndices = rowGroupIndices; + + // Create the explicit vector we need to fill later. + std::vector explicitVector(mutableRowGroupIndices.back()); + + // Next, we split the matrix into one for each group. Note that this only works if the group variables are at the very top. + std::vector, InternalAdd>> internalAddGroups = matrix.internalAdd.splitIntoGroups(*this, ddGroupVariableIndices); + std::vector, Add>> groups; + for (auto const& internalAdd : internalAddGroups) { + groups.push_back(std::make_pair(Add(this->getDdManager(), internalAdd.first, rowAndColumnMetaVariables), Add(this->getDdManager(), internalAdd.second, rowMetaVariables))); + } + + std::vector> statesWithGroupEnabled(groups.size()); + for (uint_fast64_t i = 0; i < groups.size(); ++i) { + std::pair, Add> const& ddPair = groups[i]; + Bdd matrixDdNotZero = ddPair.first.notZero(); + Bdd vectorDdNotZero = ddPair.second.notZero(); + + ddPair.second.internalAdd.composeWithExplicitVector(rowOdd, ddRowVariableIndices, mutableRowGroupIndices, explicitVector, std::plus()); + + InternalAdd statesWithGroupEnabled = (matrixDdNotZero.existsAbstract(columnMetaVariables) || vectorDdNotZero).template toAdd(); + statesWithGroupEnabled.composeWithExplicitVector(rowOdd, ddRowVariableIndices, mutableRowGroupIndices, std::plus()); + } + + return explicitVector; + } template storm::storage::SparseMatrix Add::toMatrix() const { @@ -697,7 +758,7 @@ namespace storm { } template - std::pair, std::vector> Add::toMatrixVector(storm::dd::Add const& vector, std::vector&& rowGroupSizes, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const { + std::pair, std::vector> Add::toMatrixVector(storm::dd::Add const& vector, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const { std::set rowMetaVariables; std::set columnMetaVariables; @@ -715,11 +776,11 @@ namespace storm { } // Create the canonical row group sizes and build the matrix. - return toMatrixVector(vector, std::move(rowGroupSizes), rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd); + return toMatrixVector(vector, rowMetaVariables, columnMetaVariables, groupMetaVariables, rowOdd, columnOdd); } template - std::pair, std::vector> Add::toMatrixVector(storm::dd::Add const& vector, std::vector&& rowGroupIndices, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const { + std::pair, std::vector> Add::toMatrixVector(storm::dd::Add const& vector, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const { std::vector ddRowVariableIndices; std::vector ddColumnVariableIndices; std::vector ddGroupVariableIndices; @@ -751,6 +812,9 @@ namespace storm { Bdd columnVariableCube = Bdd::getCube(this->getDdManager(), columnMetaVariables); + // Count how many choices each row group has. + std::vector rowGroupIndices = (this->notZero().existsAbstract(columnMetaVariables) || vector.notZero()).template toAdd().sumAbstract(groupMetaVariables).toVector(rowOdd); + // Transform the row group sizes to the actual row group indices. rowGroupIndices.resize(rowGroupIndices.size() + 1); uint_fast64_t tmp = 0; diff --git a/src/storm/storage/dd/Add.h b/src/storm/storage/dd/Add.h index ebd0cb9b3..bbcfb099b 100644 --- a/src/storm/storage/dd/Add.h +++ b/src/storm/storage/dd/Add.h @@ -556,7 +556,26 @@ namespace storm { std::vector toVector(storm::dd::Odd const& rowOdd) const; /*! - * Converts the ADD to a (sparse) double matrix. All contained non-primed variables are assumed to encode the + * Converts the ADD to a row-grouped vector while respecting the row group sizes of the provided matrix. + * That is, if the vector has a zero entry for some row in a row group for which the matrix has a non-zero + * row, the value at the vector will be correctly set to zero. Note: this function assumes that the meta + * variables used to distinguish different row groups are at the very top of the ADD. + * + * @param matrix The symbolic matrix whose row group sizes to respect. + * @param rowGroupSizes A vector specifying the sizes of the row groups. + * @param rowMetaVariables The meta variables that encode the rows of the matrix. + * @param columnMetaVariables The meta variables that encode the columns of the matrix. + * @param groupMetaVariables The meta variables that are used to distinguish different row groups. + * @param rowOdd The ODD used for determining the correct row. + * @param columnOdd The ODD used for determining the correct column. + * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector + * (if it was given). + * @return The vector that is represented by this ADD. + */ + std::vector toVector(storm::dd::Add const& matrix, std::vector const& rowGroupSizes, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; + + /*! + * Converts the ADD to a (sparse) matrix. All contained non-primed variables are assumed to encode the * row, whereas all primed variables are assumed to encode the column. * * @return The matrix that is represented by this ADD. @@ -564,7 +583,7 @@ namespace storm { storm::storage::SparseMatrix toMatrix() const; /*! - * Converts the ADD to a (sparse) double matrix. All contained non-primed variables are assumed to encode the + * Converts the ADD to a (sparse) matrix. All contained non-primed variables are assumed to encode the * row, whereas all primed variables are assumed to encode the column. The given offset-labeled DDs are used * to determine the correct row and column, respectively, for each entry. * @@ -575,7 +594,7 @@ namespace storm { storm::storage::SparseMatrix toMatrix(storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; /*! - * Converts the ADD to a (sparse) double matrix. The given offset-labeled DDs are used to determine the + * Converts the ADD to a (sparse) matrix. The given offset-labeled DDs are used to determine the * correct row and column, respectively, for each entry. * * @param rowMetaVariables The meta variables that encode the rows of the matrix. @@ -587,7 +606,7 @@ namespace storm { storm::storage::SparseMatrix toMatrix(std::set const& rowMetaVariables, std::set const& columnMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; /*! - * Converts the ADD to a row-grouped (sparse) double matrix. The given offset-labeled DDs are used to + * Converts the ADD to a row-grouped (sparse) matrix. The given offset-labeled DDs are used to * determine the correct row and column, respectively, for each entry. Note: this function assumes that * the meta variables used to distinguish different row groups are at the very top of the ADD. * @@ -599,19 +618,18 @@ namespace storm { storm::storage::SparseMatrix toMatrix(std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; /*! - * Converts the ADD to a row-grouped (sparse) double matrix and the given vector to a row-grouped vector. + * Converts the ADD to a row-grouped (sparse) matrix and the given vector to a row-grouped vector. * The given offset-labeled DDs are used to determine the correct row and column, respectively, for each * entry. Note: this function assumes that the meta variables used to distinguish different row groups are * at the very top of the ADD. * * @param vector The symbolic vector to convert. - * @param rowGroupSizes A vector specifying the sizes of the row groups. * @param groupMetaVariables The meta variables that are used to distinguish different row groups. * @param rowOdd The ODD used for determining the correct row. * @param columnOdd The ODD used for determining the correct column. * @return The matrix that is represented by this ADD. */ - std::pair, std::vector> toMatrixVector(storm::dd::Add const& vector, std::vector&& rowGroupSizes, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; + std::pair, std::vector> toMatrixVector(storm::dd::Add const& vector, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; /*! * Exports the DD to the given file in the dot format. @@ -702,7 +720,6 @@ namespace storm { * different row groups are at the very top of the ADD. * * @param vector The vector that is to be transformed to an equally grouped explicit vector. - * @param rowGroupSizes A vector specifying the sizes of the row groups. * @param rowMetaVariables The meta variables that encode the rows of the matrix. * @param columnMetaVariables The meta variables that encode the columns of the matrix. * @param groupMetaVariables The meta variables that are used to distinguish different row groups. @@ -711,7 +728,7 @@ namespace storm { * @return The matrix that is represented by this ADD and and a vector corresponding to the symbolic vector * (if it was given). */ - std::pair, std::vector> toMatrixVector(storm::dd::Add const& vector, std::vector&& rowGroupSizes, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; + std::pair, std::vector> toMatrixVector(storm::dd::Add const& vector, std::set const& rowMetaVariables, std::set const& columnMetaVariables, std::set const& groupMetaVariables, storm::dd::Odd const& rowOdd, storm::dd::Odd const& columnOdd) const; // The internal ADD that depends on the chosen library. InternalAdd internalAdd; diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 3dbfa2e80..3ca8d2115 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -622,25 +622,26 @@ namespace storm { choiceIt = choices->begin() + startRow; } - for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { - *targetIt = *sourceIt; - ++sourceIt; - localChoice = 1; - if (choices != nullptr) { - *choiceIt = 0; - } - - for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { - if (f(*sourceIt, *targetIt)) { - *targetIt = *sourceIt; - if (choices != nullptr) { - *choiceIt = localChoice; + for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) { + // Only traverse elements if the row group is non-empty. + if (*rowGroupingIt != *(rowGroupingIt + 1)) { + *targetIt = *sourceIt; + ++sourceIt; + localChoice = 1; + if (choices != nullptr) { + *choiceIt = 0; + } + + for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { + if (f(*sourceIt, *targetIt)) { + *targetIt = *sourceIt; + if (choices != nullptr) { + *choiceIt = localChoice; + } } } - } - - if (choices != nullptr) { - ++choiceIt; + } else { + *targetIt = storm::utility::zero(); } } } @@ -678,23 +679,25 @@ namespace storm { choiceIt = choices->begin(); } - for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt) { - *targetIt = *sourceIt; - ++sourceIt; - localChoice = 1; - if (choices != nullptr) { - *choiceIt = 0; - } - for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { - if (f(*sourceIt, *targetIt)) { - *targetIt = *sourceIt; - if (choices != nullptr) { - *choiceIt = localChoice; + for (; targetIt != targetIte; ++targetIt, ++rowGroupingIt, ++choiceIt) { + // Only traverse elements if the row group is non-empty. + if (*rowGroupingIt != *(rowGroupingIt + 1)) { + *targetIt = *sourceIt; + ++sourceIt; + localChoice = 1; + if (choices != nullptr) { + *choiceIt = 0; + } + for (sourceIte = source.begin() + *(rowGroupingIt + 1); sourceIt != sourceIte; ++sourceIt, ++localChoice) { + if (f(*sourceIt, *targetIt)) { + *targetIt = *sourceIt; + if (choices != nullptr) { + *choiceIt = localChoice; + } } } - } - if (choices != nullptr) { - ++choiceIt; + } else { + *targetIt = storm::utility::zero(); } } }