MDP-LTL model checking for explicit (still needs fixing) + correction to DTMC-DRA product construction.

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

prism/src/explicit/LTLModelChecker.java

@@ -141,18 +141,26 @@ public class LTLModelChecker extends PrismComponent
 		int numAPs = dra.getAPList().size();
 		int modelNumStates = dtmc.getNumStates();
 		int prodNumStates = modelNumStates * draSize;
+		int s_1, s_2, q_1, q_2;
+		BitSet s_labels = new BitSet(numAPs);
 
 		// Encoding: 
 		// each state s' = <s, q> = s * draSize + q
 		// s(s') = s' / draSize
 		// q(s') = s' % draSize
 
-		// Initial states
 		LinkedList<Point> queue = new LinkedList<Point>();
 		int map[] = new int[prodNumStates];
 		Arrays.fill(map, -1);
-		int q_0 = dra.getStartState();
+		
+		// Initial states
 		for (int s_0 : dtmc.getInitialStates()) {
+			// Get BitSet representing APs (labels) satisfied by initial state s_0
+			for (int k = 0; k < numAPs; k++) {
+				s_labels.set(k, labelBS.get(k).get(s_0));
+			}
+			// Find corresponding initial state in DRA
+			int q_0 = dra.getEdgeDestByLabel(dra.getStartState(), s_labels);
 			queue.add(new Point(s_0, q_0));
 			prodModel.addState();
 			prodModel.addInitialState(prodModel.getNumStates() - 1);
@@ -161,8 +169,6 @@ public class LTLModelChecker extends PrismComponent
 
 		// Product states
 		BitSet visited = new BitSet(prodNumStates);
-		int q_1 = 0, q_2 = 0, s_1 = 0, s_2 = 0;
-		BitSet s_2_labels = new BitSet(numAPs);
 		while (!queue.isEmpty()) {
 			Point p = queue.pop();
 			s_1 = p.x;
@@ -177,10 +183,10 @@ public class LTLModelChecker extends PrismComponent
 				double prob = e.getValue();
 				// Get BitSet representing APs (labels) satisfied by successor state s_2
 				for (int k = 0; k < numAPs; k++) {
-					s_2_labels.set(k, labelBS.get(k).get(s_2));
+					s_labels.set(k, labelBS.get(k).get(s_2));
 				}
 				// Find corresponding successor in DRA
-				q_2 = dra.getEdgeDestByLabel(q_1, s_2_labels);
+				q_2 = dra.getEdgeDestByLabel(q_1, s_labels);
 				// Add state/transition to model
 				if (!visited.get(s_2 * draSize + q_2) && map[s_2 * draSize + q_2] == -1) {
 					queue.add(new Point(s_2, q_2));
@@ -203,13 +209,100 @@ public class LTLModelChecker extends PrismComponent
 		return new Pair<Model, int[]>(prodModel, invMap);
 	}
 
+	/**
+	 * Construct the product of a DRA and an MDP.
+	 * @param dra The DRA
+	 * @param mdp The MDP
+	 * @param labelBS BitSets giving the set of states for each AP in the DRA
+	 * @return a Pair consisting of the product DTMC and a map from 
+	 *   (s_i * draSize + q_j) to the right state in the DRA product 
+	 */
+	public Pair<NondetModel, int[]> constructProductMDP(DRA<BitSet> dra, MDP mdp, Vector<BitSet> labelBS) throws PrismException
+	{
+		MDPSimple prodModel = new MDPSimple();
+
+		int draSize = dra.size();
+		int numAPs = dra.getAPList().size();
+		int modelNumStates = mdp.getNumStates();
+		int prodNumStates = modelNumStates * draSize;
+		int s_1, s_2, q_1, q_2;
+		BitSet s_labels = new BitSet(numAPs);
+
+		// Encoding: 
+		// each state s' = <s, q> = s * draSize + q
+		// s(s') = s' / draSize
+		// q(s') = s' % draSize
+
+		LinkedList<Point> queue = new LinkedList<Point>();
+		int map[] = new int[prodNumStates];
+		Arrays.fill(map, -1);
+		
+		// Initial states
+		for (int s_0 : mdp.getInitialStates()) {
+			// Get BitSet representing APs (labels) satisfied by initial state s_0
+			for (int k = 0; k < numAPs; k++) {
+				s_labels.set(k, labelBS.get(k).get(s_0));
+			}
+			// Find corresponding initial state in DRA
+			int q_0 = dra.getEdgeDestByLabel(dra.getStartState(), s_labels);
+			queue.add(new Point(s_0, q_0));
+			prodModel.addState();
+			prodModel.addInitialState(prodModel.getNumStates() - 1);
+			map[s_0 * draSize + q_0] = prodModel.getNumStates() - 1;
+		}
+
+		// Product states
+		BitSet visited = new BitSet(prodNumStates);
+		while (!queue.isEmpty()) {
+			Point p = queue.pop();
+			s_1 = p.x;
+			q_1 = p.y;
+			visited.set(s_1 * draSize + q_1);
+
+			// Go through transitions from state s_1 in original DTMC
+			int numChoices = mdp.getNumChoices(s_1);
+			for (int j = 0; j < numChoices; j++) {
+				Distribution prodDistr = new Distribution();
+				Iterator<Map.Entry<Integer, Double>> iter = mdp.getTransitionsIterator(s_1, j);
+				while (iter.hasNext()) {
+					Map.Entry<Integer, Double> e = iter.next();
+					s_2 = e.getKey();
+					double prob = e.getValue();
+					// Get BitSet representing APs (labels) satisfied by successor state s_2
+					for (int k = 0; k < numAPs; k++) {
+						s_labels.set(k, labelBS.get(k).get(s_2));
+					}
+					// Find corresponding successor in DRA
+					q_2 = dra.getEdgeDestByLabel(q_1, s_labels);
+					// Add state/transition to model
+					if (!visited.get(s_2 * draSize + q_2) && map[s_2 * draSize + q_2] == -1) {
+						queue.add(new Point(s_2, q_2));
+						prodModel.addState();
+						map[s_2 * draSize + q_2] = prodModel.getNumStates() - 1;
+					}
+					prodDistr.set(map[s_2 * draSize + q_2], prob);
+				}
+				prodModel.addActionLabelledChoice(map[s_1 * draSize + q_1], prodDistr, mdp.getAction(s_1, j));
+			}
+		}
+
+		int invMap[] = new int[prodModel.getNumStates()];
+		for (int i = 0; i < map.length; i++) {
+			if (map[i] != -1) {
+				invMap[map[i]] = i;
+			}
+		}
+
+		prodModel.findDeadlocks(false);
+
+		return new Pair<NondetModel, int[]>(prodModel, invMap);
+	}
+
 	/**
 	 * Find the set of states belong to accepting BSCCs in a model, according to a DRA {@code dra}.
 	 * @param dra The DRA
 	 * @param model The model
 	 * @param invMap The map returned by the constructProduct method(s)
-	 * @param sccMethod The method to use for SCC detection
-	 * @return
 	 */
 	public BitSet findAcceptingBSCCs(DRA<BitSet> dra, Model model, int invMap[]) throws PrismException
 	{
@@ -243,4 +336,44 @@ public class LTLModelChecker extends PrismComponent
 
 		return result;
 	}
+	
+	/**
+	 * Find the set of states belong to accepting MECs in a model, according to a DRA {@code dra}.
+	 * @param dra The DRA
+	 * @param model The model
+	 * @param invMap The map returned by the constructProduct method(s)
+	 */
+	public BitSet findAcceptingMECStates(DRA<BitSet> dra, NondetModel model, int invMap[]) throws PrismException
+	{
+		// Compute maximum end components (MECs)
+		ECComputer ecComputer = ECComputer.createECComputer(this, model);
+		ecComputer.computeMECStates();
+		List<BitSet> mecs = ecComputer.getMECStates();
+		mainLog.println(mecs);
+
+		int draSize = dra.size();
+		int numAcceptancePairs = dra.getNumAcceptancePairs();
+		BitSet result = new BitSet();
+
+		for (BitSet mec : mecs) {
+			boolean isLEmpty = true;
+			boolean isKEmpty = true;
+			for (int acceptancePair = 0; acceptancePair < numAcceptancePairs && isLEmpty && isKEmpty; acceptancePair++) {
+				BitSet L = dra.getAcceptanceL(acceptancePair);
+				BitSet K = dra.getAcceptanceK(acceptancePair);
+				for (int state = mec.nextSetBit(0); state != -1; state = mec.nextSetBit(state + 1)) {
+					int draState = invMap[state] % draSize;
+					isLEmpty &= !L.get(draState);
+					isKEmpty &= !K.get(draState);
+				}
+				// Stop as soon as we find the first acceptance pair that is satisfied
+				if (isLEmpty && !isKEmpty) {
+					result.or(mec);
+					break;
+				}
+			}
+		}
+
+		return result;
+	}
 }

prism/src/explicit/MDPModelChecker.java

@@ -28,16 +28,22 @@ package explicit;
 
 import java.util.BitSet;
 import java.util.Iterator;
+import java.util.LinkedList;
 import java.util.List;
 import java.util.Map;
+import java.util.Vector;
 
 import parser.ast.Expression;
 import parser.ast.ExpressionTemporal;
 import parser.ast.ExpressionUnaryOp;
+import parser.visitor.ASTTraverse;
+import prism.DRA;
+import prism.Pair;
 import prism.PrismComponent;
 import prism.PrismDevNullLog;
 import prism.PrismException;
 import prism.PrismFileLog;
+import prism.PrismLangException;
 import prism.PrismLog;
 import prism.PrismUtils;
 import strat.MDStrategyArray;
@@ -61,21 +67,21 @@ public class MDPModelChecker extends ProbModelChecker
 	/**
 	 * Compute probabilities for the contents of a P operator.
 	 */
-	protected StateValues checkProbPathFormula(Model model, Expression expr, boolean min) throws PrismException
+	protected StateValues checkProbPathFormula(NondetModel model, Expression expr, boolean min) throws PrismException
 	{
 		// Test whether this is a simple path formula (i.e. PCTL)
 		// and then pass control to appropriate method. 
 		if (expr.isSimplePathFormula()) {
 			return checkProbPathFormulaSimple(model, expr, min);
 		} else {
-			throw new PrismException("Explicit engine does not yet handle LTL-style path formulas");
+			return checkProbPathFormulaLTL(model, expr, min);
 		}
 	}
 
 	/**
 	 * Compute probabilities for a simple, non-LTL path operator.
 	 */
-	protected StateValues checkProbPathFormulaSimple(Model model, Expression expr, boolean min) throws PrismException
+	protected StateValues checkProbPathFormulaSimple(NondetModel model, Expression expr, boolean min) throws PrismException
 	{
 		StateValues probs = null;
 
@@ -125,7 +131,7 @@ public class MDPModelChecker extends ProbModelChecker
 	/**
 	 * Compute probabilities for a next operator.
 	 */
-	protected StateValues checkProbNext(Model model, ExpressionTemporal expr, boolean min) throws PrismException
+	protected StateValues checkProbNext(NondetModel model, ExpressionTemporal expr, boolean min) throws PrismException
 	{
 		BitSet target = null;
 		ModelCheckerResult res = null;
@@ -140,7 +146,7 @@ public class MDPModelChecker extends ProbModelChecker
 	/**
 	 * Compute probabilities for a bounded until operator.
 	 */
-	protected StateValues checkProbBoundedUntil(Model model, ExpressionTemporal expr, boolean min) throws PrismException
+	protected StateValues checkProbBoundedUntil(NondetModel model, ExpressionTemporal expr, boolean min) throws PrismException
 	{
 		int time;
 		BitSet b1, b2;
@@ -183,7 +189,7 @@ public class MDPModelChecker extends ProbModelChecker
 	/**
 	 * Compute probabilities for an (unbounded) until operator.
 	 */
-	protected StateValues checkProbUntil(Model model, ExpressionTemporal expr, boolean min) throws PrismException
+	protected StateValues checkProbUntil(NondetModel model, ExpressionTemporal expr, boolean min) throws PrismException
 	{
 		BitSet b1, b2;
 		StateValues probs = null;
@@ -206,10 +212,104 @@ public class MDPModelChecker extends ProbModelChecker
 		return probs;
 	}
 
+	/**
+	 * Compute probabilities for an LTL path formula
+	 */
+	protected StateValues checkProbPathFormulaLTL(NondetModel model, Expression expr, boolean min) throws PrismException
+	{
+		LTLModelChecker mcLtl;
+		StateValues probsProduct, probs;
+		Expression ltl;
+		DRA<BitSet> dra;
+		NondetModel modelProduct;
+		MDPModelChecker mcProduct;
+		long time;
+
+		// Can't do LTL with time-bounded variants of the temporal operators
+		try {
+			expr.accept(new ASTTraverse()
+			{
+				public void visitPre(ExpressionTemporal e) throws PrismLangException
+				{
+					if (e.getLowerBound() != null)
+						throw new PrismLangException(e.getOperatorSymbol());
+					if (e.getUpperBound() != null)
+						throw new PrismLangException(e.getOperatorSymbol());
+				}
+			});
+		} catch (PrismLangException e) {
+			String s = "Temporal operators (like " + e.getMessage() + ")";
+			s += " cannot have time bounds for LTL properties";
+			throw new PrismException(s);
+		}
+
+		// For LTL model checking routines
+		mcLtl = new LTLModelChecker(this);
+
+		// Model check maximal state formulas
+		Vector<BitSet> labelBS = new Vector<BitSet>();
+		ltl = mcLtl.checkMaximalStateFormulas(this, model, expr.deepCopy(), labelBS);
+
+		// Convert LTL formula to deterministic Rabin automaton (DRA)
+		// For min probabilities, need to negate the formula
+		// (add parentheses to allow re-parsing if required)
+		if (min) {
+			ltl = Expression.Not(Expression.Parenth(ltl));
+		}
+		mainLog.println("\nBuilding deterministic Rabin automaton (for " + ltl + ")...");
+		time = System.currentTimeMillis();
+		dra = LTLModelChecker.convertLTLFormulaToDRA(ltl);
+		int draSize = dra.size();
+		mainLog.println("\nDRA has " + dra.size() + " states, " + dra.getNumAcceptancePairs() + " pairs.");
+		// dra.print(System.out);
+		time = System.currentTimeMillis() - time;
+		mainLog.println("\nTime for Rabin translation: " + time / 1000.0 + " seconds.");
+
+		// Build product of MDP and automaton
+		mainLog.println("\nConstructing MDP-DRA product...");
+		Pair<NondetModel, int[]> pair = mcLtl.constructProductMDP(dra, (MDP) model, labelBS);
+		modelProduct = pair.first;
+		int invMap[] = pair.second;
+		int modelProductSize = modelProduct.getNumStates();
+
+		// Find accepting MECs + compute reachability probabilities
+		mainLog.println("\nFinding accepting MECs...");
+		BitSet acceptingMECs = mcLtl.findAcceptingMECStates(dra, modelProduct, invMap);
+		mainLog.println("\nComputing reachability probabilities...");
+		mcProduct = new MDPModelChecker(this);
+		probsProduct = StateValues.createFromDoubleArray(mcProduct.computeReachProbs((MDP) modelProduct, acceptingMECs, false).soln, modelProduct);
+
+		// Subtract from 1 if we're model checking a negated formula for regular Pmin
+		if (min) {
+			probsProduct.timesConstant(-1.0);
+			probsProduct.plusConstant(1.0);
+		}
+
+		// Mapping probabilities in the original model
+		double[] probsProductDbl = probsProduct.getDoubleArray();
+		double[] probsDbl = new double[model.getNumStates()];
+
+		LinkedList<Integer> queue = new LinkedList<Integer>();
+		for (int s : model.getInitialStates())
+			queue.add(s);
+
+		for (int i = 0; i < invMap.length; i++) {
+			int j = invMap[i];
+			int s = j / draSize;
+			// TODO: check whether this is the right way to compute probabilities in the original model
+			probsDbl[s] = Math.max(probsDbl[s], probsProductDbl[i]);
+		}
+
+		probs = StateValues.createFromDoubleArray(probsDbl, model);
+		probsProduct.clear();
+
+		return probs;
+	}
+
 	/**
 	 * Compute rewards for the contents of an R operator.
 	 */
-	protected StateValues checkRewardFormula(Model model, MDPRewards modelRewards, Expression expr, boolean min) throws PrismException
+	protected StateValues checkRewardFormula(NondetModel model, MDPRewards modelRewards, Expression expr, boolean min) throws PrismException
 	{
 		StateValues rewards = null;
 
@@ -233,7 +333,7 @@ public class MDPModelChecker extends ProbModelChecker
 	/**
 	 * Compute rewards for a reachability reward operator.
 	 */
-	protected StateValues checkRewardReach(Model model, MDPRewards modelRewards, ExpressionTemporal expr, boolean min) throws PrismException
+	protected StateValues checkRewardReach(NondetModel model, MDPRewards modelRewards, ExpressionTemporal expr, boolean min) throws PrismException
 	{
 		BitSet b;
 		StateValues rewards = null;

prism/src/explicit/ProbModelChecker.java

@@ -484,13 +484,13 @@ public class ProbModelChecker extends StateModelChecker
 			probs = ((CTMCModelChecker) this).checkProbPathFormula(model, expr.getExpression());
 			break;
 		case CTMDP:
-			probs = ((CTMDPModelChecker) this).checkProbPathFormula(model, expr.getExpression(), min);
+			probs = ((CTMDPModelChecker) this).checkProbPathFormula((NondetModel) model, expr.getExpression(), min);
 			break;
 		case DTMC:
 			probs = ((DTMCModelChecker) this).checkProbPathFormula(model, expr.getExpression());
 			break;
 		case MDP:
-			probs = ((MDPModelChecker) this).checkProbPathFormula(model, expr.getExpression(), min);
+			probs = ((MDPModelChecker) this).checkProbPathFormula((NondetModel) model, expr.getExpression(), min);
 			break;
 		/*case STPG:
 			probs = ((STPGModelChecker) this).checkProbPathFormula(model, expr.getExpression(), min);
@@ -598,7 +598,7 @@ public class ProbModelChecker extends StateModelChecker
 			rews = ((DTMCModelChecker) this).checkRewardFormula(model, mcRewards, expr.getExpression());
 			break;
 		case MDP:
-			rews = ((MDPModelChecker) this).checkRewardFormula(model, mdpRewards, expr.getExpression(), min);
+			rews = ((MDPModelChecker) this).checkRewardFormula((NondetModel) model, mdpRewards, expr.getExpression(), min);
 			break;
 		default:
 			throw new PrismException("Cannot model check " + expr + " for " + modelType + "s");

prism/src/prism/NondetModelChecker.java

@@ -1363,11 +1363,9 @@ public class NondetModelChecker extends NonProbModelChecker
 			modelProduct.exportStates(Prism.EXPORT_PLAIN, new PrismFileLog(prism.getExportProductStatesFilename()));
 		}
 
-		// Find accepting maximum end components
+		// Find accepting MECs + compute reachability probabilities
 		mainLog.println("\nFinding accepting end components...");
 		JDDNode acc = mcLtl.findAcceptingStates(dra, modelProduct, draDDRowVars, draDDColVars, fairness);
-
-		// Compute reachability probabilities
 		mainLog.println("\nComputing reachability probabilities...");
 		mcProduct = new NondetModelChecker(prism, modelProduct, null);
 		probsProduct = mcProduct.checkProbUntil(modelProduct.getReach(), acc, qual, min && fairness);