For MDP model checking (explicit, symbolic), support -pmaxquotient option (computation in MEC quotient)

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

prism/src/explicit/MDPModelChecker.java

@@ -54,6 +54,8 @@ import acceptance.AcceptanceType;
 import automata.DA;
 import automata.LTL2WDBA;
 import common.IterableBitSet;
+import explicit.modelviews.EquivalenceRelationInteger;
+import explicit.modelviews.MDPEquiv;
 import explicit.rewards.MCRewards;
 import explicit.rewards.MCRewardsFromMDPRewards;
 import explicit.rewards.MDPRewards;
@@ -339,6 +341,8 @@ public class MDPModelChecker extends ProbModelChecker
 		// Local copy of setting
 		MDPSolnMethod mdpSolnMethod = this.mdpSolnMethod;
 
+		boolean doPmaxQuotient = this.doPmaxQuotient;
+
 		// Switch to a supported method, if necessary
 		if (mdpSolnMethod == MDPSolnMethod.LINEAR_PROGRAMMING) {
 			mdpSolnMethod = MDPSolnMethod.GAUSS_SEIDEL;
@@ -358,6 +362,11 @@ public class MDPModelChecker extends ProbModelChecker
 			}
 		}
 
+		if (doPmaxQuotient && min) {
+			// for Pmin, don't do quotient
+			doPmaxQuotient = false;
+		}
+
 		// Start probabilistic reachability
 		timer = System.currentTimeMillis();
 		mainLog.println("\nStarting probabilistic reachability (" + (min ? "min" : "max") + ")...");
@@ -439,7 +448,48 @@ public class MDPModelChecker extends ProbModelChecker
 
 		// Compute probabilities (if needed)
 		if (numYes + numNo < n) {
-			res = computeReachProbsNumeric(mdp, mdpSolnMethod, no, yes, min, init, known, strat);
+
+			if (!min && doPmaxQuotient) {
+				MDPEquiv maxQuotient = maxQuotient(mdp, yes, no);
+				// MDPEquiv retains original state space, making the states that are not used
+				// trap states.
+				// yesInQuotient is the representative for the yes equivalence class
+				BitSet yesInQuotient = new BitSet();
+				yesInQuotient.set(maxQuotient.mapStateToRestrictedModel(yes.nextSetBit(0)));
+				// noInQuotient is the representative for the no equivalence class as well
+				// as the non-representative states (the states in any equivalence class
+				// that are not the representative for the class). As the latter states
+				// are traps, we can just add them to the no set
+				BitSet noInQuotient = new BitSet();
+				noInQuotient.set(maxQuotient.mapStateToRestrictedModel(no.nextSetBit(0)));
+				noInQuotient.or(maxQuotient.getNonRepresentativeStates());
+				MDPSparse quotientModel = new MDPSparse(maxQuotient);
+
+				ModelCheckerResult res1 = computeReachProbsNumeric(quotientModel,
+				                                                   mdpSolnMethod,
+				                                                   noInQuotient,
+				                                                   yesInQuotient,
+				                                                   min,
+				                                                   init,
+				                                                   known,
+				                                                   strat);
+
+				res = new ModelCheckerResult();
+				res.numIters = res1.numIters;
+				res.timeTaken = res1.timeTaken;
+				res.soln = new double[mdp.getNumStates()];
+				for (int i = 0; i < n; i++) {
+					if (yes.get(i)) {
+						res.soln[i] = 1.0;
+					} else if (no.get(i)) {
+						res.soln[i] = 0.0;
+					} else {
+						res.soln[i] = res1.soln[maxQuotient.mapStateToRestrictedModel(i)];
+					}
+				}
+			} else {
+				res = computeReachProbsNumeric(mdp, mdpSolnMethod, no, yes, min, init, known, strat);
+			}
 		} else {
 			res = new ModelCheckerResult();
 			res.soln = Utils.bitsetToDoubleArray(yes, n);
@@ -1971,6 +2021,38 @@ public class MDPModelChecker extends ProbModelChecker
 		}
 	}
 
