(interval iteration, explicit) actually perform interval iteration using the explicit engine

Supported:
DTMC reachability probability and expected reward computations.
MDP Pmax, Pmin, Rmax, Rmin for reachability.


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

prism/src/explicit/DTMCModelChecker.java

@@ -29,28 +29,38 @@ package explicit;
 import java.io.File;
 import java.util.Arrays;
 import java.util.BitSet;
+import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
+import java.util.Map.Entry;
+import java.util.PrimitiveIterator;
 import java.util.Vector;
 
 import parser.VarList;
 import parser.ast.Declaration;
 import parser.ast.DeclarationIntUnbounded;
 import parser.ast.Expression;
+import prism.OptionsIntervalIteration;
 import prism.Prism;
 import prism.PrismComponent;
 import prism.PrismException;
 import prism.PrismFileLog;
+import prism.PrismNotSupportedException;
 import prism.PrismSettings;
 import prism.PrismUtils;
 import acceptance.AcceptanceReach;
 import acceptance.AcceptanceType;
 import automata.DA;
 import common.IntSet;
+import common.PeriodicTimer;
+import common.StopWatch;
 import common.IterableBitSet;
 import common.StopWatch;
 import explicit.LTLModelChecker.LTLProduct;
+import explicit.modelviews.DTMCAlteredDistributions;
+import explicit.modelviews.MDPFromDTMC;
 import explicit.rewards.MCRewards;
