graph_decision_tree.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "66202bae-f86c-43b6-a945-61c3bdad6bc5",
   "metadata": {},
   "source": [
    "Click here to <a href=\"https://colab.research.google.com/github/Loop3D/MinEx-UQ-VoI-workshop/blob/main/graph_decision_tree.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afa6a997-5708-40a8-876d-14db96421bbb",
   "metadata": {},
   "source": [
    "# Decision tree analysis\n",
    "\n",
    "In this notebook, you will build a decision tree like the follwoing one, as illustrated in the presentation given by Pietro Guj.\n",
    "\n",
    "![A decision tree](img/example_decision_tree.png)\n",
    "\n",
    "To build and plot the decision tree, the code relies on the $igraph$ and $pycairo$ libraries.\n",
    "\n",
    "\n",
    "### Install missing libraries\n",
    "The following cell installs the required libraries to execute the code from https://colab.research.google.com/. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a3f8c0a-533a-42f7-9790-9c65e502adfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'google.colab' in str(get_ipython()):\n",
    "    !pip install igraph\n",
    "    !apt install libcairo2-dev\n",
    "    !pip install pycairo\n",
    "else:\n",
    "    print('The required libraries shall be available in your working environment!')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0598e96f-6cec-4a8a-b7fe-72fbd83b2877",
   "metadata": {},
   "source": [
    "### Importing libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bf5adb7-933a-48b8-a54e-c93741512e0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import igraph as ig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "066172b2-92e4-4197-9588-3d75504f1539",
   "metadata": {},
   "source": [
    "Define some parameters and functions to build and plot the decision tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f3a9388-69ab-4546-8831-944d1f70b81b",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# INTERNAL PARAMETERS AND FUNCTIONS\n",
    "node_decision_shape = \"rectangle\"\n",
    "node_decision_color = \"yellow\"\n",
    "node_event_shape = \"circle\"\n",
    "node_event_color = \"green\"\n",
    "node_end_shape = \"rectangle\"\n",
    "node_end_color = \"white\"\n",
    "\n",
    "\n",
    "def plot_decision_tree(g):   \n",
    "    # g is a graph and it shall contain the following attributes :\n",
    "    # vertices attributes: name, EV, vshape, vcolor\n",
    "    # edges attributes: name, probability, val\n",
    "    ne = len(g.es[\"probability\"])\n",
    "    lbl_array = (np.vstack((\n",
    "        np.asarray(ne*['p = ']),np.asarray(g.es[\"probability\"]),np.asarray(ne*['\\n ']),np.asarray(g.es[\"name\"]),np.asarray(ne*['\\n val = ']),np.asarray(g.es[\"value\"])))).T\n",
    "    \n",
    "    lbl_e = ne*['p = ']\n",
    "    for i in range(ne):\n",
    "        lbl_e[i] = lbl_array[i,0] + str(lbl_array[i,1]) + lbl_array[i,2]+ lbl_array[i,3]+ lbl_array[i,4] + str(lbl_array[i,5])\n",
    "    \n",
    "    nv = len(g.vs[\"EV\"])\n",
    "    \n",
    "    lbl_array = (np.vstack((\n",
    "        np.asarray(g.vs[\"name\"]),np.asarray(nv*['\\n EV: ']),np.asarray(g.vs[\"EV\"])))).T\n",
    "    \n",
    "    lbl_v = nv*['p = ']\n",
    "    for i in range(nv):\n",
    "        lbl_v[i] = lbl_array[i,0] + lbl_array[i,1] + str(lbl_array[i,2])\n",
    "    \n",
    "    g.vs[\"shape\"] = g.vs[\"vshape\"]\n",
    "    g.vs[\"color\"] = g.vs[\"vcolor\"]\n",
    "    g.vs[\"label\"] = lbl_v\n",
    "    g.es[\"label\"] = lbl_e\n",
    "    return ig.plot(g,layout=\"rt\",bbox=(0, 0, 1000, 600),margin=50,vertex_size=36)\n",
    "\n",
    "def compute_EV(g,verb=False):\n",
    "    print(\"***********************\")\n",
    "    print('compute_EV()')\n",
    "    print(\"***********************\")\n",
    "    print('computing all the paths from the tree root')\n",
    "    all_paths_from_root = g.get_shortest_paths(0, output=\"vpath\")\n",
    "    npaths = len(all_paths_from_root)\n",
    "    if verb: print(all_paths_from_root)\n",
    "    if verb: print(str(npaths)+' paths') \n",
    "    \n",
    "    # For each end node/vertex, compute the cumulative sum of the values along the edges as Expected Value of the project\n",
    "    print('retrieving leaves')\n",
    "\n",
    "    # Find the leaves of the tree\n",
    "    list_of_types = [\"end\"]\n",
    "    res=g.vs.select(type_in=list_of_types)\n",
    "    nres = len(res)\n",
    "    vs_end_idx = res.indices\n",
    "    \n",
    "    # For each leave, find the shortest path to the root and sum the edges values to get the expected value at the leaf\n",
    "    print('computing EV at leaves')\n",
    "    for i in vs_end_idx: # [9]:#\n",
    "        cur_val = 0\n",
    "        path = g.get_shortest_paths(0, to=i, output=\"vpath\")[0]\n",
    "        lpath = len(path)\n",
    "        for j in range(lpath-1):\n",
    "            cur_es = g.es.select(_source=path[lpath-(j+2)], _target=path[lpath-(j+1)])[0]\n",
    "            cur_val += cur_es[\"value\"]\n",
    "        g.vs[i][\"EV\"]=cur_val\n",
    "        if verb: \n",
    "            print(path)\n",
    "            print(g.vs[i])\n",
    "            print(cur_val)# Find paths not containing leaves\n",
    "    \n",
    "    print('computing EV at other nodes')\n",
    "    # Keep path that do not contain the leaves to compute expecte value backward on remaining nodes\n",
    "    paths_2_keep = []\n",
    "    for i in range(npaths):\n",
    "        if set(all_paths_from_root[npaths-i-1]).isdisjoint(set(vs_end_idx)): paths_2_keep.append(all_paths_from_root[npaths-i-1])\n",
    "        \n",
    "    # print(all_paths_from_root)\n",
    "    # print(paths_2_keep)\n",
    "    # Find nodes on which to compute EV (last node of kept path, and order them from longest path to shortest path)\n",
    "    vs_4_EV = np.zeros(len(paths_2_keep)).astype(int)\n",
    "    for i in range(len(paths_2_keep)):\n",
    "        vs_4_EV[i]=paths_2_keep[i][-1]\n",
    "    if verb: print(\"vs_4_EV: \"+str(vs_4_EV))\n",
    "    \n",
    "    # Compute EV\n",
    "    for i in vs_4_EV: \n",
    "        cur_EV = 0\n",
    "        # find downstream branch and nodes\n",
    "        edges_downstream = g.es.select(_source=i)\n",
    "        ned = len(edges_downstream)\n",
    "        for j in range(ned):\n",
    "            if verb: print(\"p=\"+str(edges_downstream[j][\"probability\"])+\" target=\"+str(edges_downstream[j].target)+\" value=\"+str(g.vs[edges_downstream[j].target][\"EV\"]))\n",
    "            cur_EV+=edges_downstream[j][\"probability\"]*g.vs[edges_downstream[j].target][\"EV\"]\n",
    "        g.vs[i][\"EV\"]=cur_EV\n",
    "        if verb: print('vs index: ' +str(i)+\" cur_EV=\"+str(cur_EV))\n",
    "    print('finished')\n",
    "    return g\n",
    "\n",
    "def create_tree(root_name):\n",
    "    tree = ig.Graph(directed=True)\n",
    "    tree.add_vertices(1, attributes={\"name\":\"0\",\"type\":\"decision\",\"EV\":np.nan,\"vshape\":node_decision_shape,\"vcolor\":node_decision_color})\n",
    "    return tree\n",
    "\n",
    "def add_node_and_branch(tree,node_name,parent_name,ntype=\"end\",branch_name=\"abandon\",probability=0,value_gained=0):\n",
    "    if ntype == \"decision\":\n",
    "        vshape = node_decision_shape\n",
    "        vcolor = node_decision_color\n",
    "    elif ntype == \"event\":\n",
    "        vshape = node_event_shape\n",
    "        vcolor = node_event_color\n",
    "    else:\n",
    "        vshape = node_end_shape\n",
    "        vcolor = node_end_color       \n",
    "    tree.add_vertices(1, attributes={\"name\":node_name,\"type\":ntype,\"EV\":np.nan,\"vshape\":vshape,\"vcolor\":vcolor})\n",
    "    tree.add_edges([(parent_name,node_name)], attributes={\"name\":branch_name, \"probability\":probability, \"value\":value_gained})\n",
    "    return\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7098b942-d431-4006-b075-da9b5bee6beb",
   "metadata": {},
   "source": [
    "### Decision tree parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5b463a0-3a5e-4c03-b268-088534bd926d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARAMETERS\n",
    "cost_geophy = -0.25    # Cost of Geophysocal data acquisition\n",
    "cost_drill = -2.5      # Cost of drilling\n",
    "EV_discovery = 196.69  # Excpected Value of discovery (based on mean distribution of some discoveries)\n",
    "prob_anomaly = 50/100  # Probability to detect a geophysical anomaly\n",
    "prob_discovery = 5/100 # Probability of ore discovery\n",
    "prop_equity = 50/100   # retained equity in case of joint venture"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87fd021d-e138-4826-8e47-7c641ddd8c36",
   "metadata": {},
   "source": [
    "### Exercise: build a decision tree\n",
    "\n",
    "Add nodes and branches to expand the tree ...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea4292d0-61b2-4039-9437-a191b3e35d34",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree = create_tree(\"0\")\n",
    "add_node_and_branch(tree,node_name=\"1b\",parent_name=\"0\",ntype=\"end\",branch_name=\"relinquish\",probability=0,value_gained=0)\n",
    "add_node_and_branch(tree,node_name=\"1a\",parent_name=\"0\",ntype=\"event\",branch_name=\"geophysics\",probability=1,value_gained=cost_geophy)\n",
    "\n",
    "# CONTINUE TO ADD NODES AND BRANCHES TO REPRODUCE THE DECISION TREE DESCRIBED AT THE BEGINNING\n",
    "\n",
    "print('tree is a tree: '+str(tree.is_tree())+'!')\n",
    "\n",
    "plot_decision_tree(tree)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8753b288-ff8e-48fa-9398-173af4dac73e",
   "metadata": {},
   "source": [
    "### Compute the Expected Value of the project at the different stages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae3bb5bd-956c-4137-b925-bd4e459180bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "res = compute_EV(tree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de64ad29-6cae-418d-9c9d-a097a78aa669",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_decision_tree(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ac5bce-cede-4fd1-96db-ca721239c170",
   "metadata": {},
   "source": [
    "### Retrieve Expected Value at desired node"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67f7bfd-6212-4401-93c4-808539854e76",
   "metadata": {},
   "outputs": [],
   "source": [
    "res.vs(name_eq = '1b')[0][\"EV\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7556b0a-5183-4bc9-ae8d-894dd8841899",
   "metadata": {},
   "source": [
    "### Exercise\n",
    " - what happens when the equity varies between 10 and 90% ?\n",
    " - how does the project expected value varies when the probability to find a geophysical anomaly varies?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b97f4ed-87f4-4de1-b514-76d25fe93256",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d0ac6a9-f8d6-4d61-97c6-5d05c87eb0d8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cad3caed-46d1-4ca0-a0a6-0da48c8a090f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}