Vulnerability and fragility curves for GMMA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vulnerability and fragility curves for GMMA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction for calculating damage and or damage states from fragility curves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start by using the matplotlib magic to allow us to display plots directly in the notebook. You can change which backend is used, but most of the time, you'll want to embed them inside the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now import the required modules to allow plotting, loading data, generating the curves and interacting with them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import lognorm\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "import seaborn\n",
    "seaborn.set_context(\"poster\")\n",
    "\n",
    "# Import widgets for interactive notebook\n",
    "from IPython.html.widgets import interact, fixed, interactive\n",
    "from IPython.html import widgets\n",
    "from IPython.display import display\n",
    "\n",
    "from StringIO import StringIO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the parameters for the vulnerability and fragility curves. These are the parameters provided by UPD-ICE under the GMMA-RAP project. The parameters are for a cumulative lognormal distribution that describes the probability of a building being in a given damage state, or a higher damage state. The damage state descriptions are based on those from the HAZUS-MH Hurricane Model (see Vickery _et al._, 2006). The assumed damage ratio gives an indication of the likely repair cost for the given damage state (as a fraction of the total replacement cost). \n",
    "\n",
    "Rather than trying to load the data from a file (which may not be readily accessible), we create a file-like object to read the data into an array (more for demonstration than out of necessity). \n",
    "\n",
    "The values given here are directly from the UPD-ICE Technical Report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fragilitycurves = StringIO(\"\"\"Type,Mean,SD,Scale,Slight,SD,Scale,Moderate,SD,Scale,Extensive,SD,Scale,Complete,SD,Scale\n",
    "C1-L-S,477.26,0.46,1.00000,209,0.23,1.00000,241,0.2,1.00000,277,0.17,1.00000,371,0.13,1.00000\n",
    "C1-L-W,416.41,0.52,1.00000,132,0.17,1.00000,208,0.16,1.00000,240,0.14,1.00000,301,0.17,1.00000\n",
    "C1-M,221.00,0.33,0.13333,144,0.17,0.133333,171,0.17,0.133333,218,0.22,0.133333,270,0.17,0.133333\n",
    "CHB-L-S,469.00,0.44,1.00000,233,0.33,1.00000,260,0.25,1.00000,290,0.23,1.00000,323,0.17,1.00000\n",
    "CHB-L-W,455.00,0.59,1.00000,164,0.15,1,188,0.19,1,230,0.22,1,324,0.17,1\n",
    "CWS-L-S,398.00,0.16,1.00000,266,0.17,1.00000,311,0.2,1.00000,354,0.09,1.00000,379,0.07,1.00000\n",
    "CWS-L-W,321.00,0.52,1.00000,139,0.34,1,162,0.3,1,218,0.29,1,263,0.31,1\n",
    "MWS-L-S,398.00,0.16,1.00000,266,0.17,1.00000,311,0.2,1.00000,354,0.09,1.00000,379,0.07,1.00000\n",
    "MWS-L-W,321.00,0.52,1.00000,139,0.34,1,162,0.3,1,218,0.29,1,263,0.31,1\n",
    "N-L,136.00,0.20,1.00000,103,0.16,1,123,0.17,1,162,0.26,1,219,0.28,1\n",
    "S1-L,379.10,0.61,1.00000,154,0.35,1,180.45,0.4,1,202.75,0.43,1,315.1,0.46,1\n",
    "S1-M,180.00,0.24,0.13333,140,0.15,0.1333333,152,0.15,0.1333333,188,0.21,0.1333333,205,0.17,0.1333333\n",
    "S3-L,387.00,0.37,1.00000,231,0.24,1,270,0.24,1,322,0.24,1,478,0.36,1\n",
    "W1-L,176.00,0.09,1.00000,143,0.16,1,163,0.17,1,200,0.26,1,259,0.28,1\n",
    "W3-L,95.31,0.25,1.00000,62,0.16,1,82,0.17,1,120,0.26,1,180,0.28,1\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `scale` factor is applied to the loss curve for larger, engineered structures (predominantly high-rise buildings). This is based on the assumption that, for these structures, the only components that can suffer complete damage are the roof and window components. See Appendix A of _Development of Vulnerability Curves for Key Building Types in the Greater Metro Manila Area, Philippines_, UPD-Institute of Civil Engineering (2013). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Building damage states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The damage states give a qualitative indication of the level of damage sustained to a buiding (or components of a building). The damage index quantifies this as a fraction of the complete replacement cost of the structure (component). \n",
    "\n",
    "\n",
    "| Damage state | Description                                |\n",
    "|--------------|:-------------------------------------------|\n",
    "| Negligible   | Little or no visible damage from the outside. No broken windows, or failed roof deck. Minimal loss of roof over, with no or very limited water penetration. |\n",
    "| Slight       | Moderate roof cover that can be covered to prevent additional water ingress. One window, door or garage door broken. |\n",
    "| Moderate     | Major roof damage, moderate window breakage. Minor roof sheathing failure. Some water damage to interior.            |\n",
    "| Extensive    | Major window damage or roof sheathing loss. Major roof cover loss. Extensive damage to interior from water.          |\n",
    "| Complete     | Complete roof failure and/or failure of wall frame. Loss of more than 50% of roof sheathing.                         |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The damage indices are building-specific. In some cases, the damage index for complete damage is < 1.0, which would be counterintuitive. For C1-M and S1-M, only window and roof components were considered, and so the damage index is scaled to reflect the proportional costs of those elements. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Type    | Negligible | Slight  | Moderate | Extensive | Complete |\n",
    "|---------|------------|---------|----------|-----------|----------|\n",
    "| C1-L-S  | 0.0003     | 0.0134  | 0.0981   | 0.185     | 0.4      |\n",
    "| C1-L-W  | 0.0003     | 0.0134  | 0.0981   | 0.185     | 0.4      |\n",
    "| C1-M    | 0.00074    | 0.00659 | 0.03025  | 0.07615   | 0.09975  |\n",
    "| CHB-L-S | 0.0003     | 0.0134  | 0.0981   | 0.185     | 0.4      |\n",
    "| CHB-L-W | 0.0003     | 0.0134  | 0.0981   | 0.185     | 0.4      |\n",
    "| CWS-L-S | 0.000      | 0.01275  | 0.04875   | 0.1125    | 0.4625   |\n",
    "| CWS-L-W | 0.000      | 0.01275  | 0.04875   | 0.1125    | 0.4625   |\n",
    "| MWS-L-S | 0.000      | 0.01275  | 0.04875   | 0.1125    | 0.4625   |\n",
    "| MWS-L-W | 0.000      | 0.01275  | 0.04875   | 0.1125    | 0.4625   |\n",
    "| N-L     | 0.003      | 0.031   | 0.12     | 0.283     | 1.000    |\n",
    "| S1-L    | 0.000      | 0.01254 | 0.067375 | 0.165     | 0.55     |\n",
    "| S1-M    | 0.000      | 0.011   | 0.049    | 0.089     | 0.133    |\n",
    "| S3-L    | 0.0003     | 0.0161  | 0.1169   | 0.24      | 0.55     |\n",
    "| W1-L    | 0.003      | 0.031   | 0.12     | 0.283     | 1.000    |\n",
    "| W3-L    | 0.003      | 0.031   | 0.12     | 0.283     | 1.000    |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "damageindices = StringIO(\"\"\"Type,Negligible,Slight,Moderate,Extensive,Complete\n",
    "C1-L-S,0.0003,0.0134,0.0981,0.185,0.4\n",
    "C1-L-W,0.0003,0.0134,0.0981,0.185,0.4\n",
    "C1-M,0.00074,0.00659,0.03025,0.07615,0.09975\n",
    "CHB-L-S,0.0003,0.0134,0.0981,0.185,0.4\n",
    "CHB-L-W,0.0003,0.0134,0.0981,0.185,0.4\n",
    "CWS-L-S,0,0.01275,0.04875,0.1125,0.4625\n",
    "CWS-L-W,0,0.01275,0.04875,0.1125,0.4625\n",
    "MWS-L-S,0,0.01275,0.04875,0.1125,0.4625\n",
    "MWS-L-W,0,0.01275,0.04875,0.1125,0.4625\n",
    "N-L,0.003,0.031,0.12,0.283,1\n",
    "S1-L,0,0.01254,0.067375,0.165,0.55\n",
    "S1-M,0,0.011,0.049,0.089,0.133\n",
    "S3-L,0.0003,0.0161,0.1169,0.24,0.55\n",
    "W1-L,0.003,0.031,0.12,0.283,1\n",
    "W3-L,0.003,0.031,0.12,0.283,1\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to load the data that describes the damage and fragility functions. This is stored in a CSV format (above), so we'll use [numpy.genfromtxt](http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html) to read it in and create a 2-d array that can be used to plot up the information. \n",
    "\n",
    "We need to first define the structure of the data, since the present format contains duplicate field names in the file header. Once that's done, we can then load it up. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "FCOLNAMES = ['BuildingType', 'Vuln_mu', 'Vuln_sd', 'Vuln_scale',\n",
    "            'S_mu','S_sd','S_scale',\n",
    "            'M_mu','M_sd','M_scale',\n",
    "            'E_mu','E_sd','E_scale',\n",
    "            'C_mu','C_sd','C_scale']\n",
    "\n",
    "FCOLTYPES = ['|S7', 'f8', 'f8', 'f8',\n",
    "            'f8', 'f8', 'f8',\n",
    "            'f8', 'f8', 'f8',\n",
    "            'f8', 'f8', 'f8',\n",
    "            'f8', 'f8', 'f8']\n",
    "\n",
    "DCOLNAMES = [\"BuildingType\", \"Negligible\", \"Slight\", \n",
    "             \"Moderate\", \"Extensive\", \"Complete\"]\n",
    "DCOLTYPES = [\"|S7\", \"f8\", \"f8\", \"f8\", \"f8\", \"f8\"]\n",
    "\n",
    "fdtype = np.dtype({'names':FCOLNAMES, 'formats':FCOLTYPES})\n",
    "ddtype = np.dtype({'names':DCOLNAMES, 'formats':DCOLTYPES})\n",
    "\n",
    "delimiter = ','\n",
    "skip_header = 1\n",
    "usecols = tuple(range(16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice here I'm using the previously defined `StringIO` object as a file object. Saves the hassle of trying to select a file (which requires importing additional notebooks). This means that whenever you run the code snippet below, you **must** run the code snippet that defines the `StringIO` object above. Reading the `StringIO` object puts the pointer at the end of the buffer, so it'll raise a `StopIteration` exception. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "curvedata = np.genfromtxt(fragilitycurves, dtype=fdtype, delimiter=delimiter, \n",
    "                     skip_header=skip_header, usecols=usecols, \n",
    "                     autostrip=True)\n",
    "damagedata = np.genfromtxt(damageindices, dtype=ddtype, delimiter=delimiter,\n",
    "                          skip_header=skip_header, usecols=tuple(range(6)),\n",
    "                          autostrip=True)\n",
    "\n",
    "bldtypes = curvedata['BuildingType'].tolist()\n",
    "nsim = 100000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Store the parameters of the damage and fragility curves in a dict, keyed by the name of the building type. \n",
    "\n",
    "The parameters define the cumulative lognormal distribution:\n",
    "\n",
    "$F_V(v;\\mu,\\sigma) = \\Phi\\bigg(\\frac{\\ln v - \\mu}{\\sigma}\\bigg)$\n",
    "\n",
    "Note that [`scipy.stats`](https://docs.scipy.org/doc/scipy-0.15.1/reference/stats.html) uses fixed names for the parameters to all the continuous distributions. This means the mean of the distribution ($\\mu$) defines the `scale` parameter, and the standard deviation ($\\sigma$) defines the `shape` parameter. \n",
    "\n",
    "(Precisely, if $log(x)$ is normally distributed with mean $\\mu$ and variance $\\sigma^2$, then $x$ is log-normally distributed with `shape` parameter $\\sigma$ and `scale` parameter $\\exp(\\mu)$ )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "curves = dict()\n",
    "damageindex = dict()\n",
    "for i, t in enumerate(bldtypes):\n",
    "    v = dict()\n",
    "    d = dict()\n",
    "    v['vloc'] = curvedata['Vuln_mu'][i]\n",
    "    v['vshape'] = curvedata['Vuln_sd'][i]\n",
    "    v['vscale'] = curvedata['Vuln_scale'][i]\n",
    "    \n",
    "    v['sloc'] = curvedata['S_mu'][i]\n",
    "    v['sshape'] = curvedata['S_sd'][i]\n",
    "    v['sscale'] = curvedata['S_scale'][i]\n",
    "    \n",
    "    v['mloc'] = curvedata['M_mu'][i]\n",
    "    v['mshape'] = curvedata['M_sd'][i]\n",
    "    v['mscale'] = curvedata['M_scale'][i]  \n",
    "    \n",
    "    v['eloc'] = curvedata['E_mu'][i]\n",
    "    v['eshape'] = curvedata['E_sd'][i]\n",
    "    v['escale'] = curvedata['E_scale'][i]\n",
    "    \n",
    "    v['cloc'] = curvedata['C_mu'][i]\n",
    "    v['cshape'] = curvedata['C_sd'][i]\n",
    "    v['cscale'] = curvedata['C_scale'][i]\n",
    "\n",
    "    curves[t] = v\n",
    "    damageindex[t] = np.array([damagedata['Negligible'][i], damagedata['Slight'][i], damagedata['Moderate'][i],\n",
    "                              damagedata['Extensive'][i],damagedata['Complete'][i]])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fitting the derived vulnerablity model\n",
    "\n",
    "We use [least squares fitting](https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.optimize.leastsq.html#scipy.optimize.leastsq) to fit the parameters of the derived loss curve. These functions define the calculation of residuals, then execute the least squares minimisation to determine the shape & location parameters for the cumulative lognormal distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def residuals(params, data, scale):\n",
    "    x = np.arange(0., 361, 1.)\n",
    "    return data - scale * lognorm.cdf(x, params[0], scale=params[1])\n",
    "\n",
    "from scipy.optimize import leastsq\n",
    "def minimise(data, scale, shape, loc):\n",
    "    plsq = leastsq(residuals, [shape, loc], args=(data, scale))\n",
    "    return plsq[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte-Carlo sampling to calculate damage ratio\n",
    "\n",
    "Define a simple function to calculate the losses using monte carlo simulation. We create $10^5$ random samples of damage states, and assign a cost to each based on the assumed damage ratios defined above (from HAZUS-MH)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def probloss(cdfs, di):\n",
    "    \"\"\"\n",
    "    Calculate a damage index for a population of buildings, given the probability of a \n",
    "    set of damage states (`cdfs`) and the damage index (`di`) associated with each damage\n",
    "    state.\n",
    "    \n",
    "    Randomly sample from the distribution of different damage states and assign a cost\n",
    "    to each, based on the damage indices. \n",
    "    \n",
    "    :param cdfs: `numpy.ndarray` of cumulative distributions for the damage states\n",
    "    :param di: `numpy.ndarray` of damage indices for the damage states.\n",
    "    \"\"\"\n",
    "    \n",
    "    u = np.random.uniform(0., 1., nsim)\n",
    "    s, m, e, c = cdfs\n",
    "    ndi, sdi, mdi, edi, cdi = di\n",
    "    loss = np.ones(nsim) * ndi\n",
    "    np.putmask(loss, ((u<=s) & (u>m)), sdi)\n",
    "    np.putmask(loss, ((u<=m) & (u>e)), mdi)\n",
    "    np.putmask(loss, ((u<=e) & (u>c)), edi)\n",
    "    np.putmask(loss, u<=c, cdi)\n",
    "    \n",
    "    return loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function uses the method as reported by UPD-ICE to calculate the damage ratio. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def probloss_upd(cdfs, di):\n",
    "    \"\"\"\n",
    "    Determine the vulnerability as described in the UPD-ICE report for GMMA RAP. \n",
    "    We are attempting to replicate the curves reported by UPD-ICE. \n",
    "    \n",
    "    In this method, a damage multiplier is multiplied to the probability of exceedance\n",
    "    to arrive at the damage for a certain damage state. For each wind speed, the damage \n",
    "    index is computed by summing the damage for all damage states. \n",
    "    \n",
    "    :param cdfs: `numpy.ndarray` of cumulative distribution values for damage states.\n",
    "    :param di: `numpy.ndarray` of damage indices for damage states (including negligible damage).\n",
    "    \n",
    "    :returns: Damage ratio (should be between 0 and 1).\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    n = 1. \n",
    "    a = np.concatenate([[n], cdfs])\n",
    "    loss = np.dot(a, di)\n",
    "    \n",
    "    return loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to loop over all wind speeds and evaluate the damage ratio, given the probabilities of damage states and corresponding damage indices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def damage(v, curves, dmgindex):\n",
    "    \n",
    "    v = np.arange(0, 361., 1.)\n",
    "    dmg_mean = np.zeros(len(v))\n",
    "    dmg_sd = np.zeros(len(v))\n",
    "    \n",
    "    for i, x in enumerate(v):\n",
    "        scdf = lognorm.cdf(x, curves['sshape'], scale=curves['sloc'])\n",
    "        mcdf = lognorm.cdf(x, curves['mshape'], scale=curves['mloc'])\n",
    "        ecdf = lognorm.cdf(x, curves['eshape'], scale=curves['eloc'])\n",
    "        ccdf = lognorm.cdf(x, curves['cshape'], scale=curves['cloc'])\n",
    "        probs = [scdf, mcdf, ecdf, ccdf]\n",
    "        loss = probloss(probs, dmgindex)\n",
    "        dmg_mean[i] = np.mean(loss)\n",
    "        dmg_sd[i] = np.std(loss)\n",
    "        \n",
    "    return dmg_mean, dmg_sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def damage_upd(v, curves, dmgindex):\n",
    "    \n",
    "    v = np.arange(0, 361., 1.)\n",
    "    dmg_upd = np.zeros(len(v))\n",
    "    for i, x in enumerate(v):\n",
    "        scdf = lognorm.cdf(x, curves['sshape'], scale=curves['sloc'])\n",
    "        mcdf = lognorm.cdf(x, curves['mshape'], scale=curves['mloc'])\n",
    "        ecdf = lognorm.cdf(x, curves['eshape'], scale=curves['eloc'])\n",
    "        ccdf = lognorm.cdf(x, curves['cshape'], scale=curves['cloc'])\n",
    "        probs = [scdf, mcdf, ecdf, ccdf]\n",
    "        dmg_upd[i] = probloss_upd(probs, dmgindex)\n",
    "        \n",
    "    return dmg_upd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plotdamage(ax, curves, windspeed, dmgindex):\n",
    "    v = np.arange(0, 361., 1.)\n",
    "    \n",
    "    dmg_mean, dmg_sd = damage(v, curves, dmgindex)\n",
    "    \n",
    "    dmg_lower = np.max(np.array([np.zeros(len(v)), dmg_mean - dmg_sd]),axis=0)\n",
    "    dmg_upper = np.min(np.array([np.ones(len(v)), dmg_mean + dmg_sd]),axis=0)\n",
    "    widx = np.where(v==windspeed)[0]\n",
    "    dmg_mean_ws = dmg_mean[widx]\n",
    "    dmg_sd_ws = dmg_sd[widx]\n",
    "    ax.plot(v, dmg_mean, 'b', linewidth=2., \n",
    "             label=\"Damage ratio (derived): %.3f\"%dmg_mean_ws)\n",
    "    ax.fill_between(v, dmg_lower, dmg_upper, facecolor='b', alpha=0.2)   \n",
    "       \n",
    "    ax.set_xlim((0,360))\n",
    "    ax.set_xticks(np.arange(0, 361., 30))\n",
    "    ax.set_ylim((0,1))\n",
    "    ax.set_ylabel('Damage ratio')\n",
    "    ax.axvline(windspeed, color='k',linestyle=':')\n",
    "    \n",
    "    if dmg_mean_ws < 0.95:\n",
    "        ax.text(windspeed+2, dmg_mean_ws+0.05, \"%.3f\"%dmg_mean_ws, \n",
    "                 fontsize='x-small', color='b', ha='left', va='center')\n",
    "    else:\n",
    "        ax.text(windspeed+2, dmg_mean_ws-0.05, \"%.3f\"%dmg_mean_ws, \n",
    "                 fontsize='x-small', color='b', \n",
    "                 ha='left', va='center')\n",
    "        \n",
    "    xx = np.where(dmg_mean<=0.5)[0][-1]\n",
    "    loc = v[xx]\n",
    "    \n",
    "    p = minimise(dmg_mean, curves['vscale'], curves['vshape'], loc)\n",
    "    pval = lognorm.cdf(windspeed, p[0], scale=p[1])\n",
    "    \n",
    "    ax.plot(v, curves['vscale']*lognorm.cdf(v, p[0], scale=p[1]), 'r', \n",
    "             linewidth=2., linestyle=':', label=\"Fitted model: %.3f\"%pval)\n",
    "    \n",
    "    fitstr = \"\"\"Fitted parameters:\n",
    "    $\\sigma$ : %.4f\n",
    "    $\\mu$ : %.1f\"\"\"%(p[0], p[1])\n",
    "    ax.text(15, 0.45, fitstr, \n",
    "             bbox=dict(facecolor='white',alpha=0.5),\n",
    "             ha='left', va='top', color='r')\n",
    "\n",
    "    return dmg_mean, dmg_sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plotdamage_upd(ax, curves, windspeed, dmgindex):\n",
    "    v = np.arange(0, 361., 1.)\n",
    "    vcdf = curves['vscale'] * lognorm.cdf(v, curves['vshape'], \n",
    "                                          scale=curves['vloc'])\n",
    "    widx = np.where(v==windspeed)[0]\n",
    "    ax.plot(v, vcdf, '0.5', linewidth=2., label=\"Damage model (UPD-ICE): %.3f\"%vcdf[widx])\n",
    "    \n",
    "    #dmg = damage_upd(v, curves, dmgindex)\n",
    "    #ax.plot(v, dmg, linewidth=2., color='g', label=\"Damage model (UPD-ICE) - cumulative: %.3f\"%dmg[widx])\n",
    "    \n",
    "    l = ax.legend(loc=2,frameon=True)\n",
    "    \n",
    "    param_str = \"\"\"UPD-ICE parameters:\n",
    "    $\\sigma$ : %.4f\n",
    "    $\\mu$ : %.1f\"\"\"%(curves['vshape'], curves['vloc'])\n",
    "    \n",
    "    ax.text(15, 0.25, param_str, \n",
    "             bbox=dict(facecolor='white',alpha=0.5),\n",
    "             ha='left', va='top', color='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function calculates a histogram, based on the likely distribution of buildings being in different damage states. We take a random sample of 100,000 'buildings' and assign a damage state based on the fragility curves for the given wind speed $v$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def problosshist(ax, curves, v, dmgindex):\n",
    "    \"\"\"\n",
    "    Plot a histogram of the distribution of damage states. \n",
    "    \n",
    "    \"\"\"\n",
    "    scdf = lognorm.cdf(v, curves['sshape'], scale=curves['sloc'])\n",
    "    mcdf = lognorm.cdf(v, curves['mshape'], scale=curves['mloc'])\n",
    "    ecdf = lognorm.cdf(v, curves['eshape'], scale=curves['eloc'])\n",
    "    ccdf = lognorm.cdf(v, curves['cshape'], scale=curves['cloc'])\n",
    "    probs = [scdf, mcdf, ecdf, ccdf]\n",
    "    loss = probloss(probs, dmgindex)\n",
    "    bins = np.concatenate([dmgindex, [1.1]])\n",
    "    losshist, bins = np.histogram(loss, bins=bins)\n",
    "    ax.bar(np.arange(5), losshist, align='center')\n",
    "    ax.set_xlim((-0.5, 4.5))\n",
    "    \n",
    "    axyl = ax.get_yticklabels()\n",
    "    ax.set_xticklabels(['','Negligible', 'Slight', 'Moderate', 'Extensive', 'Complete'],\n",
    "                       ha='center')\n",
    "    ax.set_ylim((0,nsim))\n",
    "    for t in axyl:\n",
    "        t.set_rotation('vertical')\n",
    "        \n",
    "    ax.set_title(\"Distribution of damage state\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a function to plot the fragility curves for a selected building type, and alongside it the corresponding damage curve. \n",
    "\n",
    "We add the loss curve as defined by UPD-ICE, a derived loss curve, based on monte-carlo sampling of the fragility, and a fitted curve, based on the derived curve. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plotcurve(buildingtype, windspeed, seed):\n",
    "    np.random.seed(seed)\n",
    "    bldcurves = curves[buildingtype]\n",
    "    dmgindex = damageindex[buildingtype]\n",
    "    v = np.arange(0, 361., 1.)\n",
    "    \n",
    "    vcdf = bldcurves['vscale'] * lognorm.cdf(v, bldcurves['vshape'], \n",
    "                                              scale=bldcurves['vloc'])\n",
    "    \n",
    "    scdf = lognorm.cdf(v, bldcurves['sshape'], scale=bldcurves['sloc'])\n",
    "    mcdf = lognorm.cdf(v, bldcurves['mshape'], scale=bldcurves['mloc'])\n",
    "    ecdf = lognorm.cdf(v, bldcurves['eshape'], scale=bldcurves['eloc'])\n",
    "    ccdf = lognorm.cdf(v, bldcurves['cshape'], scale=bldcurves['cloc'])\n",
    "\n",
    "    # Set up the figure canvas:\n",
    "    fig = plt.figure(figsize=(16,8))\n",
    "    gs = gridspec.GridSpec(2,6)\n",
    "    \n",
    "    ax1 = plt.subplot(gs[0:1,0:4])\n",
    "    ax2 = plt.subplot(gs[-1,0:4])\n",
    "    ax3 = plt.subplot(gs[0:2,4:])\n",
    "\n",
    "    # Plot the fragility curves:\n",
    "    ax1.fill_between(v, 0, ccdf, facecolor='k', alpha=0.25)\n",
    "    ax1.plot(v, ccdf, 'k', linewidth=1.)\n",
    "    ax1.fill_between(v, ccdf, ecdf, facecolor='r', alpha=0.25)\n",
    "    ax1.plot(v, ecdf, 'r', linewidth=1.)\n",
    "    ax1.fill_between(v, ecdf, mcdf, facecolor='b', alpha=0.25)\n",
    "    ax1.plot(v, mcdf, 'b', linewidth=1.)\n",
    "    ax1.fill_between(v, mcdf, scdf, facecolor='g', alpha=0.25)\n",
    "    ax1.plot(v, scdf, 'g', linewidth=1.)\n",
    "    \n",
    "    # Add annotations:\n",
    "    bbox = dict(facecolor='white', alpha=0.5)\n",
    "    cypos = lognorm.cdf(windspeed, bldcurves['cshape'], scale=bldcurves['cloc']) \n",
    "    ax1.plot([windspeed, windspeed], [0,cypos], color='k', \n",
    "             label=\"Complete: %.3f\"%cypos)\n",
    "    \n",
    "    if cypos > 0.05:\n",
    "        ax1.text(windspeed+2, cypos/2., \"%.3f\"%cypos, fontsize='small',\n",
    "                 color='k', ha='left', va='center', bbox=bbox)\n",
    "        \n",
    "    eypos = lognorm.cdf(windspeed, bldcurves['eshape'], scale=bldcurves['eloc']) \n",
    "    ax1.plot([windspeed, windspeed], [cypos,eypos], color='r', \n",
    "             label=\"Extensive: %.3f\"%(eypos - cypos))\n",
    "    \n",
    "    if (eypos - cypos) > 0.05:\n",
    "        ax1.text(windspeed+2, (cypos+eypos)/2., \"%.3f\"%(eypos), \n",
    "                 fontsize='small',color='r', ha='left', va='center', bbox=bbox)\n",
    "        \n",
    "    mypos = lognorm.cdf(windspeed, bldcurves['mshape'], scale=bldcurves['mloc']) \n",
    "    ax1.plot([windspeed, windspeed], [eypos,mypos], color='b', \n",
    "             label=\"Moderate: %.3f\"%(mypos - eypos))\n",
    "    if (mypos - eypos) > 0.05:\n",
    "        ax1.text(windspeed+2, (mypos+eypos)/2., \"%.3f\"%(mypos), \n",
    "                 fontsize='small',color='b', ha='left', va='center', bbox=bbox)\n",
    "        \n",
    "    sypos = lognorm.cdf(windspeed, bldcurves['sshape'], scale=bldcurves['sloc']) \n",
    "    ax1.plot([windspeed, windspeed], [mypos,sypos], color='g', \n",
    "             label=\"Slight: %.3f\"%(sypos - mypos))\n",
    "    if (sypos - mypos) > 0.05:\n",
    "        ax1.text(windspeed+2, (sypos+mypos)/2., \"%.3f\"%(sypos), \n",
    "                 fontsize='small',color='g', ha='left', va='center', bbox=bbox)\n",
    "    \n",
    "    ax1.plot([windspeed, windspeed], [sypos,1.0], color='orange', \n",
    "             label=\"No damage: %.3f\"%(1.-sypos))\n",
    "    if sypos<0.95:\n",
    "        ax1.text(windspeed+2, (sypos+1.)/2., \"%.3f\"%(1.-sypos), \n",
    "                 fontsize='small',color='brown', ha='left', va='center', bbox=bbox)\n",
    "    \n",
    "    # Arrange the axis limits, labels, etc.\n",
    "    ax1.set_xlim((0,360))\n",
    "    ax1.set_xticks(np.arange(0, 361., 30))\n",
    "    ax1.set_ylim((0,1))\n",
    "    ax1.set_ylabel(r\"$ P(D \\geq ds) $\")\n",
    "    ax1.set_xlabel(\"Wind speed (km/h)\")\n",
    "    ax1.axvline(windspeed, color='k',linestyle=':')\n",
    "    ax1.grid(True)\n",
    "    \n",
    "    l = ax1.legend(loc=2,frameon=True,)\n",
    "        \n",
    "    ax1.set_title(\"Probability of damage states\")\n",
    "    \n",
    "    # Add the other plots:\n",
    "    problosshist(ax3, bldcurves, windspeed, dmgindex)\n",
    "    plotdamage(ax2, bldcurves, windspeed, dmgindex)\n",
    "    plotdamage_upd(ax2, bldcurves, windspeed, dmgindex)\n",
    "    \n",
    "    ax2.grid(True)\n",
    "    ax2.set_title(\"Building damage ratio - %s\"%buildingtype)\n",
    "    fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the curves for a selected building type. Users can vary the wind speed to observe the probability of that building type being in one of the damage states. \n",
    "\n",
    "The labels in the legend indicate the probability of the building type being in the given damage state. The text annotations indicate the probability of exceeding the damage state - that is, the building type is in that damage state, or a higher damage state. \n",
    "\n",
    "On the right is a bar chart counting the number of instances (out of a sample of 100,000) of each damage state. It's effectively representing the labels in the legend in a graphical manner. \n",
    "\n",
    "In the lower panel, the derived vulnerability curve (blue), the published vulnerability curve (grey) and one based on the reported method of calculating the damage ratio (green). The derived curve includes the +/- 1 standard deviation range (blue band), plus a fitted lognormal cumulative distribution (red), based on the derived curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "seed = 10574\n",
    "\n",
    "# Drop down menu to select the building type:\n",
    "bldselect = widgets.Dropdown(options=bldtypes, value=bldtypes[0],\n",
    "                             description=\"Building type\")\n",
    "\n",
    "# Slider widget to set the wind speed:\n",
    "windselect = widgets.FloatSlider(min=0.0, max=360., step=1., value=180., \n",
    "                                 description=\"Wind speed (km/h)\" )\n",
    "\n",
    "w = interactive(plotcurve, \n",
    "                buildingtype=bldselect, \n",
    "                windspeed=windselect, \n",
    "                seed=fixed(seed))\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Full suite of vulnerability curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we cycle through all available building types and fit a lognormal cdf to the derived damage curve. These can then be mapped to the corresponding building types in the full (>100) suite of building type & height combinations used in GMMA RAP to evaluate impacts to the built environment. The shape parameter should correspond to the 'sd' field in the suite of curves, and the scale corresponds to the 'mean' field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "v = np.arange(0, 361., 1.)    \n",
    "newcurves = dict()\n",
    "print \"Building type: shape  : scale  : UPD-ICE shape : UPD-ICE scale\"\n",
    "print \"-------------------------------------------------------------- \"\n",
    "for bldtype in bldtypes:\n",
    "\n",
    "    bldcurves = curves[bldtype]\n",
    "    di = damageindex[bldtype]\n",
    "    vcdf = bldcurves['vscale'] * lognorm.cdf(v, bldcurves['vshape'], \n",
    "                                              scale=bldcurves['vloc'])\n",
    "    dmg_mean, dmg_sd = damage( v, bldcurves, di)\n",
    "    xx = np.where(dmg_mean<=0.5)[0][-1]\n",
    "    loc = v[xx]\n",
    "    p = minimise(dmg_mean, bldcurves['vscale'], bldcurves['vshape'], loc)\n",
    "    print \"{0:<13s}: {1:.5f}: {2:<7.2f}: {3:14.2f}: {4:<10.2f}\".format(bldtype, p[0], p[1], \n",
    "                                                                       bldcurves['vshape'], bldcurves['vloc'])\n",
    "    newcurves[bldtype] = {'vshape':p[0], 'vloc':p[1]}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll plot all the curves together to see how the derived curves compare with the reported curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(len(bldtypes)/3, 3,sharex=True,sharey=True,figsize=(18,15))\n",
    "x = np.arange(0, 360., 1.)\n",
    "\n",
    "for ax, bldtype in zip(axes.ravel(), bldtypes):\n",
    "    bldcurves = curves[bldtype]\n",
    "    ncurves = newcurves[bldtype]\n",
    "    vcdf = bldcurves['vscale'] * lognorm.cdf(x, bldcurves['vshape'], scale=bldcurves['vloc'])\n",
    "    ncdf = bldcurves['vscale'] * lognorm.cdf(x, ncurves['vshape'], scale=ncurves['vloc'])\n",
    "    ax.plot(x, ncdf, linewidth=2, label=\"Derived\")\n",
    "    ax.plot(x, vcdf, linewidth=1.5, linestyle='--',label=\"UPD-ICE\")\n",
    "    ax.set_ylabel(bldtype)\n",
    "    \n",
    "    ax.set_xlim((0,360))\n",
    "    ax.set_xticks(np.arange(0, 361., 60))\n",
    "    ax.set_ylim((0,1))\n",
    "\n",
    "axes.ravel()[-2].set_xlabel(\"Wind speed (km/h)\")\n",
    "l = ax.legend(loc=0, frameon=True,)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(len(bldtypes)/3, 3,sharex=True,sharey=True,figsize=(18,15))\n",
    "x = np.arange(0, 360., 1.)\n",
    "seaborn.set_context(\"poster\")\n",
    "for ax, bldtype in zip(axes.ravel(), bldtypes):\n",
    "    bldcurves = curves[bldtype]\n",
    "    ncurves = newcurves[bldtype]\n",
    "    vcdf = bldcurves['vscale'] * lognorm.cdf(x, bldcurves['vshape'], scale=bldcurves['vloc'])\n",
    "    ncdf = bldcurves['vscale'] * lognorm.cdf(x, ncurves['vshape'], scale=ncurves['vloc'])\n",
    "    #ax.plot(x, ncdf, linewidth=2, label=\"Derived\")\n",
    "    ax.plot(x, vcdf, linewidth=2, label=\"UPD-ICE\")\n",
    "    ax.text(15, 0.8, bldtype,fontsize='xx-large')\n",
    "    \n",
    "    ax.set_xlim((0,360))\n",
    "    ax.set_xticks(np.arange(0, 361., 60))\n",
    "    ax.set_ylim((0,1))\n",
    "\n",
    "axes.ravel()[-2].set_xlabel(\"Wind speed (km/h)\")\n",
    "axes.ravel()[-9].set_ylabel(\"Damage index\")\n",
    "#l = ax.legend(loc=4, frameon=True,)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(18,12))\n",
    "x = np.arange(0, 360., 1.)\n",
    "seaborn.set_context(\"poster\")\n",
    "for bldtype in bldtypes:\n",
    "    bldcurves = curves[bldtype]\n",
    "    ncurves = newcurves[bldtype]\n",
    "    vcdf = bldcurves['vscale'] * lognorm.cdf(x, bldcurves['vshape'], scale=bldcurves['vloc'])\n",
    "    ncdf = bldcurves['vscale'] * lognorm.cdf(x, ncurves['vshape'], scale=ncurves['vloc'])\n",
    "    #ax.plot(x, ncdf, linewidth=2, label=\"Derived\")\n",
    "    ax.plot(x, vcdf, linewidth=2, label=bldtype)\n",
    "    #ax.text(15, 0.8, bldtype,fontsize='xx-large')\n",
    "    \n",
    "    ax.set_xlim((0,360))\n",
    "    ax.set_xticks(np.arange(0, 361., 60))\n",
    "    ax.set_ylim((0,1))\n",
    "\n",
    "ax.set_xlabel(\"Wind speed (km/h)\")\n",
    "ax.set_ylabel(\"Damage index\")\n",
    "l = ax.legend(loc=2, frameon=True,ncol=3)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the derived vulnerability models closely match the reported models. Curves where there is little damage below 240 km/h result in a poorer match, largely because there are fewer non-zero data points to fit against. There may be opportunity to review the underlying fragility curves for these building types to consider the level of damage sustained at high wind speeds. \n",
    "\n",
    "Also there are some building types where the damage index assigned to a complete damage state is not 1.0. Excluding those building types where only damage to windows and roofs are considered (C1-M and S1-M), this will lead to situations where the complete damage state is reached throughout a population of buildings, but the damage index is not 1.0. \n",
    "\n",
    "The W3-L, W1-L and N-L buildings are significantly different from the reported curves. This is because the derived curves were considered too strong and replaced with empirical curves based on discussion with local builders & carpenters. It could be recommended to fit fragility curves to these empirical damage curves (e.g. usnig Geoscience Australia's eLoss tool). Also, the fragility curves for W3-L were not reported by UPD-ICE and were inferred from the N-L building type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "This document was written in an IPython notebook. The raw notebook can be downloaded [here](https://github.com/wcarthur/notebooks/blob/master/Vulnerability%20and%20fragility%20curves%20for%20GMMA.ipynb). See also [nbviewer](http://nbviewer.ipython.org/github/wcarthur/notebooks/blob/master/Vulnerability%20and%20fragility%20curves%20for%20GMMA.ipynb) for an online static view."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}