+import explicit.rewards.MDPRewards;
 import explicit.rewards.Rewards;
 
 /**
@@ -605,6 +615,10 @@ public class DTMCModelChecker extends ProbModelChecker
 			mainLog.printWarning("Switching to linear equation solution method \"" + linEqMethod.fullName() + "\"");
 		}
 
+		if (doIntervalIteration && (!precomp || !prob0 || !prob1)) {
+			throw new PrismNotSupportedException("Interval iteration requires precomputations to be active");
+		}
+
 		// Start probabilistic reachability
 		timer = System.currentTimeMillis();
 		mainLog.println("\nStarting probabilistic reachability...");
@@ -693,7 +707,11 @@ public class DTMCModelChecker extends ProbModelChecker
 			throw new PrismException("Unknown linear equation solution method " + linEqMethod.fullName());
 		}
 
-		res = doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration());
+		if (doIntervalIteration) {
+			res = doIntervalIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration());
+		} else {
+			res = doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, getDoTopologicalValueIteration());
+		}
 
 		// Finished probabilistic reachability
 		timer = System.currentTimeMillis() - timer;
@@ -1087,6 +1105,109 @@ public class DTMCModelChecker extends ProbModelChecker
 		return doValueIterationReachProbs(dtmc, no, yes, init, known, iterationMethod, false);
 	}
 
+	/**
+	 * Compute reachability probabilities using power method (interval variant).
+	 * @param dtmc The DTMC
+	 * @param no Probability 0 states
+	 * @param yes Probability 1 states
+	 * @param init Optionally, an initial solution vector (will be overwritten), will be ignored if known == null
+	 * @param known Optionally, a set of states for which the exact answer is known
+	 * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values.
+	 * @param topological do topological interval iteration?
+	 */
+	protected ModelCheckerResult doIntervalIterationReachProbs(DTMC dtmc, BitSet no, BitSet yes, double init[], BitSet known, IterationMethod iterationMethod, boolean topological) throws PrismException
+	{
+		BitSet unknown;
+		int i, n;
+		double initBelow[], initAbove[];
+		long timer;
+
+		// Start value iteration
+		timer = System.currentTimeMillis();
+		String description = "with " + iterationMethod.getDescriptionShort();
+		mainLog.println("Starting interval iteration (" + description + ")...");
+
+		ExportIterations iterationsExport = null;
+		if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) {
+			iterationsExport = new ExportIterations("Explicit DTMC ReachProbs interval iteration  (" + description + ")");
+		}
+
+		// Store num states
+		n = dtmc.getNumStates();
+
+		// Create solution vector(s)
+		initBelow = (init == null) ? new double[n] : init;
+		initAbove = new double[n];
+
+		// Initialise solution vectors. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 1.0/0.0 if in yes/no; (3) initVal
+		// where initVal is 0.0 or 1.0, depending on whether we converge from below/above.
+		if (known != null && init != null) {
+			for (i = 0; i < n; i++) {
+				initBelow[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0;
+				initAbove[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0;
+			}
+		} else {
+			for (i = 0; i < n; i++) {
+				initBelow[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 :  0.0;
+				initAbove[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 :  1.0;
+			}
+		}
+
+		// Determine set of states actually need to compute values for
+		unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(yes);
+		unknown.andNot(no);
+		if (known != null)
+			unknown.andNot(known);
+
+		if (iterationsExport != null) {
+			iterationsExport.exportVector(initBelow, 0);
+			iterationsExport.exportVector(initAbove, 1);
+		}
+
+		IntSet unknownStates = IntSet.asIntSet(unknown);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow();
+		final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove();
+		final boolean checkMonotonic = iiOptions.isCheckMonotonicity();
+
+		if (!enforceMonotonicFromAbove) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above.");
+		}
+		if (!enforceMonotonicFromBelow) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below.");
+		}
+
+		IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultInterval(dtmc, true, enforceMonotonicFromBelow, checkMonotonic);
+		IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultInterval(dtmc, false, enforceMonotonicFromAbove, checkMonotonic);
+
+		below.init(initBelow);
+		above.init(initAbove);
+
+		if (topological) {
+			// Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering
+			// states in unknown
+			SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get);
+
+			IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> {
+				soln[s] = dtmc.mvMultJacSingle(s, soln);
+			};
+
+			// run the actual value iteration
+			return iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport);
+		} else {
+			// run the actual value iteration
+			return iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport);
+		}
+
+	}
+
+
+
 	/**
 	 * Compute bounded reachability probabilities.
 	 * i.e. compute the probability of reaching a state in {@code target} within k steps.
@@ -1199,6 +1320,381 @@ public class DTMCModelChecker extends ProbModelChecker
 		return res;
 	}
 
+	/**
+	 * Compute upper bound for maximum expected reward, using the variant specified in the settings.
+	 * @param dtmc the model
+	 * @param mcRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinity states
+	 * @return upper bound on R=?[ F target ] for all states
+	 */
+	double computeReachRewardsUpperBound(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		// inf and target states become trap states (with self-loops)
+		BitSet trapStates = (BitSet) target.clone();
+		trapStates.or(inf);
+		DTMCAlteredDistributions cleanedDTMC = DTMCAlteredDistributions.addSelfLoops(dtmc, trapStates);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		double upperBound = 0.0;
+		String method = null;
+		switch (iiOptions.getBoundMethod()) {
+		case VARIANT_1_COARSE:
+			upperBound = computeReachRewardsUpperBoundVariant1Coarse(cleanedDTMC, mcRewards, target, unknown, inf);
+			method = "variant 1, coarse";
+			break;
+		case VARIANT_1_FINE:
+			upperBound = computeReachRewardsUpperBoundVariant1Fine(cleanedDTMC, mcRewards, target, unknown, inf);
+			method = "variant 1, fine";
+			break;
+		case VARIANT_2:
+			upperBound = computeReachRewardsUpperBoundVariant2(cleanedDTMC, mcRewards, target, unknown, inf);
+			method = "variant 2";
+			break;
+		case DEFAULT:
+		case DSMPI:
+		{
+			MDP mdp = new MDPFromDTMC(cleanedDTMC);
+			MDPRewards mdpRewards = new MDPRewards() {
+
+				@Override
+				public double getStateReward(int s)
+				{
+					return mcRewards.getStateReward(s);
+				}
+
+				@Override
+				public double getTransitionReward(int s, int i)
+				{
+					return 0;
+				}
+
+				@Override
+				public MDPRewards liftFromModel(Product<? extends Model> product)
+				{
+					throw new RuntimeException("Unsupported");
+				}
+
+				@Override
+				public boolean hasTransitionRewards()
+				{
+					return false;
+				}
+			};
+			upperBound = DijkstraSweepMPI.computeUpperBound(this, mdp, mdpRewards, target, unknown);
+			method = "Dijkstra Sweep MPI";
+			break;
+		}
+		}
+
+		if (method == null) {
+			throw new PrismException("Unknown upper bound heuristic");
+		}
+
+		mainLog.println("Upper bound for expectation (" + method + "): " + upperBound);
+
+		return upperBound;
+	}
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 1, coarse),
+	 * i.e., does not compute separate q_t / p_t per SCC.
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param dtmc the model
+	 * @param mcRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinity states
+	 * @return upper bound on R=?[ F target ] for all states
+	 */
+	double computeReachRewardsUpperBoundVariant1Coarse(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()];
+		int[] Ct = new int[dtmc.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, null);
+		BitSet trivial = new BitSet();
+
+		double q = 0;
+		for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) {
+			IntSet statesForSCC = sccs.getStatesForSCC(scc);
+
+			int cardinality = statesForSCC.cardinality();
+
+			PrimitiveIterator.OfInt itSCC = statesForSCC.iterator();
+			while (itSCC.hasNext()) {
+				int s = itSCC.nextInt();
+				Ct[s] = cardinality;
+
+				if (target.get(s) || inf.get(s)) {
+					// trap states
+					assert(cardinality == 1);
+					break;  // continue with next SCC
+				}
+
+				double probRemain = 0;
+				boolean allRemain = true;  // all successors remain in the SCC?
+				boolean hasSelfloop = false;
+				for (Iterator<Entry<Integer, Double>> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) {
+					Entry<Integer, Double> t = it.next();
+					if (statesForSCC.get(t.getKey())) {
+						probRemain += t.getValue();
+						hasSelfloop = true;
+					} else {
+						allRemain = false;
+					}
+				}
+
+				if (!allRemain) { // action in the set X
+					q = Math.max(q, probRemain);
+				}
+
+				if (cardinality == 1 && !hasSelfloop) {
+					trivial.set(s);
+				}
+			}
+		}
+
+		double p = 1;
+		for (int s = 0; s < dtmc.getNumStates(); s++) {
+			for (Iterator<Entry<Integer, Double>> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) {
+				Entry<Integer, Double> t = it.next();
+				p = Math.min(p, t.getValue());
+			}
+		}
+
+		double upperBound = 0;
+		for (int s = 0; s < dtmc.getNumStates(); s++) {
+			if (target.get(s) || inf.get(s)) {
+				// inf or target states: not relevant, set visits to 0, ignore in summation
+				boundsOnExpectedVisits[s] = 0.0;
+			} else if (unknown.get(s)) {
+				if (trivial.get(s)) {
+					// s is a trivial SCC: seen at most once
+					boundsOnExpectedVisits[s] = 1.0;
+				} else {
+					boundsOnExpectedVisits[s] = 1 / (Math.pow(p, Ct[s]-1) * (1.0-q));
+				}
+
+				upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s);
+			}
+		}
+
+		timer.stop();
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 1, coarse):");
+			mainLog.println("p = " + p);
+			mainLog.println("q = " + q);
+			mainLog.println("|Ct| = " + Arrays.toString(Ct));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+		}
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 1, fine).
+	 * i.e., does compute separate q_t / p_t per SCC.
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param dtmc the model
+	 * @param mcRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinity states
+	 * @return upper bound on R=?[ F target ] for all states
+	 */
+	double computeReachRewardsUpperBoundVariant1Fine(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()];
+		double[] qt = new double[dtmc.getNumStates()];
+		double[] pt = new double[dtmc.getNumStates()];
+		int[] Ct = new int[dtmc.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, null);
+		BitSet trivial = new BitSet();
+
+		for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) {
+			IntSet statesForSCC = sccs.getStatesForSCC(scc);
+
+			double q = 0;
+			double p = 1;
+			int cardinality = statesForSCC.cardinality();
+
+			PrimitiveIterator.OfInt itSCC = statesForSCC.iterator();
+			while (itSCC.hasNext()) {
+				int s = itSCC.nextInt();
+				Ct[s] = cardinality;
+
+				double probRemain = 0;
+				boolean allRemain = true;  // all successors remain in the SCC?
+				boolean hasSelfloop = false;
+				for (Iterator<Entry<Integer, Double>> it = dtmc.getTransitionsIterator(s); it.hasNext(); ) {
+					Entry<Integer, Double> t = it.next();
+					if (statesForSCC.get(t.getKey())) {
+						probRemain += t.getValue();
+						p = Math.min(p, t.getValue());
+						hasSelfloop = true;
+					} else {
+						// outgoing edge
+						allRemain = false;
+					}
+				}
+
+				if (!allRemain) { // action in the set Xt
+					q = Math.max(q, probRemain);
+				}
+
+				if (cardinality == 1 && !hasSelfloop) {
+					trivial.set(s);
+				}
+			}
+
+			for (int s : statesForSCC) {
+				qt[s] = q;
+				pt[s] = p;
+			}
+		}
+
+		double upperBound = 0;
+		for (int s = 0; s < dtmc.getNumStates(); s++) {
+			if (target.get(s) || inf.get(s)) {
+				// inf or target states: not relevant, set visits to 0, ignore in summation
+				boundsOnExpectedVisits[s] = 0.0;
+			} else if (unknown.get(s)) {
+				if (trivial.get(s)) {
+					// s is a trivial SCC: seen at most once
+					boundsOnExpectedVisits[s] = 1.0;
+				} else {
+					if (pt[s] == 1.0) {
+						//throw new PrismException("Upper bound computation had p_t = 1 for state " + s);
+					}
+					boundsOnExpectedVisits[s] = 1 / (Math.pow(pt[s], Ct[s]-1) * (1.0-qt[s]));
+				}
+
+				upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s);
+			} else {
+				throw new PrismException("Bogus arguments: inf/target/unknown should partition state space");
+			}
+		}
+
+		timer.stop();
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 1, fine):");
+			mainLog.println("pt = " + Arrays.toString(pt));
+			mainLog.println("qt = " + Arrays.toString(qt));
+			mainLog.println("|Ct| = " + Arrays.toString(Ct));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+		}
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 2).
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param dtmc the model
+	 * @param mcRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinity states
+	 * @return upper bound on R=?[ F target ] for all states
+	 */
+	double computeReachRewardsUpperBoundVariant2(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] dt = new double[dtmc.getNumStates()];
+		double[] boundsOnExpectedVisits = new double[dtmc.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get);
+
+		BitSet T = (BitSet) target.clone();
+
+		@SuppressWarnings("unused")
+		int i = 0;
+		while (true) {
+			BitSet Si = new BitSet();
+			i++;
+
+			// TODO: might be inefficient, worst-case quadratic runtime...
+			for (PrimitiveIterator.OfInt it = IterableBitSet.getClearBits(T, dtmc.getNumStates() -1 ).iterator(); it.hasNext(); ) {
+				int s = it.nextInt();
+//				mainLog.println("Check " + s + " against " + T);
+				if (dtmc.someSuccessorsInSet(s, T)) {
+					Si.set(s);
+				}
+			}
+
+			if (Si.isEmpty()) {
+				break;
+			}
+
+			// mainLog.println("S" + i + " = " + Si);
+			// mainLog.println("T = " + T);
+
+			for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(Si).iterator(); it.hasNext(); ) {
+				final int t = it.nextInt();
+				final int sccIndexForT = sccs.getSCCIndex(t);
+				double d = dtmc.sumOverTransitions(t, (int __, int u, double prob) -> {
+					// mainLog.println("t = " + t + ", u = " + u + ", prob = " + prob);
+					if (!T.get(u))
+						return 0.0;
+
+					boolean inSameSCC = (sccs.getSCCIndex(u) == sccIndexForT);
+					double d_u_t = inSameSCC ? dt[u] : 1.0;
+					// mainLog.println("d_u_t = " + d_u_t);
+					return d_u_t * prob;
+				});
+				dt[t] = d;
+				// mainLog.println("d["+t+"] = " + d);
+			}
+
+			T.or(Si);
+		}
+
+		double upperBound = 0;
+		for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(unknown).iterator(); it.hasNext();) {
+			int s = it.nextInt();
+			boundsOnExpectedVisits[s] = 1 / dt[s];
+			upperBound += boundsOnExpectedVisits[s] * mcRewards.getStateReward(s);
+		}
+
+		timer.stop();
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 2):");
+			mainLog.println("d_t = " + Arrays.toString(dt));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+		}
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
+
 	/**
 	 * Compute expected reachability rewards.
 	 * @param dtmc The DTMC
@@ -1302,7 +1798,11 @@ public class DTMCModelChecker extends ProbModelChecker
 			throw new PrismException("Unknown linear equation solution method " + linEqMethod.fullName());
 		}
 
-		res = doValueIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration());
+		if (doIntervalIteration) {
+			res = doIntervalIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration());
+		} else {
+			res = doValueIterationReachRewards(dtmc, mcRewards, target, inf, init, known, iterationMethod, getDoTopologicalValueIteration());
+		}
 
 		// Finished expected reachability
 		timer = System.currentTimeMillis() - timer;
@@ -1498,6 +1998,136 @@ public class DTMCModelChecker extends ProbModelChecker
 		}
 	}
 
+	/**
+	 * Compute expected reachability rewards using interval iteration.
+	 * @param dtmc The DTMC
+	 * @param mcRewards The rewards
+	 * @param target Target states
+	 * @param inf States for which reward is infinite
+	 * @param init Optionally, an initial solution vector (will be overwritten)
+	 * @param known Optionally, a set of states for which the exact answer is known
+	 * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values.
+	 * @param topological do topological interval iteration?
+	 */
+	protected ModelCheckerResult doIntervalIterationReachRewards(DTMC dtmc, MCRewards mcRewards, BitSet target, BitSet inf, double init[], BitSet known, IterationMethod iterationMethod, boolean topological)
+			throws PrismException
+	{
+		BitSet unknown;
+		int i, n;
+		double init_below[], init_above[];
+		long timer;
+
+		// Store num states
+		n = dtmc.getNumStates();
+
+		// Determine set of states actually need to compute values for
+		unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(target);
+		unknown.andNot(inf);
+		if (known != null)
+			unknown.andNot(known);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		double upperBound;
+		if (iiOptions.hasManualUpperBound()) {
+			upperBound = iiOptions.getManualUpperBound();
+			getLog().printWarning("Upper bound for interval iteration manually set to " + upperBound);
+		} else {
+			upperBound = computeReachRewardsUpperBound(dtmc, mcRewards, target, unknown, inf);
+		}
+
+		double lowerBound;
+		if (iiOptions.hasManualLowerBound()) {
+			lowerBound = iiOptions.getManualLowerBound();
+			getLog().printWarning("Lower bound for interval iteration manually set to " + lowerBound);
+		} else {
+			lowerBound = 0.0;
+		}
+
+		// Start value iteration
+		timer = System.currentTimeMillis();
+		String description = (topological ? "topological, " : "" ) + "with " + iterationMethod.getDescriptionShort();
+		mainLog.println("Starting interval iteration (" + description + ") ...");
+
+		ExportIterations iterationsExport = null;
+		if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) {
+			iterationsExport = new ExportIterations("Explicit DTMC ReachRewards interval iteration (" + description + ") ...");
+		}
+
+		// Create solution vector(s)
+		init_below = (init == null) ? new double[n] : init;
+		init_above = new double[n];
+
+		// Initialise solution vector from below. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) lowerBound
+		if (init != null && known != null) {
+			for (i = 0; i < n; i++)
+				init_below[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound;
+		} else {
+			for (i = 0; i < n; i++)
+				init_below[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound;
+		}
+
+		// Initialise solution vector from above. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) upperBound
+		if (init != null && known != null) {
+			for (i = 0; i < n; i++)
+				init_above[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound;
+		} else {
+			for (i = 0; i < n; i++)
+				init_above[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound;
+		}
+
+
+		if (iterationsExport != null) {
+			iterationsExport.exportVector(init_below, 0);
+			iterationsExport.exportVector(init_above, 1);
+		}
+
+		IntSet unknownStates = IntSet.asIntSet(unknown);
+
+		final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow();
+		final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove();
+		final boolean checkMonotonic = iiOptions.isCheckMonotonicity();
+
+		if (!enforceMonotonicFromAbove) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above.");
+		}
+		if (!enforceMonotonicFromBelow) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below.");
+		}
+
+		IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultRewInterval(dtmc, mcRewards, true, enforceMonotonicFromBelow, checkMonotonic);
+		IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultRewInterval(dtmc, mcRewards, false, enforceMonotonicFromAbove, checkMonotonic);
+
+		below.init(init_below);
+		above.init(init_above);
+
+		ModelCheckerResult rv;
+		if (topological) {
+			// Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering
+			// states in unknown
+			SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, dtmc, true, unknown::get);
+
+			IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> {
+				soln[s] = dtmc.mvMultRewJacSingle(s, soln, mcRewards);
+			};
+
+			// run the actual value iteration
+			rv = iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport);
+		} else {
+			// run the actual value iteration
+			rv = iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport);
+		}
+
+		double max_v = PrismUtils.findMaxFinite(rv.soln, unknownStates.iterator());
+		mainLog.println("Maximum finite value in solution vector at end of interval iteration: " + max_v);
+
+		return rv;
+	}
+
 	/**
 	 * Compute (forwards) steady-state probabilities
 	 * i.e. compute the long-run probability of being in each state,

prism/src/explicit/DijkstraSweepMPI.java

@@ -0,0 +1,312 @@
+//==============================================================================
+//	
+//	Copyright (c) 2016-
+//	Authors:
+//	* Joachim Klein <klein@tcs.inf.tu-dresden.de> (TU Dresden)
+//	
+//------------------------------------------------------------------------------
+//	
+//	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.BitSet;
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.PriorityQueue;
+import java.util.function.IntPredicate;
+
+import common.StopWatch;
+import common.IterableBitSet;
+import explicit.IncomingChoiceRelation.Choice;
+import explicit.rewards.MDPRewards;
+import prism.PrismComponent;
+
+/**
+ * An implementation of the upper bound computation for Rmin as detailed in the
+ * paper McMahan, Likhachev, Gordon "Bounded Real-Time Dynamic Programming:
+ * RTDP with monotone upper bounds and performance guarantees" (International
+ * Conference on Machine Learning, 2005).
+ * */
+public class DijkstraSweepMPI {
+
+	private static class QueueEntry implements Comparable<QueueEntry> {
+		public int y;
+		public double p;
+		public double w;
+
+		public QueueEntry(int y, double p, double w)
+		{
+			this.y = y;
+			this.p = p;
+			this.w = w;
+		}
+
+		@Override
+		public int compareTo(QueueEntry o)
+		{
+			int r = Double.compare(p, o.p);
+			if (r == 0) {
+				return Double.compare(w, o.w);
+			} else {
+				return r;
+			}
+		}
+	}
+
+	private static class ChoiceValues {
+		public double p;
+		public double w;
+
+		public ChoiceValues(double p, double w)
+		{
+			this.p = p;
+			this.w = w;
+		}
+	}
+
+	private static boolean debug = false;
+	private MDP mdp;
+	private MDPRewards rewards;
+	private PriorityQueue<QueueEntry> queue;
+	private double[] pState;
+	private double[] wState;
+	private HashMap<Choice, ChoiceValues> choiceValues = new HashMap<Choice, ChoiceValues>();
+	private QueueEntry[] pri;
+	private int[] pi;
+	private BitSet unknown, target;
+	private BitSet fin = new BitSet();
+	private IncomingChoiceRelation incoming;
+	private double lambda;
+
+	private DijkstraSweepMPI(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown)
+	{
+		this.mdp = mdp;
+		this.unknown = unknown;
+		this.target = target;
+		this.rewards = rewards;
+
+		incoming = IncomingChoiceRelation.forModel(parent, mdp);
+
+		queue = new PriorityQueue<QueueEntry>();
+		pState = new double[mdp.getNumStates()];
+		wState = new double[mdp.getNumStates()];
+		pri = new QueueEntry[mdp.getNumStates()];
+		pi = new int[mdp.getNumStates()];
+
+		for (int s : IterableBitSet.getSetBits(unknown)) {
+			for (int choice = 0, numChoices = mdp.getNumChoices(s); choice < numChoices; choice++) {
+				Choice c = new Choice(s, choice);
+				double rew = rewards.getStateReward(s);
+				rew += rewards.getTransitionReward(s, choice);
+				choiceValues.put(c, new ChoiceValues(0.0, rew));
+			}
+		}
+
+		for (int s : IterableBitSet.getSetBits(target)) {
+			pState[s] = 1.0;
+		}
+
+		HashSet<Choice> preTarget = new HashSet<Choice>();
+		for (int t : IterableBitSet.getSetBits(target)) {
+			for (Choice c : incoming.getIncomingChoices(t)) {
+				boolean newChoice = preTarget.add(c);
+				if (newChoice) {
+					if (!unknown.get(c.getState())) {
+						continue;
+					}
+					if (!validChoice(c)) {
+						continue;
+					}
+
+					update(c, target);
+				}
+			}
+		}
+		preTarget.clear();
+
+		sweep();
+
+		computeLambda();
+	}
+
+	private void sweep()
+	{
+		while (!queue.isEmpty()) {
+			int x = queue.poll().y;
+			if (fin.get(x)) {
+				// already handled
+				continue;
+			}
+
+			fin.set(x);
+			ChoiceValues v = choiceValues.get(new Choice(x, pi[x]));
+			wState[x] = v.w;
+			pState[x] = v.p;
+
+			for (Choice c : incoming.getIncomingChoices(x)) {
+				if (fin.get(c.getState())) {
+					// already handled, skip
+					continue;
+				}
+
+				if (!unknown.get(c.getState())) {
+					// uninteresting state
+					continue;
+				}
+
+				if (!validChoice(c)) {
+					// some successor go outside unknown U target (e.g., to some infinity or undefined state)
+					// skip
+					continue;
+				}
+
+				// a relevant choice, update
+				update(c, x);
+			}
+		}
+	}
+
+	private boolean validChoice(Choice choice)
+	{
+		IntPredicate outsideRelevant = (int t) -> {
+			if (unknown.get(t) || target.get(t)) return false;
+			return true;
+		};
+
+		return !mdp.someSuccessorsMatch(choice.getState(), choice.getChoice(), outsideRelevant);
+	}
+
+	private void update(Choice choice, int x)
+	{
+		double w_x = wState[x];
+		// compute P^a_yx * w(x)
+		double Pw = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> {
+			if (t != x) return 0.0;
+			return p * w_x;
+		});
+
+		double p_x = pState[x];
+		// compute P^a_yx * p_g(x)
+		double Pp = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> {
+			if (t != x) return 0.0;
+			return p * p_x;
+		});
+
+		ChoiceValues c = choiceValues.get(choice);
+		assert(c != null);
+		c.p += Pp;
+		c.w += Pw;
+
+		QueueEntry newPri = new QueueEntry(choice.getState(), 1 - c.p, c.w);
+		if (pri[choice.getState()] == null || newPri.compareTo(pri[choice.getState()]) < 0) {
+			pri[choice.getState()] = newPri;
+			pi[choice.getState()] = choice.getChoice();
+			queue.add(newPri);
+		}
+	}
+
+	private void update(Choice choice, BitSet target)
+	{
+		// compute P^a_y->target
+		double Pp = mdp.sumOverTransitions(choice.getState(), choice.getChoice(), (int s, int t, double p) -> {
+			if (target.get(t)) return p;
+			return 0.0;
+		});
+
+		ChoiceValues c = choiceValues.get(choice);
+		c.p += Pp;
+
+		QueueEntry newPri = new QueueEntry(choice.getState(), 1 - c.p, c.w);
+		if (pri[choice.getState()] == null || newPri.compareTo(pri[choice.getState()]) < 0) {
+			pri[choice.getState()] = newPri;
+			pi[choice.getState()] = choice.getChoice();
+			queue.add(newPri);
+		}
+	}
+
+	private double computeLambda()
+	{
+		lambda = 0.0;
+
+		for (int x : IterableBitSet.getSetBits(unknown)) {
+			int a = pi[x];
+			double lambda_x_a = Double.POSITIVE_INFINITY;
+
+			// check condition (I)
+			double I_sum = mdp.sumOverTransitions(x, a, (int s, int t, double p) -> {
+				return p * pState[t];
+			});
+
+			if (pState[x] < I_sum) {
+				// condition (I) holds
+				double den = rewards.getStateReward(x) + rewards.getTransitionReward(x, a); // c(x,a)
+				den += mdp.sumOverTransitions(x, a, (int s, int t, double p) -> {
+					return p * wState[t];
+				});
+				den -= wState[x];
+
+				double num = mdp.sumOverTransitions(x, a, (int s, int t, double p) -> {
+					return p * pState[t];
+				});
+				num -= pState[x];
+
+				lambda_x_a = den / num;
+			} else {
+				// TODO: check condition (II)
+				lambda_x_a = 0;
+			}
+
+			lambda = Double.max(lambda, lambda_x_a);
+		}
+
+		return lambda;
+	}
+
+	public static double[] computeUpperBounds(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown)
+	{
+		StopWatch timer = new StopWatch(parent.getLog());
+		timer.start("computing upper bound(s) for Rmin using the DSI-MP algorithm");
+
+		parent.getLog().println("Computing upper bound(s) for Rmin using the Dijkstra Sweep for Monotone Pessimistic Initialization algorithm...");
+		double[] upperBounds = new double[mdp.getNumStates()];
+		DijkstraSweepMPI dsmpi = new DijkstraSweepMPI(parent, mdp, rewards, target, unknown);
+
+		for (int x : IterableBitSet.getSetBits(unknown)) {
+			upperBounds[x] = dsmpi.wState[x] + dsmpi.lambda*(1 - dsmpi.pState[x]);
+		}
+
+		if (debug) {
+			parent.getLog().println(upperBounds);
+		}
+
+		timer.stop();
+		return upperBounds;
+	}
+
+	public static double computeUpperBound(PrismComponent parent, MDP mdp, MDPRewards rewards, BitSet target, BitSet unknown)
+	{
+		double bound = 0.0;
+		final double[] upperBounds = computeUpperBounds(parent, mdp, rewards, target, unknown);
+		for (int s : IterableBitSet.getSetBits(unknown)) {
+			bound = Double.max(bound, upperBounds[s]);
+		}
+		return bound;
+	}
+}

