s14h
/
prism-accumulation
1
0
Fork 0
Code Releases Activity
Browse Source

Add preprocessing for steady-state in explicit engine 

To guarantee convergence, the power method requires the precomputation
P = (Q * deltaT + I)
from: William J. Stewart: Introduction to the Numerical Solution of Markov Chains p. 124.
master
Steffen Märcker 8 years ago
committed by Dave Parker
parent
c25c745ae9
commit
76ab6e2d76
3 changed files with 179 additions and 14 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 56
      prism/src/explicit/DTMC.java
  2. 103
      prism/src/explicit/DTMCModelChecker.java
  3. 34
      prism/src/explicit/DTMCSparse.java

56
prism/src/explicit/DTMC.java
View File

@ -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.
* <p>
* Computes<br>
* {@code result = vect * P} with matrix
* {@code P = (Q * deltaT + I)} where<br/>
* {@code Q} is the generator matrix,
* {@code deltaT} the preconditioning factor and
* {@code I} is the the identity matrix.<br/>
* Please refer to <em>William J. Stewart: "Introduction to the Numerical Solution of Markov Chains"</em> p.124 for details.
* </p>
* <p>
* 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).
* </p>
* @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];
}
});
}
}
}

103
prism/src/explicit/DTMCModelChecker.java
View File

@ -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.
* <p>
* To ensure convergence, we use the iteration matrix<br/>
* {@code P = (Q * deltaT + I)} where<br/>
* {@code Q} is the generator matrix,
* {@code deltaT} a preconditioning factor and
* {@code I} is the the identity matrix.<br/>
* See <em>William J. Stewart: "Introduction to the Numerical Solution of Markov Chains"</em> p.124 for details.
* </p>
* @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;

34
prism/src/explicit/DTMCSparse.java
View File

@ -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 ---

Write Preview
Loading…
Cancel
Save