{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d2770b36",
   "metadata": {},
   "source": [
    "# Next Steps\n",
    "Up until now, we showcased various analysis approaches for different model types in the\n",
    "previous steps.\n",
    "This notebook presents typical use cases for the Python API: Using symbolic model representations, using preprocessing algorithms, configuring solvers, and implementing own prototypes on top of existing routines."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b910307",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "As usual, we load the Prism Orchard model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "92186079",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:04.865305Z",
     "iopub.status.busy": "2026-03-26T10:41:04.865065Z",
     "iopub.status.idle": "2026-03-26T10:41:05.100849Z",
     "shell.execute_reply": "2026-03-26T10:41:05.100288Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import prism models from the previous step\n",
    "import stormvogel\n",
    "import stormpy\n",
    "from orchard.orchard_builder import build_prism"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "facd1f72",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:05.102936Z",
     "iopub.status.busy": "2026-03-26T10:41:05.102673Z",
     "iopub.status.idle": "2026-03-26T10:41:05.214322Z",
     "shell.execute_reply": "2026-03-26T10:41:05.213765Z"
    }
   },
   "outputs": [],
   "source": [
    "orchard_prism, prism_program = build_prism()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d9fe314",
   "metadata": {},
   "source": [
    "## Model representation\n",
    "So far, we have represented the MDP through an explicit (sparse) transition matrix. Instead, we can build the MDP symbolically and represent it through **decision diagrams (DD)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "019a4bd1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:05.216214Z",
     "iopub.status.busy": "2026-03-26T10:41:05.216037Z",
     "iopub.status.idle": "2026-03-26T10:41:05.541639Z",
     "shell.execute_reply": "2026-03-26T10:41:05.539211Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------------------------------------------------------- \n",
      "Model type: \tMDP (symbolic)\n",
      "States: \t22469 (85 nodes)\n",
      "Transitions: \t44954 (916 nodes)\n",
      "Choices: \t29354\n",
      "Reward Models:  rounds\n",
      "Variables: \trows: 7 meta variables (19 DD variables), columns: 7 meta variables (19 DD variables), nondeterminism: 12 meta variables (12 DD variables)\n",
      "Labels: \t6\n",
      "   * deadlock -> 0 state(s) (1 nodes)\n",
      "   * init -> 1 state(s) (20 nodes)\n",
      "   * AllCherriesPicked\n",
      "   * PlayersWon\n",
      "   * RavenOneAway\n",
      "   * RavenWon\n",
      "-------------------------------------------------------------- \n",
      "\n"
     ]
    }
   ],
   "source": [
    "orchard_symbolic = stormpy.build_symbolic_model(prism_program)\n",
    "print(orchard_symbolic)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d756653",
   "metadata": {},
   "source": [
    "The symbolic MDP has the same size as the sparse MDP, but instead of explicitly\n",
    "storing nearly 45,000 transitions, the function is represented through a DD with\n",
    "946 nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "098bcfa9",
   "metadata": {},
   "source": [
    "We can perform model checking on this symbolic representation and obtain the same winning probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "74d53fb0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:05.547932Z",
     "iopub.status.busy": "2026-03-26T10:41:05.547675Z",
     "iopub.status.idle": "2026-03-26T10:41:06.489030Z",
     "shell.execute_reply": "2026-03-26T10:41:06.488566Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARN  (SymbolicMinMaxLinearEquationSolver.cpp:64): Selected method is not supported for this solver, switching to value iteration.\n",
      "WARN  (SymbolicMinMaxLinearEquationSolver.cpp:64): Selected method is not supported for this solver, switching to value iteration.\n",
      "Maximal probability: 0.631357\n"
     ]
    }
   ],
   "source": [
    "formula = stormpy.parse_properties('Pmax=? [F \"PlayersWon\"]')[0]\n",
    "symbolic_result = stormpy.model_checking(\n",
    "    orchard_symbolic, formula, only_initial_states=True\n",
    ")\n",
    "filter = stormpy.create_filter_initial_states_symbolic(orchard_symbolic)\n",
    "symbolic_result.filter(filter)\n",
    "print(\"Maximal probability: {}\".format(symbolic_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee55d081",
   "metadata": {},
   "source": [
    "## Bisimulation\n",
    "Bisimulation minimization can be applied to minimize the state space without affecting the model behaviour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "775307c8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:06.491002Z",
     "iopub.status.busy": "2026-03-26T10:41:06.490802Z",
     "iopub.status.idle": "2026-03-26T10:41:16.909732Z",
     "shell.execute_reply": "2026-03-26T10:41:16.907025Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model with 22469 states and 44954 transitions\n",
      "Model with 956 states and 2446 transitions\n"
     ]
    }
   ],
   "source": [
    "formula = stormpy.parse_properties('Pmax=? [F \"PlayersWon\"]')\n",
    "print(\n",
    "    \"Model with {} states and {} transitions\".format(\n",
    "        orchard_prism.nr_states, orchard_prism.nr_transitions\n",
    "    )\n",
    ")\n",
    "orchard_bisim = stormpy.perform_bisimulation(\n",
    "    orchard_prism, formula, stormpy.BisimulationType.STRONG\n",
    ")\n",
    "print(\n",
    "    \"Model with {} states and {} transitions\".format(\n",
    "        orchard_bisim.nr_states, orchard_bisim.nr_transitions\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3b988c4",
   "metadata": {},
   "source": [
    "For the orchard model, we can reduce the state space size by $95\\%$ from $22,469$ states to $956$ states.\n",
    "This reduction also speeds up subsequent analysis on the reduced model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6165568",
   "metadata": {},
   "source": [
    "## Model checking algorithms\n",
    "The underlying model checking algorithms in Storm are configured through so-called environments.\n",
    "For example, we can set the precision requirement from the default value of $10^{−6}$ to $0.1$ and see how that\n",
    "affects the resulting probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8842ba6b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:16.913522Z",
     "iopub.status.busy": "2026-03-26T10:41:16.911640Z",
     "iopub.status.idle": "2026-03-26T10:41:16.983979Z",
     "shell.execute_reply": "2026-03-26T10:41:16.983396Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5815061686029692\n"
     ]
    }
   ],
   "source": [
    "# Change precision\n",
    "env = stormpy.Environment()\n",
    "prec = stormpy.Rational(0.1)\n",
    "env.solver_environment.minmax_solver_environment.precision = prec\n",
    "result = stormpy.model_checking(orchard_prism, formula[0], environment=env)\n",
    "print(result.at(orchard_prism.initial_states[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3b1ee92",
   "metadata": {},
   "source": [
    "We can see that the resulting probability $0.5815$ is within the precision guarantee of $0.1$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd5c23af",
   "metadata": {},
   "source": [
    "We can also change the underlying algorithm and compare different approaches such as **value iteration**, **policy iteration** and **optimistic value iteration**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "728e32f3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:16.988381Z",
     "iopub.status.busy": "2026-03-26T10:41:16.988113Z",
     "iopub.status.idle": "2026-03-26T10:41:17.397920Z",
     "shell.execute_reply": "2026-03-26T10:41:17.397389Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method: MinMaxMethod.value_iteration\n",
      "Result: 0.6313572986959963\n",
      "Time: 0.066s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method: MinMaxMethod.policy_iteration\n",
      "Result: 0.6313573066006365\n",
      "Time: 0.3s\n",
      "Method: MinMaxMethod.optimistic_value_iteration\n",
      "Result: 0.6313576699315778\n",
      "Time: 0.036s\n"
     ]
    }
   ],
   "source": [
    "# Change algorithm\n",
    "import time\n",
    "\n",
    "methods = [\n",
    "    stormpy.MinMaxMethod.value_iteration,\n",
    "    stormpy.MinMaxMethod.policy_iteration,\n",
    "    stormpy.MinMaxMethod.optimistic_value_iteration,\n",
    "]\n",
    "for m in methods:\n",
    "    env = stormpy.Environment()\n",
    "    env.solver_environment.minmax_solver_environment.method = m\n",
    "    start = time.time()\n",
    "    result = stormpy.model_checking(\n",
    "        orchard_prism, formula[0], environment=env, extract_scheduler=True\n",
    "    )\n",
    "    print(f\"Method: {m}\")\n",
    "    print(f\"Result: {result.at(orchard_prism.initial_states[0])}\")\n",
    "    print(f\"Time: {time.time() - start:.2}s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99190faa",
   "metadata": {},
   "source": [
    "We see that all methods provide the same result, but their timings differ.\n",
    "For example, policy iteration is typically slower.\n",
    "\n",
    "Running the above cell multiple times allows to get better insights into the performance of the different algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69083196",
   "metadata": {},
   "source": [
    "## Writing an LP-based MDP Model Checker\n",
    "One further step is extending stormpy with own algorithms. Extending stormpy allows to profit\n",
    "from existing data structure and efficient algorithms when integrating own ideas.\n",
    "We exemplify this with a fundamental approach to solve MDPs using a **linear program (LP)** solver, an approach that is popular due to the flexibility to add additional constraints.\n",
    "\n",
    "The LP encoding below computes reachability probabilities."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea8d38e0",
   "metadata": {},
   "source": [
    "Given a target label $B \\in AP$,\n",
    "this LP characterizes the variables $(x_s)_{s \\in S}$ as follows:\n",
    "\n",
    "Minimize $\\sum_{s \\in S} x_s$ such that:\n",
    "1. If $B \\in L(s)$, then $x_s =1$,\n",
    "2. If there is no path from $s$ to any state with label $B$, then $x_s=0$,\n",
    "3. Otherwise, $0 \\leq x_s \\leq 1$, and for all actions $\\alpha \\in \\textit{Act}(s)$:\n",
    "    $x_s \\geq \\sum_{s' \\in S} \\mathbb{P}(s, \\alpha, s') \\cdot x_{s'}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31ba05ac",
   "metadata": {},
   "source": [
    "We encode the LP by using functionality of stormpy.\n",
    "We first compute all states that satisfy the requirements in (1) and (2) using graph algorithms provided in stormpy.\n",
    "For instance, we compute the states $s$ which have to path to label $B$ using the function `compute_prob01max_states()`.\n",
    "We then construct and solve the LP by inspecting the model's transition matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5792500b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:17.401065Z",
     "iopub.status.busy": "2026-03-26T10:41:17.400225Z",
     "iopub.status.idle": "2026-03-26T10:41:18.917012Z",
     "shell.execute_reply": "2026-03-26T10:41:18.916407Z"
    }
   },
   "outputs": [],
   "source": [
    "# These are the target states\n",
    "players_won = stormpy.parse_properties('\"PlayersWon\"')[0].raw_formula\n",
    "psi_states = stormpy.model_checking(orchard_prism, players_won).get_truth_values()\n",
    "\n",
    "# These are the states that can never reach the target states\n",
    "phi_states = stormpy.BitVector(orchard_prism.nr_states, True)\n",
    "prob0max_states = stormpy.compute_prob01max_states(\n",
    "    orchard_prism, phi_states, psi_states\n",
    ")[0]\n",
    "\n",
    "# SCIP is an LP solver\n",
    "from pyscipopt import Model\n",
    "\n",
    "m = Model()\n",
    "\n",
    "# Create a variable for each state\n",
    "num_states = orchard_prism.nr_states\n",
    "state_vars = [m.addVar(f\"x_{i}\", lb=0, ub=1) for i in range(num_states)]\n",
    "\n",
    "# Encode LP\n",
    "for state in range(num_states):\n",
    "    if psi_states.get(state):  # Case 1\n",
    "        m.addCons(state_vars[state] == 1)\n",
    "    elif prob0max_states.get(state):  # Case 2\n",
    "        m.addCons(state_vars[state] == 0)\n",
    "    else:  # Case 3\n",
    "        for row in orchard_prism.transition_matrix.get_rows_for_group(state):\n",
    "            summed_prob = 0\n",
    "            for transition in orchard_prism.transition_matrix.get_row(row):\n",
    "                prob = transition.value()\n",
    "                next_state = transition.column\n",
    "                summed_prob += prob * state_vars[next_state]\n",
    "            m.addCons(state_vars[state] >= summed_prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3153c9d",
   "metadata": {},
   "source": [
    "Lastly, we can solve the LP and obtain the same winning probability as before.\n",
    "We also compare the LP result with the model checking result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5ed6455e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-03-26T10:41:18.920392Z",
     "iopub.status.busy": "2026-03-26T10:41:18.920020Z",
     "iopub.status.idle": "2026-03-26T10:41:22.323029Z",
     "shell.execute_reply": "2026-03-26T10:41:22.321466Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "presolving:\n",
      "(round 1, fast)       16318 del vars, 16318 del conss, 0 add conss, 1253 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 2, fast)       16347 del vars, 16351 del conss, 0 add conss, 1471 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 3, fast)       16383 del vars, 16487 del conss, 0 add conss, 2510 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 4, fast)       16496 del vars, 16679 del conss, 0 add conss, 5251 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 5, fast)       16624 del vars, 16979 del conss, 0 add conss, 9886 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 6, fast)       16868 del vars, 17410 del conss, 0 add conss, 14995 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 7, fast)       17133 del vars, 18067 del conss, 0 add conss, 20025 c0.6313573066006342 0.6313566764945008\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 8, fast)       17550 del vars, 18959 del conss, 0 add conss, 24516 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 9, fast)       18047 del vars, 20139 del conss, 0 add conss, 28476 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 10, fast)       18674 del vars, 21479 del conss, 0 add conss, 31530 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 11, fast)       19377 del vars, 23106 del conss, 0 add conss, 34095 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 12, fast)       20099 del vars, 24631 del conss, 0 add conss, 36007 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clqs, 0 implints\n",
      "(round 13, fast)       20784 del vars, 26084 del conss, 0 add conss, 37365 chg bounds, 0 chg sides, 0 chg coeffs, 0 upgd conss, 0 impls, 0 clq"
     ]
    }
   ],
   "source": [
    "# Solve LP\n",
    "m.setObjective(sum(state_vars), sense=\"minimize\")\n",
    "m.optimize()\n",
    "sol = m.getBestSol()\n",
    "result_lp = sol[state_vars[orchard_prism.initial_states[0]]]\n",
    "\n",
    "# Compare with default VI\n",
    "properties = stormpy.parse_properties_without_context('Pmax=? [F \"PlayersWon\"]')\n",
    "vi_result = stormpy.model_checking(orchard_prism, properties[0].raw_formula)\n",
    "result_vi = vi_result.at(orchard_prism.initial_states[0])\n",
    "\n",
    "print(result_lp, result_vi)\n",
    "assert abs(result_lp - result_vi) < 1e-6  # => True"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}