+	/**
+	 * Compute the end component quotient (for use with PMax),
+	 * each maximal end component is collapsed to a single state,
+	 * likewise the yes and no regions, respectively.
+	 */
+	private MDPEquiv maxQuotient(MDP mdp, BitSet yes, BitSet no) throws PrismException
+	{
+		BitSet maybe = new BitSet();
+		maybe.set(0, mdp.getNumStates());
+		maybe.andNot(yes);
+		maybe.andNot(no);
+
+		ECComputer ec = ECComputer.createECComputer(this, mdp);
+
+		ec.computeMECStates(maybe);
+		List<BitSet> mecs = ec.getMECStates();
+		mecs.add(yes);
+		mecs.add(no);
+
+		EquivalenceRelationInteger eq = new EquivalenceRelationInteger(mecs);
+		BasicModelTransformation<MDP, MDPEquiv> quotientTransform = MDPEquiv.transformDroppingLoops(mdp, eq);
+		MDPEquiv quotient = quotientTransform.getTransformedModel();
+
+		//mdp.exportToDotFile("original.dot");
+		//quotient.exportToDotFile("maxQuotient.dot");
+
+		int realStates = quotient.getNumStates() - quotient.getNonRepresentativeStates().cardinality();
+		mainLog.println("Max-Quotient MDP: " + realStates + " equivalence classes / non-trap states.");
+
+		return quotient;
+	}
+
 	/**
 	 * Simple test program.
 	 */

prism/src/explicit/StateModelChecker.java

@@ -116,6 +116,10 @@ public class StateModelChecker extends PrismComponent
 	// Do bisimulation minimisation before model checking?
 	protected boolean doBisim = false;
 
+
+	// For Pmax computation, collapse MECs to quotient MDP?
+	protected boolean doPmaxQuotient = false;
+
 	// Model info (for reward structures, etc.)
 	protected ModulesFile modulesFile = null;
 	protected ModelInfo modelInfo = null;
@@ -211,6 +215,7 @@ public class StateModelChecker extends PrismComponent
 		setStoreVector(other.getStoreVector());
 		setGenStrat(other.getGenStrat());
 		setDoBisim(other.getDoBisim());
+		setDoPmaxQuotient(other.getDoPmaxQuotient());
 	}
 
 	/**
@@ -295,6 +300,14 @@ public class StateModelChecker extends PrismComponent
 		this.doBisim = doBisim;
 	}
 
+	/**
+	 * Specify whether or not to perform MEC quotienting for Pmax.
+	 */
+	public void setDoPmaxQuotient(boolean doPmaxQuotient)
+	{
+		this.doPmaxQuotient = doPmaxQuotient;
+	}
+
 	// Get methods for flags/settings
 
 	public int getVerbosity()
@@ -366,6 +379,14 @@ public class StateModelChecker extends PrismComponent
 		return doBisim;
 	}
 