prism/src/explicit/MDPModelChecker.java

@@ -31,6 +31,8 @@ import java.util.BitSet;
 import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
+import java.util.Map.Entry;
+import java.util.PrimitiveIterator;
 import java.util.Vector;
 
 import common.IterableStateSet;
@@ -39,6 +41,7 @@ import parser.VarList;
 import parser.ast.Declaration;
 import parser.ast.DeclarationIntUnbounded;
 import parser.ast.Expression;
+import prism.OptionsIntervalIteration;
 import prism.Prism;
 import prism.PrismComponent;
 import prism.PrismDevNullLog;
@@ -56,6 +59,7 @@ import automata.LTL2WDBA;
 import common.IntSet;
 import common.IterableBitSet;
 import explicit.modelviews.EquivalenceRelationInteger;
+import explicit.modelviews.MDPDroppedAllChoices;
 import explicit.modelviews.MDPEquiv;
 import explicit.rewards.MCRewards;
 import explicit.rewards.MCRewardsFromMDPRewards;
@@ -357,6 +361,23 @@ public class MDPModelChecker extends ProbModelChecker
 			if (!min)
 				throw new PrismException("Value iteration from above only works for minimum probabilities");
 		}
+		if (doIntervalIteration) {
+			if (!min && (genStrat || exportAdv)) {
+				throw new PrismNotSupportedException("Currently, explicit engine does not support adversary construction for interval iteration and Pmax");
+			}
+			if (mdpSolnMethod != MDPSolnMethod.VALUE_ITERATION && mdpSolnMethod != MDPSolnMethod.GAUSS_SEIDEL) {
+				throw new PrismNotSupportedException("Currently, explicit engine only supports interval iteration with value iteration or Gauss-Seidel for MDPs");
+			}
+			if (init != null)
+				throw new PrismNotSupportedException("Interval iteration currently not supported with provided initial values");
+			if (!(precomp && prob0 && prob1)) {
+				throw new PrismNotSupportedException("Precomputations (Prob0 & Prob1) must be enabled for interval iteration");
+			}
+
+			if (!min) {
+				doPmaxQuotient = true;
+			}
+		}
 		if (mdpSolnMethod == MDPSolnMethod.POLICY_ITERATION || mdpSolnMethod == MDPSolnMethod.MODIFIED_POLICY_ITERATION) {
 			if (known != null) {
 				throw new PrismException("Policy iteration methods cannot be passed 'known' values for some states");
@@ -535,9 +556,15 @@ public class MDPModelChecker extends ProbModelChecker
 			iterationMethod = new IterationMethodGS(termCrit == TermCrit.ABSOLUTE, termCritParam, false);
 			break;
 		case POLICY_ITERATION:
+			if (doIntervalIteration) {
+				throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration");
+			}
 			res = computeReachProbsPolIter(mdp, no, yes, min, strat);
 			break;
 		case MODIFIED_POLICY_ITERATION:
+			if (doIntervalIteration) {
+				throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration");
+			}
 			res = computeReachProbsModPolIter(mdp, no, yes, min, strat);
 			break;
 		default:
@@ -545,7 +572,11 @@ public class MDPModelChecker extends ProbModelChecker
 		}
 
 		if (res == null) { // not yet computed, use iterationMethod
-			res = doValueIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat);
+			if (!doIntervalIteration) {
+				res = doValueIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat);
+			} else {
+				res = doIntervalIterationReachProbs(mdp, no, yes, min, init, known, iterationMethod, getDoTopologicalValueIteration(), strat);
+			}
 		}
 
 		return res;
