diff --git a/prism/src/explicit/DTMC.java b/prism/src/explicit/DTMC.java
index 951aa255..65cc776a 100644
--- a/prism/src/explicit/DTMC.java
+++ b/prism/src/explicit/DTMC.java
@@ -31,6 +31,7 @@ import java.util.Map.Entry;
import java.util.PrimitiveIterator.OfInt;
import common.IterableStateSet;
+import common.iterable.IterableInt;
import prism.Pair;
import prism.PrismException;
import explicit.rewards.MCRewards;
@@ -607,4 +608,59 @@ public interface DTMC extends Model
});
}
}
+
+ /**
+ * Do a vector-matrix multiplication for steady-state computation with the power method.
+ *
+ * Computes
+ * {@code result = vect * P} with matrix
+ * {@code P = (Q * deltaT + I)} where
+ * {@code Q} is the generator matrix,
+ * {@code deltaT} the preconditioning factor and
+ * {@code I} is the the identity matrix.
+ * Please refer to William J. Stewart: "Introduction to the Numerical Solution of Markov Chains" p.124 for details.
+ *
+ *
+ * If the {@code states} argument only specifies a subset of the state space,
+ * only those entries of the {@code result} vector are modified that are either
+ * states in {@code states} or their successors; other entries remain unchanged.
+ * Thus, it generally only makes sense to use this method with a state sets that consists
+ * of (the union of) bottom strongly-connected components (BSCCs).
+ *
+ * @param vect Vector to multiply by
+ * @param result Vector to store result in
+ * @param diagsQ vector of the diagonal entries of the generator matrix Q, i.e., diagsQ[s] = -sum_{s!=t} prob(s,t)
+ * @param deltaT deltaT conditioning factor
+ * @param states subset of states to consider
+ */
+ public default void vmMultPowerSteadyState(double vect[], double[] result, double[] diagsQ, double deltaT, IterableInt states)
+ {
+ // Recall that the generator matrix Q has entries
+ // Q(s,s) = -sum_{t!=s} prob(s,t)
+ // and Q(s,t) = prob(s,t) for s!=t
+ // The values Q(s,s) are passed in via the diagsQ vector, while the
+ // values Q(s,t) correspond to the normal transitions
+
+ // Initialise result for relevant states to vect[s] * (deltaT * diagsQ[s] + 1),
+ // i.e., handle the product with the diagonal entries of (deltaT * Q) + I
+ for (OfInt it = states.iterator(); it.hasNext(); ) {
+ int state = it.nextInt();
+ result[state] = vect[state] * ((deltaT * diagsQ[state]) + 1.0);
+ }
+
+ // For each relevant state...
+ for (OfInt it = states.iterator(); it.hasNext();) {
+ int state = it.nextInt();
+
+ // ... handle all Q(state,t) entries of the generator matrix
+ forEachTransition(state, (s, t, prob) -> {
+ if (s != t) {
+ // ignore self loop, diagonal entries of the generator matrix handled above
+ // update result vector entry for the *successor* state
+ result[t] += deltaT * prob * vect[s];
+ }
+ });
+ }
+ }
+
}
diff --git a/prism/src/explicit/DTMCModelChecker.java b/prism/src/explicit/DTMCModelChecker.java
index 6d872e0c..f96f0af3 100644
--- a/prism/src/explicit/DTMCModelChecker.java
+++ b/prism/src/explicit/DTMCModelChecker.java
@@ -34,12 +34,14 @@ import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
import java.util.PrimitiveIterator;
+import java.util.PrimitiveIterator.OfInt;
import java.util.Vector;
import parser.VarList;
import parser.ast.Declaration;
import parser.ast.DeclarationIntUnbounded;
import parser.ast.Expression;
+import prism.ModelType;
import prism.OptionsIntervalIteration;
import prism.Prism;
import prism.PrismComponent;
@@ -2346,13 +2348,27 @@ public class DTMCModelChecker extends ProbModelChecker
* No initial distribution is specified since it does not affect the result.
* The result will be stored in the relevant portion of a full vector,
* whose size equals the number of states in the DTMC.
- * Optionally, pass in an existing vector to be used for this purpose.
+ * Optionally, pass in an existing vector to be used for this purpose;
+ * only the entries of this vector are changed that correspond to the BSCC states.
+ *
+ * To ensure convergence, we use the iteration matrix
+ * {@code P = (Q * deltaT + I)} where
+ * {@code Q} is the generator matrix,
+ * {@code deltaT} a preconditioning factor and
+ * {@code I} is the the identity matrix.
+ * See William J. Stewart: "Introduction to the Numerical Solution of Markov Chains" p.124 for details.
+ *
* @param dtmc The DTMC
- * @param bscc The BSCC to be analysed
+ * @param states The BSCC to be analysed
* @param result Storage for result (ignored if null)
*/
- public ModelCheckerResult computeSteadyStateProbsForBSCC(DTMC dtmc, BitSet bscc, double result[]) throws PrismException
+ public ModelCheckerResult computeSteadyStateProbsForBSCC(DTMC dtmc, BitSet states, double result[]) throws PrismException
{
+ if (dtmc.getModelType() != ModelType.DTMC) {
+ throw new PrismNotSupportedException("Explicit engine currently does not support steady-state computation for " + dtmc.getModelType());
+ }
+ IterableBitSet bscc = new IterableBitSet(states);
+
// Start value iteration
mainLog.println("Starting value iteration...");
StopWatch watch = new StopWatch(mainLog).start();
@@ -2363,12 +2379,49 @@ public class DTMCModelChecker extends ProbModelChecker
// Create solution vector(s)
// Use the passed in vector, if present
double[] soln = result == null ? new double[numStates] : result;
- double[] soln2 = new double[numStates];
+ double[] diagsQ = new double[numStates];
+ double maxDiagsQ = 0.0;
+
+ // Initialise the solution vector with an equiprobable distribution
+ // over the BSCC states.
+ // Additionally, compute the diagonal entries of the generator matrix Q.
+ // Recall that the entries of the generator matrix are given by
+ // Q(s,t) = prob(s,t) for s != t
+ // and Q(s,s) = -sum_{s!=t} prob(s,t),
+ // i.e., diagsQ[s] = -sum_{s!=t} prob(s,t).
+ // Furthermore, compute max |diagsQ[s]|.
+ double equiprob = 1.0 / states.cardinality();
+ for (OfInt iter = bscc.iterator(); iter.hasNext();) {
+ int state = iter.nextInt();
+
+ // Equiprobable for BSCC states.
+ soln[state] = equiprob;
+
+ // Note: diagsQ[state] = 0.0, as it was freshly created
+ // Compute negative exit rate (ignoring a possible self-loop)
+ dtmc.forEachTransition(state, (s, t, prob) -> {
+ if (s != t) {
+ diagsQ[state] -= prob;
+ }
+ });
- // Initialise solution vectors. Equiprobable for BSCC states.
- double equiprob = 1.0 / bscc.cardinality();
- for (int i = bscc.nextSetBit(0); i >= 0; i = bscc.nextSetBit(i + 1))
- soln[i] = soln2[i] = equiprob;
+ // Note: If there are no outgoing transitions, diagsQ[state] = 0, which is fine
+
+ // Update maximal absolute diagonal entry value of Q
+ // As diagsQ[s] <= 0, Math.abs(diagsQ[s]) = -diagsQ[s]
+ maxDiagsQ = Math.max(maxDiagsQ, -diagsQ[state]);
+ }
+
+ // Compute preconditioning factor deltaT
+ // In William J. Stewart: "Introduction to the Numerical Solution of Markov Chains",
+ // deltaT = 0.99 / maxDiagsQ is proposed;
+ // in the symbolic engines deltaT is computed as 0.99 / max exit[s], i.e., where
+ // the denominator corresponds to the maximal exit rate (where self loops are included).
+ // Currently, use the same deltaT values as in the symbolic engines,
+ // so for DTMCs, as exit[s]=1 for all states, deltaT is 0.99:
+ double deltaT = 0.99;
+ // TODO: Test and switch to deltaT computed as below, should lead to faster convergence.
+ // double deltaT = 0.99 / maxDiagsQ;
ExportIterations iterationsExport = null;
if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) {
@@ -2377,15 +2430,18 @@ public class DTMCModelChecker extends ProbModelChecker
mainLog.println("Exporting iterations to " + iterationsExport.getFileName());
}
+ // create copy of the solution vector
+ double[] soln2 = soln.clone();
+
// Start iterations
int iters = 0;
boolean done = false;
while (!done && iters < maxIters) {
iters++;
- // Matrix-vector multiply
- dtmc.vmMult(soln, soln2);
+ // Do vector-matrix multiplication step in (deltaT*Q + I)
+ dtmc.vmMultPowerSteadyState(soln, soln2, diagsQ, deltaT, bscc);
// Check termination
- done = PrismUtils.doublesAreClose(soln, soln2, termCritParam, termCrit == TermCrit.ABSOLUTE);
+ done = PrismUtils.doublesAreClose(soln, soln2, bscc, termCritParam, termCrit == TermCrit.ABSOLUTE);
// Swap vectors for next iter
double[] tmpsoln = soln;
soln = soln2;
@@ -2400,20 +2456,39 @@ public class DTMCModelChecker extends ProbModelChecker
watch.stop();
mainLog.println("Power method: " + iters + " iterations in " + watch.elapsedSeconds() + " seconds.");
+ // normalise solution
+ PrismUtils.normalise(soln, bscc);
+
if (iterationsExport != null) {
+ // export the normalised vector
+ iterationsExport.exportVector(soln);
iterationsExport.close();
}
+ if (result != null && result != soln) {
+ // If result vector was passed in as method argument,
+ // it can be the case that result does not point to the current soln vector (most recent values)
+ // but to the soln2 vector.
+ // In that case, we copy the relevant values from soln to result.
+ for (OfInt iter = bscc.iterator(); iter.hasNext();) {
+ int state = iter.nextInt();
+ result[state] = soln[state];
+ }
+ }
+ // store only one result vector, free temporary vectors
+ result = soln;
+ soln = soln2 = null;
+
// Non-convergence is an error (usually)
if (!done && errorOnNonConverge) {
- String msg = "Iterative method did not converge within " + iters + " iterations.";
- msg += "\nConsider using a different numerical method or increasing the maximum number of iterations";
+ String msg = "Iterative method did not converge within " + iters + " iterations.\n" +
+ "Consider using a different numerical method or increasing the maximum number of iterations";
throw new PrismException(msg);
}
// Return results
ModelCheckerResult res = new ModelCheckerResult();
- res.soln = soln;
+ res.soln = result;
res.numIters = iters;
res.timeTaken = watch.elapsedSeconds();
return res;
diff --git a/prism/src/explicit/DTMCSparse.java b/prism/src/explicit/DTMCSparse.java
index f457debb..c4fee9c4 100644
--- a/prism/src/explicit/DTMCSparse.java
+++ b/prism/src/explicit/DTMCSparse.java
@@ -38,6 +38,7 @@ import java.util.PrimitiveIterator.OfInt;
import java.util.function.Function;
import common.IterableStateSet;
+import common.iterable.IterableInt;
import common.iterable.MappingIterator;
import explicit.rewards.MCRewards;
import prism.PrismException;
@@ -356,6 +357,39 @@ public class DTMCSparse extends DTMCExplicit
}
}
+ @Override
+ public void vmMultPowerSteadyState(double vect[], double result[], double[] diagsQ, double deltaT, IterableInt states)
+ {
+ // Recall that the generator matrix Q has entries
+ // Q(s,s) = -sum_{t!=s} prob(s,t)
+ // and Q(s,t) = prob(s,t) for s!=t
+ // The values Q(s,s) are passed in via the diagsQ vector, while the
+ // values Q(s,t) correspond to the normal transitions
+
+ // Initialise result for relevant states to vect[s] * (deltaT * diagsQ[s] + 1),
+ // i.e., handle the product with the diagonal entries of (deltaT * Q) + I
+ for (OfInt it = states.iterator(); it.hasNext(); ) {
+ int state = it.nextInt();
+ result[state] = vect[state] * ((deltaT * diagsQ[state]) + 1.0);
+ }
+
+ // For each relevant state...
+ for (OfInt it = states.iterator(); it.hasNext(); ) {
+ int state = it.nextInt();
+
+ // ... handle all Q(state,t) entries of the generator matrix
+ for (int i=rows[state], stop=rows[state+1]; i < stop; i++) {
+ int target = columns[i];
+ double prob = probabilities[i];
+ if (state != target) {
+ // ignore self loop, diagonal entries of the generator matrix handled above
+ // update result vector entry for the *successor* state
+ result[target] += deltaT * prob * vect[state];
+ }
+ }
+ }
+ }
+
//--- Object ---