+	/**
+	 * Whether or not to do MEC quotient for Pmax
+	 */
+	public boolean getDoPmaxQuotient()
+	{
+		return doPmaxQuotient;
+	}
+
 	/** Get the constant values (both from the modules file and the properties file) */
 	public Values getConstantValues()
 	{

prism/src/prism/NondetModelChecker.java

@@ -72,6 +72,7 @@ import acceptance.AcceptanceType;
 import automata.DA;
 import automata.LTL2DA;
 import automata.LTL2WDBA;
+import common.StopWatch;
 import dv.DoubleVector;
 import dv.IntegerVector;
 import explicit.MinMax;
@@ -1823,6 +1824,19 @@ public class NondetModelChecker extends NonProbModelChecker
 		DoubleVector probsDV;
 		StateValues probs = null;
 
+		boolean doPmaxQuotient = getSettings().getBoolean(PrismSettings.PRISM_PMAX_QUOTIENT);
+
+		if (doPmaxQuotient && min) {
+			// don't do pmaxQuotient for min
+			doPmaxQuotient = false;
+		}
+
+		if (doPmaxQuotient) {
+			if (!(precomp && prob0 && prob1)) {
+				throw new PrismNotSupportedException("Precomputations (Prob0 & Prob1) must be enabled for -pmaxquotient setting");
+			}
+		}
+
 		// If required, export info about target states 
 		if (prism.getExportTarget()) {
 			JDDNode labels[] = { model.getStart(), b2 };
@@ -1915,9 +1929,67 @@ public class NondetModelChecker extends NonProbModelChecker
 			mainLog.println("\nComputing remaining probabilities...");
 			mainLog.println("Engine: " + Prism.getEngineString(engine));
 			try {
+				MDPQuotient transform = null;
+				NondetModel transformed = null;
+				JDDNode yesInQuotient = null;
+				JDDNode maybeInQuotient = null;
+
+				if (doPmaxQuotient) {
+					if (!tr.equals(model.getTrans()) ||
+					    !tra.equals(model.getTransActions()) ||
+					    !tr01.equals(model.getTrans01())) {
+						throw new PrismException("Can currently not compute MEC quotient for changed functions");
+					}
+
+					mainLog.println("\nBuilding quotient MDP, collapsing maximal end components as well as yes and no states...");
+					StopWatch ecWatch = new StopWatch(mainLog);
+					ecWatch.start("computing maximal end components");
+					ECComputer ec = ECComputer.createECComputer(this, model);
+					// find MECs in the maybe states
+					ec.computeMECStates(maybe);
+					ecWatch.stop("found " + ec.getMECStates().size() + " MECs");
+
+					List<JDDNode> ecs = new ArrayList<JDDNode>(ec.getMECStates());
+
+					ecs.add(yes.copy());
+					ecs.add(no.copy());
+
+					StopWatch watchTransform = new StopWatch(mainLog);
+					watchTransform.start("building MEC quotient");
+					transform = MDPQuotient.transform(this, model, ecs, model.getReach().copy());
+					watchTransform.stop();
+
+					if (false) {
+						try {
+							model.exportToFile(Prism.EXPORT_DOT, true, new File("model.dot"));
+							transform.getTransformedModel().exportToFile(Prism.EXPORT_DOT, true, new File("quotient.dot"));
+						} catch (FileNotFoundException e) {
+						}
+					}
+					transformed = transform.getTransformedModel();
+
+					mainLog.println("\nQuotient MDP:");
+					transformed.printTransInfo(mainLog);
+
+					yesInQuotient = transform.mapStateSetToQuotient(yes.copy());
+					maybeInQuotient = transform.mapStateSetToQuotient(maybe.copy());
+				}
+
 				switch (engine) {
 				case Prism.MTBDD:
-					probsMTBDD = PrismMTBDD.NondetUntil(tr, odd, nondetMask, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min);
+					if (transform != null) {
+						probsMTBDD = PrismMTBDD.NondetUntil(transformed.getTrans(),
+						                                    transformed.getODD(),
+						                                    transformed.getNondetMask(),
+						                                    transformed.getAllDDRowVars(),
+						                                    transformed.getAllDDColVars(),
+						                                    transformed.getAllDDNondetVars(),
+						                                    yesInQuotient,
+						                                    maybeInQuotient,
+						                                    min);
+					} else {
+						probsMTBDD = PrismMTBDD.NondetUntil(tr, odd, nondetMask, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min);
+					}
 					probs = new StateValuesMTBDD(probsMTBDD, model);
 					break;
 				case Prism.SPARSE:
@@ -1927,19 +1999,55 @@ public class NondetModelChecker extends NonProbModelChecker
 						strat = new IntegerVector(ddStrat, allDDRowVars, odd);
 						JDD.Deref(ddStrat);
 					}
-					probsDV = PrismSparse.NondetUntil(tr, tra, model.getSynchs(), odd, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min, strat);
-					if (genStrat) {
+					if (transform != null) {
+						strat = null;  // strategy generation with the quotient not yet supported
+						probsDV = PrismSparse.NondetUntil(transformed.getTrans(),
+						                                  transformed.getTransActions(),
+						                                  transformed.getSynchs(),
+						                                  transformed.getODD(),
+						                                  transformed.getAllDDRowVars(),
+						                                  transformed.getAllDDColVars(),
+						                                  transformed.getAllDDNondetVars(),
+						                                  yesInQuotient,
+						                                  maybeInQuotient,
+						                                  min,
+						                                  strat);
+						probs = new StateValuesDV(probsDV, transformed);
+					} else {
+						probsDV = PrismSparse.NondetUntil(tr, tra, model.getSynchs(), odd, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min, strat);
+						probs = new StateValuesDV(probsDV, model);
+					}
+					if (genStrat && strat != null) {
 						result.setStrategy(new MDStrategyIV(model, strat));
 					}
-					probs = new StateValuesDV(probsDV, model);
 					break;
 				case Prism.HYBRID:
-					probsDV = PrismHybrid.NondetUntil(tr, odd, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min);
-					probs = new StateValuesDV(probsDV, model);
+					if (transform != null) {
+						probsDV = PrismHybrid.NondetUntil(transformed.getTrans(),
+						                                  transformed.getODD(),
+						                                  transformed.getAllDDRowVars(),
+						                                  transformed.getAllDDColVars(),
+						                                  transformed.getAllDDNondetVars(),
+						                                  yesInQuotient,
+						                                  maybeInQuotient,
+						                                  min);
+						probs = new StateValuesDV(probsDV, transformed);
+					} else {
+						probsDV = PrismHybrid.NondetUntil(tr, odd, allDDRowVars, allDDColVars, allDDNondetVars, yes, maybe, min);
+						probs = new StateValuesDV(probsDV, model);
+					}
 					break;
 				default:
 					throw new PrismException("Unknown engine");
 				}
+
+				if (transform != null) {
+					// we have to project back to the original
+					probs = transform.projectToOriginalModel(probs);
+					transform.clear();
+					JDD.Deref(yesInQuotient, maybeInQuotient);
+				}
+
 			} catch (PrismException e) {
 				JDD.Deref(yes);
 				JDD.Deref(no);

prism/src/prism/Prism.java

@@ -3730,6 +3730,7 @@ public class Prism extends PrismComponent implements PrismSettingsListener
 		mc.setStoreVector(storeVector);
 		mc.setGenStrat(genStrat);
 		mc.setDoBisim(doBisim);
+		mc.setDoPmaxQuotient(settings.getBoolean(PrismSettings.PRISM_PMAX_QUOTIENT));
 
 		return mc;
 	}

prism/src/prism/PrismSettings.java

@@ -84,6 +84,7 @@ public class PrismSettings implements Observer
 	public static final	String PRISM_COMPACT						= "prism.compact";
 	public static final	String PRISM_LIN_EQ_METHOD					= "prism.linEqMethod";//"prism.iterativeMethod";
 	public static final	String PRISM_LIN_EQ_METHOD_PARAM			= "prism.linEqMethodParam";//"prism.overRelaxation";
+	public static final	String PRISM_PMAX_QUOTIENT					= "prism.pmaxQuotient";
 	public static final	String PRISM_MDP_SOLN_METHOD				= "prism.mdpSolnMethod";
 	public static final	String PRISM_MDP_MULTI_SOLN_METHOD			= "prism.mdpMultiSolnMethod";
 	public static final	String PRISM_TERM_CRIT						= "prism.termCrit";//"prism.termination";
@@ -241,6 +242,8 @@ public class PrismSettings implements Observer
 																			"Which iterative method to use when solving linear equation systems." },
 			{ DOUBLE_TYPE,		PRISM_LIN_EQ_METHOD_PARAM,				"Over-relaxation parameter",			"2.1",			new Double(0.9),															"",																							
 																			"Over-relaxation parameter for iterative numerical methods such as JOR/SOR." },
+			{ BOOLEAN_TYPE,		PRISM_PMAX_QUOTIENT,				"For Pmax computations, compute in the MEC quotient",				"4.3.1",		false,																		"",
+																				"For Pmax computations, compute in the MEC quotient."},
 			{ CHOICE_TYPE,		PRISM_MDP_SOLN_METHOD,					"MDP solution method",				"4.0",			"Value iteration",																"Value iteration,Gauss-Seidel,Policy iteration,Modified policy iteration,Linear programming",
 																			"Which method to use when solving Markov decision processes." },
 			{ CHOICE_TYPE,		PRISM_MDP_MULTI_SOLN_METHOD,			"MDP multi-objective solution method",				"4.0.3",			"Value iteration",											"Value iteration,Gauss-Seidel,Linear programming",
@@ -986,6 +989,12 @@ public class PrismSettings implements Observer
 			set(PRISM_MDP_SOLN_METHOD, "Linear programming");
 			set(PRISM_MDP_MULTI_SOLN_METHOD, "Linear programming");
 		}
+
+		// Pmax quotient
+		else if (sw.equals("pmaxquotient")) {
+			set(PRISM_PMAX_QUOTIENT, true);
+		}
+
 		// Linear equation solver over-relaxation parameter
 		else if (sw.equals("omega")) {
 			if (i < args.length - 1) {
@@ -1678,6 +1687,7 @@ public class PrismSettings implements Observer
 		mainLog.println("-ltl2datool <exec> ............. Run executable <exec> to convert LTL formulas to deterministic automata");
 		mainLog.println("-ltl2dasyntax <x> .............. Specify output format for -ltl2datool switch (lbt, spin, spot, rabinizer)");
 		mainLog.println("-exportiterations .............. Export vectors for iteration algorithms to file");
+		mainLog.println("-pmaxquotient .................. For Pmax computations in MDPs, compute in the MEC quotient.");
 		
 		mainLog.println();
 		mainLog.println("MULTI-OBJECTIVE MODEL CHECKING:");