diff --git a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp index b35fda77a..03fdaf37d 100644 --- a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.cpp @@ -5,6 +5,7 @@ #include "src/adapters/CarlAdapter.h" #include "src/models/sparse/MarkovAutomaton.h" #include "src/models/sparse/StandardRewardModel.h" +#include "src/transformer/EndComponentEliminator.h" #include "src/utility/macros.h" #include "src/utility/vector.h" @@ -23,75 +24,123 @@ namespace storm { STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected."); this->discreteActionRewards[objIndex] = rewModel.hasStateActionRewards() ? rewModel.getStateActionRewardVector() : std::vector(this->data.preprocessedModel.getTransitionMatrix().getRowCount(), storm::utility::zero()); if(rewModel.hasStateRewards()) { + // Note that state rewards are earned over time and thus play no role for probabilistic states for(auto markovianState : this->data.getMarkovianStatesOfPreprocessedModel()) { this->discreteActionRewards[objIndex][this->data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[markovianState]] += rewModel.getStateReward(markovianState) / this->data.preprocessedModel.getExitRate(markovianState); } } - } - } template void SparseMaMultiObjectiveWeightVectorChecker::boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) { - - // Set the digitization constant + // Split the preprocessed model into transitions from/to probabilistic/Markovian states. + SubModel MS = createSubModel(true, weightedRewardVector); + SubModel PS = createSubModel(false, weightedRewardVector); + + // Apply digitization to Markovian transitions ValueType digitizationConstant = getDigitizationConstant(); - std::cout << "Got delta: " << digitizationConstant << std::endl; - - /* - // Allocate some memory so this does not need to happen for each time epoch - std::vector optimalChoicesInCurrentEpoch(this->data.preprocessedModel.getNumberOfStates()); - std::vector choiceValues(weightedRewardVector.size()); - std::vector temporaryResult(this->data.preprocessedModel.getNumberOfStates()); - // Get for each occurring timeBound the indices of the objectives with that bound. - std::map> timeBounds; - storm::storage::BitVector boundedObjectives = ~this->unboundedObjectives; - for(uint_fast64_t objIndex : boundedObjectives) { - uint_fast64_t timeBound = boost::get(this->data.objectives[objIndex].timeBounds.get()); - auto timeBoundIt = timeBounds.insert(std::make_pair(timeBound, storm::storage::BitVector(this->data.objectives.size(), false))).first; - timeBoundIt->second.set(objIndex); - // There is no error for the values of these objectives. - this->offsetsToLowerBound[objIndex] = storm::utility::zero(); - this->offsetsToUpperBound[objIndex] = storm::utility::zero(); - } - storm::storage::BitVector objectivesAtCurrentEpoch = this->unboundedObjectives; - auto timeBoundIt = timeBounds.begin(); - for(uint_fast64_t currentEpoch = timeBoundIt->first; currentEpoch > 0; --currentEpoch) { - if(timeBoundIt != timeBounds.end() && currentEpoch == timeBoundIt->first) { - objectivesAtCurrentEpoch |= timeBoundIt->second; - for(auto objIndex : timeBoundIt->second) { - storm::utility::vector::addScaledVector(weightedRewardVector, this->discreteActionRewards[objIndex], weightVector[objIndex]); - } - ++timeBoundIt; - } + digitize(MS, digitizationConstant); + + // Get for each occurring (digitized) timeBound the indices of the objectives with that bound. + TimeBoundMap lowerTimeBounds; + TimeBoundMap upperTimeBounds; + digitizeTimeBounds(lowerTimeBounds, upperTimeBounds, digitizationConstant); + + // Initialize a minMaxSolver to compute an optimal scheduler (w.r.t. PS) for each epoch + // The end components that stay in PS will be removed. + // Note that the end component elimination could be omitted if we forbid zeno behavior + std::unique_ptr minMax = initMinMaxSolverData(PS); + + // Store the optimal choices of PS as computed by the minMax solver. + std::vector optimalChoicesAtCurrentEpoch(PS.getNumberOfStates(), std::numeric_limits::max()); + + // create a linear equation solver for the model induced by the optimal choice vector. + // the solver will be updated whenever the optimal choice vector has changed. + LinEqSolverData linEq; + linEq.b.resize(PS.getNumberOfStates()); + + // Stores the objectives for which we need to compute values in the current time epoch. + storm::storage::BitVector consideredObjectives = this->unboundedObjectives; + + auto lowerTimeBoundIt = lowerTimeBounds.begin(); + auto upperTimeBoundIt = upperTimeBounds.begin(); + uint_fast64_t currentEpoch = std::max(lowerTimeBounds.empty() ? 0 : lowerTimeBoundIt->first - 1, upperTimeBounds.empty() ? 0 : upperTimeBoundIt->first); // consider lowerBound - 1 since we are interested in the first epoch that passes the bound + while(true) { + // Update the objectives that are considered at the current time epoch as well as the (weighted) reward vectors. + updateDataToCurrentEpoch(MS, PS, *minMax, consideredObjectives, currentEpoch, weightVector, lowerTimeBoundIt, lowerTimeBounds, upperTimeBoundIt, upperTimeBounds); - // Get values and scheduler for weighted sum of objectives - this->data.preprocessedModel.getTransitionMatrix().multiplyWithVector(this->weightedResult, choiceValues); - storm::utility::vector::addVectors(choiceValues, weightedRewardVector, choiceValues); - storm::utility::vector::reduceVectorMax(choiceValues, this->weightedResult, this->data.preprocessedModel.getTransitionMatrix().getRowGroupIndices(), &optimalChoicesInCurrentEpoch); + // Compute the values that can be obtained at probabilistic states in the current time epoch + performPSStep(PS, MS, *minMax, linEq, optimalChoicesAtCurrentEpoch, consideredObjectives); - // get values for individual objectives - // TODO we could compute the result for one of the objectives from the weighted result, the given weight vector, and the remaining objective results. - for(auto objIndex : objectivesAtCurrentEpoch) { - std::vector& objectiveResult = this->objectiveResults[objIndex]; - std::vector objectiveRewards = this->discreteActionRewards[objIndex]; - auto rowGroupIndexIt = this->data.preprocessedModel.getTransitionMatrix().getRowGroupIndices().begin(); - auto optimalChoiceIt = optimalChoicesInCurrentEpoch.begin(); - for(ValueType& stateValue : temporaryResult){ - uint_fast64_t row = (*rowGroupIndexIt) + (*optimalChoiceIt); - ++rowGroupIndexIt; - ++optimalChoiceIt; - stateValue = objectiveRewards[row]; - for(auto const& entry : this->data.preprocessedModel.getTransitionMatrix().getRow(row)) { - stateValue += entry.getValue() * objectiveResult[entry.getColumn()]; - } + // Compute values that can be obtained at Markovian states after letting one (digitized) time unit pass. + // Only perform such a step if there is time left. + if(currentEpoch>0) { + performMSStep(MS, PS, consideredObjectives); + if(currentEpoch % 1000 == 0) { + STORM_LOG_DEBUG(currentEpoch << " digitized time steps left. Current weighted value is " << PS.weightedSolutionVector[0]); } - objectiveResult.swap(temporaryResult); + --currentEpoch; + } else { + break; } } -*/ + + // compose the results from MS and PS + storm::utility::vector::setVectorValues(this->weightedResult, MS.states, MS.weightedSolutionVector); + storm::utility::vector::setVectorValues(this->weightedResult, PS.states, PS.weightedSolutionVector); + for(uint_fast64_t objIndex = 0; objIndex < this->data.objectives.size(); ++objIndex) { + storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], MS.states, MS.objectiveSolutionVectors[objIndex]); + storm::utility::vector::setVectorValues(this->objectiveResults[objIndex], PS.states, PS.objectiveSolutionVectors[objIndex]); + } + } + + + template + typename SparseMaMultiObjectiveWeightVectorChecker::SubModel SparseMaMultiObjectiveWeightVectorChecker::createSubModel(bool createMS, std::vector const& weightedRewardVector) const { + SubModel result; + + storm::storage::BitVector probabilisticStates = ~this->data.getMarkovianStatesOfPreprocessedModel(); + result.states = createMS ? this->data.getMarkovianStatesOfPreprocessedModel() : probabilisticStates; + result.choices = this->data.preprocessedModel.getTransitionMatrix().getRowIndicesOfRowGroups(result.states); + STORM_LOG_ASSERT(!createMS || result.states.getNumberOfSetBits() == result.choices.getNumberOfSetBits(), "row groups for Markovian states should consist of exactly one row"); + + //We need to add diagonal entries for selfloops on Markovian states. + result.toMS = this->data.preprocessedModel.getTransitionMatrix().getSubmatrix(true, result.states, this->data.getMarkovianStatesOfPreprocessedModel(), createMS); + result.toPS = this->data.preprocessedModel.getTransitionMatrix().getSubmatrix(true, result.states, probabilisticStates, false); + STORM_LOG_ASSERT(result.getNumberOfStates() == result.states.getNumberOfSetBits() && result.getNumberOfStates() == result.toMS.getRowGroupCount() && result.getNumberOfStates() == result.toPS.getRowGroupCount(), "Invalid state count for subsystem"); + STORM_LOG_ASSERT(result.getNumberOfChoices() == result.choices.getNumberOfSetBits() && result.getNumberOfChoices() == result.toMS.getRowCount() && result.getNumberOfChoices() == result.toPS.getRowCount(), "Invalid state count for subsystem"); + + result.weightedRewardVector.resize(result.getNumberOfChoices()); + storm::utility::vector::selectVectorValues(result.weightedRewardVector, result.choices, weightedRewardVector); + result.objectiveRewardVectors.resize(this->data.objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->data.objectives.size(); ++objIndex) { + std::vector& objVector = result.objectiveRewardVectors[objIndex]; + objVector = std::vector(result.weightedRewardVector.size(), storm::utility::zero()); + if(this->unboundedObjectives.get(objIndex)) { + storm::utility::vector::selectVectorValues(objVector, result.choices, this->discreteActionRewards[objIndex]); + } else { + typename SparseMaModelType::RewardModelType const& rewModel = this->data.preprocessedModel.getRewardModel(this->data.objectives[objIndex].rewardModelName); + STORM_LOG_ASSERT(!rewModel.hasTransitionRewards(), "Preprocessed Reward model has transition rewards which is not expected."); + STORM_LOG_ASSERT(!rewModel.hasStateRewards(), "State rewards for bounded objectives for MAs are not expected (bounded rewards are not supported)."); + if(rewModel.hasStateActionRewards()) { + storm::utility::vector::selectVectorValues(objVector, result.choices, rewModel.getStateActionRewardVector()); + } + } + } + + result.weightedSolutionVector.resize(result.getNumberOfStates()); + storm::utility::vector::selectVectorValues(result.weightedSolutionVector, result.states, this->weightedResult); + result.objectiveSolutionVectors.resize(this->data.objectives.size()); + for(uint_fast64_t objIndex = 0; objIndex < this->data.objectives.size(); ++objIndex) { + result.objectiveSolutionVectors[objIndex].resize(result.weightedSolutionVector.size()); + storm::utility::vector::selectVectorValues(result.objectiveSolutionVectors[objIndex], result.states, this->objectiveResults[objIndex]); + } + + result.auxChoiceValues.resize(result.getNumberOfChoices()); + + return result; } template @@ -109,7 +158,7 @@ namespace storm { for(auto const& obj : this->data.objectives) { if(obj.timeBounds) { if(obj.timeBounds->which() == 0) { - lowerUpperBounds.emplace_back(storm::utility::zero(), boost::get(*obj.timeBounds)); + lowerUpperBounds.emplace_back(storm::utility::zero(), storm::utility::convertNumber(boost::get(*obj.timeBounds))); eToPowerOfMinusMaxRateTimesBound.emplace_back(storm::utility::one(), std::exp(-maxRate * lowerUpperBounds.back().second)); } else { auto const& pair = boost::get>(*obj.timeBounds); @@ -167,18 +216,205 @@ namespace storm { return delta; } - template template ::SupportsExponential, int>::type> VT SparseMaMultiObjectiveWeightVectorChecker::getDigitizationConstant() const { STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); } + + template + template ::SupportsExponential, int>::type> + void SparseMaMultiObjectiveWeightVectorChecker::digitize(SubModel& MS, VT const& digitizationConstant) const { + std::vector rateVector(MS.getNumberOfChoices()); + storm::utility::vector::selectVectorValues(rateVector, MS.states, this->data.preprocessedModel.getExitRates()); + for(uint_fast64_t row = 0; row < rateVector.size(); ++row) { + VT const eToMinusRateTimesDelta = std::exp(-rateVector[row] * digitizationConstant); + for(auto& entry : MS.toMS.getRow(row)) { + entry.setValue((storm::utility::one() - eToMinusRateTimesDelta) * entry.getValue()); + if(entry.getColumn() == row) { + entry.setValue(entry.getValue() + eToMinusRateTimesDelta); + } + } + for(auto& entry : MS.toPS.getRow(row)) { + entry.setValue((storm::utility::one() - eToMinusRateTimesDelta) * entry.getValue()); + } + MS.weightedRewardVector[row] *= storm::utility::one() - eToMinusRateTimesDelta; + for(auto& objVector : MS.objectiveRewardVectors) { + objVector[row] *= storm::utility::one() - eToMinusRateTimesDelta; + } + } + } + + template + template ::SupportsExponential, int>::type> + void SparseMaMultiObjectiveWeightVectorChecker::digitize(SubModel& subModel, VT const& digitizationConstant) const { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); + } + template + template ::SupportsExponential, int>::type> + void SparseMaMultiObjectiveWeightVectorChecker::digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant) { + + VT const maxRate = this->data.preprocessedModel.getMaximalExitRate(); + storm::storage::BitVector boundedObjectives = ~this->unboundedObjectives; + for(uint_fast64_t objIndex : boundedObjectives) { + boost::optional objLowerBound, objUpperBound; + if(this->data.objectives[objIndex].timeBounds->which() == 0) { + objUpperBound = storm::utility::convertNumber(boost::get(this->data.objectives[objIndex].timeBounds.get())); + } else { + auto const& pair = boost::get>(this->data.objectives[objIndex].timeBounds.get()); + if(!storm::utility::isZero(pair.first)) { + objLowerBound = storm::utility::convertNumber(pair.first); + } + objUpperBound = storm::utility::convertNumber(pair.second); + } + if(objLowerBound) { + uint_fast64_t digitizedBound = storm::utility::convertNumber((*objLowerBound)/digitizationConstant); + auto timeBoundIt = lowerTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->data.objectives.size(), false))).first; + timeBoundIt->second.set(objIndex); + VT digitizationError = storm::utility::one(); + digitizationError -= std::exp(-maxRate * (*objLowerBound)) * storm::utility::pow(storm::utility::one() + maxRate * digitizationConstant, digitizedBound); + this->offsetsToLowerBound[objIndex] = -digitizationError; + + } else { + this->offsetsToLowerBound[objIndex] = storm::utility::zero(); + } + + if(objUpperBound) { + uint_fast64_t digitizedBound = storm::utility::convertNumber((*objUpperBound)/digitizationConstant); + auto timeBoundIt = upperTimeBounds.insert(std::make_pair(digitizedBound, storm::storage::BitVector(this->data.objectives.size(), false))).first; + timeBoundIt->second.set(objIndex); + VT digitizationError = storm::utility::one(); + digitizationError -= std::exp(-maxRate * (*objUpperBound)) * storm::utility::pow(storm::utility::one() + maxRate * digitizationConstant, digitizedBound); + this->offsetsToUpperBound[objIndex] = digitizationError; + } else { + this->offsetsToUpperBound[objIndex] = storm::utility::zero(); + } + STORM_LOG_ASSERT(this->offsetsToUpperBound[objIndex] - this->offsetsToLowerBound[objIndex] <= this->maximumLowerUpperBoundGap, "Precision not sufficient."); + } + } + + template + template ::SupportsExponential, int>::type> + void SparseMaMultiObjectiveWeightVectorChecker::digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant) { + STORM_LOG_THROW(false, storm::exceptions::InvalidOperationException, "Computing bounded probabilities of MAs is unsupported for this value type."); + } + + template + std::unique_ptr::MinMaxSolverData> SparseMaMultiObjectiveWeightVectorChecker::initMinMaxSolverData(SubModel const& PS) const { + std::unique_ptr result(new MinMaxSolverData()); + + storm::storage::BitVector choicesStayingInPS(PS.getNumberOfChoices(), false); + for(uint_fast64_t choice = 0; choice < PS.toPS.getRowCount(); ++choice) { + if(storm::utility::isOne(PS.toPS.getRowSum(choice))) { + choicesStayingInPS.set(choice); + } + } + auto ecEliminatorResult = storm::transformer::EndComponentEliminator::transform(PS.toPS, choicesStayingInPS & storm::utility::vector::filterZero(PS.weightedRewardVector), storm::storage::BitVector(PS.getNumberOfStates(), true)); + result->matrix = std::move(ecEliminatorResult.matrix); + result->toPSChoiceMapping = std::move(ecEliminatorResult.newToOldRowMapping); + result->fromPSStateMapping = std::move(ecEliminatorResult.oldToNewStateMapping); + + result->b.resize(result->matrix.getRowCount()); + result->x.resize(result->matrix.getRowGroupCount()); + for(uint_fast64_t state=0; state < result->fromPSStateMapping.size(); ++state) { + if(result->fromPSStateMapping[state] < result->x.size()) { + result->x[result->fromPSStateMapping[state]] = PS.weightedSolutionVector[state]; + } + } + + storm::solver::GeneralMinMaxLinearEquationSolverFactory minMaxSolverFactory; + result->solver = minMaxSolverFactory.create(result->matrix); + result->solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + result->solver->setTrackScheduler(true); + result->solver->allocateAuxMemory(storm::solver::MinMaxLinearEquationSolverOperation::SolveEquations); + + return result; + } + + template + void SparseMaMultiObjectiveWeightVectorChecker::updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector const& weightVector, TimeBoundMap::iterator& lowerTimeBoundIt, TimeBoundMap const& lowerTimeBounds, TimeBoundMap::iterator& upperTimeBoundIt, TimeBoundMap const& upperTimeBounds) { + + //For lower time bounds we need to react when the currentEpoch passed the bound + // Hence, we substract 1 from the lower time bounds. + if(lowerTimeBoundIt != lowerTimeBounds.end() && currentEpoch == lowerTimeBoundIt->first - 1) { + for(auto objIndex : lowerTimeBoundIt->second) { + // No more reward is earned for this objective. + storm::utility::vector::addScaledVector(MS.weightedRewardVector, MS.objectiveRewardVectors[objIndex], -weightVector[objIndex]); + storm::utility::vector::addScaledVector(PS.weightedRewardVector, PS.objectiveRewardVectors[objIndex], -weightVector[objIndex]); + MS.objectiveRewardVectors[objIndex] = std::vector(MS.objectiveRewardVectors[objIndex].size(), storm::utility::zero()); + PS.objectiveRewardVectors[objIndex] = std::vector(PS.objectiveRewardVectors[objIndex].size(), storm::utility::zero()); + } + ++lowerTimeBoundIt; + } + + if(upperTimeBoundIt != upperTimeBounds.end() && currentEpoch == upperTimeBoundIt->first) { + consideredObjectives |= upperTimeBoundIt->second; + for(auto objIndex : upperTimeBoundIt->second) { + // This objective now plays a role in the weighted sum + storm::utility::vector::addScaledVector(MS.weightedRewardVector, MS.objectiveRewardVectors[objIndex], weightVector[objIndex]); + storm::utility::vector::addScaledVector(PS.weightedRewardVector, PS.objectiveRewardVectors[objIndex], weightVector[objIndex]); + } + ++upperTimeBoundIt; + } + + // Update the solver data + PS.toMS.multiplyWithVector(MS.weightedSolutionVector, PS.auxChoiceValues); + storm::utility::vector::addVectors(PS.auxChoiceValues, PS.weightedRewardVector, PS.auxChoiceValues); + storm::utility::vector::selectVectorValues(minMax.b, minMax.toPSChoiceMapping, PS.auxChoiceValues); + } + + template + void SparseMaMultiObjectiveWeightVectorChecker::performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector& optimalChoicesAtCurrentEpoch, storm::storage::BitVector const& consideredObjectives) const { + + // compute a choice vector for the probabilistic states that is optimal w.r.t. the weighted reward vector + std::vector newOptimalChoices(PS.getNumberOfStates()); + minMax.solver->solveEquations(minMax.x, minMax.b); + this->transformReducedSolutionToOriginalModel(minMax.matrix, minMax.x, minMax.solver->getScheduler()->getChoices(), minMax.toPSChoiceMapping, minMax.fromPSStateMapping, PS.toPS, PS.weightedSolutionVector, newOptimalChoices); + + // check whether the linEqSolver needs to be updated, i.e., whether the scheduler has changed + if(newOptimalChoices != optimalChoicesAtCurrentEpoch) { + std::cout << "Scheduler changed!"; + optimalChoicesAtCurrentEpoch.swap(newOptimalChoices); + linEq.solver = nullptr; + storm::storage::SparseMatrix linEqMatrix = PS.toPS.selectRowsFromRowGroups(optimalChoicesAtCurrentEpoch, true); + linEqMatrix.convertToEquationSystem(); + linEq.solver = linEq.factory.create(std::move(linEqMatrix)); + } + for(auto objIndex : consideredObjectives) { + PS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], PS.auxChoiceValues); + storm::utility::vector::addVectors(PS.auxChoiceValues, PS.objectiveRewardVectors[objIndex], PS.auxChoiceValues); + storm::utility::vector::selectVectorValues(linEq.b, optimalChoicesAtCurrentEpoch, PS.toPS.getRowGroupIndices(), PS.auxChoiceValues); + linEq.solver->solveEquations(PS.objectiveSolutionVectors[objIndex], linEq.b); + } + } + + template + void SparseMaMultiObjectiveWeightVectorChecker::performMSStep(SubModel& MS, SubModel const& PS, storm::storage::BitVector const& consideredObjectives) const { + + MS.toMS.multiplyWithVector(MS.weightedSolutionVector, MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.weightedRewardVector, MS.auxChoiceValues, MS.weightedSolutionVector); + MS.toPS.multiplyWithVector(PS.weightedSolutionVector, MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.weightedSolutionVector, MS.auxChoiceValues, MS.weightedSolutionVector); + + for(auto objIndex : consideredObjectives) { + MS.toMS.multiplyWithVector(MS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.objectiveRewardVectors[objIndex], MS.auxChoiceValues, MS.objectiveSolutionVectors[objIndex]); + MS.toPS.multiplyWithVector(PS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues); + storm::utility::vector::addVectors(MS.objectiveSolutionVectors[objIndex], MS.auxChoiceValues, MS.objectiveSolutionVectors[objIndex]); + } + } + + template class SparseMaMultiObjectiveWeightVectorChecker>; template double SparseMaMultiObjectiveWeightVectorChecker>::getDigitizationConstant() const; + template void SparseMaMultiObjectiveWeightVectorChecker>::digitize(SparseMaMultiObjectiveWeightVectorChecker>::SubModel& subModel, double const& digitizationConstant) const; + template void SparseMaMultiObjectiveWeightVectorChecker>::digitizeTimeBounds(SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& upperTimeBounds, double const& digitizationConstant); #ifdef STORM_HAVE_CARL - template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker>::getDigitizationConstant() const; template class SparseMaMultiObjectiveWeightVectorChecker>; + template storm::RationalNumber SparseMaMultiObjectiveWeightVectorChecker>::getDigitizationConstant() const; + template void SparseMaMultiObjectiveWeightVectorChecker>::digitize(SparseMaMultiObjectiveWeightVectorChecker>::SubModel& subModel, storm::RationalNumber const& digitizationConstant) const; + template void SparseMaMultiObjectiveWeightVectorChecker>::digitizeTimeBounds(SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& lowerTimeBounds, SparseMaMultiObjectiveWeightVectorChecker>::TimeBoundMap& upperTimeBounds, storm::RationalNumber const& digitizationConstant); #endif } diff --git a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.h b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.h index 2670fb5e3..4371856a0 100644 --- a/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.h +++ b/src/modelchecker/multiobjective/helper/SparseMaMultiObjectiveWeightVectorChecker.h @@ -5,6 +5,8 @@ #include #include "src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h" +#include "src/solver/LinearEquationSolver.h" +#include "src/solver/MinMaxLinearEquationSolver.h" #include "src/utility/NumberTraits.h" namespace storm { @@ -22,27 +24,129 @@ namespace storm { public: typedef typename SparseMaModelType::ValueType ValueType; typedef SparseMultiObjectivePreprocessorData PreprocessorData; - + SparseMaMultiObjectiveWeightVectorChecker(PreprocessorData const& data); private: + /* + * Stores (digitized) time bounds in descending order + */ + typedef std::map> TimeBoundMap; + + /* + * Stores the ingredients of a sub model + */ + struct SubModel { + storm::storage::BitVector states; // The states that are part of this sub model + storm::storage::BitVector choices; // The choices that are part of this sub model + + storm::storage::SparseMatrix toMS; // Transitions to Markovian states + storm::storage::SparseMatrix toPS; // Transitions to probabilistic states + + std::vector weightedRewardVector; + std::vector> objectiveRewardVectors; + + std::vector weightedSolutionVector; + std::vector> objectiveSolutionVectors; + + std::vector auxChoiceValues; //stores auxiliary values for every choice + + uint_fast64_t getNumberOfStates() const { return toMS.getRowGroupCount(); }; + uint_fast64_t getNumberOfChoices() const { return toMS.getRowCount(); }; + }; + + /* + * Stores the data that is relevant to invoke the minMaxSolver and retrieve the result. + */ + struct MinMaxSolverData { + std::unique_ptr> solver; // the solver itself + + storm::storage::SparseMatrix matrix; // the considered matrix + std::vector toPSChoiceMapping; // maps rows of the considered matrix to choices of the PS SubModel + std::vector fromPSStateMapping; // maps states of the PS SubModel to row groups of the considered matrix + + std::vector b; + std::vector x; + }; + + struct LinEqSolverData { + storm::solver::GeneralLinearEquationSolverFactory factory; + std::unique_ptr> solver; + + std::vector b; + }; + /*! * * @param weightVector the weight vector of the current check - * @param weightedRewardVector the weighted rewards (initially only considering the unbounded objectives, will be extended to all objectives) + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. */ virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) override; /*! - * - * Retrieves the delta used for the digitization + * Retrieves the data for a submodel of the data->preprocessedModel + * @param createMS if true, the submodel containing the Markovian states is created. + * if false, the submodel containing the probabilistic states is created. + */ + SubModel createSubModel(bool createMS, std::vector const& weightedRewardVector) const; + + /*! + * Retrieves the delta used for digitization */ template ::SupportsExponential, int>::type = 0> VT getDigitizationConstant() const; template ::SupportsExponential, int>::type = 0> VT getDigitizationConstant() const; - + + /*! + * Digitizes the given matrix and vectors w.r.t. the given digitization constant and the given rate vector. + */ + template ::SupportsExponential, int>::type = 0> + void digitize(SubModel& subModel, VT const& digitizationConstant) const; + template ::SupportsExponential, int>::type = 0> + void digitize(SubModel& subModel, VT const& digitizationConstant) const; + + /* + * Fills the given maps with the digitized time bounds. Also sets the offsetsToLowerBound / offsetsToUpperBound values + * according to the digitization error + */ + template ::SupportsExponential, int>::type = 0> + void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); + template ::SupportsExponential, int>::type = 0> + void digitizeTimeBounds(TimeBoundMap& lowerTimeBounds, TimeBoundMap& upperTimeBounds, VT const& digitizationConstant); + + + /*! + * Computes a reduction of the PStoPS submodel which does not contain end components in which no choice with reward is taken. + * With the reduction, a MinMaxSolver is initialized. + */ + std::unique_ptr initMinMaxSolverData(SubModel const& PS) const; + + /* + * Updates the reward vectors within the split model, + * the reward vector of the reduced PStoPS model, and + * objectives that are considered at the current time epoch. + */ + void updateDataToCurrentEpoch(SubModel& MS, SubModel& PS, MinMaxSolverData& minMax, storm::storage::BitVector& consideredObjectives, uint_fast64_t const& currentEpoch, std::vector const& weightVector, TimeBoundMap::iterator& lowerTimeBoundIt, TimeBoundMap const& lowerTimeBounds, TimeBoundMap::iterator& upperTimeBoundIt, TimeBoundMap const& upperTimeBounds); + + /* + * Performs a step for the probabilistic states, that is + * * Compute an optimal scheduler for the weighted reward sum + * * Compute the values for the individual objectives w.r.t. that scheduler + * + * The resulting values represent the rewards at probabilistic states that are obtained at the current time epoch. + */ + void performPSStep(SubModel& PS, SubModel const& MS, MinMaxSolverData& minMax, LinEqSolverData& linEq, std::vector& optimalChoicesAtCurrentEpoch, storm::storage::BitVector const& consideredObjectives) const; + + /* + * Performs a step for the Markovian states, that is + * * Compute values for the weighted reward sum as well as for the individual objectives + * + * The resulting values represent the rewards at Markovian states that are obtained after one (digitized) time unit has passed. + */ + void performMSStep(SubModel& MS, SubModel const& PS, storm::storage::BitVector const& consideredObjectives) const; + }; } diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.cpp index 1f473006f..476775f60 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.cpp @@ -38,11 +38,11 @@ namespace storm { this->offsetsToLowerBound[objIndex] = storm::utility::zero(); this->offsetsToUpperBound[objIndex] = storm::utility::zero(); } - storm::storage::BitVector objectivesAtCurrentEpoch = this->unboundedObjectives; + storm::storage::BitVector consideredObjectives = this->unboundedObjectives; auto timeBoundIt = timeBounds.begin(); for(uint_fast64_t currentEpoch = timeBoundIt->first; currentEpoch > 0; --currentEpoch) { if(timeBoundIt != timeBounds.end() && currentEpoch == timeBoundIt->first) { - objectivesAtCurrentEpoch |= timeBoundIt->second; + consideredObjectives |= timeBoundIt->second; for(auto objIndex : timeBoundIt->second) { storm::utility::vector::addScaledVector(weightedRewardVector, this->discreteActionRewards[objIndex], weightVector[objIndex]); } @@ -56,7 +56,7 @@ namespace storm { // get values for individual objectives // TODO we could compute the result for one of the objectives from the weighted result, the given weight vector, and the remaining objective results. - for(auto objIndex : objectivesAtCurrentEpoch) { + for(auto objIndex : consideredObjectives) { std::vector& objectiveResult = this->objectiveResults[objIndex]; std::vector objectiveRewards = this->discreteActionRewards[objIndex]; auto rowGroupIndexIt = this->data.preprocessedModel.getTransitionMatrix().getRowGroupIndices().begin(); diff --git a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.h b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.h index 93d67c996..f7b8f6747 100644 --- a/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.h +++ b/src/modelchecker/multiobjective/helper/SparseMdpMultiObjectiveWeightVectorChecker.h @@ -32,7 +32,7 @@ namespace storm { * - computes the values of these objectives w.r.t. this scheduler * * @param weightVector the weight vector of the current check - * @param weightedRewardVector the weighted rewards (initially only considering the unbounded objectives, will be extended to all objectives) + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. */ virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) override; diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp index f60295ad4..4741f2a0a 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.cpp @@ -113,49 +113,11 @@ namespace storm { solver->solveEquations(subResult, subRewardVector); this->weightedResult = std::vector(data.preprocessedModel.getNumberOfStates()); - this->scheduler = storm::storage::TotalScheduler(data.preprocessedModel.getNumberOfStates()); - std::unique_ptr solverScheduler = solver->getScheduler(); - storm::storage::BitVector statesWithUndefinedScheduler(data.preprocessedModel.getNumberOfStates(), false); - for(uint_fast64_t state = 0; state < data.preprocessedModel.getNumberOfStates(); ++state) { - uint_fast64_t stateInReducedModel = ecEliminatorResult.oldToNewStateMapping[state]; - // Check if the state exists in the reduced model - if(stateInReducedModel < ecEliminatorResult.matrix.getRowGroupCount()) { - this->weightedResult[state] = subResult[stateInReducedModel]; - // Check if the chosen row originaly belonged to the current state - uint_fast64_t chosenRowInReducedModel = ecEliminatorResult.matrix.getRowGroupIndices()[stateInReducedModel] + solverScheduler->getChoice(stateInReducedModel); - uint_fast64_t chosenRowInOriginalModel = ecEliminatorResult.newToOldRowMapping[chosenRowInReducedModel]; - if(chosenRowInOriginalModel >= data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state] && - chosenRowInOriginalModel < data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state+1]) { - this->scheduler.setChoice(state, chosenRowInOriginalModel - data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state]); - } else { - statesWithUndefinedScheduler.set(state); - } - } else { - // if the state does not exist in the reduced model, it means that the result is always zero, independent of the scheduler. - // Hence, we don't have to set the scheduler explicitly - this->weightedResult[state] = storm::utility::zero(); - } - } - while(!statesWithUndefinedScheduler.empty()) { - for(auto state : statesWithUndefinedScheduler) { - // Try to find a choice that stays inside the EC (i.e., for which all successors are represented by the same state in the reduced model) - // And at least one successor has a defined scheduler. - // This way, a scheduler is chosen that leads (with probability one) to the state of the EC for which the scheduler is defined - uint_fast64_t stateInReducedModel = ecEliminatorResult.oldToNewStateMapping[state]; - for(uint_fast64_t row = data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state]; row < data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state+1]; ++row) { - bool rowStaysInEC = true; - bool rowLeadsToDefinedScheduler = false; - for(auto const& entry : data.preprocessedModel.getTransitionMatrix().getRow(row)) { - rowStaysInEC &= ( stateInReducedModel == ecEliminatorResult.oldToNewStateMapping[entry.getColumn()]); - rowLeadsToDefinedScheduler |= !statesWithUndefinedScheduler.get(entry.getColumn()); - } - if(rowStaysInEC && rowLeadsToDefinedScheduler) { - this->scheduler.setChoice(state, row - data.preprocessedModel.getTransitionMatrix().getRowGroupIndices()[state]); - statesWithUndefinedScheduler.set(state, false); - } - } - } - } + std::vector optimalChoices(data.preprocessedModel.getNumberOfStates()); + + transformReducedSolutionToOriginalModel(ecEliminatorResult.matrix, subResult, solver->getScheduler()->getChoices(), ecEliminatorResult.newToOldRowMapping, ecEliminatorResult.oldToNewStateMapping, this->data.preprocessedModel.getTransitionMatrix(), this->weightedResult, optimalChoices); + + this->scheduler = storm::storage::TotalScheduler(std::move(optimalChoices)); } template @@ -190,6 +152,61 @@ namespace storm { } } + template + void SparseMultiObjectiveWeightVectorChecker::transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix const& reducedMatrix, + std::vector const& reducedSolution, + std::vector const& reducedOptimalChoices, + std::vector const& reducedToOriginalChoiceMapping, + std::vector const& originalToReducedStateMapping, + storm::storage::SparseMatrix const& originalMatrix, + std::vector& originalSolution, + std::vector& originalOptimalChoices) const { + + storm::storage::BitVector statesWithUndefinedScheduler(originalMatrix.getRowGroupCount(), false); + for(uint_fast64_t state = 0; state < originalMatrix.getRowGroupCount(); ++state) { + uint_fast64_t stateInReducedModel = originalToReducedStateMapping[state]; + // Check if the state exists in the reduced model + if(stateInReducedModel < reducedMatrix.getRowGroupCount()) { + originalSolution[state] = reducedSolution[stateInReducedModel]; + // Check if the chosen row originaly belonged to the current state + uint_fast64_t chosenRowInReducedModel = reducedMatrix.getRowGroupIndices()[stateInReducedModel] + reducedOptimalChoices[stateInReducedModel]; + uint_fast64_t chosenRowInOriginalModel = reducedToOriginalChoiceMapping[chosenRowInReducedModel]; + if(chosenRowInOriginalModel >= originalMatrix.getRowGroupIndices()[state] && + chosenRowInOriginalModel < originalMatrix.getRowGroupIndices()[state+1]) { + originalOptimalChoices[state] = chosenRowInOriginalModel - originalMatrix.getRowGroupIndices()[state]; + } else { + statesWithUndefinedScheduler.set(state); + } + } else { + // if the state does not exist in the reduced model, it means that the result is always zero, independent of the scheduler. + // Hence, we don't have to set the scheduler explicitly + originalSolution[state] = storm::utility::zero(); + } + } + while(!statesWithUndefinedScheduler.empty()) { + for(auto state : statesWithUndefinedScheduler) { + // Try to find a choice that stays inside the EC (i.e., for which all successors are represented by the same state in the reduced model) + // And at least one successor has a defined scheduler. + // This way, a scheduler is chosen that leads (with probability one) to the state of the EC for which the scheduler is defined + uint_fast64_t stateInReducedModel = originalToReducedStateMapping[state]; + for(uint_fast64_t row = originalMatrix.getRowGroupIndices()[state]; row < originalMatrix.getRowGroupIndices()[state+1]; ++row) { + bool rowStaysInEC = true; + bool rowLeadsToDefinedScheduler = false; + for(auto const& entry : originalMatrix.getRow(row)) { + rowStaysInEC &= ( stateInReducedModel == originalToReducedStateMapping[entry.getColumn()]); + rowLeadsToDefinedScheduler |= !statesWithUndefinedScheduler.get(entry.getColumn()); + } + if(rowStaysInEC && rowLeadsToDefinedScheduler) { + originalOptimalChoices[state] = row - originalMatrix.getRowGroupIndices()[state]; + statesWithUndefinedScheduler.set(state, false); + } + } + } + } + } + + + template class SparseMultiObjectiveWeightVectorChecker>; template class SparseMultiObjectiveWeightVectorChecker>; #ifdef STORM_HAVE_CARL diff --git a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h index 88eebb1c4..0267ac1c0 100644 --- a/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h +++ b/src/modelchecker/multiobjective/helper/SparseMultiObjectiveWeightVectorChecker.h @@ -81,10 +81,23 @@ namespace storm { * - computes the values of these objectives w.r.t. this scheduler * * @param weightVector the weight vector of the current check - * @param weightedRewardVector the weighted rewards (initially only considering the unbounded objectives, will be extended to all objectives) + * @param weightedRewardVector the weighted rewards considering the unbounded objectives. Will be invalidated after calling this. */ virtual void boundedPhase(std::vector const& weightVector, std::vector& weightedRewardVector) = 0; + /*! + * Transforms the results of a min-max-solver that considers a reduced model (without end components) to a result for the original (unreduced) model + */ + void transformReducedSolutionToOriginalModel(storm::storage::SparseMatrix const& reducedMatrix, + std::vector const& reducedSolution, + std::vector const& reducedOptimalChoices, + std::vector const& reducedToOriginalChoiceMapping, + std::vector const& originalToReducedStateMapping, + storm::storage::SparseMatrix const& originalMatrix, + std::vector& originalSolution, + std::vector& originalOptimalChoices) const; + + // stores the considered information of the multi-objective model checking problem PreprocessorData const& data; // stores the indices of the objectives for which there is no time bound diff --git a/src/storage/SparseMatrix.cpp b/src/storage/SparseMatrix.cpp index 95f65ce99..f8d3b328f 100644 --- a/src/storage/SparseMatrix.cpp +++ b/src/storage/SparseMatrix.cpp @@ -581,6 +581,17 @@ namespace storm { return rowGroupIndices.get(); } + template + storm::storage::BitVector SparseMatrix::getRowIndicesOfRowGroups(storm::storage::BitVector const& groups) const { + storm::storage::BitVector res(this->getRowCount(), false); + for(auto group : groups) { + for(uint_fast64_t row = this->getRowGroupIndices()[group]; row < this->getRowGroupIndices()[group+1]; ++row) { + res.set(row, true); + } + } + return res; + } + template void SparseMatrix::makeRowsAbsorbing(storm::storage::BitVector const& rows) { for (auto row : rows) { diff --git a/src/storage/SparseMatrix.h b/src/storage/SparseMatrix.h index 8e386f446..e6563f47d 100644 --- a/src/storage/SparseMatrix.h +++ b/src/storage/SparseMatrix.h @@ -569,6 +569,14 @@ namespace storm { */ std::vector const& getRowGroupIndices() const; + /*! + * Returns the indices of the rows that belong to one of the selected row groups. + * + * @param groups the selected row groups + * @return a bit vector that is true at position i iff the row group of row i is selected. + */ + storm::storage::BitVector getRowIndicesOfRowGroups(storm::storage::BitVector const& groups) const; + /*! * This function makes the given rows absorbing. * diff --git a/src/utility/constants.cpp b/src/utility/constants.cpp index b2cef709d..e985c7dc3 100644 --- a/src/utility/constants.cpp +++ b/src/utility/constants.cpp @@ -109,6 +109,16 @@ namespace storm { return number; } + template<> + uint_fast64_t convertNumber(double const& number){ + return std::llround(number); + } + + template<> + double convertNumber(uint_fast64_t const& number){ + return number; + } + template ValueType sqrt(ValueType const& number) { return std::sqrt(number); @@ -162,6 +172,16 @@ namespace storm { RationalNumber convertNumber(double const& number){ return carl::rationalize(number); } + + template<> + uint_fast64_t convertNumber(RationalNumber const& number){ + return carl::toInt(carl::round(number)); + } + + template<> + RationalNumber convertNumber(uint_fast64_t const& number){ + return RationalNumber(number); + } template<> RationalFunction convertNumber(double const& number){