From 65fd5cd795bcac07644d995fc41895fbafbdf59e Mon Sep 17 00:00:00 2001 From: Joachim Klein Date: Fri, 21 Jul 2017 16:43:08 +0000 Subject: [PATCH] (interval iteration, explicit) actually perform interval iteration using the explicit engine Supported: DTMC reachability probability and expected reward computations. MDP Pmax, Pmin, Rmax, Rmin for reachability. git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@12142 bbc10eb1-c90d-0410-af57-cb519fbb1720 --- prism/src/explicit/DTMCModelChecker.java | 634 ++++++++++++++++++- prism/src/explicit/DijkstraSweepMPI.java | 312 ++++++++++ prism/src/explicit/MDPModelChecker.java | 746 ++++++++++++++++++++++- 3 files changed, 1687 insertions(+), 5 deletions(-) create mode 100644 prism/src/explicit/DijkstraSweepMPI.java diff --git a/prism/src/explicit/DTMCModelChecker.java b/prism/src/explicit/DTMCModelChecker.java index 5f5f964b..6f8d241c 100644 --- a/prism/src/explicit/DTMCModelChecker.java +++ b/prism/src/explicit/DTMCModelChecker.java @@ -29,28 +29,38 @@ package explicit; import java.io.File; import java.util.Arrays; import java.util.BitSet; +import java.util.Iterator; import java.util.List; import java.util.Map; +import java.util.Map.Entry; +import java.util.PrimitiveIterator; import java.util.Vector; import parser.VarList; import parser.ast.Declaration; import parser.ast.DeclarationIntUnbounded; import parser.ast.Expression; +import prism.OptionsIntervalIteration; import prism.Prism; import prism.PrismComponent; import prism.PrismException; import prism.PrismFileLog; +import prism.PrismNotSupportedException; import prism.PrismSettings; import prism.PrismUtils; import acceptance.AcceptanceReach; import acceptance.AcceptanceType; import automata.DA; import common.IntSet; +import common.PeriodicTimer; +import common.StopWatch; import common.IterableBitSet; import common.StopWatch; import explicit.LTLModelChecker.LTLProduct; +import explicit.modelviews.DTMCAlteredDistributions; +import explicit.modelviews.MDPFromDTMC; import explicit.rewards.MCRewards; +import explicit.rewards.MDPRewards; import explicit.rewards.Rewards; /** @@ -605,6 +615,10 @@ public class DTMCModelChecker extends ProbModelChecker mainLog.printWarning("Switching to linear equation solution method \"" + linEqMethod.fullName() + "\""); } + if (doIntervalIteration && (!precomp || !prob0 || !prob1)) { + throw new PrismNotSupportedException("Interval iteration requires precomputations to be active"); + } + // Start probabilistic reachability timer = System.currentTimeMillis(); mainLog.println("\nStarting probabilistic reachability..."); @@ -693,7 +707,11 @@ public class DTMCModelChecker extends ProbModelChecker throw new PrismException("Unknown linear equation solution method " + linEqMethod.fullName()); } - res = doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration()); + if (doIntervalIteration) { + res = doIntervalIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration()); + } else { + res = doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration()); + } // Finished probabilistic reachability timer = System.currentTimeMillis() - timer; @@ -1087,6 +1105,109 @@ public class DTMCModelChecker extends ProbModelChecker return doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, false); } + /** + * Compute reachability probabilities using power method (interval variant). + * @param dtmc The DTMC + * @param no Probability 0 states + * @param yes Probability 1 states + * @param init Optionally, an initial solution vector (will be overwritten), will be ignored if known == null + * @param known Optionally, a set of states for which the exact answer is known + * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values. + * @param topological do topological interval iteration? + */ + protected ModelCheckerResult doIntervalIterationReachProbs(DTMC dtmc, BitSet no, BitSet yes, double init[], BitSet known, IterationMethod iterationMethod, boolean topological) throws PrismException + { + BitSet unknown; + int i, n; + double initBelow[], initAbove[]; + long timer; + + // Start value iteration + timer = System.currentTimeMillis(); + String description = "with " + iterationMethod.getDescriptionShort(); + mainLog.println("Starting interval iteration (" + description + ")..."); + + ExportIterations iterationsExport = null; + if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) { + iterationsExport = new ExportIterations("Explicit DTMC ReachProbs interval iteration (" + description + ")"); + } + + // Store num states + n = dtmc.getNumStates(); + + // Create solution vector(s) + initBelow = (init == null) ? new double[n] : init; + initAbove = new double[n]; + + // Initialise solution vectors. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 1.0/0.0 if in yes/no; (3) initVal + // where initVal is 0.0 or 1.0, depending on whether we converge from below/above. + if (known != null && init != null) { + for (i = 0; i < n; i++) { + initBelow[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0; + initAbove[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0; + } + } else { + for (i = 0; i < n; i++) { + initBelow[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0; + initAbove[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0; + } + } + + // Determine set of states actually need to compute values for + unknown = new BitSet(); + unknown.set(0, n); + unknown.andNot(yes); + unknown.andNot(no); + if (known != null) + unknown.andNot(known); + + if (iterationsExport != null) { + iterationsExport.exportVector(initBelow, 0); + iterationsExport.exportVector(initAbove, 1); + } + + IntSet unknownStates = IntSet.asIntSet(unknown); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow(); + final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove(); + final boolean checkMonotonic = iiOptions.isCheckMonotonicity(); + + if (!enforceMonotonicFromAbove) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above."); + } + if (!enforceMonotonicFromBelow) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below."); + } + + IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultInterval(dtmc, true, enforceMonotonicFromBelow, checkMonotonic); + IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultInterval(dtmc, false, enforceMonotonicFromAbove, checkMonotonic); + + below.init(initBelow); + above.init(initAbove); + + if (topological) { + // Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering + // states in unknown + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get); + + IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> { + soln[s] = dtmc.mvMultJacSingle(s, soln); + }; + + // run the actual value iteration + return iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport); + } else { + // run the actual value iteration + return iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport); + } + + } + + + /** * Compute bounded reachability probabilities. * i.e. compute the probability of reaching a state in {@code target} within k steps. @@ -1199,6 +1320,381 @@ public class DTMCModelChecker extends ProbModelChecker return res; } + /** + * Compute upper bound for maximum expected reward, using the variant specified in the settings. + * @param dtmc the model + * @param mcRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinity states + * @return upper bound on R=?[ F target ] for all states + */ + double computeReachRewardsUpperBound(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + // inf and target states become trap states (with self-loops) + BitSet trapStates = (BitSet) target.clone(); + trapStates.or(inf); + DTMCAlteredDistributions cleanedDTMC = DTMCAlteredDistributions.addSelfLoops(dtmc, trapStates); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + double upperBound = 0.0; + String method = null; + switch (iiOptions.getBoundMethod()) { + case VARIANT_1_COARSE: + upperBound = computeReachRewardsUpperBoundVariant1Coarse(cleanedDTMC, mcRewards, target, unknown, inf); + method = "variant 1, coarse"; + break; + case VARIANT_1_FINE: + upperBound = computeReachRewardsUpperBoundVariant1Fine(cleanedDTMC, mcRewards, target, unknown, inf); + method = "variant 1, fine"; + break; + case VARIANT_2: + upperBound = computeReachRewardsUpperBoundVariant2(cleanedDTMC, mcRewards, target, unknown, inf); + method = "variant 2"; + break; + case DEFAULT: + case DSMPI: + { + MDP mdp = new MDPFromDTMC(cleanedDTMC); + MDPRewards mdpRewards = new MDPRewards() { + + @Override + public double getStateReward(int s) + { + return mcRewards.getStateReward(s); + } + + @Override + public double getTransitionReward(int s, int i) + { + return 0; + } + + @Override + public MDPRewards liftFromModel(Product product) + { + throw new RuntimeException("Unsupported"); + } + + @Override + public boolean hasTransitionRewards() + { + return false; + } + }; + upperBound = DijkstraSweepMPI.computeUpperBound(this, mdp, mdpRewards, target, unknown); + method = "Dijkstra Sweep MPI"; + break; + } + } + + if (method == null) { + throw new PrismException("Unknown upper bound heuristic"); + } + + mainLog.println("Upper bound for expectation (" + method + "): " + upperBound); + + return upperBound; + } + + /** + * Compute upper bound for maximum expected reward (variant 1, coarse), + * i.e., does not compute separate q_t / p_t per SCC. + * Uses Rs = S, i.e., does not take reachability into account. + * @param dtmc the model + * @param mcRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinity states + * @return upper bound on R=?[ F target ] for all states + */ + double computeReachRewardsUpperBoundVariant1Coarse(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()]; + int[] Ct = new int[dtmc.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, null); + BitSet trivial = new BitSet(); + + double q = 0; + for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) { + IntSet statesForSCC = sccs.getStatesForSCC(scc); + + int cardinality = statesForSCC.cardinality(); + + PrimitiveIterator.OfInt itSCC = statesForSCC.iterator(); + while (itSCC.hasNext()) { + int s = itSCC.nextInt(); + Ct[s] = cardinality; + + if (target.get(s) || inf.get(s)) { + // trap states + assert(cardinality == 1); + break; // continue with next SCC + } + + double probRemain = 0; + boolean allRemain = true; // all successors remain in the SCC? + boolean hasSelfloop = false; + for (Iterator> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) { + Entry t = it.next(); + if (statesForSCC.get(t.getKey())) { + probRemain += t.getValue(); + hasSelfloop = true; + } else { + allRemain = false; + } + } + + if (!allRemain) { // action in the set X + q = Math.max(q, probRemain); + } + + if (cardinality == 1 && !hasSelfloop) { + trivial.set(s); + } + } + } + + double p = 1; + for (int s = 0; s < dtmc.getNumStates(); s++) { + for (Iterator> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) { + Entry t = it.next(); + p = Math.min(p, t.getValue()); + } + } + + double upperBound = 0; + for (int s = 0; s < dtmc.getNumStates(); s++) { + if (target.get(s) || inf.get(s)) { + // inf or target states: not relevant, set visits to 0, ignore in summation + boundsOnExpectedVisits[s] = 0.0; + } else if (unknown.get(s)) { + if (trivial.get(s)) { + // s is a trivial SCC: seen at most once + boundsOnExpectedVisits[s] = 1.0; + } else { + boundsOnExpectedVisits[s] = 1 / (Math.pow(p, Ct[s]-1) * (1.0-q)); + } + + upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s); + } + } + + timer.stop(); + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 1, coarse):"); + mainLog.println("p = " + p); + mainLog.println("q = " + q); + mainLog.println("|Ct| = " + Arrays.toString(Ct)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + } + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + + /** + * Compute upper bound for maximum expected reward (variant 1, fine). + * i.e., does compute separate q_t / p_t per SCC. + * Uses Rs = S, i.e., does not take reachability into account. + * @param dtmc the model + * @param mcRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinity states + * @return upper bound on R=?[ F target ] for all states + */ + double computeReachRewardsUpperBoundVariant1Fine(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()]; + double[] qt = new double[dtmc.getNumStates()]; + double[] pt = new double[dtmc.getNumStates()]; + int[] Ct = new int[dtmc.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, null); + BitSet trivial = new BitSet(); + + for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) { + IntSet statesForSCC = sccs.getStatesForSCC(scc); + + double q = 0; + double p = 1; + int cardinality = statesForSCC.cardinality(); + + PrimitiveIterator.OfInt itSCC = statesForSCC.iterator(); + while (itSCC.hasNext()) { + int s = itSCC.nextInt(); + Ct[s] = cardinality; + + double probRemain = 0; + boolean allRemain = true; // all successors remain in the SCC? + boolean hasSelfloop = false; + for (Iterator> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) { + Entry t = it.next(); + if (statesForSCC.get(t.getKey())) { + probRemain += t.getValue(); + p = Math.min(p, t.getValue()); + hasSelfloop = true; + } else { + // outgoing edge + allRemain = false; + } + } + + if (!allRemain) { // action in the set Xt + q = Math.max(q, probRemain); + } + + if (cardinality == 1 && !hasSelfloop) { + trivial.set(s); + } + } + + for (int s : statesForSCC) { + qt[s] = q; + pt[s] = p; + } + } + + double upperBound = 0; + for (int s = 0; s < dtmc.getNumStates(); s++) { + if (target.get(s) || inf.get(s)) { + // inf or target states: not relevant, set visits to 0, ignore in summation + boundsOnExpectedVisits[s] = 0.0; + } else if (unknown.get(s)) { + if (trivial.get(s)) { + // s is a trivial SCC: seen at most once + boundsOnExpectedVisits[s] = 1.0; + } else { + if (pt[s] == 1.0) { + //throw new PrismException("Upper bound computation had p_t = 1 for state " + s); + } + boundsOnExpectedVisits[s] = 1 / (Math.pow(pt[s], Ct[s]-1) * (1.0-qt[s])); + } + + upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s); + } else { + throw new PrismException("Bogus arguments: inf/target/unknown should partition state space"); + } + } + + timer.stop(); + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 1, fine):"); + mainLog.println("pt = " + Arrays.toString(pt)); + mainLog.println("qt = " + Arrays.toString(qt)); + mainLog.println("|Ct| = " + Arrays.toString(Ct)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + } + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + + + /** + * Compute upper bound for maximum expected reward (variant 2). + * Uses Rs = S, i.e., does not take reachability into account. + * @param dtmc the model + * @param mcRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinity states + * @return upper bound on R=?[ F target ] for all states + */ + double computeReachRewardsUpperBoundVariant2(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] dt = new double[dtmc.getNumStates()]; + double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get); + + BitSet T = (BitSet) target.clone(); + + @SuppressWarnings("unused") + int i = 0; + while (true) { + BitSet Si = new BitSet(); + i++; + + // TODO: might be inefficient, worst-case quadratic runtime... + for (PrimitiveIterator.OfInt it = IterableBitSet.getClearBits(T, dtmc.getNumStates() -1 ).iterator(); it.hasNext(); ) { + int s = it.nextInt(); +// mainLog.println("Check " + s + " against " + T); + if (dtmc.someSuccessorsInSet(s, T)) { + Si.set(s); + } + } + + if (Si.isEmpty()) { + break; + } + + // mainLog.println("S" + i + " = " + Si); + // mainLog.println("T = " + T); + + for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(Si).iterator(); it.hasNext(); ) { + final int t = it.nextInt(); + final int sccIndexForT = sccs.getSCCIndex(t); + double d = dtmc.sumOverTransitions(t, (int __, int u, double prob) -> { + // mainLog.println("t = " + t + ", u = " + u + ", prob = " + prob); + if (!T.get(u)) + return 0.0; + + boolean inSameSCC = (sccs.getSCCIndex(u) == sccIndexForT); + double d_u_t = inSameSCC ? dt[u] : 1.0; + // mainLog.println("d_u_t = " + d_u_t); + return d_u_t * prob; + }); + dt[t] = d; + // mainLog.println("d["+t+"] = " + d); + } + + T.or(Si); + } + + double upperBound = 0; + for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(unknown).iterator(); it.hasNext();) { + int s = it.nextInt(); + boundsOnExpectedVisits[s] = 1 / dt[s]; + upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s); + } + + timer.stop(); + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 2):"); + mainLog.println("d_t = " + Arrays.toString(dt)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + } + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + + /** * Compute expected reachability rewards. * @param dtmc The DTMC @@ -1302,7 +1798,11 @@ public class DTMCModelChecker extends ProbModelChecker throw new PrismException("Unknown linear equation solution method " + linEqMethod.fullName()); } - res = doValueIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration()); + if (doIntervalIteration) { + res = doIntervalIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration()); + } else { + res = doValueIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration()); + } // Finished expected reachability timer = System.currentTimeMillis() - timer; @@ -1498,6 +1998,136 @@ public class DTMCModelChecker extends ProbModelChecker } } + /** + * Compute expected reachability rewards using interval iteration. + * @param dtmc The DTMC + * @param mcRewards The rewards + * @param target Target states + * @param inf States for which reward is infinite + * @param init Optionally, an initial solution vector (will be overwritten) + * @param known Optionally, a set of states for which the exact answer is known + * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values. + * @param topological do topological interval iteration? + */ + protected ModelCheckerResult doIntervalIterationReachRewards(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet inf, double init[], BitSet known, IterationMethod iterationMethod, boolean topological) + throws PrismException + { + BitSet unknown; + int i, n; + double init_below[], init_above[]; + long timer; + + // Store num states + n = dtmc.getNumStates(); + + // Determine set of states actually need to compute values for + unknown = new BitSet(); + unknown.set(0, n); + unknown.andNot(target); + unknown.andNot(inf); + if (known != null) + unknown.andNot(known); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + double upperBound; + if (iiOptions.hasManualUpperBound()) { + upperBound = iiOptions.getManualUpperBound(); + getLog().printWarning("Upper bound for interval iteration manually set to " + upperBound); + } else { + upperBound = computeReachRewardsUpperBound(dtmc, mcRewards, target, unknown, inf); + } + + double lowerBound; + if (iiOptions.hasManualLowerBound()) { + lowerBound = iiOptions.getManualLowerBound(); + getLog().printWarning("Lower bound for interval iteration manually set to " + lowerBound); + } else { + lowerBound = 0.0; + } + + // Start value iteration + timer = System.currentTimeMillis(); + String description = (topological ? "topological, " : "" ) + "with " + iterationMethod.getDescriptionShort(); + mainLog.println("Starting interval iteration (" + description + ") ..."); + + ExportIterations iterationsExport = null; + if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) { + iterationsExport = new ExportIterations("Explicit DTMC ReachRewards interval iteration (" + description + ") ..."); + } + + // Create solution vector(s) + init_below = (init == null) ? new double[n] : init; + init_above = new double[n]; + + // Initialise solution vector from below. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) lowerBound + if (init != null && known != null) { + for (i = 0; i < n; i++) + init_below[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound; + } else { + for (i = 0; i < n; i++) + init_below[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound; + } + + // Initialise solution vector from above. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) upperBound + if (init != null && known != null) { + for (i = 0; i < n; i++) + init_above[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound; + } else { + for (i = 0; i < n; i++) + init_above[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound; + } + + + if (iterationsExport != null) { + iterationsExport.exportVector(init_below, 0); + iterationsExport.exportVector(init_above, 1); + } + + IntSet unknownStates = IntSet.asIntSet(unknown); + + final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow(); + final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove(); + final boolean checkMonotonic = iiOptions.isCheckMonotonicity(); + + if (!enforceMonotonicFromAbove) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above."); + } + if (!enforceMonotonicFromBelow) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below."); + } + + IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultRewInterval(dtmc, mcRewards, true, enforceMonotonicFromBelow, checkMonotonic); + IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultRewInterval(dtmc, mcRewards, false, enforceMonotonicFromAbove, checkMonotonic); + + below.init(init_below); + above.init(init_above); + + ModelCheckerResult rv; + if (topological) { + // Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering + // states in unknown + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get); + + IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> { + soln[s] = dtmc.mvMultRewJacSingle(s, soln, mcRewards); + }; + + // run the actual value iteration + rv = iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport); + } else { + // run the actual value iteration + rv = iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport); + } + + double max_v = PrismUtils.findMaxFinite(rv.soln, unknownStates.iterator()); + mainLog.println("Maximum finite value in solution vector at end of interval iteration: " + max_v); + + return rv; + } + /** * Compute (forwards) steady-state probabilities * i.e. compute the long-run probability of being in each state, diff --git a/prism/src/explicit/DijkstraSweepMPI.java b/prism/src/explicit/DijkstraSweepMPI.java new file mode 100644 index 00000000..02832dfa --- /dev/null +++ b/prism/src/explicit/DijkstraSweepMPI.java @@ -0,0 +1,312 @@ +//============================================================================== +// +// Copyright (c) 2016- +// Authors: +// * Joachim Klein (TU Dresden) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package explicit; + +import java.util.BitSet; +import java.util.HashMap; +import java.util.HashSet; +import java.util.PriorityQueue; +import java.util.function.IntPredicate; + +import common.StopWatch; +import common.IterableBitSet; +import explicit.IncomingChoiceRelation.Choice; +import explicit.rewards.MDPRewards; +import prism.PrismComponent; + +/** + * An implementation of the upper bound computation for Rmin as detailed in the + * paper McMahan, Likhachev, Gordon "Bounded Real-Time Dynamic Programming: + * RTDP with monotone upper bounds and performance guarantees" (International + * Conference on Machine Learning, 2005). + * */ +public class DijkstraSweepMPI { + + private static class QueueEntry implements Comparable { + public int y; + public double p; + public double w; + + public QueueEntry(int y, double p, double w) + { + this.y = y; + this.p = p; + this.w = w; + } + + @Override + public int compareTo(QueueEntry o) + { + int r = Double.compare(p, o.p); + if (r == 0) { + return Double.compare(w, o.w); + } else { + return r; + } + } + } + + private static class ChoiceValues { + public double p; + public double w; + + public ChoiceValues(double p, double w) + { + this.p = p; + this.w = w; + } + } + + private static boolean debug = false; + private MDP mdp; + private MDPRewards rewards; + private PriorityQueue queue; + private double[] pState; + private double[] wState; + private HashMap choiceValues = new HashMap(); + private QueueEntry[] pri; + private int[] pi; + private BitSet unknown, target; + private BitSet fin = new BitSet(); + private IncomingChoiceRelation incoming; + private double lambda; + + private DijkstraSweepMPI(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown) + { + this.mdp = mdp; + this.unknown = unknown; + this.target = target; + this.rewards = rewards; + + incoming = IncomingChoiceRelation.forModel(parent, mdp); + + queue = new PriorityQueue(); + pState = new double[mdp.getNumStates()]; + wState = new double[mdp.getNumStates()]; + pri = new QueueEntry[mdp.getNumStates()]; + pi = new int[mdp.getNumStates()]; + + for (int s : IterableBitSet.getSetBits(unknown)) { + for (int choice = 0, numChoices = mdp.getNumChoices(s); choice < numChoices; choice++) { + Choice c = new Choice(s, choice); + double rew = rewards.getStateReward(s); + rew += rewards.getTransitionReward(s, choice); + choiceValues.put(c, new ChoiceValues(0.0, rew)); + } + } + + for (int s : IterableBitSet.getSetBits(target)) { + pState[s] = 1.0; + } + + HashSet preTarget = new HashSet(); + for (int t : IterableBitSet.getSetBits(target)) { + for (Choice c : incoming.getIncomingChoices(t)) { + boolean newChoice = preTarget.add(c); + if (newChoice) { + if (!unknown.get(c.getState())) { + continue; + } + if (!validChoice(c)) { + continue; + } + + update(c, target); + } + } + } + preTarget.clear(); + + sweep(); + + computeLambda(); + } + + private void sweep() + { + while (!queue.isEmpty()) { + int x = queue.poll().y; + if (fin.get(x)) { + // already handled + continue; + } + + fin.set(x); + ChoiceValues v = choiceValues.get(new Choice(x, pi[x])); + wState[x] = v.w; + pState[x] = v.p; + + for (Choice c : incoming.getIncomingChoices(x)) { + if (fin.get(c.getState())) { + // already handled, skip + continue; + } + + if (!unknown.get(c.getState())) { + // uninteresting state + continue; + } + + if (!validChoice(c)) { + // some successor go outside unknown U target (e.g., to some infinity or undefined state) + // skip + continue; + } + + // a relevant choice, update + update(c, x); + } + } + } + + private boolean validChoice(Choice choice) + { + IntPredicate outsideRelevant = (int t) -> { + if (unknown.get(t) || target.get(t)) return false; + return true; + }; + + return !mdp.someSuccessorsMatch(choice.getState(), choice.getChoice(), outsideRelevant); + } + + private void update(Choice choice, int x) + { + double w_x = wState[x]; + // compute P^a_yx * w(x) + double Pw = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> { + if (t != x) return 0.0; + return p * w_x; + }); + + double p_x = pState[x]; + // compute P^a_yx * p_g(x) + double Pp = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> { + if (t != x) return 0.0; + return p * p_x; + }); + + ChoiceValues c = choiceValues.get(choice); + assert(c != null); + c.p += Pp; + c.w += Pw; + + QueueEntry newPri = new QueueEntry(choice.getState(), 1 - c.p, c.w); + if (pri[choice.getState()] == null || newPri.compareTo(pri[choice.getState()]) < 0) { + pri[choice.getState()] = newPri; + pi[choice.getState()] = choice.getChoice(); + queue.add(newPri); + } + } + + private void update(Choice choice, BitSet target) + { + // compute P^a_y->target + double Pp = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> { + if (target.get(t)) return p; + return 0.0; + }); + + ChoiceValues c = choiceValues.get(choice); + c.p += Pp; + + QueueEntry newPri = new QueueEntry(choice.getState(), 1 - c.p, c.w); + if (pri[choice.getState()] == null || newPri.compareTo(pri[choice.getState()]) < 0) { + pri[choice.getState()] = newPri; + pi[choice.getState()] = choice.getChoice(); + queue.add(newPri); + } + } + + private double computeLambda() + { + lambda = 0.0; + + for (int x : IterableBitSet.getSetBits(unknown)) { + int a = pi[x]; + double lambda_x_a = Double.POSITIVE_INFINITY; + + // check condition (I) + double I_sum = mdp.sumOverTransitions(x, a, (int s, int t, double p) -> { + return p * pState[t]; + }); + + if (pState[x] < I_sum) { + // condition (I) holds + double den = rewards.getStateReward(x) + rewards.getTransitionReward(x, a); // c(x,a) + den += mdp.sumOverTransitions(x, a, (int s, int t, double p) -> { + return p * wState[t]; + }); + den -= wState[x]; + + double num = mdp.sumOverTransitions(x, a, (int s, int t, double p) -> { + return p * pState[t]; + }); + num -= pState[x]; + + lambda_x_a = den / num; + } else { + // TODO: check condition (II) + lambda_x_a = 0; + } + + lambda = Double.max(lambda, lambda_x_a); + } + + return lambda; + } + + public static double[] computeUpperBounds(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown) + { + StopWatch timer = new StopWatch(parent.getLog()); + timer.start("computing upper bound(s) for Rmin using the DSI-MP algorithm"); + + parent.getLog().println("Computing upper bound(s) for Rmin using the Dijkstra Sweep for Monotone Pessimistic Initialization algorithm..."); + double[] upperBounds = new double[mdp.getNumStates()]; + DijkstraSweepMPI dsmpi = new DijkstraSweepMPI(parent, mdp, rewards, target, unknown); + + for (int x : IterableBitSet.getSetBits(unknown)) { + upperBounds[x] = dsmpi.wState[x] + dsmpi.lambda*(1 - dsmpi.pState[x]); + } + + if (debug) { + parent.getLog().println(upperBounds); + } + + timer.stop(); + return upperBounds; + } + + public static double computeUpperBound(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown) + { + double bound = 0.0; + final double[] upperBounds = computeUpperBounds(parent, mdp, rewards, target, unknown); + for (int s : IterableBitSet.getSetBits(unknown)) { + bound = Double.max(bound, upperBounds[s]); + } + return bound; + } +} diff --git a/prism/src/explicit/MDPModelChecker.java b/prism/src/explicit/MDPModelChecker.java index 00b72e98..5c2df32e 100644 --- a/prism/src/explicit/MDPModelChecker.java +++ b/prism/src/explicit/MDPModelChecker.java @@ -31,6 +31,8 @@ import java.util.BitSet; import java.util.Iterator; import java.util.List; import java.util.Map; +import java.util.Map.Entry; +import java.util.PrimitiveIterator; import java.util.Vector; import common.IterableStateSet; @@ -39,6 +41,7 @@ import parser.VarList; import parser.ast.Declaration; import parser.ast.DeclarationIntUnbounded; import parser.ast.Expression; +import prism.OptionsIntervalIteration; import prism.Prism; import prism.PrismComponent; import prism.PrismDevNullLog; @@ -56,6 +59,7 @@ import automata.LTL2WDBA; import common.IntSet; import common.IterableBitSet; import explicit.modelviews.EquivalenceRelationInteger; +import explicit.modelviews.MDPDroppedAllChoices; import explicit.modelviews.MDPEquiv; import explicit.rewards.MCRewards; import explicit.rewards.MCRewardsFromMDPRewards; @@ -357,6 +361,23 @@ public class MDPModelChecker extends ProbModelChecker if (!min) throw new PrismException("Value iteration from above only works for minimum probabilities"); } + if (doIntervalIteration) { + if (!min && (genStrat || exportAdv)) { + throw new PrismNotSupportedException("Currently, explicit engine does not support adversary construction for interval iteration and Pmax"); + } + if (mdpSolnMethod != MDPSolnMethod.VALUE_ITERATION && mdpSolnMethod != MDPSolnMethod.GAUSS_SEIDEL) { + throw new PrismNotSupportedException("Currently, explicit engine only supports interval iteration with value iteration or Gauss-Seidel for MDPs"); + } + if (init != null) + throw new PrismNotSupportedException("Interval iteration currently not supported with provided initial values"); + if (!(precomp && prob0 && prob1)) { + throw new PrismNotSupportedException("Precomputations (Prob0 & Prob1) must be enabled for interval iteration"); + } + + if (!min) { + doPmaxQuotient = true; + } + } if (mdpSolnMethod == MDPSolnMethod.POLICY_ITERATION || mdpSolnMethod == MDPSolnMethod.MODIFIED_POLICY_ITERATION) { if (known != null) { throw new PrismException("Policy iteration methods cannot be passed 'known' values for some states"); @@ -535,9 +556,15 @@ public class MDPModelChecker extends ProbModelChecker iterationMethod = new IterationMethodGS(termCrit == TermCrit.ABSOLUTE, termCritParam, false); break; case POLICY_ITERATION: + if (doIntervalIteration) { + throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration"); + } res = computeReachProbsPolIter(mdp, no, yes, min, strat); break; case MODIFIED_POLICY_ITERATION: + if (doIntervalIteration) { + throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration"); + } res = computeReachProbsModPolIter(mdp, no, yes, min, strat); break; default: @@ -545,7 +572,11 @@ public class MDPModelChecker extends ProbModelChecker } if (res == null) { // not yet computed, use iterationMethod - res = doValueIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat); + if (!doIntervalIteration) { + res = doValueIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat); + } else { + res = doIntervalIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat); + } } return res; @@ -850,6 +881,113 @@ public class MDPModelChecker extends ProbModelChecker } } + /** + * Compute reachability probabilities using interval iteration. + * Optionally, store optimal (memoryless) strategy info. + * @param mdp The MDP + * @param no Probability 0 states + * @param yes Probability 1 states + * @param min Min or max probabilities (true=min, false=max) + * @param init Optionally, an initial solution vector (will be overwritten) + * @param known Optionally, a set of states for which the exact answer is known + * @param iterationMethod The iteration method + * @param topological Do topological value iteration? + * @param strat Storage for (memoryless) strategy choice indices (ignored if null) + * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values. + */ + protected ModelCheckerResult doIntervalIterationReachProbs(MDP mdp, BitSet no, BitSet yes, boolean min, double init[], BitSet known, IterationMethod iterationMethod, boolean topological, int strat[]) + throws PrismException + { + BitSet unknown; + int i, n; + double initBelow[], initAbove[]; + long timer; + + // Start value iteration + timer = System.currentTimeMillis(); + String description = (min ? "min" : "max") + + (topological ? ", topological": "" ) + + ", with " + iterationMethod.getDescriptionShort(); + + mainLog.println("Starting interval iteration (" + description + ")..."); + + ExportIterations iterationsExport = null; + if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) { + iterationsExport = new ExportIterations("Explicit MDP ReachProbs interval iteration (" + description + ")"); + } + + // Store num states + n = mdp.getNumStates(); + + // Create solution vector(s) + initBelow = (init == null) ? new double[n] : init; + initAbove = new double[n]; + + // Initialise solution vectors. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 1.0/0.0 if in yes/no; (3) initVal + // where initVal is 0.0 or 1.0, depending on whether we converge from below/above. + if (known != null && init != null) { + for (i = 0; i < n; i++) { + initBelow[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0; + initAbove[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0; + } + } else { + for (i = 0; i < n; i++) { + initBelow[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0; + initAbove[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0; + } + } + + // Determine set of states actually need to compute values for + unknown = new BitSet(); + unknown.set(0, n); + unknown.andNot(yes); + unknown.andNot(no); + if (known != null) + unknown.andNot(known); + + if (iterationsExport != null) { + iterationsExport.exportVector(initBelow, 0); + iterationsExport.exportVector(initAbove, 1); + } + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow(); + final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove(); + final boolean checkMonotonic = iiOptions.isCheckMonotonicity(); + + if (!enforceMonotonicFromAbove) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above."); + } + if (!enforceMonotonicFromBelow) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below."); + } + + IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultMinMaxInterval(mdp, min, strat, true, enforceMonotonicFromBelow, checkMonotonic); + IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultMinMaxInterval(mdp, min, strat, false, enforceMonotonicFromAbove, checkMonotonic); + below.init(initBelow); + above.init(initAbove); + + IntSet unknownStates = IntSet.asIntSet(unknown); + + if (topological) { + // Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering + // states in unknown + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get); + + IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> { + soln[s] = mdp.mvMultJacMinMaxSingle(s, soln, min, strat); + }; + + // run the actual value iteration + return iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport); + } else { + // run the actual value iteration + return iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport); + } + } + /** * Compute reachability probabilities using Gauss-Seidel (including Jacobi-style updates). * @param mdp The MDP @@ -1255,6 +1393,453 @@ public class MDPModelChecker extends ProbModelChecker return res; } + /** + * Compute upper bound for maximum expected reward, with the method specified in the settings. + * @param mdp the model + * @param mdpRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinite states + * @return upper bound on Rmax=?[ F target ] for all states + */ + double computeReachRewardsMaxUpperBound(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + // inf and target states become trap states (with dropped choices) + BitSet trapStates = (BitSet) target.clone(); + trapStates.or(inf); + MDP cleanedMDP = new MDPDroppedAllChoices(mdp, trapStates); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + double upperBound = 0.0; + String method = null; + switch (iiOptions.getBoundMethod()) { + case VARIANT_1_COARSE: + upperBound = computeReachRewardsMaxUpperBoundVariant1Coarse(cleanedMDP, mdpRewards, target, unknown, inf); + method = "variant 1, coarse"; + break; + case VARIANT_1_FINE: + upperBound = computeReachRewardsMaxUpperBoundVariant1Fine(cleanedMDP, mdpRewards, target, unknown, inf); + method = "variant 1, fine"; + break; + case DEFAULT: + case VARIANT_2: + upperBound = computeReachRewardsMaxUpperBoundVariant2(cleanedMDP, mdpRewards, target, unknown, inf); + method = "variant 2"; + break; + case DSMPI: + throw new PrismNotSupportedException("Dijkstra Sweep MPI upper bound heuristic can not be used for Rmax"); + } + + if (method == null) { + throw new PrismException("Unknown upper bound heuristic"); + } + + mainLog.println("Upper bound for max expectation (" + method + "): " + upperBound); + return upperBound; + } + + /** + * Compute upper bound for minimum expected reward, with the method specified in the settings. + * @param mdp the model + * @param mdpRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinite states + * @return upper bound on Rmin=?[ F target ] for all unknown states + */ + double computeReachRewardsMinUpperBound(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + // inf and target states become trap states (with dropped choices) + BitSet trapStates = (BitSet) target.clone(); + trapStates.or(inf); + MDP cleanedMDP = new MDPDroppedAllChoices(mdp, trapStates); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + double upperBound = 0.0; + String method = null; + switch (iiOptions.getBoundMethod()) { + case DEFAULT: + case DSMPI: + upperBound = DijkstraSweepMPI.computeUpperBound(this, mdp, mdpRewards, target, unknown); + method = "Dijkstra Sweep MPI"; + break; + case VARIANT_1_COARSE: + upperBound = computeReachRewardsMaxUpperBoundVariant1Coarse(cleanedMDP, mdpRewards, target, unknown, inf); + method = "using Rmax upper bound via variant 1, coarse"; + break; + case VARIANT_1_FINE: + upperBound = computeReachRewardsMaxUpperBoundVariant1Fine(cleanedMDP, mdpRewards, target, unknown, inf); + method = "using Rmax upper bound via variant 1, fine"; + break; + case VARIANT_2: + upperBound = computeReachRewardsMaxUpperBoundVariant2(cleanedMDP, mdpRewards, target, unknown, inf); + method = "using Rmax upper bound via variant 2"; + break; + } + + if (method == null) { + throw new PrismException("Unknown upper bound heuristic"); + } + + mainLog.println("Upper bound for min expectation (" + method + "): " + upperBound); + return upperBound; + } + + /** + * Return true if the MDP is contracting for all states in the 'unknown' + * set, i.e., if Pmin=1( unknown U target) holds. + */ + private boolean isContracting(MDP mdp, BitSet unknown, BitSet target) + { + // compute Pmin=1( unknown U target ) + BitSet pmin1 = prob1(mdp, unknown, target, true, null); + BitSet tmp = (BitSet) unknown.clone(); + tmp.andNot(pmin1); + if (!tmp.isEmpty()) { + // unknown is not contained in pmin1, not contracting + return false; + } + return true; + } + + /** + * Compute upper bound for maximum expected reward (variant 1, coarse), + * i.e., does not compute separate q_t / p_t per SCC. + * Uses Rs = S, i.e., does not take reachability into account. + * @param mdp the model + * @param mdpRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @return upper bound on Rmax=?[ F target ] for all states + */ + double computeReachRewardsMaxUpperBoundVariant1Coarse(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] boundsOnExpectedVisits = new double[mdp.getNumStates()]; + double[] maxRews = new double[mdp.getNumStates()]; + int[] Ct = new int[mdp.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for maximal expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, null); + BitSet trivial = new BitSet(); + + double q = 0; + for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) { + IntSet statesForSCC = sccs.getStatesForSCC(scc); + + int cardinality = statesForSCC.cardinality(); + + PrimitiveIterator.OfInt itSCC = statesForSCC.iterator(); + while (itSCC.hasNext()) { + int s = itSCC.nextInt(); + Ct[s] = cardinality; + + boolean hasSelfloop = false; + for (int ch = 0; ch < mdp.getNumChoices(s); ch++) { + double probRemain = 0; + boolean allRemain = true; // all successors remain in the SCC? + for (Iterator> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) { + Entry t = it.next(); + if (statesForSCC.get(t.getKey())) { + probRemain += t.getValue(); + hasSelfloop = true; + } else { + allRemain = false; + } + } + + if (!allRemain) { // action in the set X + q = Math.max(q, probRemain); + } + } + + if (cardinality == 1 && !hasSelfloop) { + trivial.set(s); + } + } + } + + double p = 1; + for (int s = 0; s < mdp.getNumStates(); s++) { + double maxRew = 0; + for (int ch = 0; ch < mdp.getNumChoices(s); ch++) { + for (Iterator> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) { + Entry t = it.next(); + p = Math.min(p, t.getValue()); + + double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch); + maxRew = Math.max(maxRew, rew); + } + } + maxRews[s] = maxRew; + } + + double upperBound = 0; + for (int s = 0; s < mdp.getNumStates(); s++) { + if (target.get(s) || inf.get(s)) { + // inf or target states: not relevant, set visits to 0, ignore in summation + boundsOnExpectedVisits[s] = 0.0; + } else if (unknown.get(s)) { + if (trivial.get(s)) { + // s is a trivial SCC: seen at most once + boundsOnExpectedVisits[s] = 1.0; + } else { + boundsOnExpectedVisits[s] = 1 / (Math.pow(p, Ct[s]-1) * (1.0-q)); + } + upperBound += boundsOnExpectedVisits[s] * maxRews[s]; + } + } + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 1, coarse):"); + mainLog.println("p = " + p); + mainLog.println("q = " + q); + mainLog.println("|Ct| = " + Arrays.toString(Ct)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + mainLog.println("maxRews = " + Arrays.toString(maxRews)); + } + + timer.stop(); + // mainLog.println("Upper bound for max expectation (variant 1, coarse): " + upperBound); + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + + /** + * Compute upper bound for maximum expected reward (variant 1, fine). + * i.e., does compute separate q_t / p_t per SCC. + * Uses Rs = S, i.e., does not take reachability into account. + * @param mdp the model + * @param mdpRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @return upper bound on Rmax=?[ F target ] for all states + */ + double computeReachRewardsMaxUpperBoundVariant1Fine(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] boundsOnExpectedVisits = new double[mdp.getNumStates()]; + double[] qt = new double[mdp.getNumStates()]; + double[] pt = new double[mdp.getNumStates()]; + double[] maxRews = new double[mdp.getNumStates()]; + int[] Ct = new int[mdp.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for maximal expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, null); + BitSet trivial = new BitSet(); + + for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) { + IntSet statesForSCC = sccs.getStatesForSCC(scc); + + double q = 0; + double p = 1; + + int cardinality = statesForSCC.cardinality(); + + PrimitiveIterator.OfInt itSCC = statesForSCC.iterator(); + while (itSCC.hasNext()) { + int s = itSCC.nextInt(); + + Ct[s] = cardinality; + boolean hasSelfloop = false; + + for (int ch = 0; ch < mdp.getNumChoices(s); ch++) { + + double probRemain = 0; + boolean allRemain = true; // all successors remain in the SCC? + for (Iterator> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) { + Entry t = it.next(); + if (statesForSCC.get(t.getKey())) { + probRemain += t.getValue(); + p = Math.min(p, t.getValue()); + hasSelfloop = true; + } else { + allRemain = false; + } + } + + if (!allRemain) { // action in the set Xt + q = Math.max(q, probRemain); + } + } + + if (cardinality == 1 && !hasSelfloop) { + trivial.set(s); + } + } + + for (int s : statesForSCC) { + qt[s] = q; + pt[s] = p; + } + } + + for (int s = 0; s < mdp.getNumStates(); s++) { + double maxRew = 0; + for (int ch = 0; ch < mdp.getNumChoices(s); ch++) { + double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch); + maxRew = Math.max(maxRew, rew); + } + maxRews[s] = maxRew; + } + + double upperBound = 0; + for (int s = 0; s < mdp.getNumStates(); s++) { + if (target.get(s) || inf.get(s)) { + // inf or target states: not relevant, set visits to 0, ignore in summation + boundsOnExpectedVisits[s] = 0.0; + } else if (unknown.get(s)) { + if (trivial.get(s)) { + // s is a trivial SCC: seen at most once + boundsOnExpectedVisits[s] = 1.0; + } else { + boundsOnExpectedVisits[s] = 1 / (Math.pow(pt[s], Ct[s]-1) * (1.0-qt[s])); + } + upperBound += boundsOnExpectedVisits[s] * maxRews[s]; + } + } + + timer.stop(); + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 1, fine):"); + mainLog.println("pt = " + Arrays.toString(pt)); + mainLog.println("qt = " + Arrays.toString(qt)); + mainLog.println("|Ct| = " + Arrays.toString(Ct)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + mainLog.println("maxRews = " + Arrays.toString(maxRews)); + } + + // mainLog.println("Upper bound for max expectation (variant 1, fine): " + upperBound); + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + + /** + * Compute upper bound for maximum expected reward (variant 2). + * Uses Rs = S, i.e., does not take reachability into account. + * @param dtmc the model + * @param mcRewards the rewards + * @param target the target states + * @param unknown the states that are not target or infinity states + * @param inf the infinity states + * @return upper bound on R=?[ F target ] for all states + */ + double computeReachRewardsMaxUpperBoundVariant2(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException + { + double[] dt = new double[mdp.getNumStates()]; + double[] boundsOnExpectedVisits = new double[mdp.getNumStates()]; + double[] maxRews = new double[mdp.getNumStates()]; + + StopWatch timer = new StopWatch(getLog()); + timer.start("computing an upper bound for expected reward"); + + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get); + + BitSet T = (BitSet) target.clone(); + + @SuppressWarnings("unused") + int i = 0; + while (true) { + BitSet Si = new BitSet(); + i++; + + // TODO: might be inefficient, worst-case quadratic runtime... + for (PrimitiveIterator.OfInt it = IterableBitSet.getClearBits(T, mdp.getNumStates() -1 ).iterator(); it.hasNext(); ) { + int s = it.nextInt(); + // mainLog.println("Check " + s + " against " + T); + boolean allActionsReachT = true; + for (int choice = 0, choices = mdp.getNumChoices(s); choice < choices; choice++) { + if (!mdp.someSuccessorsInSet(s, choice, T)) { + allActionsReachT = false; + break; + } + } + if (allActionsReachT) { + Si.set(s); + } + } + + if (Si.isEmpty()) { + break; + } + + // mainLog.println("S" + i + " = " + Si); + // mainLog.println("T = " + T); + + for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(Si).iterator(); it.hasNext(); ) { + final int t = it.nextInt(); + final int sccIndexForT = sccs.getSCCIndex(t); + + double min = Double.POSITIVE_INFINITY; + for (int choice = 0, choices = mdp.getNumChoices(t); choice < choices; choice++) { + // mainLog.println("State " + t + ", choice = " + choice); + double d = mdp.sumOverTransitions(t, choice, (int __, int u, double prob) -> { + // mainLog.println("t = " + t + ", u = " + u + ", prob = " + prob); + if (!T.get(u)) + return 0.0; + + boolean inSameSCC = (sccs.getSCCIndex(u) == sccIndexForT); + double d_u_t = inSameSCC ? dt[u] : 1.0; + // mainLog.println("d_u_t = " + d_u_t); + return d_u_t * prob; + }); + if (d < min) { + min = d; + } + } + dt[t] = min; + // mainLog.println("d["+t+"] = " + dt[t]); + } + + T.or(Si); + } + + for (int s = 0; s < mdp.getNumStates(); s++) { + double maxRew = 0; + for (int ch = 0; ch < mdp.getNumChoices(s); ch++) { + double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch); + maxRew = Math.max(maxRew, rew); + } + maxRews[s] = maxRew; + } + + double upperBound = 0; + for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(unknown).iterator(); it.hasNext();) { + int s = it.nextInt(); + boundsOnExpectedVisits[s] = 1 / dt[s]; + upperBound += boundsOnExpectedVisits[s] * maxRews[s]; + } + + timer.stop(); + + if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) { + mainLog.println("Upper bound for max expectation computation (variant 2):"); + mainLog.println("d_t = " + Arrays.toString(dt)); + mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits)); + } + + // mainLog.println("Upper bound for expectation (variant 2): " + upperBound); + + if (!Double.isFinite(upperBound)) { + throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result"); + } + + return upperBound; + } + /** * Compute expected instantaneous reward, * i.e. compute the min/max expected reward of the states after {@code k} steps. @@ -1477,7 +2062,12 @@ public class MDPModelChecker extends ProbModelChecker throw new PrismException("Policy iteration methods cannot be passed 'known' values for some states"); } } - + if (doIntervalIteration) { + if (mdpSolnMethod != MDPSolnMethod.VALUE_ITERATION && mdpSolnMethod != MDPSolnMethod.GAUSS_SEIDEL) { + throw new PrismNotSupportedException("Currently, explicit engine only supports interval iteration with value iteration or Gauss-Seidel for MDPs"); + } + } + // Start expected reachability timer = System.currentTimeMillis(); mainLog.println("\nStarting expected reachability (" + (min ? "min" : "max") + ")..."); @@ -1624,6 +2214,9 @@ public class MDPModelChecker extends ProbModelChecker iterationMethod = new IterationMethodGS(termCrit == TermCrit.ABSOLUTE, termCritParam, false); break; case POLICY_ITERATION: + if (doIntervalIteration) { + throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration"); + } res = computeReachRewardsPolIter(mdp, mdpRewards, target, inf, min, strat); break; default: @@ -1631,7 +2224,11 @@ public class MDPModelChecker extends ProbModelChecker } if (res == null) { // not yet computed, use iterationMethod - res = doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat); + if (!doIntervalIteration) { + res = doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat); + } else { + res = doIntervalIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat); + } } return res; @@ -1760,6 +2357,149 @@ public class MDPModelChecker extends ProbModelChecker return doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, false, strat); } + /** + * Compute expected reachability rewards using interval iteration + * Optionally, store optimal (memoryless) strategy info. + * @param mdp The MDP + * @param mdpRewards The rewards + * @param target Target states + * @param inf States for which reward is infinite + * @param min Min or max rewards (true=min, false=max) + * @param init Optionally, an initial solution vector (will be overwritten) + * @param known Optionally, a set of states for which the exact answer is known + * @param topological do topological interval iteration + * @param strat Storage for (memoryless) strategy choice indices (ignored if null) + * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values. + */ + protected ModelCheckerResult doIntervalIterationReachRewards(MDP mdp, MDPRewards mdpRewards, IterationMethod iterationMethod, BitSet target, BitSet inf, boolean min, double init[], BitSet known, boolean topological, int strat[]) + throws PrismException + { + BitSet unknown; + int i, n; + double initBelow[], initAbove[]; + long timer; + + // Store num states + n = mdp.getNumStates(); + + // Determine set of states actually need to compute values for + unknown = new BitSet(); + unknown.set(0, n); + unknown.andNot(target); + unknown.andNot(inf); + if (known != null) + unknown.andNot(known); + + OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this); + + double upperBound; + if (iiOptions.hasManualUpperBound()) { + upperBound = iiOptions.getManualUpperBound(); + getLog().printWarning("Upper bound for interval iteration manually set to " + upperBound); + } else { + if (min) { + upperBound = computeReachRewardsMinUpperBound(mdp, mdpRewards, target, unknown, inf); + } else { + upperBound = computeReachRewardsMaxUpperBound(mdp, mdpRewards, target, unknown, inf); + } + } + + double lowerBound; + if (iiOptions.hasManualLowerBound()) { + lowerBound = iiOptions.getManualLowerBound(); + getLog().printWarning("Lower bound for interval iteration manually set to " + lowerBound); + } else { + lowerBound = 0.0; + } + + if (min) { + if (!isContracting(mdp, unknown, target)) { + throw new PrismNotSupportedException("Interval iteration for Rmin and non-contracting MDP currently not supported"); + } else { + mainLog.println("Relevant sub-MDP is contracting, proceed..."); + } + } + + // Start value iteration + timer = System.currentTimeMillis(); + String description = (min ? "min" : "max") + (topological ? ", topological" : "") + ", with " + iterationMethod.getDescriptionShort(); + mainLog.println("Starting interval iteration (" + description + ")..."); + + ExportIterations iterationsExport = null; + if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) { + iterationsExport = new ExportIterations("Explicit MDP ReachRewards interval iteration (" + description + ")"); + } + + // Create initial solution vector(s) + initBelow = (init == null) ? new double[n] : init; + initAbove = new double[n]; + + // Initialise solution vector from below. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) lowerBound + if (init != null && known != null) { + for (i = 0; i < n; i++) + initBelow[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound; + } else { + for (i = 0; i < n; i++) + initBelow[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound; + } + + // Initialise solution vector from above. Use (where available) the following in order of preference: + // (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) upperBound + if (init != null && known != null) { + for (i = 0; i < n; i++) + initAbove[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound; + } else { + for (i = 0; i < n; i++) + initAbove[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound; + } + + if (iterationsExport != null) { + iterationsExport.exportVector(initBelow, 0); + iterationsExport.exportVector(initAbove, 1); + } + + final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow(); + final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove(); + final boolean checkMonotonic = iiOptions.isCheckMonotonicity(); + + if (!enforceMonotonicFromAbove) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above."); + } + if (!enforceMonotonicFromBelow) { + getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below."); + } + + IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultRewMinMaxInterval(mdp, mdpRewards, min, strat, true, enforceMonotonicFromBelow, checkMonotonic); + IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultRewMinMaxInterval(mdp, mdpRewards, min, strat, false, enforceMonotonicFromAbove, checkMonotonic); + below.init(initBelow); + above.init(initAbove); + + IntSet unknownStates = IntSet.asIntSet(unknown); + + ModelCheckerResult rv; + if (topological) { + // Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering + // states in unknown + SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get); + + IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> { + soln[s] = mdp.mvMultRewJacMinMaxSingle(s, soln, mdpRewards, min, strat); + }; + + // run the actual value iteration + rv = iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport); + } else { + // run the actual value iteration + rv = iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport); + } + + double max_v = PrismUtils.findMaxFinite(rv.soln, unknownStates.iterator()); + mainLog.println("Maximum finite value in solution vector at end of interval iteration: " + max_v); + + return rv; + } + /** * Compute expected reachability rewards using policy iteration. * The array {@code strat} is used both to pass in the initial strategy for policy iteration,