TC YASI Vulnerabilty and Fragility Plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vulnerability Curves\n",
    "\n",
    "Script to create Damage Index vs Local Wind Speed best-fit curves with lognormal cdf fit, based on least square optimisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allow plots to be shown inline in notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import ma\n",
    "\n",
    "import csv\n",
    "\n",
    "from pylab import *\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import lognorm\n",
    "from scipy.stats import weibull_min\n",
    "from scipy.optimize import leastsq\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "#import matplotlib.gridspec as gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enter path of csv file containing [UNIQUE FIELD IDENTIFIER, BUILDING TYPE, LOCAL WIND SPEED, DAMAGE INDEX, HAZIMP DAMAGE INDEX*, HAZIMP WIND SPEED (0.2s gust)*] data:\n",
    "\n",
    "*Enter columns of 0 if HAZIMP not used"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Type | Roof            | Wall                        | Wall cladding        |            | Age       |\n",
    "|------|-----------------|-----------------------------|----------------------|------------|-----------|\n",
    "| 1    | Metal sheet     | Timber frame                | Fibro                | High set   | 1918-1946 |\n",
    "| 2    | Metal sheet     | Timber frame                | Fibro                | Low set    | 1960-1970 |\n",
    "| 3    | Metal sheet     | Concrete block/timber frame | Concrete block/fibro | Two-storey | 1970-1980 |\n",
    "| 4    | Metal sheet     | Brick veneer                | Brick                | Low set    | 1980-     |\n",
    "| 5    | Concrete tile   | Brick veneer                | Brick                | Low set    | 1980-     |\n",
    "| 6    | Industrial shed | Transverse portal frame     | Sheet metal          | NA         | NA        |\n",
    "| 7    | Metal sheeting  | Timber framing              | Cement sheeting      | High set   | 1960-1975 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "file_path = 'C:/WorkSpace/data/Derived/impact/yasi/input/damage_index_data.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import file and process data.\n",
    "\n",
    "**NOTE**: Appears incorrect conversion beween 0.2-sec and 3-sec gust wind speeds is being used. Confirm what value is contained in the `wind_haz` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "file = open(file_path)\n",
    "data = csv.reader(file)\n",
    "\n",
    "wind_pts = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "DI_pts = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "wind_haz = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "DI_haz = [[],[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "\n",
    "for row in data:\n",
    "    for i in range(10):\n",
    "        building_type = 'TYPE %i' %(i+1)\n",
    "        if row[1] == building_type and float(row[2]) > 10:\n",
    "            wind_pts[i].append(float(row[2]))\n",
    "            DI_pts[i].append(float(row[3]))\n",
    "            if i in [0,1,2,3,4,6,7,8,9] and row[5] != '0':\n",
    "                DI_haz[i].append(float(row[4]))\n",
    "                wind_haz[i].append(float(row[5])/1.15) # 0.2s to 3s gust conversion (1.132+1.164)/2 \n",
    "                                                       # from Holomes & Ginger 2012 Table 4\n",
    "            \n",
    "\n",
    "            if i in [3,4]:\n",
    "                wind_pts[10].append(float(row[2]))\n",
    "                DI_pts[10].append(float(row[3]))\n",
    "                DI_haz[10].append(float(row[4]))\n",
    "                wind_haz[10].append(float(row[5])/1.15)\n",
    "            elif i in [0,1,2,6]:\n",
    "                wind_pts[11].append(float(row[2]))\n",
    "                DI_pts[11].append(float(row[3]))\n",
    "                DI_haz[11].append(float(row[4]))\n",
    "                wind_haz[11].append(float(row[5]))\n",
    "            elif i == 5:\n",
    "                wind_pts[12].append(float(row[2]))\n",
    "                DI_pts[12].append(float(row[3]))\n",
    "                DI_haz[12].append(float(row[4]))\n",
    "                wind_haz[12].append(float(row[5])/1.15)\n",
    "            \n",
    "            elif i in [8,9]:\n",
    "                wind_pts[13].append(float(row[2]))\n",
    "                DI_pts[13].append(float(row[3]))\n",
    "                DI_haz[13].append(float(row[4]))\n",
    "                wind_haz[13].append(float(row[5])/1.15)            \n",
    "            elif i == 7:\n",
    "                wind_pts[14].append(float(row[2]))\n",
    "                DI_pts[14].append(float(row[3]))\n",
    "                DI_haz[14].append(float(row[4]))\n",
    "                wind_haz[14].append(float(row[5])/1.15)\n",
    "\n",
    "\n",
    "wind_pts[15] = wind_pts[10] + wind_pts[13]\n",
    "DI_pts[15] = DI_pts[10] + DI_pts[13]\n",
    "wind_pts[16] = wind_pts[11] + wind_pts[14]\n",
    "DI_pts[16] = DI_pts[11] + DI_pts[14]\n",
    "\n",
    "\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define minimise and residuals functions which will be used to calculate best fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def residuals(params, build_type):\n",
    "    x = wind_pts[build_type]\n",
    "    y = DI_pts[build_type]\n",
    "    output = y - lognorm.cdf(x, params[0], scale = params[1])\n",
    "    return output\n",
    "\n",
    "def minimise(params, build_type):\n",
    "    output = leastsq(residuals, params, args = build_type, full_output = 1, ftol=1E-10, xtol=1E-10, maxfev = 100000)\n",
    "    return output\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def function_exp(params,x):\n",
    "    alpha_inv = 1/params[0]\n",
    "    beta_exp = np.exp(params[1])\n",
    "    function_exp = np.exp(-(x/beta_exp)**(alpha_inv))\n",
    "    output = 1 - function_exp\n",
    "    return output\n",
    "\n",
    "def residuals_exp(params, build_type):\n",
    "    x = wind_pts[build_type]\n",
    "    y = DI_pts[build_type]\n",
    "    output = y - function_exp(params,x)\n",
    "    return output\n",
    "\n",
    "def minimise_exp(params, build_type):\n",
    "    output = leastsq(residuals_exp, params, args = build_type, full_output = 1, ftol=1E-10, xtol=1E-10, maxfev = 100000)\n",
    "    return output\n",
    "\n",
    "\n",
    "\n",
    "def frange(x, y, jump):\n",
    "    output = []\n",
    "    val = x\n",
    "    while val<=y:\n",
    "        output.append(val)\n",
    "        val = round(val+jump,1)\n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Carry out minimisation process and create plots for each building type.\n",
    "\n",
    "Save plots and plotting info log file in:\n",
    "\"H:\\\\\\My Documents\\Vulnerability Curves\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_Info.csv'\n",
    "info = open(infopath,'w')\n",
    "\n",
    "info.write('Building_Type,Log_Shape,Log_Scale,Log_R_Squared,Log_RMSD,Exp_Shape,Exp_Scale,Exp_R_Squared,Exp_RMSD')\n",
    "\n",
    "# MINIMISATION PROCESS - Iterate thorugh building types\n",
    "\n",
    "result = []\n",
    "result_exp = []\n",
    "bands = [[],[],[],[],[],[],[],[],[],[]]\n",
    "xval = np.arange(0.,100., 0.1)\n",
    "\n",
    "for b in range(10):\n",
    "    \n",
    "    bands[b]= [[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "    \n",
    "    for pt in range(0,len(wind_pts[b])):\n",
    "        section = int(wind_pts[b][pt]/5-4)\n",
    "        bands[b][section].append(DI_pts[b][pt])\n",
    "    \n",
    "    # Create initial estimates for parameters (Weibull)\n",
    "    \n",
    "    shape_guess_exp = np.std(wind_pts[b])\n",
    "    scale_guess_exp = np.log(1./shape_guess_exp)\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "    \n",
    "    \n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_exp.append(minimise_exp(params_guess_exp,b))\n",
    "\n",
    "    \n",
    "    # Create initial estimates for parameters (lognorm)\n",
    "\n",
    "    shape_guess = np.std(wind_pts[b])\n",
    "    scale_guess = np.mean(np.log(shape_guess))    \n",
    "    params_guess = [shape_guess, scale_guess]\n",
    "\n",
    "    # Carry out minimisation process (lognorm)\n",
    "    \n",
    "    result.append(minimise(params_guess,b))\n",
    "    \n",
    "    # Calcualte Error\n",
    "        \n",
    "    error_exp = residuals_exp(result_exp[b][0],b)\n",
    "    R_squared_exp = 1 - np.sum(np.square(error_exp)) / np.sum( np.square( DI_pts[b] - np.mean(DI_pts[b]) ) )\n",
    "    RMSD_exp = np.sqrt(np.sum(np.square(error_exp))/(len(error_exp)-2))\n",
    "    \n",
    "    error = residuals(result[b][0],b)\n",
    "    R_squared = 1 - np.sum(np.square(error)) / np.sum( np.square( DI_pts[b] - np.mean(DI_pts[b]) ) )\n",
    "    RMSD = np.sqrt(np.sum(np.square(error))/(len(error)-2))\n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,12), dpi=1000)\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "    \n",
    "    fig1 = subplot(2,1,1)\n",
    "    \n",
    "    if b == 0:\n",
    "        hazimp_alpha = 0.1476\n",
    "        hazimp_beta = 3.9768\n",
    "    elif b in [1,2,6]:\n",
    "        hazimp_alpha = 0.1457\n",
    "        hazimp_beta = 3.9844\n",
    "    elif b == 3:\n",
    "        hazimp_alpha = 0.1129\n",
    "        hazimp_beta = 4.3397\n",
    "    elif b == 4:\n",
    "        hazimp_alpha = 0.1057\n",
    "        hazimp_beta = 4.2734\n",
    "    \n",
    "    hazimp_domain = np.arange(0, 115., 0.1) #frange(0,115,.1)\n",
    "    hazimp_domain[:] = [x/1.15 for x in hazimp_domain] # conversion from 0.2s to 3s gust  #FIXME: Correct conversion factor\n",
    "    \n",
    "    plot(\n",
    "        xval,lognorm.cdf(xval,result[b][0][0], scale = result[b][0][1]),'g',\n",
    "        xval,function_exp(result_exp[b][0],xval),'k',\n",
    "        markersize=10,lw=2)\n",
    "    \n",
    "    plot(\n",
    "        xval,lognorm.cdf(xval,result[b][0][0], scale = result[b][0][1])+RMSD,'g:',\n",
    "        xval,function_exp(result_exp[b][0],xval)+RMSD_exp,'k:',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    if b in [0,1,2,3,4,6]: \n",
    "        \n",
    "        plot(\n",
    "        hazimp_domain,function_exp([hazimp_alpha, hazimp_beta], np.arange(0, 115., 0.1)),'m',\n",
    "        markersize=10,lw=2)       \n",
    "        \n",
    "        boxplot(bands[b], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "        widths= 5, whis = 2)\n",
    "        \n",
    "        fig1.legend([\"lognorm cdf\", \"weibull cdf\", \"1 RMSD (lognorm)\", \"1 RMSD (weibull)\", \"HAZIMP CURVE\", \"Observations\"],loc=2)\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        boxplot(bands[b], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "        widths= 5, whis = 2)\n",
    "        \n",
    "        fig1.legend([\"lognorm cdf\", \"weibull cdf\",\"1 RMSD (lognorm)\", \"1 RMSD (weibull)\", \"Observations\"],loc=2)\n",
    "    \n",
    "    \n",
    "    plot(\n",
    "        xval,lognorm.cdf(xval,result[b][0][0], scale = result[b][0][1])-RMSD,'g:',\n",
    "        xval,function_exp(result_exp[b][0],xval)-RMSD_exp,'k:',\n",
    "        lw=1)\n",
    "    \n",
    "    if result[b][0][1]>1e10:\n",
    "        print 'Error - Building Type %i: lognorm.cdf least squares minimisation process failed to converge.' %(b+1)\n",
    "    \n",
    "    title = 'Vulnerability Curve (Least Squares Fit) - Buiding Type %i' %(b+1)\n",
    "    fig1.set_title(title)\n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Damage Index')\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    # subplot 2 - Residuals Plot\n",
    "    \n",
    "    fig2 =subplot(2,1,2)\n",
    "    \n",
    "    title = 'Residuals Plot - Buiding Type %i' %(b+1)\n",
    "    fig2.set_title(title)\n",
    "    fig2.set_xlabel('Wind Speed (3s Gust [m/s]')\n",
    "    fig2.set_ylabel('Residual Damage Index')\n",
    "    fig2.set_xlim((0,100))\n",
    "    fig2.set_ylim((-1,1))\n",
    "    \n",
    "    plot(wind_pts[b],error,'.g', markersize = 3)\n",
    "    plot(wind_pts[b],error_exp,'.k', markersize = 3)\n",
    "    plot([0,100],[RMSD,RMSD],'g-.')\n",
    "    plot([0,100],[RMSD_exp,RMSD_exp],'k:')\n",
    "    plot([0,100],[-RMSD,-RMSD],'g-.')\n",
    "    plot([0,100],[-RMSD_exp,-RMSD_exp],'k:')\n",
    "    plot([0,100],[0,0],'k')\n",
    "    \n",
    "    \n",
    "    fig2.legend([\"lognorm\",\"weibull\", \"1 RMSD (lognorm)\",\"1 RMSD (weibull)\"],loc=2)\n",
    "    \n",
    "    # Save Figure\n",
    "    \n",
    "    imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_BT%i.png' %(b+1)\n",
    "    fig.savefig(imgpath, dpi=300)\n",
    "    \n",
    "    # Enter data into info file\n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1, result[b][0][0], result[b][0][1], R_squared, RMSD, result_exp[b][0][0], result_exp[b][0][1], R_squared_exp, RMSD_exp))\n",
    "    \n",
    "info.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Plot lognorm.cdf curves for all building types on the same axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(8,6), dpi=100)\n",
    "\n",
    "fig1 = subplot(1,1,1)\n",
    "\n",
    "fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "fig1.set_ylabel('Damage Index')\n",
    "fig1.set_xlim((0,100))\n",
    "fig1.set_ylim((0,1))\n",
    "\n",
    "legend_str = []\n",
    "\n",
    "for b in range(7):\n",
    "    plot(frange(0,100,.1),function_exp(result_exp[b][0],frange(0,100,.1)))\n",
    "    legend_str.append('Building Type %i'%(b+1))\n",
    "    if result[b][0][1]>1e10:\n",
    "        print 'Error - Building Type %i: lognorm.cdf least squares minimisation process failed to converge.' %(b+1)\n",
    "    \n",
    "fig1.set_title('Vulnerability Curves (Weibull) - All Buiding Types') \n",
    "fig1.legend(legend_str, loc = 2)\n",
    "\n",
    "imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Curves_All_Types.png'\n",
    "fig.savefig(imgpath, dpi=300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(8,6), dpi=100)\n",
    "\n",
    "fig1 = subplot(1,1,1)\n",
    "\n",
    "fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "fig1.set_ylabel('Damage Index')\n",
    "fig1.set_xlim((0,100))\n",
    "fig1.set_ylim((0,1))\n",
    "\n",
    "legend_str = []\n",
    "\n",
    "for b in range(7):\n",
    "    plot(frange(0,100,.1), lognorm.cdf(frange(0,100,.1), result[b][0][0], scale = result[b][0][1]))\n",
    "    legend_str.append('Building Type %i'%(b+1))\n",
    "    if result[b][0][1]>1e10:\n",
    "        print 'Error - Building Type %i: lognorm.cdf least squares minimisation process failed to converge.' %(b+1)\n",
    "    \n",
    "fig1.set_title('Vulnerability Curves (Log-normal) - All Buiding Types') \n",
    "fig1.legend(legend_str, loc = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate and plot lognorm.cdf curves for typical code-compliant and uncompiant rediential buildings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_Info.csv'\n",
    "info = open(infopath,'a+')\n",
    "\n",
    "result_1980=[]\n",
    "\n",
    "error_1980 = [[],[],[]]\n",
    "R_squared_1980 = [[],[],[]]\n",
    "RMSD_1980 = [[],[],[]]\n",
    "\n",
    "\n",
    "# DO CALCULATIONS AND PLOTTING\n",
    "\n",
    "for c in [10,11,12]:\n",
    "    \n",
    "    # Create initial estimates for parameters\n",
    "    \n",
    "    shape_guess = np.std(wind_pts[c])\n",
    "    loc_guess = np.log(1/shape_guess)\n",
    "        \n",
    "    params_guess = [shape_guess, loc_guess]\n",
    "    \n",
    "    # Carry out minimisation process\n",
    "    \n",
    "    result_1980.append(minimise_exp(params_guess,c))\n",
    "    \n",
    "    error_1980[c-10] = residuals_exp(result_1980[c-10][0],c)\n",
    "    R_squared_1980[c-10] = 1 - np.sum(np.square(error_1980[c-10])) / np.sum( np.square( DI_pts[c] - np.mean(DI_pts[c]) ) )\n",
    "    RMSD_1980[c-10] = np.sqrt(np.sum(np.square(error_1980[c-10]))/(len(error_1980[c-10])-2))\n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(c+1,0 ,0 ,0 ,0 ,result_1980[c-10][0][0], result_1980[c-10][0][1], R_squared_1980[c-10],RMSD_1980[c-10]))\n",
    "\n",
    "# PLOT\n",
    "    \n",
    "fig,ax = subplots(figsize=(8,6), dpi=600)\n",
    "    \n",
    "title = 'Vulnerability Curve (Weibull) - 1980 Building Code Compliance'\n",
    "ax.set_title(title)\n",
    "ax.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "ax.set_ylabel('Damage Index')\n",
    "ax.set_xlim((0,100))\n",
    "ax.set_ylim((0,1))\n",
    "    \n",
    "    \n",
    "ax.plot(frange(0,100,.1),function_exp(result_1980[0][0],frange(0,100,.1)),'-k',\n",
    "        frange(0,100,.1),function_exp(result_1980[1][0],frange(0,100,.1)),'-r',\n",
    "        frange(0,100,.1),function_exp(result_1980[2][0],frange(0,100,.1)),'-g',\n",
    "        lw=2)\n",
    "\n",
    "ax.plot(frange(0,100,.1),function_exp(result_1980[0][0],frange(0,100,.1))+RMSD_1980[0],'k--',\n",
    "        frange(0,100,.1),function_exp(result_1980[1][0],frange(0,100,.1))+RMSD_1980[1],'r--',\n",
    "        frange(0,100,.1),function_exp(result_1980[2][0],frange(0,100,.1))-RMSD_1980[2],'g--',\n",
    "        frange(0,100,.1),function_exp(result_1980[0][0],frange(0,100,.1))-RMSD_1980[0],'k--',\n",
    "        frange(0,100,.1),function_exp(result_1980[1][0],frange(0,100,.1))-RMSD_1980[1],'r--',\n",
    "        frange(0,100,.1),function_exp(result_1980[2][0],frange(0,100,.1))+RMSD_1980[2],'g--',  \n",
    "        lw=1)        \n",
    "\n",
    "\n",
    "ax.plot(wind_pts[10],DI_pts[10], '.k',\n",
    "        wind_pts[11],DI_pts[11],'.r',\n",
    "        wind_pts[12],DI_pts[12], '.g',\n",
    "        markersize = 1)\n",
    "\n",
    "ax.legend([\"Code Compliant Buildings\", \"Pre-Code Buildings\", \"Industrial Sheds\", \"1 RMSD (compliant)\", \"1 RMSD (pre-code)\", \"1 RMSD (sheds)\"],loc=2)\n",
    "\n",
    "info.close()\n",
    "\n",
    "# Save Figure\n",
    "    \n",
    "imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_Compliance.png'\n",
    "fig.savefig(imgpath, dpi=300)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Functions for fragility curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def residuals_fragility0(params, i_value, curve, total_sections_value, X_sections):\n",
    "    \n",
    "    x=[]\n",
    "    x[:] = [(100/total_sections_value)/2 + m*(100/total_sections_value) for m in X_sections]\n",
    "    \n",
    "    y = []\n",
    "    y[:] = [section_probs[i_value][m][curve] for m in range(len(x))]\n",
    "    \n",
    "    output = y - function_exp(params,x)\n",
    "    return output\n",
    "\n",
    "def minimise_fragility0(params, i_value, curve, total_sections_value,X_sections):\n",
    "    output = leastsq(residuals_fragility0, params, args = (i_value, curve, total_sections_value,X_sections),\n",
    "                     full_output = 1, ftol=1E-10, xtol=1E-10, maxfev = 100000)\n",
    "    return output\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def residuals_fragility( scale, shape, i_value, curve, total_sections_value, X_sections):\n",
    "    \n",
    "    x=[]\n",
    "    x[:] = [(100/total_sections_value)/2 + m*(100/total_sections_value) for m in X_sections]\n",
    "    \n",
    "    y = []\n",
    "    y[:] = [section_probs[i_value][m][curve] for m in range(len(x))]\n",
    "    \n",
    "    output = y - function_exp([shape,scale],x)\n",
    "    return output\n",
    "\n",
    "def minimise_fragility(scale, shape, i_value, curve, total_sections_value,X_sections):\n",
    "    output = leastsq(residuals_fragility, scale, args = (shape, i_value, curve, total_sections_value,X_sections),\n",
    "                     full_output = 1, ftol=1E-10, xtol=1E-10, maxfev = 100000)\n",
    "    return output\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fragility Curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "section_bands = []\n",
    "section_probs = []\n",
    "section_x = []\n",
    "result_fragility = []\n",
    "\n",
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_Info.csv'\n",
    "info = open(infopath,'w')\n",
    "\n",
    "info.write('Building_Type,\\\n",
    "\\0.02_Shape,0.02_Scale,0.02_RMSD,\\\n",
    "\\0.12_Shape,0.12_Scale,0.12_RMSD,\\\n",
    "\\0.50_Shape,0.50_Scale,0.50_RMSD,\\\n",
    "\\0.90_Shape,0.90_Scale,0.90_RMSD')\n",
    "\n",
    "\n",
    "error_frag = [[],[],[],[]]\n",
    "R_squared_frag = [[],[],[],[]]\n",
    "RMSD_frag = [[],[],[],[]]\n",
    "\n",
    "## post 1980\n",
    "\n",
    "for i in [0]:\n",
    "    \n",
    "    total_sections = 20  # must be a factor of 100\n",
    "    curve_shape_multiplier = [.04,.035] # use to alter shaopes of DI=0.12, DI=0.5 and DI=0.9 curves\n",
    "    \n",
    "    section_bands.append([])\n",
    "    section_probs.append([])\n",
    "    section_x.append([])\n",
    "    result_fragility.append([[],[],[],[]])\n",
    "    \n",
    "    b=i+10\n",
    "    \n",
    "    wind = wind_pts[b]\n",
    "    DI = DI_pts[b]\n",
    "    \n",
    "    section_bands[i][:] = [ [0,0,0,0,0] for element in range(total_sections) ]\n",
    "    \n",
    "    for pt in range(len(wind)):\n",
    "        \n",
    "        section = int(wind[pt]/(100/total_sections))        \n",
    "        if section == total_sections:\n",
    "            section -= 1\n",
    "        \n",
    "        section_bands[i][section][0] += 1\n",
    "        \n",
    "        if DI[pt] >= 0.02:\n",
    "            section_bands[i][section][1] += 1\n",
    "        if DI[pt] >= 0.12:\n",
    "            section_bands[i][section][2] += 1\n",
    "        if DI[pt] >= 0.50:\n",
    "            section_bands[i][section][3] += 1\n",
    "        if DI[pt] >= 0.90:\n",
    "            section_bands[i][section][4] += 1\n",
    "    \n",
    "    for k in range(len(section_bands[i])): \n",
    "        if section_bands[i][k][0]!=0:\n",
    "            section_probs[i].append([])\n",
    "            section_x[i].append(k)\n",
    "            section_probs[i][-1][:] = [ float(section_bands[i][k][j])/float(section_bands[i][k][0]) for j in range(1,5)]\n",
    "    \n",
    "    for k in range(100):\n",
    "        section_x[i].append(total_sections*1.5) #1.2\n",
    "        section_probs[i].append([1,1,1,1])\n",
    "        section_x[i].append(total_sections*.3) #.3\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "\n",
    "    shape_guess_exp = .2\n",
    "    scale_guess_exp = 4\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "\n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_fragility[i][0] = minimise_fragility0( params_guess_exp , i, 0, total_sections, section_x[i])\n",
    "    shape_true = result_fragility[i][0][0][0]\n",
    "###    \n",
    "    result_fragility[i][1:] = [ minimise_fragility( scale_guess_exp, shape_true-multiplier*curve_shape_multiplier[i], i, curve, total_sections, section_x[i] )\n",
    "                              for (multiplier,curve) in [(1,1),(3.1,2),(3,3)] ]\n",
    "    \n",
    "    \n",
    "    # Error\n",
    "    error_frag[0] = residuals_fragility0(result_fragility[i][0][0],i,0, total_sections, section_x[i])\n",
    "    RMSD_frag[0] = np.sqrt(np.sum(np.square(error_frag[0][:-200]))/(len(error_frag[0][:-200])-2))\n",
    "    \n",
    "###    \n",
    "    for (multiplier,curve) in [(1,1),(3.1,2),(3,3)]:\n",
    "        error_frag[curve] = residuals_fragility(result_fragility[i][curve][0],shape_true-multiplier*curve_shape_multiplier[i],i,curve, total_sections, section_x[i])\n",
    "        RMSD_frag[curve] = np.sqrt(np.sum(np.square(error_frag[curve][:-200]))/(len(error_frag[curve][:-200])-1))\n",
    "    \n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1,result_fragility[i][0][0][0] ,result_fragility[i][0][0][1] ,RMSD_frag[0],shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0],RMSD_frag[1],shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0],RMSD_frag[2],shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0],RMSD_frag[3]))        \n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,6), dpi=1000)\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "        \n",
    "    fig1 = subplot(1,1,1)\n",
    "    \n",
    "    xx=[]\n",
    "    xx[:] = [(100/total_sections)/2 + m*(100/total_sections) for m in section_x[i]]\n",
    "    yy = []\n",
    "  \n",
    "    for curve in range(4):\n",
    "        \n",
    "        yy.append([])\n",
    "        yy[curve][:] = [section_probs[i][n][curve] for n in range(len(xx))]\n",
    "###\n",
    "    multiplier = [1,3.1,3]\n",
    "        \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1)), '-g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1)), '-b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1)), '-r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1)), '-k',\n",
    "        xx[:-200],yy[3][:-200],'.k', [-1],[-1], '--k',\n",
    "        xx[:-200],yy[0][:-200],'.g',\n",
    "        xx[:-200],yy[1][:-200],'.b',\n",
    "        xx[:-200],yy[2][:-200],'.r',   \n",
    "        lw=2,markersize=15)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))+RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))+RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))+RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))+RMSD_frag[3], '--k',\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))-RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))-RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))-RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))-RMSD_frag[3], '--k',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Probability of Damage')\n",
    "    fig1.legend(['DI = 0.02', 'DI = 0.12','DI = 0.50','DI = 0.90', 'data point (5m/s section)', '1 RMSD'],loc=2)\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    if i==0:\n",
    "        title = 'Fragility Plot (Weibull) - Post 1980 Buildings (Code Compliant)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_post1980.png'\n",
    "    elif i==1:\n",
    "        title = 'Fragility Plot (Weibull) - Pre 1980 Buildings'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_pre1980.png'\n",
    "    \n",
    "    fig1.set_title(title)\n",
    "    fig.savefig(imgpath, dpi=600)\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "## pre 1980 \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for i in [1]:\n",
    "    \n",
    "    total_sections = 20  # must be a factor of 100\n",
    "    curve_shape_multiplier = [.01,.036] # use to alter shaopes of DI=0.12, DI=0.5 and DI=0.9 curves\n",
    "    \n",
    "    section_bands.append([])\n",
    "    section_probs.append([])\n",
    "    section_x.append([])\n",
    "    result_fragility.append([[],[],[],[]])\n",
    "    \n",
    "    b=i+10\n",
    "    \n",
    "    wind = wind_pts[b]\n",
    "    DI = DI_pts[b]\n",
    "    \n",
    "    section_bands[i][:] = [ [0,0,0,0,0] for element in range(total_sections) ]\n",
    "    \n",
    "    for pt in range(len(wind)):\n",
    "        \n",
    "        section = int(wind[pt]/(100/total_sections))        \n",
    "        if section == total_sections:\n",
    "            section -= 1\n",
    "        \n",
    "        section_bands[i][section][0] += 1\n",
    "        \n",
    "        if DI[pt] >= 0.02:\n",
    "            section_bands[i][section][1] += 1\n",
    "        if DI[pt] >= 0.12:\n",
    "            section_bands[i][section][2] += 1\n",
    "        if DI[pt] >= 0.50:\n",
    "            section_bands[i][section][3] += 1\n",
    "        if DI[pt] >= 0.90:\n",
    "            section_bands[i][section][4] += 1\n",
    "    \n",
    "    for k in range(len(section_bands[i])): \n",
    "        if section_bands[i][k][0]!=0:\n",
    "            section_probs[i].append([])\n",
    "            section_x[i].append(k)\n",
    "            section_probs[i][-1][:] = [ float(section_bands[i][k][j])/float(section_bands[i][k][0]) for j in range(1,5)]\n",
    "    \n",
    "    for k in range(100):\n",
    "        section_x[i].append(total_sections*1.5) #1.2\n",
    "        section_probs[i].append([1,1,1,1])\n",
    "        section_x[i].append(total_sections*.3) #.3\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "\n",
    "    shape_guess_exp = .2\n",
    "    scale_guess_exp = 4\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "\n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_fragility[i][0] = minimise_fragility0( params_guess_exp , i, 0, total_sections, section_x[i])\n",
    "    shape_true = result_fragility[i][0][0][0]\n",
    "###    \n",
    "    result_fragility[i][1:] = [ minimise_fragility( scale_guess_exp, shape_true-multiplier*curve_shape_multiplier[i], i, curve, total_sections, section_x[i] )\n",
    "                              for (multiplier,curve) in [(1,1),(3,2),(3.2,3)] ]\n",
    "    \n",
    "    \n",
    "    # Error\n",
    "    error_frag[0] = residuals_fragility0(result_fragility[i][0][0],i,0, total_sections, section_x[i])\n",
    "    RMSD_frag[0] = np.sqrt(np.sum(np.square(error_frag[0][:-200]))/(len(error_frag[0][:-200])-2))\n",
    "    \n",
    "###    \n",
    "    for (multiplier,curve) in [(1,1),(3,2),(3.2,3)]:\n",
    "        error_frag[curve] = residuals_fragility(result_fragility[i][curve][0],shape_true-multiplier*curve_shape_multiplier[i],i,curve, total_sections, section_x[i])\n",
    "        RMSD_frag[curve] = np.sqrt(np.sum(np.square(error_frag[curve][:-200]))/(len(error_frag[curve][:-200])-1))\n",
    "    \n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1,result_fragility[i][0][0][0] ,result_fragility[i][0][0][1] ,RMSD_frag[0],shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0],RMSD_frag[1],shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0],RMSD_frag[2],shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0],RMSD_frag[3]))        \n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,6), dpi=1000)\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "        \n",
    "    fig1 = subplot(1,1,1)\n",
    "    \n",
    "    xx=[]\n",
    "    xx[:] = [(100/total_sections)/2 + m*(100/total_sections) for m in section_x[i]]\n",
    "    yy = []\n",
    "\n",
    "    for curve in range(4):\n",
    "        \n",
    "        yy.append([])\n",
    "        yy[curve][:] = [section_probs[i][n][curve] for n in range(len(xx))]\n",
    "###\n",
    "    multiplier = [1,3,3.2]\n",
    "        \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1)), '-g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1)), '-b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1)), '-r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1)), '-k',\n",
    "        xx[:-200],yy[3][:-200],'.k', [-1],[-1], '--k',\n",
    "        xx[:-200],yy[0][:-200],'.g',\n",
    "        xx[:-200],yy[1][:-200],'.b',\n",
    "        xx[:-200],yy[2][:-200],'.r',   \n",
    "        lw=2,markersize=15)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))+RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))+RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))+RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))+RMSD_frag[3], '--k',\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))-RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))-RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))-RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))-RMSD_frag[3], '--k',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Probability of Damage')\n",
    "    fig1.legend(['DI = 0.02', 'DI = 0.12','DI = 0.50','DI = 0.90', 'data point (5m/s section)', '1 RMSD'],loc=2)\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    if i==0:\n",
    "        title = 'Fragility Plot (Weibull) - Post 1980 Buildings (Code Compliant)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_post1980.png'\n",
    "    elif i==1:\n",
    "        title = 'Fragility Plot (Weibull) - Pre 1980 Buildings'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_pre1980.png'\n",
    "    \n",
    "    fig1.set_title(title)\n",
    "    fig.savefig(imgpath, dpi=600)\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "info.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Combine TC Larry and TC Yasi data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_Info.csv'\n",
    "info = open(infopath,'a+')\n",
    "\n",
    "result_LY = []\n",
    "bands_LY = [[],[],[],[]]\n",
    "\n",
    "error_LY = [[],[]]\n",
    "R_squared_LY = [[],[]]\n",
    "RMSD_LY = [[],[]]\n",
    "\n",
    "\n",
    "\n",
    "for i, c in [(0,10),(1,11),(2,13),(3,14)]:\n",
    "    bands_LY[i]= [[],[],[],[],[],[],[],[],[],[],[],[],[],[]]\n",
    "    for pt in range(0,len(wind_pts[c])):\n",
    "        section = int(wind_pts[c][pt]/5-4)\n",
    "        bands_LY[i][section].append(DI_pts[c][pt])\n",
    "\n",
    "        \n",
    "for i,c in [(0,15),(1,16)]:\n",
    "    \n",
    "    # Create initial estimates for parameters\n",
    "    \n",
    "    shape_guess = np.std(wind_pts[c])\n",
    "    loc_guess = np.log(1/shape_guess)\n",
    "        \n",
    "    params_guess = [shape_guess, loc_guess]\n",
    "\n",
    "    # Carry out minimisation process\n",
    "    \n",
    "    result_LY.append(minimise_exp(params_guess,c))\n",
    "    \n",
    "    # Error\n",
    "    \n",
    "    error_LY[i] = residuals_exp(result_LY[i][0],c)\n",
    "    R_squared_LY[i] = 1 - np.sum(np.square(error_LY[i])) / np.sum( np.square( DI_pts[c] - np.mean(DI_pts[c]) ) )\n",
    "    RMSD_LY[i] = np.sqrt(np.sum(np.square(error_LY[i]))/(len(error_LY[i])-2))\n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(c+1,0 ,0 ,0 ,0 ,result_LY[i][0][0], result_LY[i][0][1], R_squared_LY[i],RMSD_LY[i]))\n",
    "\n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig,ax = subplots(figsize=(8,6), dpi=600)\n",
    "    \n",
    "    \n",
    "    ax.plot([],[],'-k',lw=2)\n",
    "    ax.plot([],[],'-g',lw=2)\n",
    "    ax.plot([],[],'-b',lw=2)\n",
    "    ax.plot([],[],'-r',lw=2)\n",
    "    \n",
    "    \n",
    "    if c == 15:\n",
    "        bp1 = boxplot(bands_LY[0], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "            widths= 5, whis = 2, patch_artist = True)\n",
    "        \n",
    "        bp2 = boxplot(bands_LY[2], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "            widths= 5, whis = 2)\n",
    "      \n",
    "    if c == 16:\n",
    "        bp1 =boxplot(bands_LY[1], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "            widths= 5, whis = 2, patch_artist = True)\n",
    "        \n",
    "        bp2 = boxplot(bands_LY[3], positions = (22.5,27.5,32.5,37.5,42.5,47.5,52.5,57.5,62.5,67.5,72.5,77.5,82.5,87.5),\n",
    "            widths= 5, whis = 2)\n",
    "        \n",
    "         \n",
    "    \n",
    "    setp(bp1['boxes'], color='blue')\n",
    "    setp(bp1['whiskers'], color='blue') \n",
    "    setp(bp1['fliers'], color='blue')\n",
    "    setp(bp1['medians'], color='white',lw=2)\n",
    "    setp(bp1['caps'], color='blue')\n",
    "    \n",
    "    setp(bp2['boxes'], color='red')\n",
    "    setp(bp2['whiskers'], color='red')\n",
    "    setp(bp2['fliers'], color='red')\n",
    "    setp(bp2['medians'], color='red',lw=2)\n",
    "    setp(bp2['caps'], color='red')\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_LY[i][0],frange(0,100,.1)),'k',\n",
    "        markersize=10,lw=2)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_1980[i][0],frange(0,100,.1)),'g',\n",
    "        markersize=10,lw=2)\n",
    "    \n",
    "    ax.plot(frange(0,100,.1),function_exp(result_LY[i][0],frange(0,100,.1))+RMSD_LY[i],'k--',\n",
    "        frange(0,100,.1),function_exp(result_LY[i][0],frange(0,100,.1))-RMSD_LY[i],'k--',\n",
    "        frange(0,100,.1),function_exp(result_1980[i][0],frange(0,100,.1))+RMSD_1980[i],'g--',\n",
    "        frange(0,100,.1),function_exp(result_1980[i][0],frange(0,100,.1))-RMSD_1980[i],'g--',\n",
    "        lw=1)  \n",
    "    \n",
    "    \n",
    "    \n",
    "    ax.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    ax.set_ylabel('Damage Index')\n",
    "    ax.set_xlim((0,100))\n",
    "    ax.set_ylim((0,1))\n",
    "    ax.set_xticks([20, 40, 60, 80, 100])\n",
    "    ax.legend([\"Weibull - Combined\", \"Weibull - TC Yasi\", \"Boxplot - TC Yasi\", \"Boxplot - TC Larry\"],loc=2)\n",
    "    \n",
    "    \n",
    "    # Save Figure\n",
    "\n",
    "    if c == 15:\n",
    "        title = 'Vulnerability Curve - TC Yasi and TC Larry - Code Compliant Post 1980'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_YasiLarry_post1980.png'\n",
    "    elif c ==16:\n",
    "        title = 'Vulnerability Curve - TC Yasi and TC Larry - Pre 1980'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Vulnerability_Plot_YasiLarry_pre1980.png'\n",
    "\n",
    "    ax.set_title(title)\n",
    "    \n",
    "    fig.savefig(imgpath, dpi=600)\n",
    "\n",
    "info.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fragility plots - TC Yasi and TC Larry Combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "section_bands = []\n",
    "section_probs = []\n",
    "section_x = []\n",
    "result_fragility = []\n",
    "\n",
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_Info.csv'\n",
    "info = open(infopath,'a+')\n",
    "\n",
    "\n",
    "error_frag = [[],[],[],[]]\n",
    "R_squared_frag = [[],[],[],[]]\n",
    "RMSD_frag = [[],[],[],[]]\n",
    "\n",
    "\n",
    "### post-1980\n",
    "\n",
    "for i in [0]:\n",
    "    \n",
    "    total_sections = 20  # must be a factor of 100\n",
    "    curve_shape_multiplier = [.033,.037] # use to alter shaopes of DI=0.12, DI=0.5 and DI=0.9 curves\n",
    "    \n",
    "    section_bands.append([])\n",
    "    section_probs.append([])\n",
    "    section_x.append([])\n",
    "    result_fragility.append([[],[],[],[]])\n",
    "    \n",
    "    b=i+15\n",
    "    \n",
    "    wind = wind_pts[b]\n",
    "    DI = DI_pts[b]\n",
    "    \n",
    "    section_bands[i][:] = [ [0,0,0,0,0] for element in range(total_sections) ]\n",
    "    \n",
    "    for pt in range(len(wind)):\n",
    "        \n",
    "        section = int(wind[pt]/(100/total_sections))        \n",
    "        if section == total_sections:\n",
    "            section -= 1\n",
    "        \n",
    "        section_bands[i][section][0] += 1\n",
    "        \n",
    "        if DI[pt] >= 0.02:\n",
    "            section_bands[i][section][1] += 1\n",
    "        if DI[pt] >= 0.12:\n",
    "            section_bands[i][section][2] += 1\n",
    "        if DI[pt] >= 0.50:\n",
    "            section_bands[i][section][3] += 1\n",
    "        if DI[pt] >= 0.90:\n",
    "            section_bands[i][section][4] += 1\n",
    "    \n",
    "    for k in range(len(section_bands[i])): \n",
    "        if section_bands[i][k][0]!=0:\n",
    "            section_probs[i].append([])\n",
    "            section_x[i].append(k)\n",
    "            section_probs[i][-1][:] = [ float(section_bands[i][k][j])/float(section_bands[i][k][0]) for j in range(1,5)]\n",
    "    \n",
    "    for k in range(100):\n",
    "        section_x[i].append(total_sections*1.5) #1.2\n",
    "        section_probs[i].append([1,1,1,1])\n",
    "        section_x[i].append(total_sections*.3) #.3\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "\n",
    "    shape_guess_exp = .2\n",
    "    scale_guess_exp = 4\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "\n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_fragility[i][0] = minimise_fragility0( params_guess_exp , i, 0, total_sections, section_x[i])\n",
    "    shape_true = result_fragility[i][0][0][0]\n",
    "###    \n",
    "    result_fragility[i][1:] = [ minimise_fragility( scale_guess_exp, shape_true-multiplier*curve_shape_multiplier[i], i, curve, total_sections, section_x[i] )\n",
    "                              for (multiplier,curve) in [(1,1),(3.1,2),(3,3)] ]\n",
    "    \n",
    "    \n",
    "    # Error\n",
    "    error_frag[0] = residuals_fragility0(result_fragility[i][0][0],i,0, total_sections, section_x[i])\n",
    "    RMSD_frag[0] = np.sqrt(np.sum(np.square(error_frag[0][:-200]))/(len(error_frag[0][:-200])-2))\n",
    "    \n",
    "###    \n",
    "    for (multiplier,curve) in [(1,1),(3.1,2),(3,3)]:\n",
    "        error_frag[curve] = residuals_fragility(result_fragility[i][curve][0],shape_true-multiplier*curve_shape_multiplier[i],i,curve, total_sections, section_x[i])\n",
    "        RMSD_frag[curve] = np.sqrt(np.sum(np.square(error_frag[curve][:-200]))/(len(error_frag[curve][:-200])-1))\n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1,result_fragility[i][0][0][0] ,result_fragility[i][0][0][1] ,RMSD_frag[0],shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0],RMSD_frag[1],shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0],RMSD_frag[2],shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0],RMSD_frag[3]))        \n",
    "    \n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,6), dpi=1000)\n",
    "\n",
    "###\n",
    "    multiplier =[1,3.1,3]\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "        \n",
    "    fig1 = subplot(1,1,1)\n",
    "    \n",
    "    xx=[]\n",
    "    xx[:] = [(100/total_sections)/2 + m*(100/total_sections) for m in section_x[i]]\n",
    "    yy = []\n",
    "    \n",
    "    for curve in range(4):\n",
    "        \n",
    "        yy.append([])\n",
    "        yy[curve][:] = [section_probs[i][n][curve] for n in range(len(xx))]\n",
    "        \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1)), '-g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1)), '-b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1)), '-r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1)), '-k',\n",
    "        xx[:-200],yy[3][:-200],'.k', [-1],[-1], '--k',\n",
    "        xx[:-200],yy[0][:-200],'.g',\n",
    "        xx[:-200],yy[1][:-200],'.b',\n",
    "        xx[:-200],yy[2][:-200],'.r',   \n",
    "        lw=2,markersize=15)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))+RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))+RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))+RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))+RMSD_frag[3], '--k',\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))-RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))-RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))-RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))-RMSD_frag[3], '--k',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Probability of Damage')\n",
    "    fig1.legend(['DI = 0.02', 'DI = 0.12','DI = 0.50','DI = 0.90', 'data point (5m/s section)', '1 RMSD'],loc=2)\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    if i==0:\n",
    "        title = 'Fragility Plot (Weibull) - Post 1980 (TC Yasi & TC Larry Combined)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_post1980_combined.png'\n",
    "    elif i==1:\n",
    "        title = 'Fragility Plot (Weibull) - Pre 1980 Buildings (TC Yasi & TC Larry Combined)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_pre1980_combined.png'\n",
    "    \n",
    "    fig1.set_title(title)\n",
    "    fig.savefig(imgpath, dpi=600)\n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "### pre-1980    \n",
    "    \n",
    "    \n",
    "for i in [1]:\n",
    "    \n",
    "    total_sections = 20  # must be a factor of 100\n",
    "    curve_shape_multiplier = [.01,.036] # use to alter shaopes of DI=0.12, DI=0.5 and DI=0.9 curves\n",
    "    \n",
    "    section_bands.append([])\n",
    "    section_probs.append([])\n",
    "    section_x.append([])\n",
    "    result_fragility.append([[],[],[],[]])\n",
    "    \n",
    "    b=i+15\n",
    "    \n",
    "    wind = wind_pts[b]\n",
    "    DI = DI_pts[b]\n",
    "    \n",
    "    section_bands[i][:] = [ [0,0,0,0,0] for element in range(total_sections) ]\n",
    "    \n",
    "    for pt in range(len(wind)):\n",
    "        \n",
    "        section = int(wind[pt]/(100/total_sections))        \n",
    "        if section == total_sections:\n",
    "            section -= 1\n",
    "        \n",
    "        section_bands[i][section][0] += 1\n",
    "        \n",
    "        if DI[pt] >= 0.02:\n",
    "            section_bands[i][section][1] += 1\n",
    "        if DI[pt] >= 0.12:\n",
    "            section_bands[i][section][2] += 1\n",
    "        if DI[pt] >= 0.50:\n",
    "            section_bands[i][section][3] += 1\n",
    "        if DI[pt] >= 0.90:\n",
    "            section_bands[i][section][4] += 1\n",
    "    \n",
    "    for k in range(len(section_bands[i])): \n",
    "        if section_bands[i][k][0]!=0:\n",
    "            section_probs[i].append([])\n",
    "            section_x[i].append(k)\n",
    "            section_probs[i][-1][:] = [ float(section_bands[i][k][j])/float(section_bands[i][k][0]) for j in range(1,5)]\n",
    "    \n",
    "    for k in range(100):\n",
    "        section_x[i].append(total_sections*1.5) #1.2\n",
    "        section_probs[i].append([1,1,1,1])\n",
    "        section_x[i].append(total_sections*.3) #.3\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "\n",
    "    shape_guess_exp = .2\n",
    "    scale_guess_exp = 4\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "\n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_fragility[i][0] = minimise_fragility0( params_guess_exp , i, 0, total_sections, section_x[i])\n",
    "    shape_true = result_fragility[i][0][0][0]\n",
    "###    \n",
    "    result_fragility[i][1:] = [ minimise_fragility( scale_guess_exp, shape_true-multiplier*curve_shape_multiplier[i], i, curve, total_sections, section_x[i] )\n",
    "                              for (multiplier,curve) in [(1,1),(3,2),(3.2,3)] ]\n",
    "    \n",
    "    \n",
    "    # Error\n",
    "    error_frag[0] = residuals_fragility0(result_fragility[i][0][0],i,0, total_sections, section_x[i])\n",
    "    RMSD_frag[0] = np.sqrt(np.sum(np.square(error_frag[0][:-200]))/(len(error_frag[0][:-200])-2))\n",
    "    \n",
    "###    \n",
    "    for (multiplier,curve) in [(1,1),(3,2),(3.2,3)]:\n",
    "        error_frag[curve] = residuals_fragility(result_fragility[i][curve][0],shape_true-multiplier*curve_shape_multiplier[i],i,curve, total_sections, section_x[i])\n",
    "        RMSD_frag[curve] = np.sqrt(np.sum(np.square(error_frag[curve][:-200]))/(len(error_frag[curve][:-200])-1))\n",
    "    \n",
    "    info.write('\\n%i, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1,result_fragility[i][0][0][0] ,result_fragility[i][0][0][1] ,RMSD_frag[0],shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0],RMSD_frag[1],shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0],RMSD_frag[2],shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0],RMSD_frag[3]))        \n",
    "    \n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,6), dpi=1000)\n",
    "\n",
    "###\n",
    "    multiplier =[1,3,3.2]\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "        \n",
    "    fig1 = subplot(1,1,1)\n",
    "    \n",
    "    xx=[]\n",
    "    xx[:] = [(100/total_sections)/2 + m*(100/total_sections) for m in section_x[i]]\n",
    "    yy = []\n",
    "    \n",
    "    for curve in range(4):\n",
    "        \n",
    "        yy.append([])\n",
    "        yy[curve][:] = [section_probs[i][n][curve] for n in range(len(xx))]\n",
    "        \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1)), '-g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1)), '-b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1)), '-r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1)), '-k',\n",
    "        xx[:-200],yy[3][:-200],'.k', [-1],[-1], '--k',\n",
    "        xx[:-200],yy[0][:-200],'.g',\n",
    "        xx[:-200],yy[1][:-200],'.b',\n",
    "        xx[:-200],yy[2][:-200],'.r',   \n",
    "        lw=2,markersize=15)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))+RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))+RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))+RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))+RMSD_frag[3], '--k',\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))-RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[0]*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))-RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[1]*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))-RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-multiplier[2]*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))-RMSD_frag[3], '--k',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Probability of Damage')\n",
    "    fig1.legend(['DI = 0.02', 'DI = 0.12','DI = 0.50','DI = 0.90', 'data point (5m/s section)', '1 RMSD'],loc=2)\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    if i==0:\n",
    "        title = 'Fragility Plot (Weibull) - Post 1980 (TC Yasi & TC Larry Combined)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_post1980_combined.png'\n",
    "    elif i==1:\n",
    "        title = 'Fragility Plot (Weibull) - Pre 1980 Buildings (TC Yasi & TC Larry Combined)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_pre1980_combined.png'\n",
    "    \n",
    "    fig1.set_title(title)\n",
    "    fig.savefig(imgpath, dpi=600)    \n",
    "    \n",
    "    \n",
    "    \n",
    "info.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fragility Curves - TC Yasi only (original unmodified shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "section_bands = []\n",
    "section_probs = []\n",
    "section_x = []\n",
    "result_fragility = []\n",
    "\n",
    "infopath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_Info.csv'\n",
    "info = open(infopath,'a+')\n",
    "\n",
    "\n",
    "error_frag = [[],[],[],[]]\n",
    "R_squared_frag = [[],[],[],[]]\n",
    "RMSD_frag = [[],[],[],[]]\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "for i in range(2):\n",
    "    \n",
    "    total_sections = 20  # must be a factor of 100\n",
    "    curve_shape_multiplier = [.0,.0] # use to alter shaopes of DI=0.12, DI=0.5 and DI=0.9 curves\n",
    "    \n",
    "    section_bands.append([])\n",
    "    section_probs.append([])\n",
    "    section_x.append([])\n",
    "    result_fragility.append([[],[],[],[]])\n",
    "    \n",
    "    b=i+10\n",
    "    \n",
    "    wind = wind_pts[b]\n",
    "    DI = DI_pts[b]\n",
    "    \n",
    "    section_bands[i][:] = [ [0,0,0,0,0] for element in range(total_sections) ]\n",
    "    \n",
    "    for pt in range(len(wind)):\n",
    "        \n",
    "        section = int(wind[pt]/(100/total_sections))        \n",
    "        if section == total_sections:\n",
    "            section -= 1\n",
    "        \n",
    "        section_bands[i][section][0] += 1\n",
    "        \n",
    "        if DI[pt] >= 0.02:\n",
    "            section_bands[i][section][1] += 1\n",
    "        if DI[pt] >= 0.12:\n",
    "            section_bands[i][section][2] += 1\n",
    "        if DI[pt] >= 0.50:\n",
    "            section_bands[i][section][3] += 1\n",
    "        if DI[pt] >= 0.90:\n",
    "            section_bands[i][section][4] += 1\n",
    "    \n",
    "    for k in range(len(section_bands[i])): \n",
    "        if section_bands[i][k][0]!=0:\n",
    "            section_probs[i].append([])\n",
    "            section_x[i].append(k)\n",
    "            section_probs[i][-1][:] = [ float(section_bands[i][k][j])/float(section_bands[i][k][0]) for j in range(1,5)]\n",
    "    \n",
    "#    if i ==0:\n",
    "    for k in range(100):\n",
    "        section_x[i].append(total_sections*0) #1.2\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "        section_x[i].append(total_sections*0) #.3\n",
    "        section_probs[i].append([0,0,0,0])\n",
    "        \n",
    "    shape_guess_exp = .2\n",
    "    scale_guess_exp = 4\n",
    "    params_guess_exp = [shape_guess_exp, scale_guess_exp]\n",
    "\n",
    "    # Carry out minimisation process (Weibull)\n",
    "    \n",
    "    result_fragility[i][0] = minimise_fragility0( params_guess_exp , i, 0, total_sections, section_x[i])\n",
    "    shape_true = result_fragility[i][0][0][0]\n",
    "    \n",
    "    result_fragility[i][1:] = [ minimise_fragility( scale_guess_exp, shape_true-curve*curve_shape_multiplier[i], i, curve, total_sections, section_x[i] )\n",
    "                              for (multiplier,curve) in [(1,1),(2,2),(3,3)] ]\n",
    "    \n",
    "    \n",
    "    # Error\n",
    "    error_frag[0] = residuals_fragility0(result_fragility[i][0][0],i,0, total_sections, section_x[i])\n",
    "    RMSD_frag[0] = np.sqrt(np.sum(np.square(error_frag[0][:-200]))/(len(error_frag[0][:-200])-2))\n",
    "    \n",
    "    \n",
    "    for u in range(1,4):\n",
    "        error_frag[u] = residuals_fragility(result_fragility[i][u][0],shape_true-u*curve_shape_multiplier[i],i,u, total_sections, section_x[i])\n",
    "        RMSD_frag[u] = np.sqrt(np.sum(np.square(error_frag[u][:-200]))/(len(error_frag[u][:-200])-1))\n",
    "    \n",
    "    info.write('\\n%i (original unmodified shape), %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f, %f'\n",
    "               %(b+1,result_fragility[i][0][0][0] ,result_fragility[i][0][0][1] ,RMSD_frag[0],shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0],RMSD_frag[1],shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0],RMSD_frag[2],shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0],RMSD_frag[3]))        \n",
    "    \n",
    "    \n",
    "    # PLOT RESULTS\n",
    "    \n",
    "    fig = plt.figure(figsize=(8,6), dpi=1000)\n",
    "    \n",
    "    # Subplot 1 - Vulnerability Curve\n",
    "        \n",
    "    fig1 = subplot(1,1,1)\n",
    "    \n",
    "    xx=[]\n",
    "    xx[:] = [(100/total_sections)/2 + m*(100/total_sections) for m in section_x[i]]\n",
    "    yy = []\n",
    "    \n",
    "    for curve in range(4):\n",
    "        \n",
    "        yy.append([])\n",
    "        yy[curve][:] = [section_probs[i][n][curve] for n in range(len(xx))]\n",
    "        \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1)), '-g',\n",
    "        frange(0,100,.1),function_exp([shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1)), '-b',\n",
    "        frange(0,100,.1),function_exp([shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1)), '-r',\n",
    "        frange(0,100,.1),function_exp([shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1)), '-k',\n",
    "        xx[:-200],yy[3][:-200],'.k', [-1],[-1], '--k',\n",
    "        xx[:-200],yy[0][:-200],'.g',\n",
    "        xx[:-200],yy[1][:-200],'.b',\n",
    "        xx[:-200],yy[2][:-200],'.r',   \n",
    "        lw=2,markersize=15)\n",
    "    \n",
    "    plot(\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))+RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))+RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))+RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))+RMSD_frag[3], '--k',\n",
    "        frange(0,100,.1),function_exp(result_fragility[i][0][0],frange(0,100,.1))-RMSD_frag[0], '--g',\n",
    "        frange(0,100,.1),function_exp([shape_true-1*curve_shape_multiplier[i],result_fragility[i][1][0]],frange(0,100,.1))-RMSD_frag[1], '--b',\n",
    "        frange(0,100,.1),function_exp([shape_true-2*curve_shape_multiplier[i],result_fragility[i][2][0]],frange(0,100,.1))-RMSD_frag[2], '--r',\n",
    "        frange(0,100,.1),function_exp([shape_true-3*curve_shape_multiplier[i],result_fragility[i][3][0]],frange(0,100,.1))-RMSD_frag[3], '--k',\n",
    "        lw=1)\n",
    "    \n",
    "    \n",
    "    \n",
    "    fig1.set_xlabel('Wind Speed (3s Gust) [m/s]')\n",
    "    fig1.set_ylabel('Probability of Damage')\n",
    "    fig1.legend(['DI = 0.02', 'DI = 0.12','DI = 0.50','DI = 0.90', 'data point (5m/s section)', '1 RMSD'],loc=2)\n",
    "    fig1.set_xlim((0,100))\n",
    "    fig1.set_ylim((0,1))\n",
    "    fig1.set_xticks([20, 40, 60, 80, 100])\n",
    "    \n",
    "    if i==0:\n",
    "        title = 'Fragility Plot (Weibull) - Post 1980 (TC Yasi only, Original Unmodified Shape)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_post1980_unmodified.png'\n",
    "    elif i==1:\n",
    "        title = 'Fragility Plot (Weibull) - Pre 1980 Buildings (TC Yasi only, Original Unmodified Shape)'\n",
    "        imgpath = 'C:/WorkSpace/data/Derived/impact/yasi/output/Fragility_Plot_pre1980_unmodified.png'\n",
    "    \n",
    "    fig1.set_title(title)\n",
    "    fig.savefig(imgpath, dpi=600)\n",
    "\n",
    "info.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}