diff --git a/prism/src/explicit/BirthProcess.java b/prism/src/explicit/BirthProcess.java new file mode 100644 index 00000000..fc7bc648 --- /dev/null +++ b/prism/src/explicit/BirthProcess.java @@ -0,0 +1,264 @@ +//============================================================================== +// +// Copyright (c) 2013- +// Authors: +// * Dave Parker (University of Oxford) +// * Frits Dannenberg (University of Oxford) +// * Ernst Moritz Hahn (University of Oxford) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package explicit; + +import java.util.ArrayList; +import prism.PrismException; + +/** + * Class to efficiently compute transient probabilities of a birth process. + * + * @author Dave Parker (University of Oxford) + * @author Frits Dannenberg (University of Oxford) + * @author Ernst Moritz Hahn (University of Oxford) + */ +public class BirthProcess +{ + /* uniformisation rate to compute probabilities. + * Must be at least as large as largest rate. */ + double unifRate; + /* precision of computations using Fox-Glynn algorithm */ + double termCritParam; + /* time to compute probabilities for */ + double time; + /* values used to compute probability to be in a given stage */ + double[] probs; + /* same as above */ + double[] newProbs; + /* whether rates added by calculateNextRate will be stored */ + boolean withRateArray; + ArrayList jumpRates; + boolean initialising; + /* current birth process stage */ + int stageNr; + /* whether to try to avoid birth process computations in + * case all rates are the same + */ + boolean avoidBirthComputation; + + /** + * Construct birth process. + */ + public BirthProcess() + { + stageNr = 0; + termCritParam = 1E-7; + withRateArray = true; + initialising = true; + avoidBirthComputation = true; + } + + /** + * Sets whether rates added by calculateNextRate will be stored. + * The default is to do so. Setting it to false allows to save space, + * but will not allow rates larger than @a unifRate to be added. + * + * @param withRateArray whether to store rates in array + */ + public void setWithRateArray(boolean withRateArray) + { + if (!initialising) { + throw new IllegalArgumentException("this method might not be called after calculateNextRate"); + } + this.withRateArray = withRateArray; + } + + /** + * Sets the time to compute probabilities for. + * + * @param time time to compute probabilities for + */ + public void setTime(double time) + { + if (!initialising) { + throw new IllegalArgumentException("this method might not be called after calculateNextRate"); + } + if (time < 0.0) { + throw new IllegalArgumentException("time must be nonnegative"); + } + this.time = time; + } + + /** + * Sets precision to be used to compute probabilities. + * + * @param termCritParam precision to be used to compute probabilities + */ + public void setTermCritParam(double termCritParam) + { + if (!initialising) { + throw new IllegalArgumentException("this method might not be called after calculateNextRate"); + } + + this.termCritParam = termCritParam; + } + + /** + * Chooses whether to try to avoid birth process construction. + * In case all rates provided to calculateNextProb are the same, it is + * possible to use the Fox-Glynn algorithm to compute probabilities in + * the birth process, which is then a Poisson process. If delayBirthComputation + * is true, the more expensive computations will only be computed in case + * this is necessary. Thus, if calculateNextProb is called with the same + * rate all the time, computations are significantly faster. + * + * @param avoidBirthComputation true iff birth process construction shall be delayed + */ + public void setAvoidBirthComputation(boolean avoidBirthComputation) + { + if (!initialising) { + throw new IllegalArgumentException("this method might not be called after calculateNextRate"); + } + this.avoidBirthComputation = avoidBirthComputation; + } + + /** + * Computes probability to reside in next process stage at given time. + * If this is the nth call to the process, computes the probability + * to reside in the nth stage of the birth process. The rates of the + * birth process up to the n-1th stage have been set by previous calls + * to the function, the nth rate is set by the current call. The time + * probabilities are computed for must have been set previously by + * setTime. + * + * @param rate leaving rate of current state + * @return probability to reside in current stage at given time + * @throws PrismException + */ + public double calculateNextProb(double rate) throws PrismException + { + if (initialising && 0.0 == unifRate && !withRateArray) { + throw new IllegalArgumentException("unifRate must be set if withRateArray is false"); + } + if (withRateArray && initialising) { + jumpRates = new ArrayList(); + } + initialising = false; + if (!withRateArray && rate > unifRate) { + throw new IllegalArgumentException("cannot use rates larger than initial rate if withRateArray is false"); + } + if (withRateArray) { + jumpRates.add(rate); + } + boolean recompute = false; + if (rate > unifRate) { + if (!avoidBirthComputation) { + recompute = true; + unifRate = rate * 1.25 * 1.02; + } else { + unifRate = rate; + } + } + if ((jumpRates.size() != 1) && (Math.abs(rate - jumpRates.get(jumpRates.size() - 2)) > 1E-100)) { + if (avoidBirthComputation) { + recompute = true; + } + avoidBirthComputation = false; + } + + if (null == probs || recompute) { + initPoisson(); + } + double result = 0.0; + if (recompute) { + for (stageNr = 0; stageNr < jumpRates.size(); stageNr++) { + result = compNextStageProb(jumpRates.get(stageNr)); + } + } else { + if (avoidBirthComputation) { + result = (stageNr < probs.length) ? probs[stageNr] : 0.0; + stageNr++; + } else { + result = compNextStageProb(rate); + stageNr++; + } + } + + return result; + } + + private double compNextStageProb(double rate) + { + assert(rate > 0.0); + double prob = rate / unifRate; // p = r / q + double omprob = 1.0 - prob; // 1-p + double result = 0.0; + + double omprobtti = 1.0; // (1-p)^i + for (int i = 0; i < probs.length; i++) { + result += omprobtti * probs[i]; + omprobtti *= omprob; + } + + newProbs[newProbs.length - 1] = 0.0; + for (int i = probs.length - 1; i >= 1; i--) { + newProbs[i - 1] = newProbs[i] * omprob + probs[i] * prob; + } + double[] temp = probs; + probs = newProbs; + newProbs = temp; + return result; + } + + /** + * Initialises probability vectors by Poisson probabilities. + */ + private void initPoisson() throws PrismException + { + long left, right; + double qt = unifRate * time; + double acc = termCritParam / 8.0; + double[] weights; + double totalWeight; + if (unifRate * time == 0.0) { + left = 0; + right = 0; + totalWeight = 1.0; + weights = new double[1]; + weights[0] = 1.0; + } else { + FoxGlynn fg = new FoxGlynn(qt, 1e-300, 1e+300, acc); + left = fg.getLeftTruncationPoint(); + right = fg.getRightTruncationPoint(); + if (right < 0) { + throw new PrismException("Overflow in Fox-Glynn computation (time bound too big?)"); + } + weights = fg.getWeights(); + totalWeight = fg.getTotalWeight(); + } + for (long i = left; i <= right; i++) { + weights[(int) (i - left)] /= totalWeight; + } + probs = new double[(int) (right + 1)]; + newProbs = new double[(int) (right + 1)]; + for (long entryNr = left; entryNr <= right; entryNr++) { + probs[(int) entryNr] = weights[(int) (entryNr - left)]; + } + } +} diff --git a/prism/src/explicit/FastAdaptiveUniformisation.java b/prism/src/explicit/FastAdaptiveUniformisation.java new file mode 100644 index 00000000..b47c1f59 --- /dev/null +++ b/prism/src/explicit/FastAdaptiveUniformisation.java @@ -0,0 +1,1206 @@ +//============================================================================== +// +// Copyright (c) 2013- +// Authors: +// * Dave Parker (University of Oxford) +// * Frits Dannenberg (University of Oxford) +// * Ernst Moritz Hahn (University of Oxford) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package explicit; + +import java.io.File; +import java.util.ArrayList; +import java.util.Arrays; +import java.util.HashMap; +import java.util.HashSet; +import java.util.LinkedHashMap; +import java.util.ListIterator; +import java.util.Map; + +import parser.ast.Expression; +import parser.ast.ExpressionIdent; +import parser.ast.LabelList; +import parser.ast.RewardStruct; +import parser.type.TypeDouble; +import parser.Values; + +import parser.State; +import prism.*; + +/* + * TODO + * - add options for state removal, e.g. + * - by delta (as current) + * - by max probability loss per iteration (requires sort by prob) + * - by max number of states (requires sort by prob) + * - compress states to Bitset (memory waste currently excessive e.g. for mapk) + * - do not delete states immediately but only after they have been below + * delta for a specified number of iterations to avoid deleting and exploring + * the same states over and over again + * - dynamic adaption of interval width (with of 1 seems to work best, however) + * - improve birth process - only worth it if we find case study where it makes a + * difference, but could contribute to publicability then + * - plot number of active states over time for mapk to get idea of model behaviour + * - in gen-dat.pl, mark runs as dead if we can derive that they cannot succeed + * - discuss public interface with Dave + * - make stop of deletion after half of Birth threshold reached optional + * - check whether it's worth storing incoming transitions rather than outgoing + * would be faster, but edges are more costly to remove + * - if we reach a point where we only delete states and don't add new ones, it + * might make sense to switch to array representation and ignore the fact that + * we could delete further states + */ + +/** + * Implementation of fast adaptive uniformisation + * @author Dave Parker (University of Oxford) + * @author Frits Dannenberg (University of Oxford) + * @author Ernst Moritz Hahn (University of Oxford) + */ +public class FastAdaptiveUniformisation +{ + /** + * Stores properties of states needed for fast adaptive method. + * This includes the current-step probability, next-state probability, + * and the transient probability (sum of step probabilities weighted + * with birth process distributions). It also contains the list of successor + * states and the rates to them, the number of incoming transitions + * (references) and a flag whether the state has a significant probability + * mass (alive). + * + * @author Ernst Moritz Hahn (University of Oxford) + */ + private final class StateProp + { + /** current-step probability. + * should contain initial probability before actual analysis is started. + * will contain transient probability after analysis. */ + private double prob; + /** next-state probability */ + private double nextProb; + /** sum probability weighted with birth process distribution */ + private double sum; + /** reward of this state */ + private double reward; + /** rates to successor states */ + private double[] succRates; + /** successor states */ + private StateProp[] succStates; + /** number of incoming transitions of relevant states */ + private int references; + /** true if and only if state probability above relevance threshold */ + private boolean alive; + + /** + * Constructs a new state property object. + */ + StateProp() + { + prob = 0.0; + nextProb = 0.0; + prob = 0.0; + reward = 0.0; + references = 0; + alive = true; + succRates = null; + succStates = null; + totalProbLoss = 0.0; + } + + /** + * Set current state probability. + * + * @param prob current state probability to set + */ + void setProb(double prob) + { + this.prob = prob; + } + + /** + * Gets current state probability. + * + * @return current state probability + */ + double getProb() + { + return prob; + } + + /** + * Sets next state probability. + * + * @param nextProb next state probability to set + */ + void setNextProb(double nextProb) + { + this.nextProb = nextProb; + } + + /** + * Adds value to next state probability. + * + * @param add value to add to next state probability + */ + void addToNextProb(double add) + { + this.nextProb += add; + } + + /** + * Sets weighted sum probability. + * + * @param sum weighted sum probability to set. + */ + void setSum(double sum) + { + this.sum = sum; + } + + /** + * Adds current probability times {@code poisson} to weighted sum probability. + * + * @param poisson this value times current probability will be added to sum probability + */ + void addToSum(double poisson) + { + sum += poisson * prob; + } + + /** + * Gets weighted sum probability. + * + * @return weighted sum probability + */ + double getSum() + { + return sum; + } + + /** + * Prepares next iteration step. + * Sets current probability to next probability, and sets next + * probability to zero. + */ + void prepareNextIteration() + { + prob = nextProb; + nextProb = 0.0; + } + + /** + * Set state reward. + * + * @param reward state reward to set + */ + void setReward(double reward) + { + this.reward = reward; + } + + /** + * Get state reward. + * + * @return state reward + */ + double getReward() + { + return reward; + } + + /** + * Sets rates to successor states. + * Expects an array of rates, so that the rate to the successor + * state set by {@setSuccStates} is the one given by the corresponding + * index. The value {@code null} is allowed here. + * + * @param succRates rates to successor states. + */ + void setSuccRates(double[] succRates) + { + this.succRates = succRates; + } + + /** + * Sets successor states. + * Expects an array of successor states, so that the rate set by + * {@setSuccRates} is the one given by the corresponding index. + * The value {@code null} is allowed here. + * + * @param succStates successor states + */ + void setSuccStates(StateProp[] succStates) + { + this.succStates = succStates; + if (succStates != null) { + for (int succNr = 0; succNr < succStates.length; succNr++) { + succStates[succNr].incReferences(); + } + } + } + + /** + * Returns number of successor states of this state. + * + * @return number of successor states + */ + int getNumSuccs() + { + if (succRates == null) { + return 0; + } else { + return succRates.length; + } + } + + /** + * Gets successor rates. + * + * @return successor rates + */ + double[] getSuccRates() + { + return succRates; + } + + /** + * Gets successor states. + * + * @return successor states + */ + StateProp[] getSuccStates() + { + return succStates; + } + + /** + * Sets whether state is alive. + * + * @param alive whether state should be set to being alive + */ + void setAlive(boolean alive) + { + this.alive = alive; + } + + /** + * Checks whether state is alive. + * + * @return true iff state is alive + */ + boolean isAlive() + { + return alive; + } + + /** + * Increments the number of references of this state. + * The number of references should correspond to the number of alive + * states which have this state as successor state. + */ + void incReferences() + { + references++; + } + + /** + * Decrements the number of references of this state. + * The number of references should correspond to the number of alive + * states which have this state as successor state. + */ + void decReferences() + { + references--; + } + + /** + * Deletes this state. + * This means basically removing all of its successors. Beforehand, + * their reference counter is decreased, because this state does no + * longer count as a model state. It is left in the model however, + * because it might still be the successor state of some alive state. + */ + void delete() + { + if (null != succStates) { + for (int succNr = 0; succNr < succStates.length; succNr++) { + succStates[succNr].decReferences(); + } + } + succStates = null; + succRates = null; + alive = false; + prob = 0.0; + nextProb = 0.0; + } + + /** + * Checks whether this state can be removed. + * This is only the case if its probability is below the threshold + * specified, and then only if there are no transitions from alive + * states into this state. + * + * @return true if and only if this state can be removed + */ + boolean canRemove() + { + return !alive && (0 == references); + } + + /** + * Checks whether this state has successors or not. + * Will be true if and only if successor state array is nonnull. + * + * @return whether this state has successors or not + */ + boolean hasSuccs() + { + return succStates != null; + } + + /** + * Returns the sum of all rates leaving to successor states. + * + * @return sum of all rates leaving to successor states + */ + double sumRates() + { + if (null == succRates) { + return 0.0; + } + double sumRates = 0.0; + for (int rateNr = 0; rateNr < succRates.length; rateNr++) { + sumRates += succRates[rateNr]; + } + return sumRates; + } + } + + /** + * Enum to store type of analysis to perform. + */ + public enum AnalysisType { + /** transient probability distribution */ + TRANSIENT, + /** reachability, for "F" or "U" PCTL operators */ + REACH, + /** instantaneous rewards */ + REW_INST, + /** cumulative rewards */ + REW_CUMUL + } + + /** PRISM settings to read analysis parameters from */ + private PrismSettings settings = null; + /** model exploration component to generate new states */ + private ModelExplorer modelExplorer; + /** probability allowed to drop birth process */ + private double termCritParam; + /** probability threshold when to drop states in discrete-time process */ + private double delta; + /** number of intervals to divide time into */ + private int numIntervals; + /** iterations after which switch to sparse matrix if no new/dropped states */ + private int arrayThreshold; + + /** reward structure to use for analysis */ + private RewardStruct rewStruct = null; + /** result value of analysis */ + private double value; + /** model constants */ + private Values constantValues = null; + /** maps from state (assignment of variable values) to property object */ + private LinkedHashMap states; + /** states for which successor rates are to be computed */ + private ArrayList addDistr; + /** states which are to be deleted */ + private ArrayList deleteStates; + /** initial size of state hash map */ + private final int initSize = 3000; + /** maximal total leaving rate of all states alive */ + private double maxRate = 0.0; + /** target state set - used for reachability (until or finally properties) */ + private Expression target; + /** number of consecutive iterations without new states are state drops */ + private int itersUnchanged; + /** sum of probabilities in stages of birth process seen so far */ + private double birthProbSum; + /** birth process used for time discretisation */ + private BirthProcess birthProc; + /** states which fulfill this will be made absorbing - for until props */ + private Expression sink; + /** if true, don't drop further states. + * Used to avoid excessive probability loss in some cases. */ + private boolean keepSumProb; + /** maximal number of states ever stored during analysis */ + private int maxNumStates; + /** list of special labels we need to maintain, like "init", "deadlock", etc. */ + private LabelList specialLabels; + /** set of initial states of the model */ + private HashSet initStates; + /** type of analysis to perform */ + private AnalysisType analysisType; + /** total loss of probability in discrete-time process */ + private double totalProbLoss; + + /** + * Constructor. + */ + public FastAdaptiveUniformisation(PrismSettings settings, ModelExplorer modelExplorer) throws PrismException + { + maxNumStates = 0; + this.settings = settings; + this.modelExplorer = modelExplorer; + + termCritParam = settings.getDouble(PrismSettings.PRISM_TERM_CRIT_PARAM); + delta = settings.getDouble(PrismSettings.PRISM_FAU_DELTA); + numIntervals = settings.getInteger(PrismSettings.PRISM_FAU_INTERVALS); + arrayThreshold = settings.getInteger(PrismSettings.PRISM_FAU_ARRAYTHRESHOLD); + analysisType = AnalysisType.TRANSIENT; + rewStruct = null; + target = Expression.False(); + sink = Expression.False(); + specialLabels = new LabelList(); + specialLabels.addLabel(new ExpressionIdent("deadlock"), new ExpressionIdent("deadlock")); + specialLabels.addLabel(new ExpressionIdent("init"), new ExpressionIdent("init")); + } + + /** + * Sets analysis type to perform. + * + * @param analysisType analysis type to perform + */ + public void setAnalysisType(AnalysisType analysisType) + { + this.analysisType = analysisType; + } + + /** + * Sets values for model constants. + * + * @param constantValues values for model constants + */ + public void setConstantValues(Values constantValues) + { + this.constantValues = constantValues; + } + + /** + * Sets reward structure to use. + * + * @param rewStruct reward structure to use + */ + public void setRewardStruct(RewardStruct rewStruct) + { + this.rewStruct = rewStruct; + } + + /** + * + * @param target + */ + public void setTarget(Expression target) + { + this.target = target; + } + + /** + * Returns maximal number of states used during analysis. + * + * @return maximal number of states used during analysis + */ + public int getMaxNumStates() + { + return maxNumStates; + } + + /** + * Returns the value of the analysis. + * For reachability analyses, this is the probability to reach state in + * the reach set, for instantaneous reward properties this is the + * instantaneous reward and for cumulative reward analysis it is the + * cumulative reward. For the computation of transient probabilities + * without doing model checking, this value is not significant. + * + * @return value of the analysis + */ + public double getValue() + { + return value; + } + + /** + * Sets which states shall be treated as sink states. + * To be used for properties like "a U<=T b" where states "b || !a" have + * to be made absorbing. + * + * @param sink expressing stating which states are sink states + * @throws PrismException thrown if problems in underlying function occurs + */ + public void setSink(Expression sink) throws PrismException + { + this.sink = sink; + if (states != null) { + for (Map.Entry statePair : states.entrySet()) { + State state = statePair.getKey(); + StateProp prop = statePair.getValue(); + modelExplorer.queryState(state); + specialLabels.setLabel(0, modelExplorer.getNumTransitions() == 0 ? Expression.True() : Expression.False()); + specialLabels.setLabel(1, initStates.contains(state) ? Expression.True() : Expression.False()); + Expression evSink = sink.deepCopy(); + evSink = (Expression) evSink.expandLabels(specialLabels); + if (evSink.evaluateBoolean(constantValues, state)) { + double[] succRates = new double[1]; + StateProp[] succStates = new StateProp[1]; + succRates[0] = 1.0; + succStates[0] = states.get(state); + prop.setSuccRates(succRates); + prop.setSuccStates(succStates); + } + } + } + } + + /** + * Get the number of states in the current window. + */ + public int getNumStates() + { + return states.size(); + } + + /** + * Compute transient probability distribution (forwards). + * Start from initial state (or uniform distribution over multiple initial states). + */ + public StateValues doTransient(double time) throws PrismException + { + return doTransient(time, (StateValues) null); + } + + /** + * Compute transient 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). + * @param t Time point + * @param initDistFile File containing initial distribution + */ + public StateValues doTransient(double t, File initDistFile) throws PrismException + { + StateValues initDist = null; + return doTransient(t, initDist); + } + + /** + * Compute transient probability distribution (forwards). + * Use the passed in vector initDist as the initial probability distribution (time 0). + * In case initDist is null starts at the default initial state with prob 1. + * + * @param time Time point + * @param initDist Initial distribution + */ + public StateValues doTransient(double time, StateValues initDist) throws PrismException + { + if (initDist == null) { + initDist = new StateValues(); + initDist.type = TypeDouble.getInstance(); + initDist.size = 1; + initDist.valuesD = new double[1]; + initDist.statesList = new ArrayList(); + initDist.valuesD[0] = 1.0; + initDist.statesList.add(modelExplorer.getDefaultInitialState()); + } + + /* prepare fast adaptive uniformisation */ + addDistr = new ArrayList(); + deleteStates = new ArrayList(); + states = new LinkedHashMap(initSize); + value = 0.0; + initStates = new HashSet(); + ListIterator it = initDist.statesList.listIterator(); + double[] values = initDist.getDoubleArray(); + maxRate = 0.0; + for (int stateNr = 0; stateNr < initDist.size; stateNr++) { + State initState = it.next(); + addToModel(initState); + computeStateRatesAndRewards(initState); + states.get(initState).setProb(values[stateNr]); + maxRate = Math.max(maxRate, states.get(initState).sumRates() * 1.02); + } + + /* run fast adaptive uniformisation */ + computeTransientProbsAdaptive(time); + + /* prepare and return results */ + ArrayList statesList = new ArrayList(states.size()); + double[] probsArr = new double[states.size()]; + int probsArrEntry = 0; + for (Map.Entry statePair : states.entrySet()) { + statesList.add(statePair.getKey()); + probsArr[probsArrEntry] = statePair.getValue().getProb(); + probsArrEntry++; + } + StateValues probs = new StateValues(); + probs.type = TypeDouble.getInstance(); + probs.size = probsArr.length; + probs.valuesD = probsArr; + probs.statesList = statesList; + return probs; + } + + /** + * Compute transient probabilities using fast adaptive uniformisation + * Compute the probability of being in each state at time {@code t}. + * If corresponding options are set, also computes cumulative rewards. + * For space efficiency, the initial distribution vector will be modified and values over-written, + * so if you wanted it, take a copy. + * @param time time point + */ + public void computeTransientProbsAdaptive(double time) throws PrismException + { + if (addDistr == null) { + addDistr = new ArrayList(); + deleteStates = new ArrayList(); + states = new LinkedHashMap(initSize); + value = 0.0; + prepareInitialDistribution(); + } + + double initIval = settings.getDouble(PrismSettings.PRISM_FAU_INITIVAL); + if (time - initIval < 0.0) { + initIval = 0.0; + } + if (initIval != 0.0) { + iterateAdaptiveInterval(initIval); + for (StateProp prop : states.values()) { + prop.setProb(prop.getSum()); + prop.setSum(0.0); + prop.setNextProb(0.0); + } + updateStates(); + } + + for (int ivalNr = 0; ivalNr < numIntervals; ivalNr++) { + double interval = (time - initIval) / numIntervals; + iterateAdaptiveInterval(interval); + for (StateProp prop : states.values()) { + prop.setProb(prop.getSum()); + prop.setSum(0.0); + prop.setNextProb(0.0); + } + updateStates(); + } + if (AnalysisType.REW_INST == analysisType) { + for (StateProp prop : states.values()) { + value += prop.getProb() * prop.getReward(); + } + } else { + for (Map.Entry statePair : states.entrySet()) { + State state = statePair.getKey(); + StateProp prop = statePair.getValue(); + modelExplorer.queryState(state); + specialLabels.setLabel(0, modelExplorer.getNumTransitions() == 0 ? Expression.True() : Expression.False()); + specialLabels.setLabel(1, initStates.contains(state) ? Expression.True() : Expression.False()); + Expression evTarget = target.deepCopy(); + evTarget = (Expression) evTarget.expandLabels(specialLabels); + if (AnalysisType.REACH == analysisType) { + value += prop.getProb() * (evTarget.evaluateBoolean(constantValues, state) ? 1.0 : 0.0); + } + } + } + } + + /** + * Performs fast adaptive uniformisation for a single time interval. + * + * @param interval duration of time interval + * @throws PrismException + */ + private void iterateAdaptiveInterval(double interval) throws PrismException + { + birthProc = new BirthProcess(); + birthProc.setTime(interval); + birthProc.setTermCritParam(termCritParam); + + int iters = 0; + birthProbSum = 0.0; + itersUnchanged = 0; + keepSumProb = false; + while (birthProbSum < (1 - termCritParam)) { + if (birthProbSum >= termCritParam/2) { + keepSumProb = true; + } + if ((itersUnchanged == arrayThreshold)) { + iters = arrayIterate(iters); + } else { + long birthProcTimer = System.currentTimeMillis(); + double prob = birthProc.calculateNextProb(maxRate); + birthProcTimer = System.currentTimeMillis() - birthProcTimer; + birthProbSum += prob; + collectValuePostIter(prob, birthProbSum); + for (StateProp prop : states.values()) { + prop.addToSum(prob); + } + + mvMult(maxRate); + updateStates(); + iters++; + } + } + + computeTotalDiscreteLoss(); + } + + /** + * Transforms the current submodel to array form. + * In case there are no further changes in the states discovered, or + * further states only become relevant after a large number of + * iterations, this allows the analysis to be performed much faster. + * After the analysis has finished or after it has to be terminated as + * formerly irrelevant states become relevant, results are mapped back + * to the original data structure. The method returns the current + * iteration. + * + * In case border states become + * relevant, this data structure can + * + * @param iters current iteration number + * @return current iteration after termination of this method + * @throws PrismException thrown if problems in underlying methods occur + */ + private int arrayIterate(int iters) throws PrismException + { + /* build backwards matrix and map values */ + int numStates = states.size(); + int numTransitions = 0; + for (StateProp prop : states.values()) { + numTransitions += prop.getNumSuccs() + 1; + } + int stateNr = 0; + HashMap stateToNumber = new HashMap(numStates); + StateProp[] numberToState = new StateProp[numStates]; + for (StateProp prop : states.values()) { + if (prop.isAlive()) { + stateToNumber.put(prop, stateNr); + numberToState[stateNr] = prop; + stateNr++; + } + } + int numAlive = stateNr; + for (StateProp prop : states.values()) { + if (!prop.isAlive()) { + stateToNumber.put(prop, stateNr); + numberToState[stateNr] = prop; + stateNr++; + } + } + + double[] inProbs = new double[numTransitions]; + int[] rows = new int[numStates + 1]; + int[] cols = new int[numTransitions]; + double[] outRates = new double[numStates]; + for (StateProp prop : states.values()) { + StateProp[] succStates = prop.getSuccStates(); + if (succStates != null) { + for (StateProp succ : succStates) { + rows[stateToNumber.get(succ) + 1]++; + } + } + rows[stateToNumber.get(prop) + 1]++; + } + for (stateNr = 0; stateNr < numStates; stateNr++) { + rows[stateNr + 1] += rows[stateNr]; + } + + for (StateProp prop : states.values()) { + int stateNumber = stateToNumber.get(prop); + StateProp[] succStates = prop.getSuccStates(); + double[] succRates = prop.getSuccRates(); + if (succStates != null) { + for (int i = 0; i < succStates.length; i++) { + StateProp succState = succStates[i]; + int succStateNumber = stateToNumber.get(succState); + double succRate = succRates[i]; + cols[rows[succStateNumber]] = stateNumber; + inProbs[rows[succStateNumber]] = succRate / maxRate; + rows[succStateNumber]++; + outRates[stateNumber] += succRate; + } + } + } + + for (stateNr = 0; stateNr < numStates; stateNr++) { + cols[rows[stateNr]] = stateNr; + inProbs[rows[stateNr]] = (maxRate - outRates[stateNr]) / maxRate; + } + + Arrays.fill(rows, 0); + for (StateProp prop : states.values()) { + StateProp[] succStates = prop.getSuccStates(); + if (succStates != null) { + for (StateProp succ : succStates) { + rows[stateToNumber.get(succ) + 1]++; + } + } + rows[stateToNumber.get(prop) + 1]++; + } + for (stateNr = 0; stateNr < numStates; stateNr++) { + rows[stateNr + 1] += rows[stateNr]; + } + + double[] rewards = new double[numStates]; + double[] probs = new double[numStates]; + double[] nextProbs = new double[numStates]; + double[] sum = new double[numStates]; + for (stateNr = 0; stateNr < numberToState.length; stateNr++) { + StateProp prop = numberToState[stateNr]; + if (analysisType == AnalysisType.REW_CUMUL) { + rewards[stateNr] = prop.getReward(); + } + probs[stateNr] = prop.getProb(); + sum[stateNr] = prop.getSum(); + } + + /* iterate using matrix */ + boolean canArray = true; + while (birthProbSum < (1 - 1E-9) && canArray) { + // timer2 = System.currentTimeMillis(); + double prob = birthProc.calculateNextProb(maxRate); + birthProbSum += prob; + double mixed = (1.0 - birthProbSum) / maxRate; + for (stateNr = 0; stateNr < numStates; stateNr++) { + value += probs[stateNr] * mixed * rewards[stateNr]; + sum[stateNr] += prob * probs[stateNr]; + nextProbs[stateNr] = 0.0; + for (int succNr = rows[stateNr]; succNr < rows[stateNr+1]; succNr++) { + nextProbs[stateNr] += inProbs[succNr] * probs[cols[succNr]]; + } + if ((stateNr < numAlive) != (nextProbs[stateNr] > delta)) { + canArray = false; + } else if (stateNr >= numAlive) { + nextProbs[stateNr] = 0.0; + } + } + double[] swap = probs; + probs = nextProbs; + nextProbs = swap; + + iters++; + } + + /* map back, update states and return current iteration */ + for (stateNr = 0; stateNr < numberToState.length; stateNr++) { + StateProp prop = numberToState[stateNr]; + prop.setProb(probs[stateNr]); + prop.setSum(sum[stateNr]); + } + updateStates(); + return iters; + } + + /** + * Update analysis value after iteration. + * For certain analyses (currently cumulative rewards) we have to modify + * the analysis value after each iteration. + * + * @param prob + * @param probSum + */ + private void collectValuePostIter(double prob, double probSum) + { + switch (analysisType) { + case TRANSIENT: + // nothing to do here, we're just computing distributions + break; + case REACH: + // nothing to do here, we're collecting values later on + break; + case REW_INST: + // nothing to do here, we're collecting rewards later on + break; + case REW_CUMUL: + double mixed = (1.0 - probSum) / maxRate; + for (StateProp prop : states.values()) { + value += prop.getProb() * mixed * prop.getReward(); + } + break; + } + + } + + /** + * Updates state values once a transient analysis of time interval finished. + * Deletes states which can be deleted according to their current + * probability and the threshold. Computes new maximal rate for remaining + * states. Computes transitions to successors of states which have become + * alive to to probability threshold only after a transient analysis has + * finished. + * + * @throws PrismException thrown if something goes wrong + */ + private void updateStates() throws PrismException + { + maxRate = 0.0; + addDistr.clear(); + for (Map.Entry statePair : states.entrySet()) { + State state = statePair.getKey(); + StateProp prop = statePair.getValue(); + if (prop.getProb() > delta) { + prop.setAlive(true); + if (!prop.hasSuccs()) { + itersUnchanged = 0; + addDistr.add(state); + } else { + maxRate = Math.max(maxRate, prop.sumRates()); + } + } else { + prop.delete(); + } + } + for (int stateNr = 0; stateNr < addDistr.size(); stateNr++) { + computeStateRatesAndRewards(addDistr.get(stateNr)); + maxRate = Math.max(maxRate, states.get(addDistr.get(stateNr)).sumRates()); + } + maxRate *= 1.02; + + removeDeletedStates(); + } + + /** + * Removes all states subject to removal. + * This affects states which both have a present-state probability below + * the given threshold, and do not have incoming transitions from states + * with a relevant probability mass. + */ + private void removeDeletedStates() + { + boolean unchanged = true; + for (Map.Entry statePair : states.entrySet()) { + State state = statePair.getKey(); + StateProp prop = statePair.getValue(); + if (prop.canRemove()) { + deleteStates.add(state); + unchanged = false; + } + } + if (!keepSumProb) { + for (int i = 0; i < deleteStates.size(); i++) { + states.remove(deleteStates.get(i)); + } + } + if (unchanged) { + itersUnchanged++; + } else { + itersUnchanged = 0; + } + deleteStates.clear(); + } + + /** + * Prepares initial distribution for the case of a single initial state. + * + * @throws PrismException + */ + private void prepareInitialDistribution() throws PrismException + { + initStates = new HashSet(); + State initState = modelExplorer.getDefaultInitialState(); + initStates.add(initState); + addToModel(initState); + computeStateRatesAndRewards(initState); + states.get(initState).setProb(1.0); + maxRate = states.get(initState).sumRates() * 1.02; + } + + /** + * Computes total sum of lost probabilities. + * + * @return total probability still in model + */ + public void computeTotalDiscreteLoss() + { + double totalProb = 0; + for (StateProp prop : states.values()) { + totalProb += prop.getSum(); + } + totalProbLoss = 1.0 - totalProb; + } + + /** + * Returns the total probability loss. + * + * @return + */ + public double getTotalDiscreteLoss() + { + return totalProbLoss; + } + + /** + * Adds @a state to model. + * Computes reward for this states, creates entry in map of states, + * and updates number of states + * + * @param state state to add + * @throws PrismException thrown if something wrong happens in underlying methods + */ + private void addToModel(State state) throws PrismException + { + StateProp prop = new StateProp(); + prop.setReward(computeRewards(state)); + states.put(state, prop); + maxNumStates = Math.max(maxNumStates, states.size()); + } + + /** + * Computes successor rates and rewards for a given state. + * Rewards computed depend on the reward structure set by + * {@code setRewardStruct}. + * + * @param state state to compute successor rates and rewards for + * @throws PrismException thrown if something goes wrong + */ + private void computeStateRatesAndRewards(State state) throws PrismException + { + double[] succRates; + StateProp[] succStates; + modelExplorer.queryState(state); + specialLabels.setLabel(0, modelExplorer.getNumTransitions() == 0 ? Expression.True() : Expression.False()); + specialLabels.setLabel(1, initStates.contains(state) ? Expression.True() : Expression.False()); + Expression evSink = sink.deepCopy(); + evSink = (Expression) evSink.expandLabels(specialLabels); + if (evSink.evaluateBoolean(constantValues, state)) { + succRates = new double[1]; + succStates = new StateProp[1]; + succRates[0] = 1.0; + succStates[0] = states.get(state); + } else { + int nt = modelExplorer.getNumTransitions(); + succRates = new double[nt]; + succStates = new StateProp[nt]; + for (int i = 0; i < nt; i++) { + State succState = modelExplorer.computeTransitionTarget(i); + StateProp succProp = states.get(succState); + if (null == succProp) { + addToModel(succState); + modelExplorer.queryState(state); + succProp = states.get(succState); + } + succRates[i] = modelExplorer.getTransitionProbability(i); + succStates[i] = succProp; + } + if (nt == 0) { + succRates = new double[1]; + succStates = new StateProp[1]; + succRates[0] = 1.0; + succStates[0] = states.get(state); + } + } + states.get(state).setSuccRates(succRates); + states.get(state).setSuccStates(succStates); + } + + /** + * Perform a single matrix-vector multiplication. + * + * @param maxRate maximal total leaving rate sum in living states + */ + private void mvMult(double maxRate) + { + for (StateProp prop : states.values()) { + double[] succRates = prop.getSuccRates(); + StateProp[] succStates = prop.getSuccStates(); + double stateProb = prop.getProb(); + if (null != succStates) { + double sumRates = 0.0; + for (int succ = 0; succ < succStates.length; succ++) { + double rate = succRates[succ]; + sumRates += rate; + succStates[succ].addToNextProb((rate / maxRate) * stateProb); + } + prop.addToNextProb(((maxRate - sumRates) / maxRate) * prop.getProb()); + } + } + for (StateProp prop : states.values()) { + prop.prepareNextIteration(); + } + } + + /** + * Checks if rewards are needed for analysis. + * + * @return true if and only if rewards are needed + */ + private boolean isRewardAnalysis() + { + return (analysisType == AnalysisType.REW_INST) + || (analysisType == AnalysisType.REW_CUMUL); + } + + /** + * Computes the reward for a given state. + * In case a cumulative reward analysis is to be performed, transition + * rewards are transformed into equivalent state rewards. + * + * @param state the state to compute the reward of + * @return the reward for state @a state + * @throws PrismException thrown if problems occur in PRISM functions called + */ + private double computeRewards(State state) throws PrismException + { + if (!isRewardAnalysis()) { + return 0.0; + } + int numTransitions = 0; + if (AnalysisType.REW_CUMUL == analysisType) { + modelExplorer.queryState(state); + numTransitions = modelExplorer.getNumTransitions(); + } + double sumReward = 0.0; + int numStateItems = rewStruct.getNumItems(); + for (int i = 0; i < numStateItems; i++) { + Expression guard = rewStruct.getStates(i); + if (guard.evaluateBoolean(constantValues, state)) { + double reward = rewStruct.getReward(i).evaluateDouble(constantValues, state); + String action = rewStruct.getSynch(i); + if (action != null) { + if (AnalysisType.REW_CUMUL == analysisType) { + for (int j = 0; j < numTransitions; j++) { + String tAction = modelExplorer.getTransitionAction(j); + if (tAction == null) { + tAction = ""; + } + if (tAction.equals(action)) { + sumReward += reward * modelExplorer.getTransitionProbability(j); + } + } + } + } else { + sumReward += reward; + } + } + } + + return sumReward; + } +} diff --git a/prism/src/explicit/FastAdaptiveUniformisationModelChecker.java b/prism/src/explicit/FastAdaptiveUniformisationModelChecker.java new file mode 100644 index 00000000..e821e0f5 --- /dev/null +++ b/prism/src/explicit/FastAdaptiveUniformisationModelChecker.java @@ -0,0 +1,275 @@ +//============================================================================== +// +// Copyright (c) 2013- +// Authors: +// * Dave Parker (University of Oxford) +// * Ernst Moritz Hahn (University of Oxford) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package explicit; + +import parser.Values; +import parser.ast.Expression; +import parser.ast.ExpressionProb; +import parser.ast.ExpressionReward; +import parser.ast.ExpressionTemporal; +import parser.ast.ModulesFile; +import parser.ast.PropertiesFile; +import parser.ast.LabelList; +import prism.Prism; +import prism.PrismException; +import prism.PrismLog; +import prism.Result; +import simulator.PrismModelExplorer; +import parser.ast.RewardStruct; + +/** + * CTMC model checker based on fast adaptive uniformisation. + */ +public class FastAdaptiveUniformisationModelChecker +{ + // Prism object + private Prism prism; + // Log + private PrismLog mainLog; + // Model file + private ModulesFile modulesFile; + // Properties file + private PropertiesFile propertiesFile; + // Constants from model + private Values constantValues; + // Labels from the model + private LabelList labelListModel; + // Labels from the property file + private LabelList labelListProp; + + /** + * Constructor. + */ + public FastAdaptiveUniformisationModelChecker(Prism prism, ModulesFile modulesFile, PropertiesFile propertiesFile) throws PrismException + { + this.prism = prism; + mainLog = prism.getMainLog(); + this.modulesFile = modulesFile; + this.propertiesFile = propertiesFile; + + // Get combined constant values from model/properties + constantValues = new Values(); + constantValues.addValues(modulesFile.getConstantValues()); + if (propertiesFile != null) + constantValues.addValues(propertiesFile.getConstantValues()); + this.labelListModel = modulesFile.getLabelList(); + this.labelListProp = propertiesFile.getLabelList(); + } + + /** + * Model check a property. + */ + public Result check(Expression expr) throws PrismException + { + Result res; + String resultString; + long timer; + + // Starting model checking + timer = System.currentTimeMillis(); + + // Do model checking + res = checkExpression(expr); + + // Model checking complete + timer = System.currentTimeMillis() - timer; + mainLog.println("\nModel checking completed in " + (timer / 1000.0) + " secs."); + + // Print result to log + resultString = "Result"; + if (!("Result".equals(expr.getResultName()))) + resultString += " (" + expr.getResultName().toLowerCase() + ")"; + resultString += ": " + res; + mainLog.print("\n" + resultString + "\n"); + + // Return result + return res; + } + + /** + * Model check an expression (used recursively). + */ + private Result checkExpression(Expression expr) throws PrismException + { + Result res; + + // Current range of supported properties is quite limited... + if (expr instanceof ExpressionProb) + res = checkExpressionProb((ExpressionProb) expr); + else if (expr instanceof ExpressionReward) + res = checkExpressionReward((ExpressionReward) expr); + else + throw new PrismException("Fast adaptive uniformisation not yet supported for this operator"); + + return res; + } + + /** + * Model check a P operator. + */ + private Result checkExpressionProb(ExpressionProb expr) throws PrismException + { + // Check whether P=? (only case allowed) + if (expr.getProb() != null) { + throw new PrismException("Fast adaptive uniformisation model checking currently only supports P=? properties"); + } + + if (!(expr.getExpression() instanceof ExpressionTemporal)) { + throw new PrismException("Fast adaptive uniformisation model checking currently only supports simple path operators"); + } + ExpressionTemporal exprTemp = (ExpressionTemporal) expr.getExpression(); + if (!exprTemp.isSimplePathFormula()) { + throw new PrismException("Fast adaptive uniformisation window model checking currently only supports simple until operators"); + } + + double timeLower = 0.0; + if (exprTemp.getLowerBound() != null) { + timeLower = exprTemp.getLowerBound().evaluateDouble(constantValues); + } + if (exprTemp.getUpperBound() == null) { + throw new PrismException("Fast adaptive uniformisation window model checking currently requires an upper time bound"); + } + double timeUpper = exprTemp.getUpperBound().evaluateDouble(constantValues); + + if (!exprTemp.hasBounds()) { + throw new PrismException("Fast adaptive uniformisation window model checking currently only supports timed properties"); + } + + mainLog.println("Starting transient probability computation using fast adaptive uniformisation..."); + PrismModelExplorer modelExplorer = new PrismModelExplorer(prism.getSimulator(), modulesFile); + FastAdaptiveUniformisation fau = new FastAdaptiveUniformisation(prism.getSettings(), modelExplorer); + fau.setConstantValues(constantValues); + + Expression op1 = exprTemp.getOperand1(); + if (op1 == null) { + op1 = Expression.True(); + } + Expression op2 = exprTemp.getOperand2(); + op1 = (Expression) op1.expandPropRefsAndLabels(propertiesFile, labelListModel); + op1 = (Expression) op1.expandPropRefsAndLabels(propertiesFile, labelListProp); + op2 = (Expression) op2.expandPropRefsAndLabels(propertiesFile, labelListModel); + op2 = (Expression) op2.expandPropRefsAndLabels(propertiesFile, labelListProp); + int operator = exprTemp.getOperator(); + + Expression sink = null; + Expression target = null; + switch (operator) { + case ExpressionTemporal.P_U: + case ExpressionTemporal.P_F: + sink = Expression.Not(op1); + break; + case ExpressionTemporal.P_G: + sink = Expression.False(); + break; + case ExpressionTemporal.P_W: + case ExpressionTemporal.P_R: + default: + throw new PrismException("operator currently not supported for fast adaptive uniformisation"); + } + + fau.setSink(sink); + fau.computeTransientProbsAdaptive(timeLower); + switch (operator) { + case ExpressionTemporal.P_U: + case ExpressionTemporal.P_F: + sink = Expression.Or(Expression.Not(op1), op2); + target = op2; + break; + case ExpressionTemporal.P_G: + sink = Expression.Not(op2); + target = op2; + break; + case ExpressionTemporal.P_W: + case ExpressionTemporal.P_R: + default: + throw new PrismException("operator currently not supported for fast adaptive uniformisation"); + } + Values varValues = new Values(); + varValues.addValue("deadlock", "true"); + sink.replaceVars(varValues); + fau.setAnalysisType(FastAdaptiveUniformisation.AnalysisType.REACH); + fau.setSink(sink); + fau.setTarget(target); + fau.computeTransientProbsAdaptive(timeUpper - timeLower); + mainLog.println("\nTotal probability lost is : " + fau.getTotalDiscreteLoss()); + mainLog.println("Maximal number of states stored during analysis : " + fau.getMaxNumStates()); + + return new Result(new Double(fau.getValue())); + } + + private RewardStruct findRewardStruct(ExpressionReward expr) throws PrismException + { + RewardStruct rewStruct = null; + Object rs = expr.getRewardStructIndex(); + if (modulesFile == null) + throw new PrismException("No model file to obtain reward structures"); + if (modulesFile.getNumRewardStructs() == 0) + throw new PrismException("Model has no rewards specified"); + if (rs == null) { + rewStruct = modulesFile.getRewardStruct(0); + } else if (rs instanceof Expression) { + int i = ((Expression) rs).evaluateInt(constantValues); + rs = new Integer(i); // for better error reporting below + rewStruct = modulesFile.getRewardStruct(i - 1); + } else if (rs instanceof String) { + rewStruct = modulesFile.getRewardStructByName((String) rs); + } + if (rewStruct == null) + throw new PrismException("Invalid reward structure index \"" + rs + "\""); + return rewStruct; + } + + /** + * Model check an R operator. + */ + private Result checkExpressionReward(ExpressionReward expr) throws PrismException + { + mainLog.println("Starting transient probability computation using fast adaptive uniformisation..."); + PrismModelExplorer modelExplorer = new PrismModelExplorer(prism.getSimulator(), modulesFile); + FastAdaptiveUniformisation fau = new FastAdaptiveUniformisation(prism.getSettings(), modelExplorer); + ExpressionTemporal temporal = (ExpressionTemporal) expr.getExpression(); + switch (temporal.getOperator()) { + case ExpressionTemporal.R_I: + fau.setAnalysisType(FastAdaptiveUniformisation.AnalysisType.REW_INST); + break; + case ExpressionTemporal.R_C: + fau.setAnalysisType(FastAdaptiveUniformisation.AnalysisType.REW_CUMUL); + break; + default: + throw new PrismException("Currently only instantaneous or cumulative rewards are allowed."); + } + double time = temporal.getUpperBound().evaluateDouble(constantValues); + RewardStruct rewStruct = findRewardStruct(expr); + fau.setRewardStruct(rewStruct); + fau.setConstantValues(constantValues); + fau.computeTransientProbsAdaptive(time); + mainLog.println("\nTotal probability lost is : " + fau.getTotalDiscreteLoss()); + mainLog.println("Maximal number of states stored during analysis : " + fau.getMaxNumStates()); + return new Result(new Double(fau.getValue())); + } +} diff --git a/prism/src/explicit/ModelExplorer.java b/prism/src/explicit/ModelExplorer.java new file mode 100644 index 00000000..28ae099e --- /dev/null +++ b/prism/src/explicit/ModelExplorer.java @@ -0,0 +1,99 @@ +//============================================================================== +// +// Authors: +// * Dave Parker (University of Oxford) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package explicit; + +import parser.State; +import prism.PrismException; + +public interface ModelExplorer +{ + public State getDefaultInitialState() throws PrismException; + + public void queryState(State state) throws PrismException; + + public void queryState(State state, double time) throws PrismException; + + /** + * Returns the current number of available choices. + * @throws PrismException + */ + public int getNumChoices() throws PrismException; + + /** + * Returns the current (total) number of available transitions. + * @throws PrismException + */ + public int getNumTransitions() throws PrismException; + + /** + * Returns the current number of available transitions in choice i. + * @throws PrismException + */ + public int getNumTransitions(int i) throws PrismException; + + /** + * Get the probability/rate of a transition within a choice, specified by its index/offset. + */ + public double getTransitionProbability(int i, int offset) throws PrismException; + + /** + * Get the action label of a transition as a string, specified by its index/offset. + * (null for asynchronous/independent transitions) + * (see also {@link #getTransitionModuleOrAction(int, int)} and {@link #getTransitionModuleOrActionIndex(int, int)}) + * TODO: change return type to Object? + * @throws PrismException + */ + public String getTransitionAction(int i, int offset) throws PrismException; + + /** + * Get the action label of a transition as a string, specified by its index. + * (null for asynchronous/independent transitions) + * (see also {@link #getTransitionModuleOrAction(int)} and {@link #getTransitionModuleOrActionIndex(int)}) + * TODO: change return type to Object? + * @throws PrismException + */ + public String getTransitionAction(int index) throws PrismException; + + /** + * Get the probability/rate of a transition, specified by its index. + */ + public double getTransitionProbability(int i) throws PrismException; + + /** + * Get the sum of probabilities/rates for transitions. + */ +// public double getTransitionProbabilitySum() throws PrismException; + + /** + * Get the target (as a new State object) of a transition within a choice, specified by its index/offset. + */ + public State computeTransitionTarget(int i, int offset) throws PrismException; + + /** + * Get the target of a transition (as a new State object), specified by its index. + */ + public State computeTransitionTarget(int i) throws PrismException; +} diff --git a/prism/src/prism/Prism.java b/prism/src/prism/Prism.java index 4bb82716..add346a9 100644 --- a/prism/src/prism/Prism.java +++ b/prism/src/prism/Prism.java @@ -2395,8 +2395,14 @@ public class Prism implements PrismSettingsListener return modelCheckPTA(propertiesFile, prop.getExpression(), definedPFConstants); } + // For fast adaptive uniformisation + if (currentModelType == ModelType.CTMC && settings.getString(PrismSettings.PRISM_TRANSIENT_METHOD).equals("Fast adaptive uniformisation")) { + FastAdaptiveUniformisationModelChecker fauMC; + fauMC = new FastAdaptiveUniformisationModelChecker(this, currentModulesFile, propertiesFile); + return fauMC.check(prop.getExpression()); + } // Auto-switch engine if required - if (currentModelType == ModelType.MDP && !Expression.containsMultiObjective(prop.getExpression())) { + else if (currentModelType == ModelType.MDP && !Expression.containsMultiObjective(prop.getExpression())) { if (getMDPSolnMethod() != Prism.MDP_VALITER && !getExplicit()) { mainLog.printWarning("Switching to explicit engine to allow use of chosen MDP solution method."); engineSwitch = true; @@ -2782,7 +2788,16 @@ public class Prism implements PrismSettingsListener l = System.currentTimeMillis(); - if (!getExplicit()) { + if (currentModelType == ModelType.CTMC && settings.getString(PrismSettings.PRISM_TRANSIENT_METHOD).equals("Fast adaptive uniformisation")) { + PrismModelExplorer modelExplorer = new PrismModelExplorer(getSimulator(), currentModulesFile); + FastAdaptiveUniformisation fau = new FastAdaptiveUniformisation(settings, modelExplorer); + fau.setConstantValues(currentModulesFile.getConstantValues()); + mainLog.println("Starting transient probability computation using fast adaptive uniformisation..."); + probsExpl = fau.doTransient(time, fileIn); + mainLog.println("\nTotal probability lost is : " + fau.getTotalDiscreteLoss()); + mainLog.println("Maximal number of states stored during analysis : " + fau.getMaxNumStates()); + } + else if (!getExplicit()) { if (currentModelType == ModelType.DTMC) { mc = new ProbModelChecker(this, currentModel, null); probs = ((ProbModelChecker) mc).doTransient((int) time, fileIn); diff --git a/prism/src/prism/PrismSettings.java b/prism/src/prism/PrismSettings.java index 6ec11223..475a3ec7 100644 --- a/prism/src/prism/PrismSettings.java +++ b/prism/src/prism/PrismSettings.java @@ -97,6 +97,7 @@ public class PrismSettings implements Observer public static final String PRISM_SCC_METHOD = "prism.sccMethod"; public static final String PRISM_SYMM_RED_PARAMS = "prism.symmRedParams"; public static final String PRISM_PTA_METHOD = "prism.ptaMethod"; + public static final String PRISM_TRANSIENT_METHOD = "prism.transientMethod"; public static final String PRISM_AR_OPTIONS = "prism.arOptions"; public static final String PRISM_EXPORT_ADV = "prism.exportAdv"; public static final String PRISM_EXPORT_ADV_FILENAME = "prism.exportAdvFilename"; @@ -156,6 +157,11 @@ public class PrismSettings implements Observer public static final String LOG_BG_COLOUR = "log.bgColour"; public static final String LOG_BUFFER_LENGTH = "log.bufferLength"; + //FAU + public static final String PRISM_FAU_DELTA = "prism.faudelta"; + public static final String PRISM_FAU_INTERVALS = "prism.fauintervals"; + public static final String PRISM_FAU_INITIVAL = "prism.fauinitival"; + public static final String PRISM_FAU_ARRAYTHRESHOLD = "prism.fauarraythreshold"; //Defaults, types and constaints @@ -165,7 +171,8 @@ public class PrismSettings implements Observer "Simulator", "Model", "Properties", - "Log" + "Log", + "FAU" }; public static final int[] propertyOwnerIDs = { @@ -173,7 +180,8 @@ public class PrismSettings implements Observer PropertyConstants.SIMULATOR, PropertyConstants.MODEL, PropertyConstants.PROPERTIES, - PropertyConstants.LOG + PropertyConstants.LOG, + PropertyConstants.FAU }; @@ -198,6 +206,8 @@ public class PrismSettings implements Observer "Which engine (hybrid, sparse, MTBDD, explicit) should be used for model checking." }, { CHOICE_TYPE, PRISM_PTA_METHOD, "PTA model checking method", "3.3", "Stochastic games", "Digital clocks,Stochastic games", "Which method to use for model checking of PTAs." }, + { CHOICE_TYPE, PRISM_TRANSIENT_METHOD, "Transient probability computation method", "3.3", "Uniformisation", "Uniformisation,Fast adaptive uniformisation", + "Which method to use for computing transient probabilities in CTMCs." }, // NUMERICAL SOLUTION OPTIONS: { CHOICE_TYPE, PRISM_LIN_EQ_METHOD, "Linear equations method", "2.1", "Jacobi", "Power,Jacobi,Gauss-Seidel,Backwards Gauss-Seidel,Pseudo-Gauss-Seidel,Backwards Pseudo-Gauss-Seidel,JOR,SOR,Backwards SOR,Pseudo-SOR,Backwards Pseudo-SOR", "Which iterative method to use when solving linear equation systems." }, @@ -327,6 +337,12 @@ public class PrismSettings implements Observer { FONT_COLOUR_TYPE, LOG_FONT, "Display font", "2.1", new FontColorPair(new Font("monospaced", Font.PLAIN, 12), Color.black), "", "Font used for the log display." }, { COLOUR_TYPE, LOG_BG_COLOUR, "Background colour", "2.1", new Color(255,255,255), "", "Background colour for the log display." }, { INTEGER_TYPE, LOG_BUFFER_LENGTH, "Buffer length", "2.1", new Integer(10000), "1,", "Length of the buffer for the log display." } + }, + { + { DOUBLE_TYPE, PRISM_FAU_DELTA, "Cut off delta", "4.0.1", new Double(10E-12), "", "States which get a probability below this number during the fast adaptive analysis will be removed." }, + { INTEGER_TYPE, PRISM_FAU_ARRAYTHRESHOLD, "Threshold to swap to array mode", "4.0.1", new Integer(100), "", "If this number of iterations happened during fast adaptive uniformisation without changes to the state space, assume that further changes are unlikely." }, + { INTEGER_TYPE, PRISM_FAU_INTERVALS, "Number of time intervals", "4.0.1", new Integer(1), "", "Splits the time of a time-bounded property into the specified number of intervals." }, + { DOUBLE_TYPE, PRISM_FAU_INITIVAL, "Length of initial time interval", "4.0.1", new Double(1.0), "", "Length of initial time interval in addition to regular time intervals." }, } }; @@ -781,6 +797,20 @@ public class PrismSettings implements Observer throw new PrismException("No parameter specified for -" + sw + " switch"); } } + // Transient methods + else if (sw.equals("transientmethod")) { + if (i < args.length - 1) { + s = args[++i]; + if (s.equals("unif")) + set(PRISM_TRANSIENT_METHOD, "Uniformisation"); + else if (s.equals("fau")) + set(PRISM_TRANSIENT_METHOD, "Fast adaptive uniformisation"); + else + throw new PrismException("Unrecognised option for -" + sw + " switch (options are: unif, fau)"); + } else { + throw new PrismException("No parameter specified for -" + sw + " switch"); + } + } // NUMERICAL SOLUTION OPTIONS: @@ -1114,6 +1144,68 @@ public class PrismSettings implements Observer } } + // Fast Adaptive Uniformisation + + // Delta for fast adaptive uniformisation + else if (sw.equals("faudelta")) { + if (i < args.length - 1) { + try { + d = Double.parseDouble(args[++i]); + if (d < 0) + throw new NumberFormatException(""); + set(PRISM_FAU_DELTA, d); + } catch (NumberFormatException e) { + throw new PrismException("Invalid value for -" + sw + " switch"); + } + } else { + throw new PrismException("No value specified for -" + sw + " switch"); + } + } + else if (sw.equals("fauarraythreshold")) { + if (i < args.length - 1) { + try { + j = Integer.parseInt(args[++i]); + if (j < 0) + throw new NumberFormatException(""); + set(PRISM_FAU_ARRAYTHRESHOLD, j); + } catch (NumberFormatException e) { + throw new PrismException("Invalid value for -" + sw + " switch"); + } + } else { + throw new PrismException("No value specified for -" + sw + " switch"); + } + } + + // number of intervals for fast adaptive uniformisation + else if (sw.equals("fauintervals")) { + if (i < args.length - 1) { + try { + j = Integer.parseInt(args[++i]); + if (j < 0) + throw new NumberFormatException(""); + set(PRISM_FAU_INTERVALS, j); + } catch (NumberFormatException e) { + throw new PrismException("Invalid value for -" + sw + " switch"); + } + } else { + throw new PrismException("No value specified for -" + sw + " switch"); + } + } + else if (sw.equals("fauinitival")) { + if (i < args.length - 1) { + try { + d = Double.parseDouble(args[++i]); + if (d < 0.0) + throw new NumberFormatException(""); + set(PRISM_FAU_INITIVAL, d); + } catch (NumberFormatException e) { + throw new PrismException("Invalid value for -" + sw + " switch"); + } + } else { + throw new PrismException("No value specified for -" + sw + " switch"); + } + } + // unknown switch - error else { throw new PrismException("Invalid switch -" + sw + " (type \"prism -help\" for full list)"); @@ -1135,6 +1227,7 @@ public class PrismSettings implements Observer mainLog.println("-hybrid (or -h) ................ Use the Hybrid engine [default]"); mainLog.println("-explicit (or -ex) ............. Use the explicit engine"); mainLog.println("-ptamethod .............. Specify PTA engine (games, digital) [default: games]"); + mainLog.println("-transientmethod ........ CTMC transient analysis methof (unif, fau) [default: unif]"); mainLog.println(); mainLog.println("SOLUTION METHODS (LINEAR EQUATIONS):"); mainLog.println("-power (or -pow, -pwr) ......... Use the Power method for numerical computation"); @@ -1198,6 +1291,12 @@ public class PrismSettings implements Observer mainLog.println("-gsmax (or sormax ) ..... Set memory limit (KB) for hybrid GS/SOR [default: 1024]"); mainLog.println("-cuddmaxmem ................ Set max memory for CUDD package (KB) [default: 200x1024]"); mainLog.println("-cuddepsilon ............... Set epsilon value for CUDD package [default: 1e-15]"); + mainLog.println(); + mainLog.println("FAST ADAPTIVE UNIFORMISATION OPTIONS:"); + mainLog.println("-faudelta .................. Set probability threshold for irrelevant states [default: 1e-12]"); + mainLog.println("-fauarraythreshold ......... Set threshold when to switch to sparse matrix [default: 100]"); + mainLog.println("-fauintervals .............. Set number of intervals to divide time intervals to [default: 1]"); + mainLog.println("-fauinitival ............... Set length of additional initial time interval [default: 1.0]"); } /** diff --git a/prism/src/prism/PropertyConstants.java b/prism/src/prism/PropertyConstants.java index 5aa3f054..27f2643d 100644 --- a/prism/src/prism/PropertyConstants.java +++ b/prism/src/prism/PropertyConstants.java @@ -43,8 +43,9 @@ public final class PropertyConstants public static final int PROPERTIES = 11; public static final int SIMULATOR = 12; public static final int LOG = 13; - public static final int SSHHOST = 14; - public static final int FILESYSTEM = 15; - public static final int NETWORK_PROFILE = 16; - public static final int GRAPH_DISPLAY = 17; + public static final int FAU = 14; + public static final int SSHHOST = 15; + public static final int FILESYSTEM = 16; + public static final int NETWORK_PROFILE = 17; + public static final int GRAPH_DISPLAY = 18; } diff --git a/prism/src/simulator/PrismModelExplorer.java b/prism/src/simulator/PrismModelExplorer.java new file mode 100644 index 00000000..4e1ba1c3 --- /dev/null +++ b/prism/src/simulator/PrismModelExplorer.java @@ -0,0 +1,125 @@ +//============================================================================== +// +// Authors: +// * Dave Parker (University of Oxford) +// +//------------------------------------------------------------------------------ +// +// This file is part of PRISM. +// +// PRISM is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// PRISM is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with PRISM; if not, write to the Free Software Foundation, +// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +//============================================================================== + +package simulator; + +import parser.State; +import parser.ast.ModulesFile; +import prism.PrismException; +import explicit.ModelExplorer; + +public class PrismModelExplorer implements ModelExplorer +{ + private SimulatorEngine simEngine; + private ModulesFile modulesFile; + + public PrismModelExplorer(SimulatorEngine simEngine, ModulesFile modulesFile) throws PrismException + { + this.simEngine = simEngine; + this.modulesFile = modulesFile; + simEngine.createNewOnTheFlyPath(modulesFile); + + } + + @Override + public State getDefaultInitialState() throws PrismException + { + return modulesFile.getDefaultInitialState(); + } + + @Override + public void queryState(State state) throws PrismException + { + simEngine.initialisePath(state); + } + + @Override + public void queryState(State state, double time) throws PrismException + { + queryState(state); + } + + @Override + public int getNumChoices() throws PrismException + { + return simEngine.getNumChoices(); + } + + @Override + public int getNumTransitions() throws PrismException + { + return simEngine.getNumTransitions(); + } + + @Override + public int getNumTransitions(int i) throws PrismException + { + return simEngine.getNumTransitions(i); + } + + @Override + public String getTransitionAction(int i, int offset) throws PrismException + { + return simEngine.getTransitionAction(i, offset); + } + + @Override + public String getTransitionAction(int i) throws PrismException + { + return simEngine.getTransitionAction(i); + } + + @Override + public double getTransitionProbability(int i, int offset) throws PrismException + { + return simEngine.getTransitionProbability(i, offset); + } + + @Override + public double getTransitionProbability(int i) throws PrismException + { + return simEngine.getTransitionProbability(i); + } + + /* + @Override + public double getTransitionProbabilitySum() throws PrismException + { + return simEngine.getTransitionProbabilitySum(); + } + */ + + @Override + public State computeTransitionTarget(int i, int offset) throws PrismException + { + return simEngine.computeTransitionTarget(i, offset); + } + + @Override + public State computeTransitionTarget(int i) throws PrismException + { + return simEngine.computeTransitionTarget(i); + } +}