From 4ec5f0f9aeab3b4c627f04d4824e2842a13a3985 Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Fri, 17 Jun 2011 15:54:58 +0000 Subject: [PATCH] Transient probability computation in explicit engine + some connection to CL. git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@3110 bbc10eb1-c90d-0410-af57-cb519fbb1720 --- prism/src/explicit/CTMCModelChecker.java | 142 +++++++++++++++++- prism/src/explicit/DTMC.java | 30 ++-- prism/src/explicit/DTMCEmbeddedSimple.java | 30 ++++ .../DTMCFromMDPMemorylessAdversary.java | 6 + prism/src/explicit/DTMCSimple.java | 29 +++- prism/src/explicit/DTMCUniformisedSimple.java | 29 ++++ prism/src/explicit/PrismExplicit.java | 74 +++++++++ prism/src/explicit/StateValues.java | 31 ++++ prism/src/prism/PrismCL.java | 34 ++++- 9 files changed, 382 insertions(+), 23 deletions(-) diff --git a/prism/src/explicit/CTMCModelChecker.java b/prism/src/explicit/CTMCModelChecker.java index 938c13ae..082712d4 100644 --- a/prism/src/explicit/CTMCModelChecker.java +++ b/prism/src/explicit/CTMCModelChecker.java @@ -28,9 +28,10 @@ package explicit; import java.util.*; +import jdd.JDD; +import jdd.JDDNode; import explicit.rewards.MCRewards; import explicit.rewards.MCRewardsStateArray; - import parser.ast.ExpressionTemporal; import prism.*; @@ -77,6 +78,43 @@ public class CTMCModelChecker extends DTMCModelChecker return probs; } + // Transient/steady-state probability computation + + /** + * Compute transient 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. + */ + public StateValues doTransient(CTMC ctmc, double time, double initDist[]) throws PrismException + { + ModelCheckerResult res = null; + int n; + double initDistNew[] = null; + StateValues probs = null; + + // Store num states + n = ctmc.getNumStates(); + + // Build initial distribution (if not specified) + if (initDist == null) { + initDistNew = new double[n]; + for (int in : ctmc.getInitialStates()) { + initDistNew[in] = 1 / ctmc.getNumInitialStates(); + } + } else { + initDistNew = initDist; + throw new PrismException("Not implemented yet"); // TODO + } + + // Compute transient probabilities + res = computeTransientProbs(ctmc, initDistNew, time); + probs = StateValues.createFromDoubleArray(res.soln); + + return probs; + } + // Numerical computation functions /** @@ -119,8 +157,7 @@ public class CTMCModelChecker extends DTMCModelChecker * @param t Time bound * @param init Initial solution vector - pass null for default */ - public ModelCheckerResult computeTimeBoundedUntilProbs(CTMC ctmc, BitSet target, double t, double init[]) - throws PrismException + public ModelCheckerResult computeTimeBoundedUntilProbs(CTMC ctmc, BitSet target, double t, double init[]) throws PrismException { ModelCheckerResult res = null; int i, n, iters; @@ -159,7 +196,7 @@ public class CTMCModelChecker extends DTMCModelChecker // Build (implicit) uniformised DTMC dtmc = ctmc.buildImplicitUniformisedDTMC(q); - + // Create solution vector(s) soln = new double[n]; soln2 = (init == null) ? new double[n] : init; @@ -235,6 +272,103 @@ public class CTMCModelChecker extends DTMCModelChecker return super.computeReachRewards(dtmcEmb, rewEmb, target); } + /** + * Compute transient probabilities. + * i.e. compute the probability of being in each state at time {@code t}, + * assuming the initial distribution . + * @param ctmc The CTMC + * @param initDist Initial distribution + * @param t Time point + */ + public ModelCheckerResult computeTransientProbs(CTMC ctmc, double initDist[], double t) throws PrismException + { + ModelCheckerResult res = null; + int i, n, iters; + double soln[], soln2[], tmpsoln[], sum[]; + DTMC dtmc; + long timer; + // Fox-Glynn stuff + FoxGlynn fg; + int left, right; + double q, qt, acc, weights[], totalWeight; + + // Start bounded probabilistic reachability + timer = System.currentTimeMillis(); + mainLog.println("Starting transient probability computation..."); + + // Store num states + n = ctmc.getNumStates(); + + // Get uniformisation rate; do Fox-Glynn + q = ctmc.getDefaultUniformisationRate(); + qt = q * t; + mainLog.println("\nUniformisation: q.t = " + q + " x " + t + " = " + qt); + termCritParam = 1e-6; + acc = termCritParam / 8.0; + fg = new FoxGlynn(qt, 1e-300, 1e+300, acc); + left = fg.getLeftTruncationPoint(); + right = fg.getRightTruncationPoint(); + if (right < 0) { + throw new PrismException("Overflow in Fox-Glynn computation (time bound too big?)"); + } + weights = fg.getWeights(); + totalWeight = fg.getTotalWeight(); + for (i = left; i <= right; i++) { + weights[i - left] /= totalWeight; + } + mainLog.println("Fox-Glynn (" + acc + "): left = " + left + ", right = " + right); + + // Build (implicit) uniformised DTMC + dtmc = ctmc.buildImplicitUniformisedDTMC(q); + + // Create solution vector(s) + // For soln, we just use init (since we are free to modify this vector) + soln = initDist; + soln2 = new double[n]; + sum = new double[n]; + + // Initialise solution vectors + // (don't need to do soln2 since will be immediately overwritten) + for (i = 0; i < n; i++) + sum[i] = 0.0; + + // If necessary, do 0th element of summation (doesn't require any matrix powers) + if (left == 0) + for (i = 0; i < n; i++) + sum[i] += weights[0] * soln[i]; + + // Start iterations + iters = 1; + while (iters <= right) { + // Matrix-vector multiply + dtmc.vmMult(soln, soln2); + // Swap vectors for next iter + tmpsoln = soln; + soln = soln2; + soln2 = tmpsoln; + // Add to sum + if (iters >= left) { + for (i = 0; i < n; i++) + sum[i] += weights[iters - left] * soln[i]; + } + iters++; + } + + // Finished bounded probabilistic reachability + timer = System.currentTimeMillis() - timer; + mainLog.print("Transient probability computation"); + mainLog.println(" took " + iters + " iters and " + timer / 1000.0 + " seconds."); + + // Return results + res = new ModelCheckerResult(); + res.soln = sum; + res.lastSoln = soln2; + res.numIters = iters; + res.timeTaken = timer / 1000.0; + res.timePre = 0.0; + return res; + } + /** * Simple test program. */ diff --git a/prism/src/explicit/DTMC.java b/prism/src/explicit/DTMC.java index 12a78f1c..feb94b84 100644 --- a/prism/src/explicit/DTMC.java +++ b/prism/src/explicit/DTMC.java @@ -51,17 +51,6 @@ public interface DTMC extends Model */ public double getTransitionReward(int s); - /** - * Do a matrix-vector multiplication for - * the DTMC's transition probability matrix P and the vector {@code vect} passed in. - * i.e. for all s: result[s] = sum_j P(s,j)*vect[j] - * @param vect Vector to multiply by - * @param result Vector to store result in - * @param subset Only do multiplication for these rows (ignored if null) - * @param complement If true, {@code subset} is taken to be its complement (ignored if {@code subset} is null) - */ - public void mvMult(double vect[], double result[], BitSet subset, boolean complement); - /** * Perform a single step of precomputation algorithm Prob0, i.e., for states i in {@code subset}, * set bit i of {@code result} iff there is a transition to a state in {@code u}. @@ -81,6 +70,17 @@ public interface DTMC extends Model */ public void prob1step(BitSet subset, BitSet u, BitSet v, BitSet result); + /** + * Do a matrix-vector multiplication for + * the DTMC's transition probability matrix P and the vector {@code vect} passed in. + * i.e. for all s: result[s] = sum_j P(s,j)*vect[j] + * @param vect Vector to multiply by + * @param result Vector to store result in + * @param subset Only do multiplication for these rows (ignored if null) + * @param complement If true, {@code subset} is taken to be its complement (ignored if {@code subset} is null) + */ + public void mvMult(double vect[], double result[], BitSet subset, boolean complement); + /** * Do a single row of matrix-vector multiplication for * the DTMC's transition probability matrix P and the vector {@code vect} passed in. @@ -132,4 +132,12 @@ public interface DTMC extends Model */ public double mvMultRewSingle(int s, double vect[], MCRewards mcRewards); + /** + * Do a vector-matrix multiplication for + * the DTMC's transition probability matrix P and the vector {@code vect} passed in. + * i.e. for all s: result[s] = sum_i P(i,s)*vect[i] + * @param vect Vector to multiply by + * @param result Vector to store result in + */ + public void vmMult(double vect[], double result[]); } diff --git a/prism/src/explicit/DTMCEmbeddedSimple.java b/prism/src/explicit/DTMCEmbeddedSimple.java index 4e88ae30..d2e58c16 100644 --- a/prism/src/explicit/DTMCEmbeddedSimple.java +++ b/prism/src/explicit/DTMCEmbeddedSimple.java @@ -382,6 +382,36 @@ public class DTMCEmbeddedSimple implements DTMC return d; } + @Override + public void vmMult(double vect[], double result[]) + { + int i, j; + double prob, er; + Distribution distr; + + // Initialise result to 0 + for (j = 0; j < numStates; j++) { + result[j] = 0; + } + // Go through matrix elements (by row) + for (i = 0; i < numStates; i++) { + distr = ctmc.getTransitions(i); + er = exitRates[i]; + // Exit rate 0: prob 1 self-loop + if (er == 0) { + result[i] += vect[i]; + } + // Exit rate > 0 + else { + for (Map.Entry e : distr) { + j = (Integer) e.getKey(); + prob = (Double) e.getValue(); + result[j] += (prob / er) * vect[i]; + } + } + } + } + @Override public String toString() { diff --git a/prism/src/explicit/DTMCFromMDPMemorylessAdversary.java b/prism/src/explicit/DTMCFromMDPMemorylessAdversary.java index f4ba8bea..78c63037 100644 --- a/prism/src/explicit/DTMCFromMDPMemorylessAdversary.java +++ b/prism/src/explicit/DTMCFromMDPMemorylessAdversary.java @@ -290,6 +290,12 @@ public class DTMCFromMDPMemorylessAdversary implements DTMC //return mdp.mvMultRewSingle(s, adv[s], vect); } + @Override + public void vmMult(double vect[], double result[]) + { + throw new RuntimeException("Not implemented yet"); // TODO + } + @Override public String toString() { diff --git a/prism/src/explicit/DTMCSimple.java b/prism/src/explicit/DTMCSimple.java index 28ee980c..9aa29760 100644 --- a/prism/src/explicit/DTMCSimple.java +++ b/prism/src/explicit/DTMCSimple.java @@ -312,7 +312,7 @@ public class DTMCSimple extends ModelSimple implements DTMC } if (fix) { for (i = deadlocks.nextSetBit(0); i >= 0; i = deadlocks.nextSetBit(i + 1)) { - setProbability(i, i, 1.0); + setProbability(i, i, 1.0); } } return deadlocks; @@ -387,7 +387,7 @@ public class DTMCSimple extends ModelSimple implements DTMC // Output transitions to PRISM language file out = new FileWriter(filename); out.write(getModelType().keyword() + "\n"); - out.write("module M\nx : [0.." + (numStates-1) + "];\n"); + out.write("module M\nx : [0.." + (numStates - 1) + "];\n"); sorted = new TreeMap(); for (i = 0; i < numStates; i++) { // Extract transitions and sort by destination state index (to match PRISM-exported files) @@ -442,7 +442,7 @@ public class DTMCSimple extends ModelSimple implements DTMC } @Override - public Iterator> getTransitionsIterator(int s) + public Iterator> getTransitionsIterator(int s) { return trans.get(s).iterator(); } @@ -609,6 +609,29 @@ public class DTMCSimple extends ModelSimple implements DTMC return d; } + @Override + public void vmMult(double vect[], double result[]) + { + int i, j; + double prob; + Distribution distr; + + // Initialise result to 0 + for (j = 0; j < numStates; j++) { + result[j] = 0; + } + // Go through matrix elements (by row) + for (i = 0; i < numStates; i++) { + distr = trans.get(i); + for (Map.Entry e : distr) { + j = (Integer) e.getKey(); + prob = (Double) e.getValue(); + result[j] += prob * vect[i]; + } + + } + } + // Accessors (other) /** diff --git a/prism/src/explicit/DTMCUniformisedSimple.java b/prism/src/explicit/DTMCUniformisedSimple.java index 99b6f7f0..569f75f7 100644 --- a/prism/src/explicit/DTMCUniformisedSimple.java +++ b/prism/src/explicit/DTMCUniformisedSimple.java @@ -351,6 +351,35 @@ public class DTMCUniformisedSimple implements DTMC throw new Error("Not yet supported"); } + @Override + public void vmMult(double vect[], double result[]) + { + int i, j; + double prob, sum; + Distribution distr; + + // Initialise result to 0 + for (j = 0; j < numStates; j++) { + result[j] = 0; + } + // Go through matrix elements (by row) + for (i = 0; i < numStates; i++) { + distr = ctmc.getTransitions(i); + sum = 0.0; + for (Map.Entry e : distr) { + j = (Integer) e.getKey(); + prob = (Double) e.getValue(); + // Non-diagonal entries only + if (j != i) { + sum += prob; + result[j] += (prob / q) * vect[i]; + } + } + // Diagonal entry is 1 - sum/q + result[i] += (1 - sum / q) * vect[i]; + } + } + @Override public String toString() { diff --git a/prism/src/explicit/PrismExplicit.java b/prism/src/explicit/PrismExplicit.java index 4953959b..ba06420f 100644 --- a/prism/src/explicit/PrismExplicit.java +++ b/prism/src/explicit/PrismExplicit.java @@ -168,6 +168,80 @@ public class PrismExplicit return result; } + /** + * Compute transient probabilities (for a DTMC or CTMC). + * Output probability distribution to log. + */ + public void doTransient(Model model, double time) throws PrismException + { + doTransient(model, time, Prism.EXPORT_PLAIN, null, null); + } + + /** + * Compute transient probabilities (for a DTMC or CTMC). + * Output probability distribution to a file (or, if file is null, to log). + * The exportType should be EXPORT_PLAIN or EXPORT_MATLAB. + * Optionally (if non-null), read in the initial probability distribution from a file. + */ + public void doTransient(Model model, double time, int exportType, File fileOut, File fileIn) throws PrismException + { + long l = 0; // timer + StateValues probs = null; + PrismLog tmpLog; + + if (time < 0) throw new PrismException("Cannot compute transient probabilities for negative time value"); + + // no specific states format for MRMC + if (exportType == Prism.EXPORT_MRMC) exportType = Prism.EXPORT_PLAIN; + // rows format does not apply to states output + if (exportType == Prism.EXPORT_ROWS) exportType = Prism.EXPORT_PLAIN; + + l = System.currentTimeMillis(); + + if (model.getModelType() == ModelType.DTMC) { + throw new PrismException("Not implemented yet"); // TODO + } + else if (model.getModelType() == ModelType.CTMC) { + mainLog.println("\nComputing transient probabilities (time = " + time + ")..."); + CTMCModelChecker mcCTMC = new CTMCModelChecker(); + // TODO: pass settings + probs = mcCTMC.doTransient((CTMC) model, time, null); + } + else { + throw new PrismException("Transient probabilities only computed for DTMCs/CTMCs"); + } + + l = System.currentTimeMillis() - l; + + // print message + mainLog.print("\nPrinting transient probabilities "); + switch (exportType) { + case Prism.EXPORT_PLAIN: mainLog.print("in plain text format "); break; + case Prism.EXPORT_MATLAB: mainLog.print("in Matlab format "); break; + } + if (fileOut != null) mainLog.println("to file \"" + fileOut + "\"..."); else mainLog.println("below:"); + + // create new file log or use main log + if (fileOut != null) { + tmpLog = new PrismFileLog(fileOut.getPath()); + if (!tmpLog.ready()) { + throw new PrismException("Could not open file \"" + fileOut + "\" for output"); + } + } else { + tmpLog = mainLog; + } + + // print out or export probabilities + probs.print(tmpLog, fileOut == null, exportType == Prism.EXPORT_MATLAB, fileOut == null); + + // print out computation time + mainLog.println("\nTime for transient probability computation: " + l/1000.0 + " seconds."); + + // tidy up + probs.clear(); + if (fileOut != null) tmpLog.close(); + } + /** * Simple test program. */ diff --git a/prism/src/explicit/StateValues.java b/prism/src/explicit/StateValues.java index de5065fe..1c7d077f 100644 --- a/prism/src/explicit/StateValues.java +++ b/prism/src/explicit/StateValues.java @@ -32,6 +32,7 @@ import dv.DoubleVector; import parser.type.*; import prism.PrismException; +import prism.PrismLog; /** * Class for explicit-state storage of a state-indexed vector of (integer or double) values. @@ -158,6 +159,18 @@ public class StateValues // ... + // clear (free memory) + + /** + * Clear the vector, i.e. free any used memory. + * (Well, actually, just set pointer to null and wait for later garbage collection.) + */ + public void clear() + { + valuesI = null; + valuesD = null; + } + // METHODS TO ACCESS VECTOR DATA /** @@ -188,6 +201,24 @@ public class StateValues } */ + // PRINTING STUFF + + /** + * Print vector to a log/file (non-zero entries only) + */ + public void print(PrismLog log) throws PrismException + { + int i; + for (i = 0; i < size; i++) { + log.println(getValue(i)); + } + } + + public void print(PrismLog log, boolean printSparse, boolean printMatlab, boolean printStates) throws PrismException + { + print(log); //TODO + } + /** * Make a (deep) copy of this vector */ diff --git a/prism/src/prism/PrismCL.java b/prism/src/prism/PrismCL.java index 3036d11e..f6c0b2b0 100644 --- a/prism/src/prism/PrismCL.java +++ b/prism/src/prism/PrismCL.java @@ -912,6 +912,7 @@ public class PrismCL private void doSteadyState() throws PrismException { + ModelType modelType; File exportSteadyStateFile = null; // choose destination for output (file or log) @@ -920,8 +921,15 @@ public class PrismCL else exportSteadyStateFile = new File(exportSteadyStateFilename); + // Determine model type + if (explicit) { + modelType = modelExpl.getModelType(); + } else { + modelType = model.getModelType(); + } + // compute steady-state probabilities - if (model.getModelType() == ModelType.CTMC || model.getModelType() == ModelType.DTMC) { + if (modelType == ModelType.CTMC || modelType == ModelType.DTMC) { prism.doSteadyState(model, exportType, exportSteadyStateFile); } else { mainLog.println("\nWarning: Steady-state probabilities only computed for DTMCs/CTMCs."); @@ -932,6 +940,7 @@ public class PrismCL private void doTransient() throws PrismException { + ModelType modelType; double d; int i; File exportTransientFile = null; @@ -942,21 +951,36 @@ public class PrismCL else exportTransientFile = new File(exportTransientFilename); + // Determine model type + if (explicit) { + modelType = modelExpl.getModelType(); + } else { + modelType = model.getModelType(); + } + // compute transient probabilities - if (model.getModelType() == ModelType.CTMC) { + if (modelType == ModelType.CTMC) { try { d = Double.parseDouble(transientTime); } catch (NumberFormatException e) { throw new PrismException("Invalid value \"" + transientTime + "\" for transient probability computation"); } - prism.doTransient(model, d, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); - } else if (model.getModelType() == ModelType.DTMC) { + if (explicit) { + prismExpl.doTransient(modelExpl, d, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); + } else { + prism.doTransient(model, d, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); + } + } else if (modelType == ModelType.DTMC) { try { i = Integer.parseInt(transientTime); } catch (NumberFormatException e) { throw new PrismException("Invalid value \"" + transientTime + "\" for transient probability computation"); } - prism.doTransient(model, i, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); + if (explicit) { + prismExpl.doTransient(modelExpl, i, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); + } else { + prism.doTransient(model, i, exportType, exportTransientFile, importinitdist ? new File(importInitDistFilename) : null); + } } else { mainLog.println("\nWarning: Transient probabilities only computed for DTMCs/CTMCs."); }