Fixes in explicit CTMC solving + some CTMDP stuff.

git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@1824 bbc10eb1-c90d-0410-af57-cb519fbb1720

prism/src/explicit/CTMC.java

@@ -37,7 +37,7 @@ public interface CTMC extends DTMC
 	public static ModelType modelType = ModelType.CTMC;
 
 	/**
-	 * Compute the maximum exit rate (ignoring self-loops).
+	 * Compute the maximum exit rate.
 	 */
 	public double getMaxExitRate();
 	

prism/src/explicit/CTMCModelChecker.java

@@ -91,11 +91,12 @@ public class CTMCModelChecker extends DTMCModelChecker
 		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, weights[], totalWeight;
+		double q, qt, acc, weights[], totalWeight;
 
 		// Start bounded probabilistic reachability
 		timer = System.currentTimeMillis();
@@ -105,10 +106,11 @@ public class CTMCModelChecker extends DTMCModelChecker
 		n = ctmc.getNumStates();
 
 		// Get uniformisation rate; do Fox-Glynn
-		q = ctmc.getMaxExitRate();
+		q = ctmc.getDefaultUniformisationRate();
 		qt = q * t;
 		mainLog.println("\nUniformisation: q.t = " + q + " x " + t + " = " + qt);
-		fg = new FoxGlynn(qt, 1e-300, 1e+300, termCritParam / 8.0);
+		acc = termCritParam / 8.0;
+		fg = new FoxGlynn(qt, 1e-300, 1e+300, acc);
 		left = fg.getLeftTruncationPoint();
 		right = fg.getRightTruncationPoint();
 		if (right < 0) {
@@ -119,8 +121,11 @@ public class CTMCModelChecker extends DTMCModelChecker
 		for (i = left; i <= right; i++) {
 			weights[i - left] /= totalWeight;
 		}
-		mainLog.println("Fox-Glynn: left = " + left + ", right = " + right);
+		mainLog.println("Fox-Glynn (" + acc + "): left = " + left + ", right = " + right);
 
+		// Build (implicit) uniformised DTMC
+		dtmc = ctmc.buildImplicitUniformisedDTMC(q);
+		
 		// Create solution vector(s)
 		soln = new double[n];
 		soln2 = (init == null) ? new double[n] : init;
@@ -146,7 +151,7 @@ public class CTMCModelChecker extends DTMCModelChecker
 		iters = 1;
 		while (iters < right) {
 			// Matrix-vector multiply
-			ctmc.mvMult(soln, soln2, target, true);
+			dtmc.mvMult(soln, soln2, target, true);
 			// Swap vectors for next iter
 			tmpsoln = soln;
 			soln = soln2;
@@ -224,7 +229,7 @@ public class CTMCModelChecker extends DTMCModelChecker
 			target = labels.get(args[2]);
 			if (target == null)
 				throw new PrismException("Unknown label \"" + args[2] + "\"");
-			for (int i = 3; i < args.length; i++) {
+			for (int i = 4; i < args.length; i++) {
 				if (args[i].equals("-nopre"))
 					mc.setPrecomp(false);
 			}

prism/src/explicit/CTMCSimple.java

@@ -67,7 +67,7 @@ public class CTMCSimple extends DTMCSimple implements CTMC
 		int i;
 		double d, max = Double.NEGATIVE_INFINITY;
 		for (i = 0; i < numStates; i++) {
-			d = trans.get(i).sumAllBut(i);
+			d = trans.get(i).sum();
 			if (d > max)
 				max = d;
 		}
@@ -125,8 +125,7 @@ public class CTMCSimple extends DTMCSimple implements CTMC
 	@Override
 	public DTMC buildImplicitUniformisedDTMC(double q)
 	{
-		// TODO
-		return null;
+		return new DTMCUniformisedSimple(this, q);
 	}
 	
 	@Override

prism/src/explicit/CTMDP.java

@@ -38,6 +38,11 @@ public interface CTMDP extends MDP
 
 	// TODO: copy/modify functions from CTMC
 	
+	/**
+	 * Compute the maximum exit rate.
+	 */
+	public double getMaxExitRate();
+	
 	/**
 	 * Build the discretised (DT)MDP for this CTMDP, in implicit form
 	 * (i.e. where the details are computed on the fly from this one).

prism/src/explicit/CTMDPModelChecker.java

@@ -51,11 +51,17 @@ public class CTMDPModelChecker extends MDPModelChecker
 		MDPModelChecker mc;
 		ModelCheckerResult res;
 		
-		mdp = ctmdp.buildDiscretisedMDP(0.0002);
+		// TODO: check locally uniform
+		double epsilon = 1e-3;
+		double q = ctmdp.getMaxExitRate();
+		int k = (int)Math.ceil((q * t * q * t) / (2 * epsilon));
+		double tau = t / k;
+		mainLog.println("q = " + q + ", k = " + k + ", tau = " + tau);
+		mdp = ctmdp.buildDiscretisedMDP(tau);
 		mainLog.println(mdp);
 		mc = new MDPModelChecker();
 		mc.inheritSettings(this);
-		res = mc.probReachBounded(mdp, target, 4500, min);
+		res = mc.probReachBounded(mdp, target, k, min);
 		
 		return res;
 	}
@@ -185,7 +191,7 @@ public class CTMDPModelChecker extends MDPModelChecker
 			target = labels.get(args[2]);
 			if (target == null)
 				throw new PrismException("Unknown label \"" + args[2] + "\"");
-			for (int i = 3; i < args.length; i++) {
+			for (int i = 4; i < args.length; i++) {
 				if (args[i].equals("-min"))
 					min = true;
 				else if (args[i].equals("-max"))
@@ -193,7 +199,7 @@ public class CTMDPModelChecker extends MDPModelChecker
 				else if (args[i].equals("-nopre"))
 					mc.setPrecomp(false);
 			}
-			res = mc.probReachBounded2(ctmdp, target, 1.0, min, null, null);
+			res = mc.probReachBounded2(ctmdp, target, Double.parseDouble(args[3]), min, null, null);
 			System.out.println(res.soln[0]);
 		} catch (PrismException e) {
 			System.out.println(e);

prism/src/explicit/CTMDPSimple.java

@@ -61,6 +61,21 @@ public class CTMDPSimple extends MDPSimple implements CTMDP
 
 	// Accessors (for CTMDP)
 
+	@Override
+	public double getMaxExitRate()
+	{
+		int i;
+		double d, max = Double.NEGATIVE_INFINITY;
+		for (i = 0; i < numStates; i++) {
+			for (Distribution distr : trans.get(i)) {
+				d = distr.sum();
+				if (d > max)
+					max = d;
+			}
+		}
+		return max;
+	}
+
 	@Override
 	public MDP buildImplicitDiscretisedMDP(double tau)
 	{
@@ -74,7 +89,7 @@ public class CTMDPSimple extends MDPSimple implements CTMDP
 		MDPSimple mdp;
 		Distribution distrNew;
 		int i;
-		double d;
+		double sum, d;
 		mdp = new MDPSimple(numStates);
 		for (int in : getInitialStates()) {
 			mdp.addInitialState(in);
@@ -82,9 +97,10 @@ public class CTMDPSimple extends MDPSimple implements CTMDP
 		for (i = 0; i < numStates; i++) {
 			for (Distribution distr : trans.get(i)) {
 				distrNew = new Distribution();
-				d = Math.exp(-distr.sum() * tau);
+				sum = distr.sum();
+				d = Math.exp(-sum * tau);
 				for (Map.Entry<Integer, Double> e : distr) {
-					distrNew.add(e.getKey(), (1 - d) * e.getValue());
+					distrNew.add(e.getKey(), (1 - d) * (e.getValue() / sum));
 				}
 				distrNew.add(i, d);
 				mdp.addChoice(i, distrNew);

prism/src/explicit/DTMCEmbeddedSimple.java

@@ -40,7 +40,7 @@ public class DTMCEmbeddedSimple implements DTMC
 	protected CTMCSimple ctmc;
 	// Exit rates vector
 	protected double exitRates[];
-	// Number of extra transitions added
+	// Number of extra transitions added (just for stats)
 	protected int numExtraTransitions;
 
 	/**

prism/src/explicit/DTMCUniformisedSimple.java

@@ -0,0 +1,252 @@
+//==============================================================================
+//
+//Copyright (c) 2002-
+//Authors:
+//* Dave Parker <david.parker@comlab.ox.ac.uk> (University of Oxford)
+//
+//------------------------------------------------------------------------------
+//
+//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.*;
+
+import prism.PrismException;
+
+/**
+* Simple explicit-state representation of a DTMC, constructed (implicitly) as the uniformised DTMC of a CTMC.
+* This class is read-only: most of data is pointers to other model info.
+*/
+public class DTMCUniformisedSimple implements DTMC
+{
+	// Parent CTMC
+	protected CTMCSimple ctmc;
+	// Uniformisation rate
+	protected double q;
+	// Number of extra transitions added (just for stats)
+	protected int numExtraTransitions;
+
+	/**
+	 * Constructor: create from CTMC and uniformisation rate q.
+	 */
+	public DTMCUniformisedSimple(CTMCSimple ctmc, double q)
+	{
+		this.ctmc = ctmc;
+		this.q = q;
+		int numStates = ctmc.getNumStates();
+		numExtraTransitions = 0;
+		for (int i = 0; i < numStates; i++) {
+			if (ctmc.getTransitions(i).get(i) == 0 && ctmc.getTransitions(i).sumAllBut(i) < q) {
+				numExtraTransitions++;
+			}
+		}
+	}
+
+	/**
+	 * Constructor: create from CTMC and its default uniformisation rate.
+	 */
+	public DTMCUniformisedSimple(CTMCSimple ctmc)
+	{
+		this(ctmc, ctmc.getDefaultUniformisationRate());
+	}
+
+	// Accessors (for Model)
+
+	public int getNumStates()
+	{
+		return ctmc.getNumStates();
+	}
+
+	public int getNumInitialStates()
+	{
+		return ctmc.getNumInitialStates();
+	}
+
+	public Iterable<Integer> getInitialStates()
+	{
+		return ctmc.getInitialStates();
+	}
+
+	public int getFirstInitialState()
+	{
+		return ctmc.getFirstInitialState();
+	}
+
+	public boolean isInitialState(int i)
+	{
+		return ctmc.isInitialState(i);
+	}
+
+	public int getNumTransitions()
+	{
+		return ctmc.getNumTransitions() + numExtraTransitions;
+	}
+
+	public boolean isSuccessor(int s1, int s2)
+	{
+		// TODO
+		throw new Error("Not yet supported");
+	}
+
+	public boolean allSuccessorsInSet(int s, BitSet set)
+	{
+		// TODO
+		throw new Error("Not yet supported");
+	}
+
+	public boolean someSuccessorsInSet(int s, BitSet set)
+	{
+		// TODO
+		throw new Error("Not yet supported");
+	}
+
+	public int getNumChoices(int s)
+	{
+		// Always 1 for a DTMC
+		return 1;
+	}
+
+	public void checkForDeadlocks() throws PrismException
+	{
+		// No deadlocks by definition
+	}
+
+	public void checkForDeadlocks(BitSet except) throws PrismException
+	{
+		// No deadlocks by definition
+	}
+
+	public void exportToPrismExplicit(String baseFilename) throws PrismException
+	{
+		throw new PrismException("Export not yet supported");
+	}
+
+	public void exportToDotFile(String filename) throws PrismException
+	{
+		throw new PrismException("Export not yet supported");
+	}
+
+	public void exportToDotFile(String filename, BitSet mark) throws PrismException
+	{
+		throw new PrismException("Export not yet supported");
+	}
+
+	public String infoString()
+	{
+		return ctmc.infoString() + " + " + numExtraTransitions + " self-loops";
+	}
+
+	// Accessors (for DTMC)
+
+	public double getTransitionReward(int s)
+	{
+		// TODO
+		throw new Error("Not yet supported");
+	}
+
+	public void mvMult(double vect[], double result[], BitSet subset, boolean complement)
+	{
+		int s, numStates;
+		numStates = ctmc.getNumStates();
+		// Loop depends on subset/complement arguments
+		if (subset == null) {
+			for (s = 0; s < numStates; s++)
+				result[s] = mvMultSingle(s, vect);
+		} else if (complement) {
+			for (s = subset.nextClearBit(0); s < numStates; s = subset.nextClearBit(s + 1))
+				result[s] = mvMultSingle(s, vect);
+		} else {
+			for (s = subset.nextSetBit(0); s >= 0; s = subset.nextSetBit(s + 1))
+				result[s] = mvMultSingle(s, vect);
+		}
+	}
+
+	public double mvMultSingle(int s, double vect[])
+	{
+		int k;
+		double sum, d, prob;
+		Distribution distr;
+
+		distr = ctmc.getTransitions(s);
+		sum = d = 0.0;
+		for (Map.Entry<Integer, Double> e : distr) {
+			k = (Integer) e.getKey();
+			prob = (Double) e.getValue();
+			// Non-diagonal entries
+			if (k != s) {
+				sum += prob;
+				d += (prob / q) * vect[k];
+			}
+		}
+		// Diagonal entry
+		if (sum < q) {
+			d += (1 - sum/q) * vect[s];
+		}
+
+		return d;
+	}
+
+	public void mvMultRew(double vect[], double result[], BitSet subset, boolean complement)
+	{
+		int s, numStates;
+		numStates = ctmc.getNumStates();
+		// Loop depends on subset/complement arguments
+		if (subset == null) {
+			for (s = 0; s < numStates; s++)
+				result[s] = mvMultRewSingle(s, vect);
+		} else if (complement) {
+			for (s = subset.nextClearBit(0); s < numStates; s = subset.nextClearBit(s + 1))
+				result[s] = mvMultRewSingle(s, vect);
+		} else {
+			for (s = subset.nextSetBit(0); s >= 0; s = subset.nextSetBit(s + 1))
+				result[s] = mvMultRewSingle(s, vect);
+		}
+	}
+
+	public double mvMultRewSingle(int s, double vect[])
+	{
+		// TODO
+		throw new Error("Not yet supported");
+	}
+
+	@Override
+	public String toString()
+	{
+		String s = "";
+		s += "ctmc: " + ctmc;
+		s = ", q: " + q;
+		return s;
+	}
+
+	@Override
+	public boolean equals(Object o)
+	{
+		if (o == null || !(o instanceof DTMCUniformisedSimple))
+			return false;
+		DTMCUniformisedSimple dtmc = (DTMCUniformisedSimple) o;
+		if (!ctmc.equals(dtmc.ctmc))
+			return false;
+		if (q != dtmc.q)
+			return false;
+		if (numExtraTransitions != dtmc.numExtraTransitions)
+			return false;
+		return true;
+	}
+}