diff --git a/prism/src/explicit/DTMCModelChecker.java b/prism/src/explicit/DTMCModelChecker.java index 33c985c2..85980d85 100644 --- a/prism/src/explicit/DTMCModelChecker.java +++ b/prism/src/explicit/DTMCModelChecker.java @@ -154,7 +154,7 @@ public class DTMCModelChecker extends ProbModelChecker protected StateValues checkRewardFormula(Model model, MCRewards modelRewards, Expression expr) throws PrismException { StateValues rewards = null; - + if (expr instanceof ExpressionTemporal) { ExpressionTemporal exprTemp = (ExpressionTemporal) expr; switch (exprTemp.getOperator()) { @@ -165,7 +165,7 @@ public class DTMCModelChecker extends ProbModelChecker throw new PrismException("Explicit engine does not yet handle the " + exprTemp.getOperatorSymbol() + " operator in the R operator"); } } - + if (rewards == null) throw new PrismException("Unrecognised operator in R operator"); @@ -198,43 +198,40 @@ public class DTMCModelChecker extends ProbModelChecker /** * Compute steady-state probability distribution (forwards). - * Optionally, use the passed in vector initDist as the initial probability distribution (time step 0). + * Start from initial state (or uniform distribution over multiple initial states). + */ + public StateValues doSteadyState(DTMC dtmc) throws PrismException + { + return doSteadyState(dtmc, (StateValues) null); + } + + /** + * Compute steady-state probability distribution (forwards). + * Optionally, use the passed in file initDistFile to give the initial probability distribution (time 0). + * If null, start from initial state (or uniform distribution over multiple initial states). + */ + public StateValues doSteadyState(DTMC dtmc, File initDistFile) throws PrismException + { + StateValues initDist = readDistributionFromFile(initDistFile, dtmc); + return doSteadyState(dtmc, initDist); + } + + /** + * Compute steady-state probability distribution (forwards). + * Optionally, use the passed in vector initDist as the initial probability distribution (time 0). * If null, start from initial state (or uniform distribution over multiple initial states). * For reasons of efficiency, when a vector is passed in, it will be trampled over, * so if you wanted it, take a copy. * @param dtmc The DTMC * @param initDist Initial distribution (will be overwritten) */ - public StateValues doSteadyState(DTMC dtmc, double initDist[]) throws PrismException + public StateValues doSteadyState(DTMC dtmc, StateValues initDist) throws PrismException { - ModelCheckerResult res = null; - int n; - double initDistNew[] = null; - StateValues probs = null; - - // TODO: BSCC computation - - // Store num states - n = dtmc.getNumStates(); - - // Build initial distribution (if not specified) - if (initDist == null) { - initDistNew = new double[n]; - for (int in : dtmc.getInitialStates()) { - initDistNew[in] = 1 / dtmc.getNumInitialStates(); - } - } else { - initDistNew = initDist; - throw new PrismException("Not implemented yet"); // TODO - } - - // Compute transient probabilities - res = computeSteadyStateProbs(dtmc, initDistNew); - probs = StateValues.createFromDoubleArray(res.soln, dtmc); - - return probs; + StateValues initDistNew = (initDist == null) ? buildInitialDistribution(dtmc) : initDist; + ModelCheckerResult res = computeSteadyStateProbs(dtmc, initDistNew.getDoubleArray()); + return StateValues.createFromDoubleArray(res.soln, dtmc); } - + /** * Compute transient probability distribution (forwards). * Optionally, use the passed in vector initDist as the initial probability distribution (time step 0). @@ -249,7 +246,9 @@ public class DTMCModelChecker extends ProbModelChecker { throw new PrismException("Not implemented yet"); } - + + // Utility methods for probability distributions + /** * Generate a probability distribution, stored as a StateValues object, from a file. * If {@code distFile} is null, so is the return value. @@ -269,6 +268,27 @@ public class DTMCModelChecker extends ProbModelChecker return dist; } + /** + * Build a probability distribution, stored as a StateValues object, + * from the initial states info of the current model: either probability 1 for + * the (single) initial state or equiprobable over multiple initial states. + * The type of storage (MTBDD or double vector) matches the current engine. + */ + public StateValues buildInitialDistribution(Model model) throws PrismException + { + StateValues dist = null; + + // Build an empty vector + dist = new StateValues(TypeDouble.getInstance(), model); + // Populate vector (equiprobable over initial states) + double d = 1.0 / model.getNumInitialStates(); + for (int in : model.getInitialStates()) { + dist.setDoubleValue(in, d); + } + + return dist; + } + // Numerical computation functions /** @@ -320,7 +340,7 @@ public class DTMCModelChecker extends ProbModelChecker linEqMethod = LinEqMethod.GAUSS_SEIDEL; mainLog.printWarning("Switching to linear equation solution method \"" + linEqMethod.fullName() + "\""); } - + // Start probabilistic reachability timer = System.currentTimeMillis(); mainLog.println("Starting probabilistic reachability..."); @@ -607,7 +627,7 @@ public class DTMCModelChecker extends ProbModelChecker msg += "\nConsider using a different numerical method or increasing the maximum number of iterations"; throw new PrismException(msg); } - + // Return results res = new ModelCheckerResult(); res.soln = soln; @@ -691,7 +711,7 @@ public class DTMCModelChecker extends ProbModelChecker msg += "\nConsider using a different numerical method or increasing the maximum number of iterations"; throw new PrismException(msg); } - + // Return results res = new ModelCheckerResult(); res.soln = soln; @@ -740,7 +760,7 @@ public class DTMCModelChecker extends ProbModelChecker public ModelCheckerResult computeBoundedReachProbs(DTMC dtmc, BitSet remain, BitSet target, int k, double init[], double results[]) throws PrismException { // TODO: implement until - + ModelCheckerResult res = null; int i, n, iters; double soln[], soln2[], tmpsoln[]; @@ -901,7 +921,8 @@ public class DTMCModelChecker extends ProbModelChecker * @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. */ - protected ModelCheckerResult computeReachRewardsValIter(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet inf, double init[], BitSet known) throws PrismException + protected ModelCheckerResult computeReachRewardsValIter(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet inf, double init[], BitSet known) + throws PrismException { ModelCheckerResult res; BitSet unknown; @@ -971,7 +992,7 @@ public class DTMCModelChecker extends ProbModelChecker msg += "\nConsider using a different numerical method or increasing the maximum number of iterations"; throw new PrismException(msg); } - + // Return results res = new ModelCheckerResult(); res.soln = soln; @@ -990,6 +1011,98 @@ public class DTMCModelChecker extends ProbModelChecker * @param initDist Initial distribution (will be overwritten) */ public ModelCheckerResult computeSteadyStateProbs(DTMC dtmc, double initDist[]) throws PrismException + { + ModelCheckerResult res; + BitSet startNot, bscc; + double probBSCCs[], solnProbs[], probsReach[]; + int n, numBSCCs = 0, allInOneBSCC; + + // Store num states + n = dtmc.getNumStates(); + // Create results vector + solnProbs = new double[n]; + + // Compute bottom strongly connected components (BSCCs) + SCCComputer sccComputer = SCCComputer.createSCCComputer(sccMethod, dtmc); + sccComputer.computeBSCCs(); + List bsccs = sccComputer.getBSCCs(); + BitSet notInBSCCs = sccComputer.getNotInBSCCs(); + numBSCCs = bsccs.size(); + + // See which states in the initial distribution do *not* have non-zero prob + startNot = new BitSet(); + for (int i = 0; i < n; i++) { + if (initDist[i] == 0) + startNot.set(i); + } + // Determine whether initial states are all in a single BSCC + allInOneBSCC = -1; + for (int b = 0; b < numBSCCs; b++) { + if (!bsccs.get(b).intersects(startNot)) { + allInOneBSCC = b; + break; + } + } + + // If all initial states are in a single BSCC, it's easy... + // Just compute steady-state probabilities for the BSCC + if (allInOneBSCC != -1) { + mainLog.println("\nInitial states all in one BSCC (so no reachability probabilities computed)"); + bscc = bsccs.get(allInOneBSCC); + computeSteadyStateProbsForBSCC(dtmc, bscc, solnProbs); + } + + // Otherwise, have to consider all the BSCCs + else { + + // Compute probability of reaching each BSCC from initial distribution + probBSCCs = new double[numBSCCs]; + for (int b = 0; b < numBSCCs; b++) { + mainLog.println("\nComputing probability of reaching BSCC " + (b + 1)); + bscc = bsccs.get(b); + // Compute probabilities + probsReach = computeReachProbs(dtmc, bscc).soln; + // Compute probability of reaching BSCC, which is dot product of + // vectors for initial distribution and probabilities of reaching it + probBSCCs[b] = 0.0; + for (int i = 0; i < n; i++) { + probBSCCs[b] += initDist[i] * probsReach[i]; + } + mainLog.print("\nProbability of reaching BSCC " + (b + 1) + ": " + probBSCCs[b] + "\n"); + } + + // Compute steady-state probabilities for each BSCC + for (int b = 0; b < numBSCCs; b++) { + mainLog.println("\nComputing steady-state probabilities for BSCC " + (b + 1)); + bscc = bsccs.get(b); + // Compute steady-state probabilities for the BSCC + computeSteadyStateProbsForBSCC(dtmc, bscc, solnProbs); + // Multiply by BSCC reach prob + for (int i = bscc.nextSetBit(0); i >= 0; i = bscc.nextSetBit(i + 1)) + solnProbs[i] *= probBSCCs[b]; + } + } + + // Return results + res = new ModelCheckerResult(); + res.soln = solnProbs; + // TODO + //res.timeTaken = timer / 1000.0; + return res; + } + + /** + * Compute steady-state probabilities for a BSCC + * i.e. compute the long-run probability of being in each state of the BSCC. + * No initial distribution is specified since it does not affect the result. + * The result will be stored in the relevant portion of a full vector, + * whose size equals the number of states in the DTMC. + * Optionally, pass in an existing vector to be used for this purpose. + * @param dtmc The DTMC + * @param bscc The BSCC to be analysed + * @param result Storage for result (ignored if null) + */ + public ModelCheckerResult computeSteadyStateProbsForBSCC(DTMC dtmc, BitSet bscc, double result[]) throws PrismException { ModelCheckerResult res; int n, iters; @@ -1005,12 +1118,14 @@ public class DTMCModelChecker extends ProbModelChecker n = dtmc.getNumStates(); // Create solution vector(s) - // For soln, we just use init (since we are free to modify this vector) - soln = initDist; + // Use the passed in vector, if present + soln = result == null ? new double[n] : result; soln2 = new double[n]; - - // No need to initialise solution vectors - // (soln is done, soln2 will be immediately overwritten) + + // Initialise solution vectors. Equiprobable for BSCC states. + double equiprob = 1.0 / bscc.cardinality(); + for (int i = bscc.nextSetBit(0); i >= 0; i = bscc.nextSetBit(i + 1)) + soln[i] = soln2[i] = equiprob; // Start iterations iters = 0; @@ -1038,7 +1153,7 @@ public class DTMCModelChecker extends ProbModelChecker msg += "\nConsider using a different numerical method or increasing the maximum number of iterations"; throw new PrismException(msg); } - + // Return results res = new ModelCheckerResult(); res.soln = soln; @@ -1046,7 +1161,7 @@ public class DTMCModelChecker extends ProbModelChecker res.timeTaken = timer / 1000.0; return res; } - + /** * Compute transient probabilities * i.e. compute the probability of being in each state at time step {@code k}, @@ -1061,7 +1176,7 @@ public class DTMCModelChecker extends ProbModelChecker { throw new PrismException("Not implemented yet"); } - + /** * Simple test program. */ diff --git a/prism/src/explicit/PrismExplicit.java b/prism/src/explicit/PrismExplicit.java index ed7a2e20..4ed2acb0 100644 --- a/prism/src/explicit/PrismExplicit.java +++ b/prism/src/explicit/PrismExplicit.java @@ -268,7 +268,7 @@ public class PrismExplicit DTMCModelChecker mcDTMC = new DTMCModelChecker(); mcDTMC.setLog(mainLog); mcDTMC.setSettings(settings); - probs = mcDTMC.doSteadyState((DTMC) model, null); + probs = mcDTMC.doSteadyState((DTMC) model); } else if (model.getModelType() == ModelType.CTMC) { throw new PrismException("Not implemented yet"); // TODO diff --git a/prism/src/explicit/StateModelChecker.java b/prism/src/explicit/StateModelChecker.java index fb1052c6..62165015 100644 --- a/prism/src/explicit/StateModelChecker.java +++ b/prism/src/explicit/StateModelChecker.java @@ -36,6 +36,8 @@ import java.util.HashMap; import java.util.List; import java.util.Map; +import explicit.SCCComputer.SCCMethod; + import parser.State; import parser.Values; import parser.ast.Expression; @@ -78,6 +80,12 @@ public class StateModelChecker // PRISM settings object //protected PrismSettings settings = new PrismSettings(); + + // Flags/settings + // (NB: defaults do not necessarily coincide with PRISM) + + // Method used for numerical solution + protected SCCMethod sccMethod = SCCMethod.TARJAN; // Model file (for reward structures, etc.) protected ModulesFile modulesFile = null; diff --git a/prism/src/prism/Prism.java b/prism/src/prism/Prism.java index 34549959..f5f8c4a9 100644 --- a/prism/src/prism/Prism.java +++ b/prism/src/prism/Prism.java @@ -2628,7 +2628,7 @@ public class Prism implements PrismSettingsListener DTMCModelChecker mcDTMC = new DTMCModelChecker(); mcDTMC.setLog(mainLog); mcDTMC.setSettings(settings); - probsExpl = mcDTMC.doSteadyState((DTMC) currentModelExpl, (double[]) null); + probsExpl = mcDTMC.doSteadyState((DTMC) currentModelExpl, (File) null); //TODO: probsExpl = mcDTMC.doSteadyState((DTMC) currentModelExpl, fileIn); } else if (currentModelType == ModelType.CTMC) { throw new PrismException("Not implemented yet");