explicit: add support for steady-state probability computation for CTMCs

First version to get basic support. Uses the implicit embedded DTMC,
might be worth it to convert to a DTMCSparse for performance.

prism/src/explicit/CTMCModelChecker.java

@@ -29,6 +29,8 @@ package explicit;
 import java.io.File;
 import java.util.*;
 
+import common.IterableBitSet;
+import explicit.StateValues;
 import explicit.rewards.MCRewards;
 import explicit.rewards.Rewards;
 import explicit.rewards.StateRewardsArray;
@@ -218,6 +220,42 @@ public class CTMCModelChecker extends ProbModelChecker
 
 	// Steady-state/transient probability computation
 
+	/**
+	 * Compute steady-state probability distribution (forwards).
+	 * Start from initial state (or uniform distribution over multiple initial states).
+	 */
+	public StateValues doSteadyState(CTMC ctmc) throws PrismException
+	{
+		return doSteadyState(ctmc, (StateValues) null);
+	}
+
+	/**
+	 * Compute steady-state probability distribution (forwards).
+	 * Optionally, use the passed in file initDistFile to give the initial probability distribution (time 0).
+	 * If null, start from initial state (or uniform distribution over multiple initial states).
+	 */
+	public StateValues doSteadyState(CTMC ctmc, File initDistFile) throws PrismException
+	{
+		StateValues initDist = readDistributionFromFile(initDistFile, ctmc);
+		return doSteadyState(ctmc, initDist);
+	}
+
+	/**
+	 * Compute steady-state 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.
+	 * @param ctmc The CTMC
+	 * @param initDist Initial distribution (will be overwritten)
+	 */
+	public StateValues doSteadyState(CTMC ctmc, StateValues initDist) throws PrismException
+	{
+		StateValues initDistNew = (initDist == null) ? buildInitialDistribution(ctmc) : initDist;
+		ModelCheckerResult res = computeSteadyStateProbs(ctmc, initDistNew.getDoubleArray());
+		return StateValues.createFromDoubleArray(res.soln, ctmc);
+	}
+
 	/**
 	 * Compute transient probability distribution (forwards).
 	 * Start from initial state (or uniform distribution over multiple initial states).
@@ -739,6 +777,70 @@ public class CTMCModelChecker extends ProbModelChecker
 		return createDTMCModelChecker().computeReachRewards(dtmcEmb, rewEmb, target);
 	}
 
+	/**
+	 * We compute steady-state probabilities in the embedded DTMC.
+	 * To take the exit rates into account, we have to weight the
+	 * steady-state probabilities in each BSCC, using this post-processor.
+	 * See: Baier et al, "Approximate Symbolic Model Checking of Continuous-Time Markov Chains"
+	 * CONCUR'99, p. 151
+	 */
+	private static class SteadyStateBSCCPostProcessor implements DTMCModelChecker.BSCCPostProcessor {
+		private CTMC ctmc;
+
+		public SteadyStateBSCCPostProcessor(CTMC ctmc)
+		{
+			this.ctmc = ctmc;
+		}
+
+		@Override
+		public void apply(double[] soln, BitSet bscc)
+		{
+			// compute sum_{s in BSCC} pi'[s] / E[S]
+			// where pi' are the steady-state probabilities in the BSCC of the embedded DTMC
+			double sum = 0.0;
+			for (int s : IterableBitSet.getSetBits(bscc)) {
+				double E = ctmc.getExitRate(s);
+				if (E == 0.0) // corner case: no outgoing transitions -> self-loop with rate 1
+					E = 1.0;
+
+				sum += soln[s] / E;
+			}
+
+			// set pi[s] = pi'[s] / sum for each state s in BSCC, where again
+			// pi' are the steady-state probabilities in the BSCC of the embedded DTMC
+			// and pi are the steady-state probabilities in the BSCC of the original CTMC
+			for (int s : IterableBitSet.getSetBits(bscc)) {
+				double E = ctmc.getExitRate(s);
+				if (E == 0.0) // corner case: no outgoing transitions -> self-loop with rate 1
+					E = 1.0;
+
+				soln[s] /= E;
+				soln[s] /= sum;
+			}
+		}
+
+	}
+
+	/**
+	 * Compute (forwards) steady-state probabilities
+	 * i.e. compute the long-run probability of being in each state,
+	 * assuming the initial distribution {@code initDist}.
+	 * For space efficiency, the initial distribution vector will be modified and values over-written,
+	 * so if you wanted it, take a copy.
+	 * @param ctmc The CTMC
+	 * @param initDist Initial distribution (will be overwritten)
+	 */
+	public ModelCheckerResult computeSteadyStateProbs(CTMC ctmc, double initDist[]) throws PrismException
+	{
+		// We construct the embedded DTMC and do the steady-state computation there
+		mainLog.println("Building embedded DTMC...");
+		DTMC dtmcEmb = ctmc.getImplicitEmbeddedDTMC();
+
+		// compute the steady-state probabilities in the embedded DTMC, applying the BSCC value post-processing
+		mainLog.println("Doing steady-state computation in embedded DTMC (with exit-rate weighting for BSCC probabilities)...");
+		return createDTMCModelChecker().computeSteadyStateProbs(dtmcEmb, initDist, new SteadyStateBSCCPostProcessor(ctmc));
+	}
+
 	/**
 	 * Compute transient probabilities.
 	 * i.e. compute the probability of being in each state at time {@code t},

prism/src/prism/Prism.java

@@ -3438,14 +3438,19 @@ public class Prism extends PrismComponent implements PrismSettingsListener
 	 */
 	protected explicit.StateValues computeSteadyStateProbabilitiesExplicit(explicit.Model model, File fileIn) throws PrismException
 	{
-		DTMCModelChecker mcDTMC;
 		explicit.StateValues probs;
-		if (model.getModelType() == ModelType.DTMC) {
-			mcDTMC = new DTMCModelChecker(this);
+		switch (model.getModelType()) {
+		case DTMC: {
+			DTMCModelChecker mcDTMC = new DTMCModelChecker(this);
 			probs = mcDTMC.doSteadyState((DTMC) model, fileIn);
-		} else if (model.getModelType() == ModelType.CTMC) {
-			throw new PrismException("Not implemented yet");
-		} else {
+			break;
+		}
+		case CTMC: {
+			CTMCModelChecker mcCTMC = new CTMCModelChecker(this);
+			probs = mcCTMC.doSteadyState((CTMC) model, fileIn);
+			break;
+		}
+		default:
 			throw new PrismException("Steady-state probabilities only computed for DTMCs/CTMCs");
 		}
 		return probs;