From f1536369be981480d26cbf1c8e93f2fdc81a1b9f Mon Sep 17 00:00:00 2001 From: Tim Quatmann Date: Fri, 5 Jul 2019 15:44:48 +0200 Subject: [PATCH] DeterministicSchedsLpChecker: Bug fixes, further experiments with upper reward bounds --- .../DeterministicSchedsLpChecker.cpp | 148 +++++++++------- .../DeterministicSchedsObjectiveHelper.cpp | 163 +++++++++++++++--- src/storm/solver/GurobiLpSolver.cpp | 4 +- src/storm/storage/geometry/PolytopeTree.h | 17 ++ 4 files changed, 252 insertions(+), 80 deletions(-) diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp index d3d741377..44aabf04e 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsLpChecker.cpp @@ -136,7 +136,7 @@ namespace storm { template std::pair>, std::vector>>> DeterministicSchedsLpChecker::check(storm::Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps) { - std::cout << "Checking " << polytopeTree.toString() << std::endl << "\t"; + STORM_LOG_INFO("Checking " << polytopeTree.toString()); swCheck.start(); STORM_LOG_ASSERT(!currentWeightVector.empty(), "Checking invoked before specifying a weight vector."); if (polytopeTree.isEmpty()) { @@ -159,7 +159,6 @@ namespace storm { checkRecursive(env, polytopeTree, eps, foundPoints, infeasableAreas, 0); swCheck.stop(); - std::cout << " done!" << std::endl; return {foundPoints, infeasableAreas}; } @@ -210,28 +209,39 @@ namespace storm { // Compute an upper bound on the expected number of visits of the states in this ec. // First get a lower bound l on the probability of a path that leaves this MEC. 1-l is an upper bound on Pr_s(X F s). // The desired upper bound is thus 1/(1-(1-l)) = 1/l. See Baier et al., CAV'17 - ValueType expVisitsUpperBound; - uint64_t numStates = std::distance(ec.begin(), ec.end()); - { - std::vector leastPathProbs(transitions.getRowGroupCount(), storm::utility::one()); - std::vector prevLeastPathProbs(transitions.getRowGroupCount(), storm::utility::one()); - for (uint64_t i = 0; i < numStates; ++i) { - for (auto const& stateChoices : ec) { - auto state = stateChoices.first; - auto currVal = prevLeastPathProbs[state]; - for (auto const& transition : transitions.getRowGroup(state)) { + + // To compute l, we multiply the smallest transition probabilities occurring at each state and MEC-Choice + // as well as the smallest 'exit' probability + ValueType lpath = storm::utility::one(); + ValueType minExitProbability = storm::utility::one(); + for (auto const& stateChoices : ec) { + auto state = stateChoices.first; + ValueType minProb = storm::utility::one(); + for (uint64_t choice = transitions.getRowGroupIndices()[state]; choice < transitions.getRowGroupIndices()[state + 1]; ++state) { + if (stateChoices.second.count(choice) == 0) { + // The choice leaves the EC, so we take the sum over the exiting probabilities + ValueType exitProbabilitySum = storm::utility::zero(); + for (auto const& transition : transitions.getRow(choice)) { + if (!ec.containsState(transition.getColumn())) { + exitProbabilitySum += transition.getValue(); + } + } + minExitProbability = std::min(minExitProbability, exitProbabilitySum); + } else { + // Get the minimum over all transition probabilities + for (auto const& transition : transitions.getRow(choice)) { if (!storm::utility::isZero(transition.getValue())) { - currVal = std::min(currVal, transition.getValue() * prevLeastPathProbs[transition.getColumn()]); + minProb = std::min(minProb, transition.getValue()); } } - leastPathProbs[state] = currVal; + } - leastPathProbs.swap(prevLeastPathProbs); } - ValueType l = *std::min_element(leastPathProbs.begin(), leastPathProbs.end()); - expVisitsUpperBound = storm::utility::one() / l; + lpath *= minProb; } - expVisitsUpperBound = storm::utility::one() / expVisitsUpperBound; + lpath *= minExitProbability; + ValueType expVisitsUpperBound = storm::utility::one() / lpath; + STORM_LOG_WARN_COND(expVisitsUpperBound <= storm::utility::convertNumber(1000.0), "Large upper bound for expected visiting times: " << expVisitsUpperBound); // create variables for (auto const& stateChoices : ec) { ecStateVars.emplace(stateChoices.first, lpModel.addBoundedIntegerVariable("e" + std::to_string(stateChoices.first) + varNameSuffix, storm::utility::zero(), storm::utility::one()).getExpression()); @@ -271,6 +281,7 @@ namespace storm { } outs.emplace(stateChoices.first, out); } + uint64_t numStates = std::distance(ec.begin(), ec.end()); for (auto const& stateChoices : ec) { auto in = ins.find(stateChoices.first); STORM_LOG_ASSERT(in != ins.end(), "ec state does not seem to have an incoming transition."); @@ -301,52 +312,59 @@ namespace storm { storm::storage::MaximalEndComponentDecomposition mecs(model.getTransitionMatrix(), backwardTransitions, storm::storage::BitVector(model.getNumberOfStates(), true), choicesWithValueZero); for (auto const& mec : mecs) { + // For each objective we might need to split this mec into several subECs, if the objective yields a non-zero scheduler-independent state value for some states of this ec. std::map, std::vector> excludedStatesToObjIndex; + bool mecContainsSchedulerDependentValue = false; for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { std::set excludedStates; for (auto const& stateChoices : mec) { auto schedIndValueIt = objectiveHelper[objIndex].getSchedulerIndependentStateValues().find(stateChoices.first); - if (schedIndValueIt != objectiveHelper[objIndex].getSchedulerIndependentStateValues().end() && !storm::utility::isZero(schedIndValueIt->second)) { + if (schedIndValueIt == objectiveHelper[objIndex].getSchedulerIndependentStateValues().end()) { + mecContainsSchedulerDependentValue = true; + } else if (!storm::utility::isZero(schedIndValueIt->second)) { excludedStates.insert(stateChoices.first); } } excludedStatesToObjIndex[excludedStates].push_back(objIndex); } - for (auto const& exclStates : excludedStatesToObjIndex) { - if (exclStates.first.empty()) { - auto ecVars = processEc(mec, model.getTransitionMatrix(), "", choiceVariables, *lpModel); - ++ecCounter; - for (auto const& stateVar : ecVars) { - for (auto const& objIndex : exclStates.second) { - result[objIndex][stateVar.first] = stateVar.second; - } - } - } else { - // Compute sub-end components - storm::storage::BitVector subEcStates(model.getNumberOfStates(), false), subEcChoices(model.getNumberOfChoices(), false); - for (auto const& stateChoices : mec) { - if (exclStates.first.count(stateChoices.first) == 0) { - subEcStates.set(stateChoices.first, true); - for (auto const& choice : stateChoices.second) { - subEcChoices.set(choice, true); - } - } - } - storm::storage::MaximalEndComponentDecomposition subEcs(model.getTransitionMatrix(), backwardTransitions, subEcStates, subEcChoices); - for (auto const& subEc : subEcs) { - auto ecVars = processEc(subEc, model.getTransitionMatrix(), "o" + std::to_string(*exclStates.second.begin()), choiceVariables, *lpModel); + // Skip this mec if all state values are independent of the scheduler (e.g. no reward is reachable from here). + if (mecContainsSchedulerDependentValue) { + for (auto const& exclStates : excludedStatesToObjIndex) { + if (exclStates.first.empty()) { + auto ecVars = processEc(mec, model.getTransitionMatrix(), "", choiceVariables, *lpModel); ++ecCounter; for (auto const& stateVar : ecVars) { for (auto const& objIndex : exclStates.second) { result[objIndex][stateVar.first] = stateVar.second; } } + } else { + // Compute sub-end components + storm::storage::BitVector subEcStates(model.getNumberOfStates(), false), subEcChoices(model.getNumberOfChoices(), false); + for (auto const& stateChoices : mec) { + if (exclStates.first.count(stateChoices.first) == 0) { + subEcStates.set(stateChoices.first, true); + for (auto const& choice : stateChoices.second) { + subEcChoices.set(choice, true); + } + } + } + storm::storage::MaximalEndComponentDecomposition subEcs(model.getTransitionMatrix(), backwardTransitions, subEcStates, subEcChoices); + for (auto const& subEc : subEcs) { + auto ecVars = processEc(subEc, model.getTransitionMatrix(), "o" + std::to_string(*exclStates.second.begin()), choiceVariables, *lpModel); + ++ecCounter; + for (auto const& stateVar : ecVars) { + for (auto const& objIndex : exclStates.second) { + result[objIndex][stateVar.first] = stateVar.second; + } + } + } } } } } - std::cout << "found " << ecCounter << "many ECs" << std::endl; + STORM_LOG_WARN_COND(ecCounter == 0, "Processed " << ecCounter << " End components."); return result; } @@ -437,7 +455,8 @@ namespace storm { std::vector ecValVars(model.getNumberOfStates()); if (hasEndComponents) { for (auto const& state : nonBottomStates) { - // For the in-out-encoding, all objectives have the same ECs. Hence, we only care for the variables of the first objective. + // For the in-out-encoding, all objectives have the same ECs (because there are no non-zero scheduler independend state values). + // Hence, we only care for the variables of the first objective. if (ecVars.front()[state].isInitialized()) { ecValVars[state] = lpModel->addBoundedContinuousVariable("z" + std::to_string(state), storm::utility::zero(), visitingTimesUpperBounds[state]).getExpression(); lpModel->addConstraint("", ecValVars[state] <= lpModel->getConstant(visitingTimesUpperBounds[state]) * ecVars.front()[state]); @@ -490,10 +509,13 @@ namespace storm { } } auto objValueVariable = lpModel->addBoundedContinuousVariable("x" + std::to_string(objIndex), objectiveHelper[objIndex].getLowerValueBoundAtState(env, initialState), objectiveHelper[objIndex].getUpperValueBoundAtState(env, initialState)); - lpModel->addConstraint("", objValueVariable == storm::expressions::sum(objValue)); + if (objValue.empty()) { + lpModel->addConstraint("", objValueVariable == zero); + } else { + lpModel->addConstraint("", objValueVariable == storm::expressions::sum(objValue)); + } initialStateResults.push_back(objValueVariable); } - } else { // 'classic' backward encoding. for (uint64_t objIndex = 0; objIndex < objectiveHelper.size(); ++objIndex) { @@ -622,7 +644,7 @@ namespace storm { } template - void DeterministicSchedsLpChecker::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree & polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth) { + void DeterministicSchedsLpChecker::checkRecursive(Environment const& env, storm::storage::geometry::PolytopeTree& polytopeTree, GeometryValueType const& eps, std::vector& foundPoints, std::vector& infeasableAreas, uint64_t const& depth) { std::cout << "."; std::cout.flush(); STORM_LOG_ASSERT(!polytopeTree.isEmpty(), "Tree node is empty"); @@ -646,16 +668,19 @@ namespace storm { lpModel->optimize(); STORM_LOG_TRACE("\tDone solving MILP..."); swLpSolve.stop(); - - //std::cout << "writing model to file out.lp"<< std::endl; - //lpModel->writeModelToFile("out.lp"); + +#ifndef NDEBUG + STORM_PRINT_AND_LOG("Writing model to file '" << polytopeTree.toId() << ".lp'" << std::endl;); + lpModel->writeModelToFile(polytopeTree.toId() + ".lp"); +#endif if (lpModel->isInfeasible()) { infeasableAreas.push_back(polytopeTree.getPolytope()); polytopeTree.clear(); } else { STORM_LOG_ASSERT(!lpModel->isUnbounded(), "LP result is unbounded."); - // TODO: only for debugging +#ifndef NDEBUG validateCurrentModel(env); +#endif Point newPoint; for (auto const& objVar : currentObjectiveVariables) { newPoint.push_back(storm::utility::convertNumber(lpModel->getContinuousValue(objVar))); @@ -688,20 +713,21 @@ namespace storm { if (infeasableAreas.back()->isEmpty()) { infeasableAreas.pop_back(); } - swAux.start(); polytopeTree.setMinus(storm::storage::geometry::Polytope::create({halfspace})); for (auto const& p : newPoints) { foundPoints.push_back(p); polytopeTree.substractDownwardClosure(p, eps); } - swAux.stop(); if (!polytopeTree.isEmpty()) { checkRecursive(env, polytopeTree, eps, foundPoints, infeasableAreas, depth); } } } else { - // Traverse all the children. + // Traverse all the non-empty children. for (uint64_t childId = 0; childId < polytopeTree.getChildren().size(); ++childId) { + if (polytopeTree.getChildren()[childId].isEmpty()) { + continue; + } uint64_t newPointIndex = foundPoints.size(); checkRecursive(env, polytopeTree.getChildren()[childId], eps, foundPoints, infeasableAreas, depth + 1); STORM_LOG_ASSERT(polytopeTree.getChildren()[childId].isEmpty(), "expected empty children."); @@ -734,7 +760,13 @@ namespace storm { uint64_t globalChoiceOffset = model.getTransitionMatrix().getRowGroupIndices()[state]; bool choiceFound = false; for (uint64_t localChoice = 0; localChoice < numChoices; ++localChoice) { - if (lpModel->getIntegerValue(choiceVariables[globalChoiceOffset + localChoice].getBaseExpression().asVariableExpression().getVariable()) == 1) { + bool localChoiceEnabled = false; + if (choiceVarReduction() && localChoice + 1 == numChoices) { + localChoiceEnabled = !choiceFound; + } else { + localChoiceEnabled = (lpModel->getIntegerValue(choiceVariables[globalChoiceOffset + localChoice].getBaseExpression().asVariableExpression().getVariable()) == 1); + } + if (localChoiceEnabled) { STORM_LOG_THROW(!choiceFound, storm::exceptions::UnexpectedException, "Multiple choices selected at state " << state); scheduler.setChoice(localChoice, state); choiceFound = true; @@ -750,9 +782,9 @@ namespace storm { expectedValue = -expectedValue; } ValueType actualValue = objectiveHelper[objIndex].evaluateOnModel(env, *inducedModel); - std::cout << "obj" << objIndex << ": LpSolver: " << storm::utility::convertNumber(expectedValue) << " (" << expectedValue << ")" << std::endl; - std::cout << "obj" << objIndex << ": model checker: " << storm::utility::convertNumber(actualValue) << " (" << actualValue << ")" << std::endl; - STORM_LOG_THROW(storm::utility::convertNumber(storm::utility::abs(actualValue - expectedValue)) <= 1e-4, storm::exceptions::UnexpectedException, "Invalid value for objective " << objIndex << ": expected " << expectedValue << " but got " << actualValue); + double diff = storm::utility::convertNumber(storm::utility::abs(expectedValue - actualValue)); + STORM_PRINT_AND_LOG("Validating Lp solution for objective " << objIndex << ": LP" << storm::utility::convertNumber(expectedValue) << " InducedScheduler=" << storm::utility::convertNumber(actualValue) << " (difference is " << diff << ")" << std::endl); + STORM_LOG_WARN_COND(diff <= 1e-4 * std::abs(storm::utility::convertNumber(actualValue)), "Invalid value for objective " << objIndex << ": expected " << expectedValue << " but got " << actualValue << " (difference is " << diff << ")"); } std::cout << std::endl; diff --git a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp index 7b9fe26c8..deeb71b68 100644 --- a/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp +++ b/src/storm/modelchecker/multiobjective/deterministicScheds/DeterministicSchedsObjectiveHelper.cpp @@ -18,6 +18,12 @@ #include "storm/logic/CloneVisitor.h" #include "storm/environment/solver/MinMaxSolverEnvironment.h" #include "storm/transformer/EndComponentEliminator.h" +#include "storm/utility/solver.h" +#include "storm/solver/LpSolver.h" +#include "storm/solver/GurobiLpSolver.h" +#include "storm/solver/SmtSolver.h" +#include "storm/storage/expressions/Expressions.h" + #include "storm/exceptions/UnexpectedException.h" namespace storm { @@ -260,6 +266,7 @@ namespace storm { return objective.formula->isRewardOperatorFormula() && objective.formula->getSubformula().isTotalRewardFormula(); } + /* template ValueType getLpathDfs(uint64_t currentState, ValueType currentValue, storm::storage::BitVector& visited, storm::storage::BitVector const& mecStates, storm::storage::SparseMatrix const& transitions, uint64_t& counter) { // exhausive dfs @@ -279,6 +286,82 @@ namespace storm { } } + void oneOfManyEncoding(std::vector& operands, storm::solver::SmtSolver& solver) { + assert(!operands.empty()); + uint64_t currentOperand = 0; + while (currentOperand + 2 < operands.size()) { + solver.add(!(operands[currentOperand] && operands[currentOperand + 1])); + auto auxVar = solver.getManager().declareFreshBooleanVariable(true).getExpression(); + solver.add(storm::expressions::iff(auxVar, (operands[currentOperand] || operands[currentOperand + 1]))); + operands.push_back(auxVar); + currentOperand += 2; + } + if (currentOperand + 1 == operands.size()) { + solver.add(operands[currentOperand]); + } else if (currentOperand + 2 == operands.size()) { + solver.add(storm::expressions::xclusiveor(operands[currentOperand], operands[currentOperand + 1])); + } + } + */ + + /* + template + ValueType getExpVisitsUpperBoundForMec(storm::storage::BitVector const& mecStates, storm::storage::SparseMatrix const& transitions) { + auto generalSolver = storm::utility::solver::getLpSolver("mecBounds", storm::solver::LpSolverTypeSelection::Gurobi); + auto solver = dynamic_cast*>(generalSolver.get()); + solver->setOptimizationDirection(storm::solver::OptimizationDirection::Maximize); + + auto one = solver->getConstant(storm::utility::one()); + auto zero = solver->getConstant(storm::utility::zero()); + + std::vector> ins(mecStates.size()), outs(mecStates.size()); + std::vector choiceVars(transitions.getRowCount()), ecVars(mecStates.size()); + std::vector choicesDisjunction; + std::vector leavingSum; + for (auto const& mecState : mecStates) { + std::vector choiceIndicatorVars; + ecVars[mecState] = solver->addLowerBoundedContinuousVariable("e" + std::to_string(mecState), storm::utility::zero()).getExpression(); + ins[mecState].push_back(one); + outs[mecState].push_back(ecVars[mecState]); + leavingSum.push_back(ecVars[mecState]); + choiceIndicatorVars.push_back(solver->addBoundedIntegerVariable("c_mec" + std::to_string(mecState), storm::utility::zero(), storm::utility::one())); + solver->addIndicatorConstraint("", choiceIndicatorVars.back().getBaseExpression().asVariableExpression().getVariable(), 0, ecVars[mecState] <= zero); + for (uint64_t choice = transitions.getRowGroupIndices()[mecState]; choice < transitions.getRowGroupIndices()[mecState + 1]; ++choice) { + choiceVars[choice] = solver->addLowerBoundedContinuousVariable("y" + std::to_string(choice), storm::utility::zero()).getExpression(); + choiceIndicatorVars.push_back(solver->addBoundedIntegerVariable("c" + std::to_string(choice), storm::utility::zero(), storm::utility::one()).getExpression()); + solver->addIndicatorConstraint("", choiceIndicatorVars.back().getBaseExpression().asVariableExpression().getVariable(), 0, choiceVars[choice] <= zero); + choicesDisjunction.push_back(choiceVars[choice]); + outs[mecState].push_back(choiceVars[choice]); + for (auto const& entry : transitions.getRow(choice)) { + if (!storm::utility::isZero(entry.getValue())) { + if (mecStates.get(entry.getColumn())) { + ins[entry.getColumn()].push_back(choiceVars[choice] * solver->getConstant(entry.getValue())); + } else { + leavingSum.push_back(choiceVars[choice] * solver->getConstant(entry.getValue())); + } + } + } + } + //oneOfManyEncoding(choicesSum, *solver); + solver->addConstraint("", storm::expressions::sum(choiceIndicatorVars) == one); + } + for (auto const& mecState : mecStates) { + STORM_LOG_ASSERT(!ins[mecState].empty(), "empty in set at a state"); + STORM_LOG_ASSERT(!outs[mecState].empty(), "empty out set at a state"); + solver->addConstraint("", storm::expressions::sum(ins[mecState]) == storm::expressions::sum(outs[mecState])); + } + STORM_LOG_ASSERT(!leavingSum.empty(), "empty leaving sum at a mec"); + solver->addConstraint("", storm::expressions::sum(leavingSum) == solver->getConstant(storm::utility::convertNumber(mecStates.getNumberOfSetBits()))); + choicesDisjunction.push_back(one); + auto boundVar = solver->addUnboundedContinuousVariable("bound", storm::utility::one()); + solver->addGeneralConstraint("", boundVar, storm::solver::GurobiLpSolver::GeneralConstraintOperator::Max, choicesDisjunction); + solver->optimize(); + STORM_LOG_THROW(!solver->isInfeasible(), storm::exceptions::UnexpectedException, "MEC LP has infeasable solution"); + STORM_LOG_THROW(!solver->isUnbounded(), storm::exceptions::UnexpectedException, "MEC LP has unbounded solution"); + return solver->getObjectiveValue(); + } + */ + template std::vector DeterministicSchedsObjectiveHelper::computeUpperBoundOnExpectedVisitingTimes(storm::storage::SparseMatrix const& modelTransitions, storm::storage::BitVector const& bottomStates, storm::storage::BitVector const& nonBottomStates, bool hasEndComponents) { storm::storage::SparseMatrix transitions; @@ -306,31 +389,71 @@ namespace storm { // slow down mec states by adding self loop probability 1-lpath for (auto const& mec : mecs) { ValueType lpath = storm::utility::one(); - if (true) { + /* //smt/lp + storm::storage::BitVector mecStates(modelTransitions.getRowGroupCount(), false); + uint64_t numTransitions = 0; + uint64_t numChoices = 0; for (auto const& stateChoices : mec) { - auto state = stateChoices.first; - ValueType minProb = storm::utility::one(); - for (uint64_t choice = modelTransitions.getRowGroupIndices()[state]; choice < modelTransitions.getRowGroupIndices()[state + 1]; ++state) { - if (stateChoices.second.count(choice) > 0) { - for (auto const& transition : modelTransitions.getRow(choice)) { - if (!storm::utility::isZero(transition.getValue())) { - minProb = std::min(minProb, transition.getValue()); - } + mecStates.set(stateChoices.first); + numTransitions += modelTransitions.getRowGroup(stateChoices.first).getNumberOfEntries(); + numChoices += modelTransitions.getRowGroupSize(stateChoices.first); + } + std::cout << "Checking a mec with " << mecStates.getNumberOfSetBits() << " states " << numChoices << " choices and " << numTransitions << " transitions." << std::endl; + lpath = storm::utility::one() / getExpVisitsUpperBoundForMec(mecStates, modelTransitions); + }*/ + // Multiply the smallest probability occurring at each state. + for (auto const& stateChoices : mec) { + auto state = stateChoices.first; + ValueType minProb = storm::utility::one(); + for (uint64_t choice = modelTransitions.getRowGroupIndices()[state]; choice < modelTransitions.getRowGroupIndices()[state + 1]; ++state) { + if (stateChoices.second.count(choice) > 0) { + for (auto const& transition : modelTransitions.getRow(choice)) { + if (!storm::utility::isZero(transition.getValue())) { + minProb = std::min(minProb, transition.getValue()); } - } else { - ValueType sum = storm::utility::zero(); - for (auto const& transition : modelTransitions.getRow(choice)) { - if (!mec.containsState(transition.getColumn())) { - sum += transition.getValue(); - } + } + } else { + ValueType sum = storm::utility::zero(); + for (auto const& transition : modelTransitions.getRow(choice)) { + if (!mec.containsState(transition.getColumn())) { + sum += transition.getValue(); } - minProb = std::min(minProb, sum); } + minProb = std::min(minProb, sum); } - lpath *= minProb; } + lpath *= minProb; } - if (false) { + // We multiply the smallest transition probabilities occurring at each state and MEC-Choice + // as well as the smallest 'exit' probability + ValueType minExitProbability = storm::utility::one(); + for (auto const& stateChoices : mec) { + auto state = stateChoices.first; + ValueType minProb = storm::utility::one(); + for (uint64_t choice = transitions.getRowGroupIndices()[state]; choice < transitions.getRowGroupIndices()[state + 1]; ++state) { + if (stateChoices.second.count(choice) == 0) { + // The choice leaves the EC, so we take the sum over the exiting probabilities + ValueType exitProbabilitySum = storm::utility::zero(); + for (auto const& transition : transitions.getRow(choice)) { + if (!mec.containsState(transition.getColumn())) { + exitProbabilitySum += transition.getValue(); + } + } + minExitProbability = std::min(minExitProbability, exitProbabilitySum); + } else { + // Get the minimum over all transition probabilities + for (auto const& transition : transitions.getRow(choice)) { + if (!storm::utility::isZero(transition.getValue())) { + minProb = std::min(minProb, transition.getValue()); + } + } + + } + } + lpath *= minProb; + } + lpath *= minExitProbability; + /* other ideas... std::vector leastPathProbs(modelTransitions.getRowGroupCount(), storm::utility::one()); std::vector prevLeastPathProbs(modelTransitions.getRowGroupCount(), storm::utility::one()); uint64_t mecSize = std::distance(mec.begin(), mec.end()); @@ -370,9 +493,9 @@ namespace storm { std::cout << "\tnew lpath is " << storm::utility::convertNumber(lpath) << ". checked " << counter << " paths." << std::endl; } } - } + }*/ STORM_LOG_ASSERT(!storm::utility::isZero(lpath), "unexpected value of lpath"); - STORM_LOG_WARN_COND(lpath >= storm::utility::convertNumber(0.01), "Small lower bound for the probability to leave a mec: " << storm::utility::convertNumber(lpath) << ". Numerical issues might occur."); + STORM_LOG_WARN_COND(lpath >= storm::utility::convertNumber(0.001), "Small lower bound for the probability to leave a mec: " << storm::utility::convertNumber(lpath) << ". Numerical issues might occur."); uint64_t mecState = modelToSubsystemStateMapping.get()[mec.begin()->first]; // scale all the probabilities at this state with lpath for (uint64_t mecChoice = transitions.getRowGroupIndices()[mecState]; mecChoice < transitions.getRowGroupIndices()[mecState + 1]; ++mecChoice) { diff --git a/src/storm/solver/GurobiLpSolver.cpp b/src/storm/solver/GurobiLpSolver.cpp index 7056b5093..445491fc3 100644 --- a/src/storm/solver/GurobiLpSolver.cpp +++ b/src/storm/solver/GurobiLpSolver.cpp @@ -184,7 +184,7 @@ namespace storm { template void GurobiLpSolver::addConstraint(std::string const& name, storm::expressions::Expression const& constraint) { - STORM_LOG_TRACE("Adding constraint " << name << " to GurobiLpSolver:" << std::endl << "\t" << constraint); + STORM_LOG_TRACE("Adding constraint " << (name == "" ? std::to_string(nextConstraintIndex) : name) << " to GurobiLpSolver:" << std::endl << "\t" << constraint); STORM_LOG_THROW(constraint.isRelationalExpression(), storm::exceptions::InvalidArgumentException, "Illegal constraint is not a relational expression."); STORM_LOG_THROW(constraint.getOperator() != storm::expressions::OperatorType::NotEqual, storm::exceptions::InvalidArgumentException, "Illegal constraint uses inequality operator."); @@ -569,7 +569,7 @@ namespace storm { if (relative) { return result; } else { - return result * getObjectiveValue(); + return storm::utility::abs(result * getObjectiveValue()); } } diff --git a/src/storm/storage/geometry/PolytopeTree.h b/src/storm/storage/geometry/PolytopeTree.h index 54ddb23d9..34a5b16a6 100644 --- a/src/storm/storage/geometry/PolytopeTree.h +++ b/src/storm/storage/geometry/PolytopeTree.h @@ -114,6 +114,23 @@ namespace storm { return children; } + std::string toId() { + if (isEmpty()) { + return "empty"; + } + std::stringstream s; + s << "p"; + auto vertices = getPolytope()->getVertices(); + for (auto const& v : vertices) { + s << "_"; + for (auto const& vi : v) { + s << storm::utility::convertNumber(vi) << "-"; + } + } + s << "_id" << children.data(); + return s.str(); + } + /*! * Returns a string representation of this node (for debugging purposes) */