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

@@ -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.
 	 */

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[]);
 }

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<Integer, Double> e : distr) {
+					j = (Integer) e.getKey();
+					prob = (Double) e.getValue();
+					result[j] += (prob / er) * vect[i];
+				}
+			}
+		}
+	}
+
 	@Override
 	public String toString()
 	{

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()
 	{

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<Integer, Double>();
 			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<Entry<Integer,Double>> getTransitionsIterator(int s)
+	public Iterator<Entry<Integer, Double>> 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<Integer, Double> e : distr) {
+				j = (Integer) e.getKey();
+				prob = (Double) e.getValue();
+				result[j] += prob * vect[i];
+			}
+
+		}
+	}
+
 	// Accessors (other)
 
 	/**

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<Integer, Double> 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()
 	{

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.
 	 */

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
 	 */

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.");
 		}