@@ -850,6 +881,113 @@ public class MDPModelChecker extends ProbModelChecker
 		}
 	}
 
+	/**
+	 * Compute reachability probabilities using interval iteration.
+	 * Optionally, store optimal (memoryless) strategy info.
+	 * @param mdp The MDP
+	 * @param no Probability 0 states
+	 * @param yes Probability 1 states
+	 * @param min Min or max probabilities (true=min, false=max)
+	 * @param init Optionally, an initial solution vector (will be overwritten)
+	 * @param known Optionally, a set of states for which the exact answer is known
+	 * @param iterationMethod The iteration method
+	 * @param topological Do topological value iteration?
+	 * @param strat Storage for (memoryless) strategy choice indices (ignored if null)
+	 * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values.
+	 */
+	protected ModelCheckerResult doIntervalIterationReachProbs(MDP mdp, BitSet no, BitSet yes, boolean min, double init[], BitSet known, IterationMethod iterationMethod, boolean topological, int strat[])
+			throws PrismException
+	{
+		BitSet unknown;
+		int i, n;
+		double initBelow[], initAbove[];
+		long timer;
+
+		// Start value iteration
+		timer = System.currentTimeMillis();
+		String description = (min ? "min" : "max")
+				+ (topological ? ", topological": "" )
+				+ ", with " + iterationMethod.getDescriptionShort();
+
+		mainLog.println("Starting interval iteration (" + description + ")...");
+
+		ExportIterations iterationsExport = null;
+		if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) {
+			iterationsExport = new ExportIterations("Explicit MDP ReachProbs interval iteration (" + description + ")");
+		}
+
+		// Store num states
+		n = mdp.getNumStates();
+
+		// Create solution vector(s)
+		initBelow = (init == null) ? new double[n] : init;
+		initAbove = new double[n];
+
+		// Initialise solution vectors. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 1.0/0.0 if in yes/no; (3) initVal
+		// where initVal is 0.0 or 1.0, depending on whether we converge from below/above.
+		if (known != null && init != null) {
+			for (i = 0; i < n; i++) {
+				initBelow[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 0.0;
+				initAbove[i] = known.get(i) ? init[i] : yes.get(i) ? 1.0 : no.get(i) ? 0.0 : 1.0;
+			}
+		} else {
+			for (i = 0; i < n; i++) {
+				initBelow[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 :  0.0;
+				initAbove[i] = yes.get(i) ? 1.0 : no.get(i) ? 0.0 :  1.0;
+			}
+		}
+
+		// Determine set of states actually need to compute values for
+		unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(yes);
+		unknown.andNot(no);
+		if (known != null)
+			unknown.andNot(known);
+
+		if (iterationsExport != null) {
+			iterationsExport.exportVector(initBelow, 0);
+			iterationsExport.exportVector(initAbove, 1);
+		}
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow();
+		final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove();
+		final boolean checkMonotonic = iiOptions.isCheckMonotonicity();
+
+		if (!enforceMonotonicFromAbove) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above.");
+		}
+		if (!enforceMonotonicFromBelow) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below.");
+		}
+
+		IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultMinMaxInterval(mdp, min, strat, true, enforceMonotonicFromBelow, checkMonotonic);
+		IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultMinMaxInterval(mdp, min, strat, false, enforceMonotonicFromAbove, checkMonotonic);
+		below.init(initBelow);
+		above.init(initAbove);
+
+		IntSet unknownStates = IntSet.asIntSet(unknown);
+
+		if (topological) {
+			// Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering
+			// states in unknown
+			SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get);
+
+			IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> {
+				soln[s] = mdp.mvMultJacMinMaxSingle(s, soln, min, strat);
+			};
+
+			// run the actual value iteration
+			return iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport);
+		} else {
+			// run the actual value iteration
+			return iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport);
+		}
+	}
+
 	/**
 	 * Compute reachability probabilities using Gauss-Seidel (including Jacobi-style updates).
 	 * @param mdp The MDP
@@ -1255,6 +1393,453 @@ public class MDPModelChecker extends ProbModelChecker
 		return res;
 	}
 
+	/**
+	 * Compute upper bound for maximum expected reward, with the method specified in the settings.
+	 * @param mdp the model
+	 * @param mdpRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinite states
+	 * @return upper bound on Rmax=?[ F target ] for all states
+	 */
+	double computeReachRewardsMaxUpperBound(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		// inf and target states become trap states (with dropped choices)
+		BitSet trapStates = (BitSet) target.clone();
+		trapStates.or(inf);
+		MDP cleanedMDP = new MDPDroppedAllChoices(mdp, trapStates);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		double upperBound = 0.0;
+		String method = null;
+		switch (iiOptions.getBoundMethod()) {
+		case VARIANT_1_COARSE:
+			upperBound = computeReachRewardsMaxUpperBoundVariant1Coarse(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "variant 1, coarse";
+			break;
+		case VARIANT_1_FINE:
+			upperBound = computeReachRewardsMaxUpperBoundVariant1Fine(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "variant 1, fine";
+			break;
+		case DEFAULT:
+		case VARIANT_2:
+			upperBound = computeReachRewardsMaxUpperBoundVariant2(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "variant 2";
+			break;
+		case DSMPI:
+			throw new PrismNotSupportedException("Dijkstra Sweep MPI upper bound heuristic can not be used for Rmax");
+		}
+
+		if (method == null) {
+			throw new PrismException("Unknown upper bound heuristic");
+		}
+
+		mainLog.println("Upper bound for max expectation (" + method + "): " + upperBound);
+		return upperBound;
+	}
+
+	/**
+	 * Compute upper bound for minimum expected reward, with the method specified in the settings.
+	 * @param mdp the model
+	 * @param mdpRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinite states
+	 * @return upper bound on Rmin=?[ F target ] for all unknown states
+	 */
+	double computeReachRewardsMinUpperBound(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		// inf and target states become trap states (with dropped choices)
+		BitSet trapStates = (BitSet) target.clone();
+		trapStates.or(inf);
+		MDP cleanedMDP = new MDPDroppedAllChoices(mdp, trapStates);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		double upperBound = 0.0;
+		String method = null;
+		switch (iiOptions.getBoundMethod()) {
+		case DEFAULT:
+		case DSMPI:
+			upperBound = DijkstraSweepMPI.computeUpperBound(this, mdp, mdpRewards, target, unknown);
+			method = "Dijkstra Sweep MPI";
+			break;
+		case VARIANT_1_COARSE:
+			upperBound = computeReachRewardsMaxUpperBoundVariant1Coarse(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "using Rmax upper bound via variant 1, coarse";
+			break;
+		case VARIANT_1_FINE:
+			upperBound = computeReachRewardsMaxUpperBoundVariant1Fine(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "using Rmax upper bound via variant 1, fine";
+			break;
+		case VARIANT_2:
+			upperBound = computeReachRewardsMaxUpperBoundVariant2(cleanedMDP, mdpRewards, target, unknown, inf);
+			method = "using Rmax upper bound via variant 2";
+			break;
+		}
+
+		if (method == null) {
+			throw new PrismException("Unknown upper bound heuristic");
+		}
+
+		mainLog.println("Upper bound for min expectation (" + method + "): " + upperBound);
+		return upperBound;
+	}
+
+	/**
+	 * Return true if the MDP is contracting for all states in the 'unknown'
+	 * set, i.e., if Pmin=1( unknown U target) holds.
+	 */
+	private boolean isContracting(MDP mdp, BitSet unknown, BitSet target)
+	{
+		// compute Pmin=1( unknown U target )
+		BitSet pmin1 = prob1(mdp, unknown, target, true, null);
+		BitSet tmp = (BitSet) unknown.clone();
+		tmp.andNot(pmin1);
+		if (!tmp.isEmpty()) {
+			// unknown is not contained in pmin1, not contracting
+			return false;
+		}
+		return true;
+	}
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 1, coarse),
+	 * i.e., does not compute separate q_t / p_t per SCC.
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param mdp the model
+	 * @param mdpRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @return upper bound on Rmax=?[ F target ] for all states
+	 */
+	double computeReachRewardsMaxUpperBoundVariant1Coarse(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] boundsOnExpectedVisits = new double[mdp.getNumStates()];
+		double[] maxRews = new double[mdp.getNumStates()];
+		int[] Ct = new int[mdp.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for maximal expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, null);
+		BitSet trivial = new BitSet();
+
+		double q = 0;
+		for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) {
+			IntSet statesForSCC = sccs.getStatesForSCC(scc);
+
+			int cardinality = statesForSCC.cardinality();
+
+			PrimitiveIterator.OfInt itSCC = statesForSCC.iterator();
+			while (itSCC.hasNext()) {
+				int s = itSCC.nextInt();
+				Ct[s] = cardinality;
+
+				boolean hasSelfloop = false;
+				for (int ch = 0; ch < mdp.getNumChoices(s); ch++) {
+					double probRemain = 0;
+					boolean allRemain = true;  // all successors remain in the SCC?
+					for (Iterator<Entry<Integer, Double>> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) {
+						Entry<Integer, Double> t = it.next();
+						if (statesForSCC.get(t.getKey())) {
+							probRemain += t.getValue();
+							hasSelfloop = true;
+						} else {
+							allRemain = false;
+						}
+					}
+
+					if (!allRemain) { // action in the set X
+						q = Math.max(q, probRemain);
+					}
+				}
+
+				if (cardinality == 1 && !hasSelfloop) {
+					trivial.set(s);
+				}
+			}
+		}
+
+		double p = 1;
+		for (int s = 0; s < mdp.getNumStates(); s++) {
+			double maxRew = 0;
+			for (int ch = 0; ch < mdp.getNumChoices(s); ch++) {
+				for (Iterator<Entry<Integer, Double>> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) {
+					Entry<Integer, Double> t = it.next();
+					p = Math.min(p, t.getValue());
+
+					double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch);
+					maxRew = Math.max(maxRew, rew);
+				}
+			}
+			maxRews[s] = maxRew;
+		}
+
+		double upperBound = 0;
+		for (int s = 0; s < mdp.getNumStates(); s++) {
+			if (target.get(s) || inf.get(s)) {
+				// inf or target states: not relevant, set visits to 0, ignore in summation
+				boundsOnExpectedVisits[s] = 0.0;
+			} else if (unknown.get(s)) {
+				if (trivial.get(s)) {
+					// s is a trivial SCC: seen at most once
+					boundsOnExpectedVisits[s] = 1.0;
+				} else {
+					boundsOnExpectedVisits[s] = 1 / (Math.pow(p, Ct[s]-1) * (1.0-q));
+				}
+				upperBound += boundsOnExpectedVisits[s] * maxRews[s];
+			}
+		}
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 1, coarse):");
+			mainLog.println("p = " + p);
+			mainLog.println("q = " + q);
+			mainLog.println("|Ct| = " + Arrays.toString(Ct));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+			mainLog.println("maxRews = " + Arrays.toString(maxRews));
+		}
+
+		timer.stop();
+		// mainLog.println("Upper bound for max expectation (variant 1, coarse): " + upperBound);
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 1, fine).
+	 * i.e., does compute separate q_t / p_t per SCC.
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param mdp the model
+	 * @param mdpRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @return upper bound on Rmax=?[ F target ] for all states
+	 */
+	double computeReachRewardsMaxUpperBoundVariant1Fine(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] boundsOnExpectedVisits = new double[mdp.getNumStates()];
+		double[] qt = new double[mdp.getNumStates()];
+		double[] pt = new double[mdp.getNumStates()];
+		double[] maxRews = new double[mdp.getNumStates()];
+		int[] Ct = new int[mdp.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for maximal expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, null);
+		BitSet trivial = new BitSet();
+
+		for (int scc = 0, numSCCs = sccs.getNumSCCs(); scc < numSCCs; scc++) {
+			IntSet statesForSCC = sccs.getStatesForSCC(scc);
+
+			double q = 0;
+			double p = 1;
+
+			int cardinality = statesForSCC.cardinality();
+
+			PrimitiveIterator.OfInt itSCC = statesForSCC.iterator();
+			while (itSCC.hasNext()) {
+				int s = itSCC.nextInt();
+
+				Ct[s] = cardinality;
+				boolean hasSelfloop = false;
+
+				for (int ch = 0; ch < mdp.getNumChoices(s); ch++) {
+
+					double probRemain = 0;
+					boolean allRemain = true;  // all successors remain in the SCC?
+					for (Iterator<Entry<Integer, Double>> it = mdp.getTransitionsIterator(s, ch); it.hasNext(); ) {
+						Entry<Integer, Double> t = it.next();
+						if (statesForSCC.get(t.getKey())) {
+							probRemain += t.getValue();
+							p = Math.min(p, t.getValue());
+							hasSelfloop = true;
+						} else {
+							allRemain = false;
+						}
+					}
+
+					if (!allRemain) { // action in the set Xt
+						q = Math.max(q, probRemain);
+					}
+				}
+
+				if (cardinality == 1 && !hasSelfloop) {
+					trivial.set(s);
+				}
+			}
+
+			for (int s : statesForSCC) {
+				qt[s] = q;
+				pt[s] = p;
+			}
+		}
+
+		for (int s = 0; s < mdp.getNumStates(); s++) {
+			double maxRew = 0;
+			for (int ch = 0; ch < mdp.getNumChoices(s); ch++) {
+				double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch);
+				maxRew = Math.max(maxRew, rew);
+			}
+			maxRews[s] = maxRew;
+		}
+
+		double upperBound = 0;
+		for (int s = 0; s < mdp.getNumStates(); s++) {
+			if (target.get(s) || inf.get(s)) {
+				// inf or target states: not relevant, set visits to 0, ignore in summation
+				boundsOnExpectedVisits[s] = 0.0;
+			} else if (unknown.get(s)) {
+				if (trivial.get(s)) {
+					// s is a trivial SCC: seen at most once
+					boundsOnExpectedVisits[s] = 1.0;
+				} else {
+					boundsOnExpectedVisits[s] = 1 / (Math.pow(pt[s], Ct[s]-1) * (1.0-qt[s]));
+				}
+				upperBound += boundsOnExpectedVisits[s] * maxRews[s];
+			}
+		}
+
+		timer.stop();
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 1, fine):");
+			mainLog.println("pt = " + Arrays.toString(pt));
+			mainLog.println("qt = " + Arrays.toString(qt));
+			mainLog.println("|Ct| = " + Arrays.toString(Ct));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+			mainLog.println("maxRews = " + Arrays.toString(maxRews));
+		}
+
+		// mainLog.println("Upper bound for max expectation (variant 1, fine): " + upperBound);
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
+	/**
+	 * Compute upper bound for maximum expected reward (variant 2).
+	 * Uses Rs = S, i.e., does not take reachability into account.
+	 * @param dtmc the model
+	 * @param mcRewards the rewards
+	 * @param target the target states
+	 * @param unknown the states that are not target or infinity states
+	 * @param inf the infinity states
+	 * @return upper bound on R=?[ F target ] for all states
+	 */
+	double computeReachRewardsMaxUpperBoundVariant2(MDP mdp, MDPRewards mdpRewards, BitSet target, BitSet unknown, BitSet inf) throws PrismException
+	{
+		double[] dt = new double[mdp.getNumStates()];
+		double[] boundsOnExpectedVisits = new double[mdp.getNumStates()];
+		double[] maxRews = new double[mdp.getNumStates()];
+
+		StopWatch timer = new StopWatch(getLog());
+		timer.start("computing an upper bound for expected reward");
+
+		SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get);
+
+		BitSet T = (BitSet) target.clone();
+
+		@SuppressWarnings("unused")
+		int i = 0;
+		while (true) {
+			BitSet Si = new BitSet();
+			i++;
+
+			// TODO: might be inefficient, worst-case quadratic runtime...
+			for (PrimitiveIterator.OfInt it = IterableBitSet.getClearBits(T, mdp.getNumStates() -1 ).iterator(); it.hasNext(); ) {
+				int s = it.nextInt();
+				// mainLog.println("Check " + s + " against " + T);
+				boolean allActionsReachT = true;
+				for (int choice = 0, choices = mdp.getNumChoices(s); choice < choices; choice++) {
+					if (!mdp.someSuccessorsInSet(s, choice, T)) {
+						allActionsReachT = false;
+						break;
+					}
+				}
+				if (allActionsReachT) {
+					Si.set(s);
+				}
+			}
+
+			if (Si.isEmpty()) {
+				break;
+			}
+
+			// mainLog.println("S" + i + " = " + Si);
+			// mainLog.println("T = " + T);
+
+			for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(Si).iterator(); it.hasNext(); ) {
+				final int t = it.nextInt();
+				final int sccIndexForT = sccs.getSCCIndex(t);
+
+				double min = Double.POSITIVE_INFINITY;
+				for (int choice = 0, choices = mdp.getNumChoices(t); choice < choices; choice++) {
+					// mainLog.println("State " + t + ", choice = " + choice);
+					double d = mdp.sumOverTransitions(t, choice, (int __, int u, double prob) -> {
+						// mainLog.println("t = " + t + ", u = " + u + ", prob = " + prob);
+						if (!T.get(u))
+							return 0.0;
+
+						boolean inSameSCC = (sccs.getSCCIndex(u) == sccIndexForT);
+						double d_u_t = inSameSCC ? dt[u] : 1.0;
+						// mainLog.println("d_u_t = " + d_u_t);
+						return d_u_t * prob;
+					});
+					if (d < min) {
+						min = d;
+					}
+				}
+				dt[t] = min;
+				// mainLog.println("d["+t+"] = " + dt[t]);
+			}
+
+			T.or(Si);
+		}
+
+		for (int s = 0; s < mdp.getNumStates(); s++) {
+			double maxRew = 0;
+			for (int ch = 0; ch < mdp.getNumChoices(s); ch++) {
+				double rew = mdpRewards.getStateReward(s) + mdpRewards.getTransitionReward(s, ch);
+				maxRew = Math.max(maxRew, rew);
+			}
+			maxRews[s] = maxRew;
+		}
+
+		double upperBound = 0;
+		for (PrimitiveIterator.OfInt it = IterableBitSet.getSetBits(unknown).iterator(); it.hasNext();) {
+			int s = it.nextInt();
+			boundsOnExpectedVisits[s] = 1 / dt[s];
+			upperBound += boundsOnExpectedVisits[s] * maxRews[s];
+		}
+
+		timer.stop();
+
+		if (OptionsIntervalIteration.from(this).isBoundComputationVerbose()) {
+			mainLog.println("Upper bound for max expectation computation (variant 2):");
+			mainLog.println("d_t = " + Arrays.toString(dt));
+			mainLog.println("ζ* = " + Arrays.toString(boundsOnExpectedVisits));
+		}
+
+		// mainLog.println("Upper bound for expectation (variant 2): " + upperBound);
+
+		if (!Double.isFinite(upperBound)) {
+			throw new PrismException("Problem computing an upper bound for the expectation, did not get finite result");
+		}
+
+		return upperBound;
+	}
+
 	/**
 	 * Compute expected instantaneous reward,
 	 * i.e. compute the min/max expected reward of the states after {@code k} steps.
@@ -1477,7 +2062,12 @@ public class MDPModelChecker extends ProbModelChecker
 				throw new PrismException("Policy iteration methods cannot be passed 'known' values for some states");
 			}
 		}
-		
+		if (doIntervalIteration) {
+			if (mdpSolnMethod != MDPSolnMethod.VALUE_ITERATION && mdpSolnMethod != MDPSolnMethod.GAUSS_SEIDEL) {
+				throw new PrismNotSupportedException("Currently, explicit engine only supports interval iteration with value iteration or Gauss-Seidel for MDPs");
+			}
+		}
+
 		// Start expected reachability
 		timer = System.currentTimeMillis();
 		mainLog.println("\nStarting expected reachability (" + (min ? "min" : "max") + ")...");
@@ -1624,6 +2214,9 @@ public class MDPModelChecker extends ProbModelChecker
 			iterationMethod = new IterationMethodGS(termCrit == TermCrit.ABSOLUTE, termCritParam, false);
 			break;
 		case POLICY_ITERATION:
+			if (doIntervalIteration) {
+				throw new PrismNotSupportedException("Interval iteration currently not supported for policy iteration");
+			}
 			res = computeReachRewardsPolIter(mdp, mdpRewards, target, inf, min, strat);
 			break;
 		default:
@@ -1631,7 +2224,11 @@ public class MDPModelChecker extends ProbModelChecker
 		}
 
 		if (res == null) { // not yet computed, use iterationMethod
-			res = doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat);
+			if (!doIntervalIteration) {
+				res = doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat);
+			} else {
+				res = doIntervalIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, getDoTopologicalValueIteration(), strat);
+			}
 		}
 
 		return res;
@@ -1760,6 +2357,149 @@ public class MDPModelChecker extends ProbModelChecker
 		return doValueIterationReachRewards(mdp, mdpRewards, iterationMethod, target, inf, min, init, known, false, strat);
 	}
 
+	/**
+	 * Compute expected reachability rewards using interval iteration
+	 * Optionally, store optimal (memoryless) strategy info.
+	 * @param mdp The MDP
+	 * @param mdpRewards The rewards
+	 * @param target Target states
+	 * @param inf States for which reward is infinite
+	 * @param min Min or max rewards (true=min, false=max)
+	 * @param init Optionally, an initial solution vector (will be overwritten)
+	 * @param known Optionally, a set of states for which the exact answer is known
+	 * @param topological do topological interval iteration
+	 * @param strat Storage for (memoryless) strategy choice indices (ignored if null)
+	 * Note: if 'known' is specified (i.e. is non-null, 'init' must also be given and is used for the exact values.
+	 */
+	protected ModelCheckerResult doIntervalIterationReachRewards(MDP mdp, MDPRewards mdpRewards, IterationMethod iterationMethod, BitSet target, BitSet inf, boolean min, double init[], BitSet known, boolean topological, int strat[])
+			throws PrismException
+	{
+		BitSet unknown;
+		int i, n;
+		double initBelow[], initAbove[];
+		long timer;
+
+		// Store num states
+		n = mdp.getNumStates();
+
+		// Determine set of states actually need to compute values for
+		unknown = new BitSet();
+		unknown.set(0, n);
+		unknown.andNot(target);
+		unknown.andNot(inf);
+		if (known != null)
+			unknown.andNot(known);
+
+		OptionsIntervalIteration iiOptions = OptionsIntervalIteration.from(this);
+
+		double upperBound;
+		if (iiOptions.hasManualUpperBound()) {
+			upperBound = iiOptions.getManualUpperBound();
+			getLog().printWarning("Upper bound for interval iteration manually set to " + upperBound);
+		} else {
+			if (min) {
+				upperBound = computeReachRewardsMinUpperBound(mdp, mdpRewards, target, unknown, inf);
+			} else {
+				upperBound = computeReachRewardsMaxUpperBound(mdp, mdpRewards, target, unknown, inf);
+			}
+		}
+
+		double lowerBound;
+		if (iiOptions.hasManualLowerBound()) {
+			lowerBound = iiOptions.getManualLowerBound();
+			getLog().printWarning("Lower bound for interval iteration manually set to " + lowerBound);
+		} else {
+			lowerBound = 0.0;
+		}
+
+		if (min) {
+			if (!isContracting(mdp, unknown, target)) {
+				throw new PrismNotSupportedException("Interval iteration for Rmin and non-contracting MDP currently not supported");
+			} else {
+				mainLog.println("Relevant sub-MDP is contracting, proceed...");
+			}
+		}
+
+		// Start value iteration
+		timer = System.currentTimeMillis();
+		String description = (min ? "min" : "max") + (topological ? ", topological" : "") + ", with " + iterationMethod.getDescriptionShort();
+		mainLog.println("Starting interval iteration (" + description + ")...");
+
+		ExportIterations iterationsExport = null;
+		if (settings.getBoolean(PrismSettings.PRISM_EXPORT_ITERATIONS)) {
+			iterationsExport = new ExportIterations("Explicit MDP ReachRewards interval iteration (" + description + ")");
+		}
+
+		// Create initial solution vector(s)
+		initBelow = (init == null) ? new double[n] : init;
+		initAbove = new double[n];
+
+		// Initialise solution vector from below. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) lowerBound
+		if (init != null && known != null) {
+			for (i = 0; i < n; i++)
+				initBelow[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound;
+		} else {
+			for (i = 0; i < n; i++)
+				initBelow[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : lowerBound;
+		}
+
+		// Initialise solution vector from above. Use (where available) the following in order of preference:
+		// (1) exact answer, if already known; (2) 0.0/infinity if in target/inf; (3) upperBound
+		if (init != null && known != null) {
+			for (i = 0; i < n; i++)
+				initAbove[i] = known.get(i) ? init[i] : target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound;
+		} else {
+			for (i = 0; i < n; i++)
+				initAbove[i] = target.get(i) ? 0.0 : inf.get(i) ? Double.POSITIVE_INFINITY : upperBound;
+		}
+
+		if (iterationsExport != null) {
+			iterationsExport.exportVector(initBelow, 0);
+			iterationsExport.exportVector(initAbove, 1);
+		}
+
+		final boolean enforceMonotonicFromBelow = iiOptions.isEnforceMonotonicityFromBelow();
+		final boolean enforceMonotonicFromAbove = iiOptions.isEnforceMonotonicityFromAbove();
+		final boolean checkMonotonic = iiOptions.isCheckMonotonicity();
+
+		if (!enforceMonotonicFromAbove) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from above.");
+		}
+		if (!enforceMonotonicFromBelow) {
+			getLog().println("Note: Interval iteration is configured to not enforce monotonicity from below.");
+		}
+
+		IterationMethod.IterationIntervalIter below = iterationMethod.forMvMultRewMinMaxInterval(mdp, mdpRewards, min, strat, true, enforceMonotonicFromBelow, checkMonotonic);
+		IterationMethod.IterationIntervalIter above = iterationMethod.forMvMultRewMinMaxInterval(mdp, mdpRewards, min, strat, false, enforceMonotonicFromAbove, checkMonotonic);
+		below.init(initBelow);
+		above.init(initAbove);
+
+		IntSet unknownStates = IntSet.asIntSet(unknown);
+
+		ModelCheckerResult rv;
+		if (topological) {
+			// Compute SCCInfo, including trivial SCCs in the subgraph obtained when only considering
+			// states in unknown
+			SCCInfo sccs = SCCComputer.computeTopologicalOrdering(this, mdp, true, unknown::get);
+
+			IterationMethod.SingletonSCCSolver singletonSCCSolver = (int s, double[] soln) -> {
+				soln[s] = mdp.mvMultRewJacMinMaxSingle(s, soln, mdpRewards, min, strat);
+			};
+
+			// run the actual value iteration
+			rv = iterationMethod.doTopologicalIntervalIteration(this, description, sccs, below, above, singletonSCCSolver, timer, iterationsExport);
+		} else {
+			// run the actual value iteration
+			rv = iterationMethod.doIntervalIteration(this, description, below, above, unknownStates, timer, iterationsExport);
+		}
+
+		double max_v = PrismUtils.findMaxFinite(rv.soln, unknownStates.iterator());
+		mainLog.println("Maximum finite value in solution vector at end of interval iteration: " + max_v);
+
+		return rv;
+	}
+
 	/**
 	 * Compute expected reachability rewards using policy iteration.
 	 * The array {@code strat} is used both to pass in the initial strategy for policy iteration,