{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Statistical constraints on observing low-angle normal fault seismicity"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Introduction\n",
      "\n",
      "Low-angle normal faults (LANFs), with dips less than 30$^\\circ$, are\n",
      "well described in the geologic record. They are hypothesized to have\n",
      "an important role in accommodating large-magnitude continental\n",
      "extension \\citep{howard1987crustal} and crustal thinning\n",
      "\\citep{lister1986detachment}, and their recognition has been one of\n",
      "the most important developments in tectonics over the past several\n",
      "decades \\citep{wernicke2009detachment}. However, despite widespread\n",
      "field observations of inactive LANFs, and their central role in modern\n",
      "extensional tectonic theory, they remain enigmatic and contentious\n",
      "structures, and it is not clear if they are seismically active at low dip \n",
      "angles in the upper crust. This is for two reasons: because brittle faulting\n",
      "on LANFs is in apparent conflict with standard Andersonian rock mechanical\n",
      "theory as typically applied to the upper crust\n",
      "\\citep{axen2004lanfmech}, and because observations of active faulting\n",
      "on LANFs is sparse and at times ambiguous \\citep{wernicke1995seis}. A\n",
      "considerable amount of research has been performed to address the\n",
      "former concern, reconciling LANF slip with rock mechanics \\citep [e.g.,]\n",
      "[]{axenbartley1997, collettini2011lanfmech}. The latter issue, the paucity of \n",
      "definitive LANF earthquakes in focal mechanism catalogs, has inhibited\n",
      "hypothesis testing of LANF fault theory, and has also contributed to a mode\n",
      "of thought where the  absence of evidence for LANF seismicity is taken\n",
      "as evidence that LANFs are inactive or aseismic \\citep{jackson1987,\n",
      "collettinisibson2001}. However, the lack of observed seismic slip on\n",
      "continental LANFs may be simply due to the rarity of earthquakes on\n",
      "LANFs, coupled with the small number of potential active structures.\n",
      "\n",
      "In this work, we choose to directly address the question of whether\n",
      "the lack of observed large earthquakes on LANFs may be \n",
      "interpreted as an indication that\n",
      "LANFs do not slip seismically, or if it is better explained as an\n",
      "effect of the small number of active seismogenic LANFs with low moment\n",
      "accumulation rates (and hence long recurrence intervals). \n",
      "We then calculate the maximum likelihood of observing a significant\n",
      "earthquake on a LANF by treating all potentially active LANFs described in the \n",
      "literature as seismically active at their surface dip angles \n",
      "throughout the upper crust. %, and displaying typical seismic behavior.\n",
      "Under these assumptions, we create synthetic earthquake catalogs with\n",
      "Gutenberg-Richter and `characteristic' frequency-magnitude distributions\n",
      "based on each fault's geometry and slip rate, and then we calculate the\n",
      "probability of observing earthquakes on each, and any, fault over different\n",
      "durations of observation.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### LANF Slip, Mohr-Coulomb Failure Theory, and Seismicity"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Areas of the crust undergoing active extension are generally assumed\n",
      "to have a subvertical maximum compressive stress.  Mohr-Coulomb\n",
      "theory, as applied to the crust, predicts that a fault with a typical\n",
      "coefficient of friction for rocks (0.6--0.8) should lock up if it is\n",
      "oriented at an angle greater than 60$^\\circ$ to the maximum compressive stress,\n",
      "and new, optimally oriented faults should form \\citep{sibson1985}.  Therefore,\n",
      "for normal faults with dips less than 30$^\\circ$, either much lower\n",
      "fault friction or elevated pore fluid pressure is required for fault slip.\n",
      "\n",
      "Evidence for seismic slip on LANFs is sparse.  This is partly due to the\n",
      "ambiguity of the rupture plane in earthquake focal mechanisms, which have\n",
      "nodal planes 90$^\\circ$ apart, and therefore have a\n",
      "high-angle nodal plane if they have a low-angle nodal plane.  Without\n",
      "ancillary information indicating which plane slipped, searches of earthquake\n",
      "catalogs cannot yield unique results as to whether they contain LANF events.\n",
      "Several collections of normal fault earthquakes\n",
      "with known surface breaks \\citep{jackson1987, collettinisibson2001}, thereby\n",
      "resolving dip ambiguity, contain no low-angle events, though the collections\n",
      "are small ($\\le$ 25 events).  Some candidate events exist, but they are\n",
      "undersea \\citep[e.g.,][]{abers2001} or difficult to verify \\citep[e.g.,]\n",
      "[]{doser1987ancash}."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Potentiall Active LANFs"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Over the past decade or so, many field studies have found evidence for LANF \n",
      "activity in orogens throughout the world. These studies typically find arrays of \n",
      "Quaternary normal fault scarps on the fault traces and/or in the hanging walls \n",
      "of mapped or inferred low-angle detachment faults \\citep [e.g.,]\n",
      "[]{axen1999baja}. Some studies also have bedrock thermochronology data from the \n",
      "exhumed footwalls of the detachments that is suggestive of ongoing rapid \n",
      "exhumation \\citep [e.g.,][]{sundell2013lunggar}, although this data does not \n",
      "preclude a recent cessation of faulting. In some cases, additional evidence for \n",
      "LANF activity comes from geophysical data such as GPS geodesy \\citep [e.g.,]\n",
      "[]{hreinsdottir2009altotib} and seismic waves \\citep [e.g.,][]{doser1987ancash}.\n",
      "\n",
      "We have compiled all potentially active LANFs with known subareal\n",
      "fault traces from a thorough review of the literature, finding twenty\n",
      "total (Figure~\\ref{fig:lanf_map}).  We have then mapped the approximate fault\n",
      "traces into a GIS file (available at https://github.com/cossatot/LANF\\_gis), \n",
      "with metadata such as slip rate and source. Though the fault traces of many\n",
      "LANFs considered here are obscured by vegetation or agriculture, others\n",
      "display large fault scarps in Quaternary sediments, particularly those in the\n",
      "dry climates of Tibet \\citep[e.g.,][]{styron2013slr, kapp2005nqtl} and the\n",
      "western US \\citep[e.g.,][]{axen1999baja, hayman2003dv}, which are commonly\n",
      "interpreted as evidence for past seismic slip.  About half are in Tibet,\n",
      "consistent with hypotheses that LANFs and metamorphic core complexes\n",
      "form in areas of hot, thick crust \\citep [e.g.,][]{buck1991mcc}.  The\n",
      "rest are distributed through other areas of active continental\n",
      "extension: the North American Basin and Range, the Malay Archipelago,\n",
      "western Turkey, Italy, and Peru. Several of the most-commonly cited\n",
      "candidates for seismically active LANFs were not included because they\n",
      "do not have a clearly-defined, mappable fault trace, which is\n",
      "necessary for our earthquake likelihood calculations.  These include the \n",
      "submarine core complexes in the Woodlark Basin \\citep{abers2001}, the fault\n",
      "responsible for the 1995 Aigion, Greece earthquake \\citep{bernard1997}\n",
      "and other potential LANFs underneath the Gulf of Corinth, and the\n",
      "fault responsible for the 1952 Ancash, Peru earthquake\n",
      "\\citep{doser1987ancash}."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# show map and figure caption"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Likelihood of observing a LANF event\n",
      "### Earthquake likelihood on an individual LANF"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To estimate the likelihood of observing a significant earthquake on an\n",
      "individual LANF over some contiguous time window of length $t$ (in\n",
      "years), we perform a Monte Carlo simulation in which we create 4000\n",
      "synthetic time series of earthquakes, with unique values for fault\n",
      "geometry and slip rate for each time series. Then, for each time\n",
      "series we calculate the fraction of unique time windows of length $t$\n",
      "in which an earthquake as large or larger than a given magnitude\n",
      "occurs.  We take this value as the probability of observing an\n",
      "earthquake greater than or equal to moment magnitude $M$ over\n",
      "time period $t$, which we will refer to in general as $P(M,t)$.\n",
      "\n",
      "The geometry for each fault is estimated based on the length of the\n",
      "fault trace, the dip of the fault, and the estimated seismogenic\n",
      "thickness or fault locking depth in the area.  The fault is treated as\n",
      "planar for simplicity of calculations, even though the exposed\n",
      "footwalls of many detachment faults are highly corrugated.  We\n",
      "determine the fault length by measuring the approximate length of the\n",
      "mapped fault trace perpendicular to the assumed extension direction;\n",
      "for faults that change dip significantly along strike (e.g., the Dixie\n",
      "Valley fault), we only consider the low-angle segments of the fault.\n",
      "Values for the dip are taken from the literature in most cases, and\n",
      "measured of the dip of footwall triangular facets (interpreted as the exhumed\n",
      "fault plane) from SRTM data otherwise. In all cases, ranges of fault\n",
      "geometries are considered, encompassing the degree to which the values are \n",
      "known. The fault locking depth is assumed to be 10 km in the absence of other\n",
      "evidence (such as a geodetic study, \\citep[e.g.,][]{hreinsdottir2009altotib}).\n",
      "\n",
      "Slip rates of the 20 LANFs are similarly gathered from the\n",
      "literature if possible, or given broad ranges if not documented\n",
      "(e.g., 1--10 mm yr$^{-1}$).  In the Monte Carlo simulation, samples\n",
      "for slip rate and dip are drawn from uniform distributions defined\n",
      "by the maximum and minimum values.  Based on field observations,\n",
      "some faults have dip ranges that go above 30$^\\circ$; this study\n",
      "only considers slip on faults shallower than this, so for these\n",
      "faults, dip values are sampled from the minimum to 30$^\\circ$ and\n",
      "the resulting probabilities are then multiplied by the fraction of\n",
      "the dip range that is under 30$^\\circ$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Each earthquake synthetic earthquake sequence is generated by\n",
      "randomly sampling either 50,000 events from a tapered Gutenberg-Richter (GR)\n",
      "distribution with corner magnitude $M_c = 7.64$ and $\\beta = 0.65$\n",
      "(from values estimated by \\citet{birdkagan2004f_m} for continental\n",
      "rifts), or a 25,000 events from `characteristic' distribution. It is not\n",
      "certain which distribution more appropriately describes seismicity on a single\n",
      "LANF, though studies of many individual fault rupture histories suggests\n",
      "that the characteristic distribution is more accurate \\citep{hecker2013eqdist}.\n",
      "The smaller number of samples from the characteristic distribution is due to \n",
      "the increased computation time associated with a higher proportion of large \n",
      "events.  The samples are taken \n",
      "from an interval $M = [5.0, \\, M_{max}]$, where $M_{max}$ is calculated as the\n",
      "moment magnitude $M$ resulting from fault slip $D$ = 15 m over a fault of\n",
      "length $L$ cutting through a seismogenic thickness $z$ at dip $\\delta$,\n",
      "given the relations\n",
      "\n",
      "\\begin{equation}\n",
      " M_o = \\mu L z D \\,/ \\, \\sin \\delta \n",
      " \\end{equation}\n",
      "\n",
      "and\n",
      "\n",
      "\\begin{equation}\n",
      "M = 2/3 \\; \\log_{10} (M_o) - 6\n",
      "\\end{equation}\n",
      "\n",
      "where $\\mu = 30$ GPa is the shear modulus and $M_o$ is the seismic\n",
      "moment in N m \\citep{kagan2003pepi}.  The characteristic distribution has a\n",
      "large-magnitude mode corresponding to $D$ = 1.5 m on the fault, a typical\n",
      "slip distance for normal fault events \\citep[e.g.][]{wesnousky2008displacement}.\n",
      "The distributions are shown in Figure~\\ref{fig:fms}.\n",
      "\n",
      "These calculations rely on two important assumptions that warrant some\n",
      "discussion.  The first is that each earthquake ruptures the entire\n",
      "fault patch uniformly.  Though this is highly improbable fault behavior,\n",
      "the long-term statistical distribution of earthquake recurrence is \n",
      "insensitive to assumptions about slip distribution in individual events:\n",
      "if $n$ different, equal fault patches rupture independently, each \n",
      "requires $n$ times the interseismic strain accumulation time to rupture with\n",
      "some earthquake of magnitude $M$ versus a single fault rupturing uniformly; \n",
      "but on the whole, magnitude $M$ events would happen with the same long-term\n",
      "frequency.  The next assumption is that earthquakes are ordered randomly and\n",
      "separated by the time necessary for sufficient strain to accumulate for each\n",
      "earthquake to occur.  This means that foreshock and aftershock sequences\n",
      "and other types of event clustering are not taken into account.  However,\n",
      "the modal inter-event times for earthquakes $\\ge M \\,6$ or so are greater than\n",
      "a hundred years for many LANFs [supplemental figure xx], so the ordering of\n",
      "events does not impact the results, as this is longer than our maximum\n",
      "observation window.  Furthermore, any clustering resulting in\n",
      "event spacing less than the observation window would decrease $P(M,t)$, and\n",
      "we are concerned with calculating the maximum $P(M,t)$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Setting up the problem"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Import necessary modules"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#%pylab inline\n",
      "#%config InlineBackend.figure_format = 'retina'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib.pyplot import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "sys.path.append('../eq_stats')\n",
      "\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "import eq_stats as eqs\n",
      "import time\n",
      "from joblib import Parallel, delayed\n",
      "from itertools import chain"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Read in fault data table\n",
      "\n",
      "Makes a Pandas dataframe of fault data (length, slip rates, etc.)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = pd.read_csv('../data/lanf_stats.csv', index_col=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Define some variables to be used later"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n_cores = -3\n",
      "n_eq_samp = 3e4 # number of earthquakes in time series\n",
      "time_window = np.hstack( (1, np.arange(5, 105, step=5) ) ) # observation times\n",
      "mc_iters = 2e3 # number of Monte Carlo iterations\n",
      "mc_index = np.arange(mc_iters, dtype='int')\n",
      "mc_cols = ['dip', 'Ddot'] + [t for t in time_window]\n",
      "max_eq_slip = 15 #m\n",
      "char_eq_slip= 1.5 #m\n",
      "Mc = 7.64 # Hypothetical corner magnitude for Continental Rift Boundaries"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make list of minimum search magnitude $M_{min}$, and then make MultiIndex for Pandas dataframes"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "min_M_list = [5, 5.5, 6, 6.5, 7, 7.5]\n",
      "\n",
      "df_ind_tuples = [[i, M] for i in mc_index for M in min_M_list]\n",
      "df_multi_ind = pd.MultiIndex.from_tuples(df_ind_tuples, names=['mc_iter','M'])\n",
      "\n",
      "rec_int_bins = np.logspace(1, 5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Calculate $P(M,t)$ for Gutenberg-Richter frequency-magnitude distributions"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Define a function for Joblib Parallel to calculate probabilities for each iteration.\n",
      "\n",
      "Function is defined here so it can access all variables generated by script, not just passed variables.  This makes the code cleaner even if it's not very abstracted.\n",
      "\n",
      "Here is what this function does:\n",
      "\n",
      "- Get the dip, Ddot and maximum earthquake magnitude for each iteration.\n",
      "- Take this info and make the earthquake sequence:\n",
      "    - Take the max earthquake magnitude and make a frequency-magnitude distribution based on a Gutenburg-Richter exponential model.\n",
      "    - Take 50k samples from this distribution, \n",
      "- Make an earthquake time series form the EQ sequence\n",
      "    - Calculate the interseismic strain accumulation time for each event\n",
      "    - Separate each earthquake in the sequence with the appropriate number of years with no events.\n",
      "- Calculate the probability of observation\n",
      "    - Run a rolling maximum for each $t$ in [1, 5, 10, 15, ..., 95, 100]\n",
      "    - Calculate the observation probability above $M_{min}$ in [5, 5.5, 6, 6.5, 7, 7.5]\n",
      "- Calculate inter-event times for EQs $\\ge \\, M$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def calc_iter_probs(iter):\n",
      "    df_iter = fdf.loc[iter].copy()\n",
      "    df_iter['dip'] = mc_d['dip_samp'][iter]\n",
      "    df_iter['Ddot'] = mc_d['Ddot_samp'][iter]\n",
      "\n",
      "    # Generate EQ sample/sequence from F(M) dist.\n",
      "    m_vec = np.linspace(5, mc_d['max_M'][iter], num=1000)\n",
      "    fm_vec = eqs.F(m_vec, Mc=Mc)\n",
      "    M_samp = eqs.sample_from_pdf(m_vec, fm_vec, n_eq_samp)\n",
      "    Mo_samp = eqs.calc_Mo_from_M(M_samp)\n",
      "    \n",
      "    # Make time series of earthquakes, including no eq years\n",
      "    recur_int = eqs.calc_recurrence_interval(Mo=Mo_samp, \n",
      "                                             dip=mc_d['dip_samp'][iter],\n",
      "                                             slip_rate=mc_d['Ddot_samp'][iter],\n",
      "                                             L=params['L_km'],\n",
      "                                             z=params['z_km'])\n",
      "\n",
      "    cum_yrs = eqs.calc_cumulative_yrs(recur_int)\n",
      "    eq_series = eqs.make_eq_time_series(M_samp, cum_yrs)\n",
      "    \n",
      "    # calculate probability of observing EQ in time_window\n",
      "    for t in time_window:\n",
      "        roll_max = pd.rolling_max(eq_series, t)\n",
      "        df_iter[t] = (eqs.get_probability_above_value(roll_max, min_M_list)\n",
      "                      * mc_d['dip_frac'])\n",
      "        \n",
      "    # calculate histgrams of recurrence intervals\n",
      "    rec_int_counts_df = rec_int_df.loc[iter].copy()\n",
      "    for mm in np.array(min_M_list):\n",
      "        ints = np.diff( np.where(eq_series >= mm) )\n",
      "        rec_int_counts_df.loc[mm] = np.histogram(ints, bins=rec_int_bins)[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Iterate through the faults in the fault database, doing all the calculations for each.\n",
      "\n",
      "The setup of this for loop is basically this:\n",
      "\n",
      "- Make DataFrame for each fault.\n",
      "    - Columns are dip, Ddot, and observation time windows.\n",
      "    - Rows are values for each Monte Carlo iteration.  Values for time windows are calculated probabilities.\n",
      "    \n",
      "- Calculate maximum earthquake magnitude for each MC iteration.\n",
      "\n",
      "- Run the above 'calc_iter_probs' function (parallelized over the MC iterations) and concatenate the results"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Note! Not actually running the calcs here, as it's a looong process!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for fault in list(f.index):\n",
      "    fdf = pd.DataFrame(index=df_multi_ind, columns=mc_cols, dtype='float')\n",
      "    params = f.loc[fault]\n",
      "    mc_d = {}\n",
      "    mc_d['dip_samp'], mc_d['dip_frac'] = eqs.dip_rand_samp( params['dip_deg'], \n",
      "                                                         params['dip_err_deg'], \n",
      "                                                         mc_iters)\n",
      "\n",
      "    mc_d['Ddot_samp'] = eqs.Ddot_rand_samp(params['slip_rate_mm_a'],\n",
      "                                           params['sr_err_mm_a'], mc_iters)\n",
      "\n",
      "    mc_d['max_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], \n",
      "                                                   z=params['z_km'], \n",
      "                                                   dip=mc_d['dip_samp'], \n",
      "                                                   D=max_eq_slip)\n",
      "\n",
      "    mc_d['max_M'] = eqs.calc_M_from_Mo(mc_d['max_Mo'])\n",
      "    \n",
      "    rec_int_df = pd.DataFrame(columns = rec_int_bins[1:],\n",
      "                              index=df_multi_ind, dtype='float')\n",
      "    t0 = time.time()\n",
      "    out_list = Parallel(n_jobs=-3)( delayed( calc_iter_probs)(ii) \n",
      "                                    for ii in mc_index)\n",
      "    print 'done with', fault, 'parallel calcs in {} s'.format((time.time()-t0))\n",
      "    for ii in mc_index:\n",
      "        fdf.loc[ii][:] = out_list[ii][0]\n",
      "        rec_int_df.loc[ii][:] = out_list[ii][1]\n",
      "        \n",
      "    fdf.to_csv('../results/{}_test.csv'.format(fault))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with kongur_shan parallel calcs in 51.947163105 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with leo_pargil parallel calcs in 418.76543808 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with gurla_mandhata parallel calcs in 113.975022078 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with s_lunggar parallel calcs in 241.274301052 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with n_lunggar parallel calcs in 259.45785284 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with pqx_qingdu parallel calcs in 836.985168934 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with pqx_n parallel calcs in 917.926791191 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with nqtl parallel calcs in 113.513681889 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with alto_tiberina parallel calcs in 233.532536983 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with death_valley parallel calcs in 201.206351042 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with panamint_valley parallel calcs in 102.183755159 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with laguna_salada parallel calcs in 217.666782141 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with dixie_valley parallel calcs in 1027.02227688 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with dayman_dome parallel calcs in 103.56942296 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with papangeo parallel calcs in 183.855463028 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with tokorondo parallel calcs in 78.5454099178 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with cordillera_blanca parallel calcs in 66.0141530037 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with kuzey parallel calcs in 268.01406312 s\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "done with guney parallel calcs in 158.371592999 s\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Calculate $P(M,t)$ for faults with characteristic frequency-magnitude distributions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def calc_iter_probs(ii):\n",
      "    df_iter = fdf.loc[ii].copy()\n",
      "    df_iter['dip'] = mc_d['dip_samp'][ii]\n",
      "    df_iter['Ddot'] = mc_d['Ddot_samp'][ii]\n",
      "\n",
      "    # Generate EQ sample/sequence from F(M) dist.\n",
      "    m_vec = np.linspace(5, mc_d['max_M'][ii], num=1000)\n",
      "    fm_vec = eqs.F_char(m_vec, Mc=Mc, char_M=mc_d['char_M'][ii])\n",
      "    M_samp = eqs.sample_from_pdf(m_vec, fm_vec, n_eq_samp)\n",
      "    Mo_samp = eqs.calc_Mo_from_M(M_samp)\n",
      "    \n",
      "    # Make time series of earthquakes, including no eq years\n",
      "    recur_int = eqs.calc_recurrence_interval(Mo=Mo_samp, \n",
      "                                             dip=mc_d['dip_samp'][ii],\n",
      "                                             slip_rate=mc_d['Ddot_samp'][ii],\n",
      "                                             L=params['L_km'],\n",
      "                                             z=params['z_km'])\n",
      "\n",
      "    cum_yrs = eqs.calc_cumulative_yrs(recur_int)\n",
      "    eq_series = eqs.make_eq_time_series(M_samp, cum_yrs)\n",
      "    \n",
      "    # calculate probability of observing EQ in time_window\n",
      "    for t in time_window:\n",
      "        roll_max = pd.rolling_max(eq_series, t)\n",
      "        df_iter[t] = (eqs.get_probability_above_value(roll_max, min_M_list)\n",
      "                      * mc_d['dip_frac'])\n",
      "        \n",
      "    # calculate histgrams of recurrence intervals\n",
      "    rec_int_counts_df = rec_int_df.loc[ii].copy()\n",
      "    for mm in np.array(min_M_list):\n",
      "        ints = np.diff( np.where(eq_series >= mm) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for fault in list(f.index):\n",
      "    fdf = pd.DataFrame(index=df_multi_ind, columns=mc_cols, dtype='float')\n",
      "    params = f.loc[fault]\n",
      "    mc_d = {}\n",
      "    mc_d['dip_samp'], mc_d['dip_frac'] = eqs.dip_rand_samp( params['dip_deg'], \n",
      "                                                         params['dip_err_deg'], \n",
      "                                                         mc_iters)\n",
      "\n",
      "    mc_d['Ddot_samp'] = eqs.Ddot_rand_samp(params['slip_rate_mm_a'],\n",
      "                                           params['sr_err_mm_a'], mc_iters)\n",
      "\n",
      "    mc_d['max_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], \n",
      "                                                   z=params['z_km'], \n",
      "                                                   dip=mc_d['dip_samp'], \n",
      "                                                   D=max_eq_slip)\n",
      "\n",
      "    mc_d['max_M'] = eqs.calc_M_from_Mo(mc_d['max_Mo'])\n",
      "    \n",
      "    mc_d['char_Mo'] = eqs.calc_Mo_from_fault_params(L=params['L_km'], \n",
      "                                                    z=params['z_km'], \n",
      "                                                    dip=mc_d['dip_samp'], \n",
      "                                                    D=char_eq_slip)\n",
      "    \n",
      "    mc_d['char_M'] = eqs.calc_M_from_Mo(mc_d['char_Mo'])\n",
      "    \n",
      "    rec_int_df = pd.DataFrame(columns = rec_int_bins[1:],\n",
      "                              index=df_multi_ind, dtype='float')\n",
      "    t0 = time.time()\n",
      "    out_list = Parallel(n_jobs=n_cores)( delayed( calc_iter_probs)(ii) \n",
      "                                    for ii in mc_index)\n",
      "    print 'done with', fault, 'parallel calcs in {} s'.format((time.time()-t0))\n",
      "    for ii in mc_index:\n",
      "        fdf.loc[ii][:] = out_list[ii][0]\n",
      "        rec_int_df.loc[ii][:] = out_list[ii][1]\n",
      "        \n",
      "    fdf.to_csv('../results/{}_char.csv'.format(fault))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Examining individual fault results"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Load datasets into Pandas dataframes\n",
      "\n",
      "test with one:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pv = pd.read_csv('../results/panamint_valley_gr.csv', index_col=[0])\n",
      "pv_c = pd.read_csv('../results/panamint_valley_char.csv', index_col=[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tw_xarray = np.tile(time_window, (mc_iters,1))\n",
      "\n",
      "tw_xarray.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "(2000, 21)"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tw_cols = list(time_window.astype('str'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "make some plotting functions"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plot_probs(df, lw=0.5, alpha=0.2):\n",
      "\n",
      "    df_5 = df[df.M==5]\n",
      "    df_55 = df[df.M==5.5]\n",
      "    df_6 = df[df.M==6]\n",
      "    df_65 = df[df.M==6.5]\n",
      "    df_7 = df[df.M==7]\n",
      "    df_75 = df[df.M==7.5]\n",
      "\n",
      "    aa = plot(tw_xarray.T, df_5[tw_cols].T, 'b', \n",
      "              lw=lw, alpha=alpha)\n",
      "    bb = plot(tw_xarray.T, df_55[tw_cols].T, 'c', \n",
      "              lw=lw, alpha=alpha)\n",
      "    cc = plot(tw_xarray.T, df_6[tw_cols].T, 'g', \n",
      "              lw=lw, alpha=alpha)\n",
      "    dd = plot(tw_xarray.T, df_65[tw_cols].T, 'gold', \n",
      "              lw=lw, alpha=alpha)\n",
      "    ee = plot(tw_xarray.T, df_7[tw_cols].T, 'orange', \n",
      "              lw=lw, alpha=alpha)\n",
      "    ff = plot(tw_xarray.T, df_75[tw_cols].T, 'r', \n",
      "              lw=lw, alpha=alpha)\n",
      "\n",
      "    return aa, bb, cc, dd, ee, ff"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def pmf_x_secs(df, tw, n_bins=101, hist_type='stepfilled', alpha=0.5):\n",
      "    hist_bins = np.linspace(0, 1, num=n_bins)\n",
      "    tw = str(tw)\n",
      "    \n",
      "    aa = df[df.M==5][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                              color='b', alpha=alpha)\n",
      "    bb = df[df.M==5.5][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                                color='c', alpha=alpha)\n",
      "    cc = df[df.M==6][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                              color='g', alpha=alpha)\n",
      "    dd = df[df.M==6.5][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                                color='gold', alpha=alpha)\n",
      "    ee = df[df.M==7][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                              color='orange', alpha=alpha)\n",
      "    ff = df[df.M==7.5][tw].hist(bins=hist_bins, histtype=hist_type, \n",
      "                                color='red', alpha=alpha)\n",
      "    \n",
      "    return aa, bb, cc, dd, ee, ff"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#plot_probs(pv)\n",
      "\n",
      "#xlabel('observation time, yrs', fontsize=18)\n",
      "#ylabel('probability', fontsize=18)\n",
      "#title('Figure 1: Likelihood of observing an earthquake \\n'\n",
      "#      + 'larger than M on the Kongur Shan fault', fontsize=21)\n",
      "#show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tw=35\n",
      "nbins=1001\n",
      "x_sec_lims = [0,0.4,0,400]\n",
      "\n",
      "figure()\n",
      "\n",
      "subplot2grid((3,2), (0,0), colspan=2)\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(pv, lw=0.2, alpha=0.2)\n",
      "\n",
      "subplot2grid((3,2), (1,0), colspan=2)\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(pv_c, lw=0.2, alpha=0.2)\n",
      "\n",
      "subplot2grid((3,2), (2,0))\n",
      "pmf_x_secs(pv, tw, n_bins=nbins)\n",
      "axis(x_sec_lims)\n",
      "\n",
      "subplot2grid((3,2), (2,1))\n",
      "pmf_x_secs(pv_c, tw, n_bins=nbins)\n",
      "axis([0,0.1,0,400])\n",
      "\n",
      "show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The results for faults with a GR frequency-magnitude distribution\n",
      "indicate that it is unlikely that any individual fault would have an earthquake\n",
      "greater than $M \\, 5$ in any modeled observation time (up to 100 years).  \n",
      "As an example, the results for the Panamint Valley fault are shown in\n",
      "Figure~\\ref{fig:pv}a;\n",
      "this fault has the highest $P(M,t)$ of any of the well-studied LANFs.  The\n",
      "probability of observing a $\\ge M \\, 6.0$ event on the Panamint Valley fault \n",
      "is about 0.5 for $t$ = 100 years, and about 0.15 for $t$ = 35 years, which is\n",
      "the length of the Global CMT catalog. As expected given the GR distribution,\n",
      "the $P(M,t)$ is much higher for smaller, more frequent events than for larger\n",
      "events.  \n",
      "\n",
      "The results for faults with a characteristic frequency-magnitude distribution\n",
      "yield much lower $P(M,t)$ for small to moderate events, but for large events,\n",
      "$P(M,t)$ is higher (Figure~\\ref{fig:pv}b,d); this is because the earthquake\n",
      "sequences are dominated by large, infrequent events, so the inter-event times\n",
      "for moderate events are several times greater. For the Panamint Valley fault,\n",
      "$P(M,t)$ for $ \\ge M \\, 5$ is about 0.07 (versus 0.25 for the GR distribution),\n",
      "but $P(M\\ge 7, t=35)$ is around 0.025 (versus essentially zero for the GR\n",
      "distribution).  As the characteristic distribution likely better represents\n",
      "earthquakes on an individual large fault, these results suggest that is \n",
      "very unlikely that we would expect to capture any significant seismicity\n",
      "on an single LANF in the moment tensor catalogs. A similar conclusion was \n",
      "found by \\citet{wernicke1995seis} based on a simple calculation, \n",
      "assuming perfectly repeating large earthquakes on an idealized fault. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Earthquake Likelihood on All LANFs\n",
      "\n",
      "To calculate the probability of observing an earthquake over the time\n",
      "window on any of the LANFs studied, we first assume that seismicity on\n",
      "each fault is independent and uncorrelated with seismicity on all\n",
      "other faults. This assumption is likely true for most faults, but may\n",
      "not be true for some proximal faults (such as the North and South\n",
      "Lunggar detachments, or the Papangeo and Tokorondo detachments),\n",
      "though it is unclear how to determine how these faults may interact \n",
      "such that an appropriate conditional or joint probability may be calculated. \n",
      "We determine the probability for each time window and minimum magnitude\n",
      "with the equation\n",
      "\n",
      "\\begin{equation}\n",
      "P_{AT \\, or \\, LP\\, \\ldots \\, or \\, DV} = 1 - (Q_{AT} \\cdot Q_{LP} \\cdot \\ldots \\, \\cdot Q_{DV})\n",
      "%\\label{ProbUnion}\n",
      "\\end{equation}\n",
      "\n",
      "where $P_{AT}$ is the probability of observing an earthquake on a\n",
      "single LANF (e.g., the Alto-Tiberina fault), and $Q_{AT} = 1 -\n",
      "P_{AT}$. Equation (\\ref{ProbUnion}) is the union of probabilities for\n",
      "non mutually exclusive random events."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Calculating these probabilities"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "####now load the rest:\n",
      "\n",
      "Gutenberg-Richter models"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ks = pd.read_csv('../results/kongur_shan_gr.csv', index_col=0)\n",
      "gm = pd.read_csv('../results/gurla_mandhata_gr.csv', index_col=0)\n",
      "slr = pd.read_csv('../results/s_lunggar_gr.csv', index_col=0)\n",
      "nlr = pd.read_csv('../results/n_lunggar_gr.csv', index_col=0)\n",
      "pqxn = pd.read_csv('../results/pqx_n_gr.csv', index_col=0)\n",
      "pqxq = pd.read_csv('../results/pqx_qingdu_gr.csv', index_col=0)\n",
      "lp = pd.read_csv('../results/leo_pargil_gr.csv', index_col=0)\n",
      "nqtl = pd.read_csv('../results/nqtl_gr.csv', index_col=0)\n",
      "at = pd.read_csv('../results/alto_tiberina_gr.csv', index_col=0)\n",
      "dv = pd.read_csv('../results/death_valley_gr.csv', index_col=0)\n",
      "cd = pd.read_csv('../results/canada_david_gr.csv', index_col=0)\n",
      "sd = pd.read_csv('../results/sevier_desert_gr.csv', index_col=0)\n",
      "dxv = pd.read_csv('../results/dixie_valley_gr.csv', index_col=0)\n",
      "dd = pd.read_csv('../results/dayman_dome_gr.csv', index_col=0)\n",
      "pp = pd.read_csv('../results/papangeo_gr.csv', index_col=0)\n",
      "tk = pd.read_csv('../results/tokorondo_gr.csv', index_col=0)\n",
      "cb = pd.read_csv('../results/cordillera_blanca_gr.csv', index_col=0)\n",
      "kz = pd.read_csv('../results/kuzey_gr.csv', index_col=0)\n",
      "gn = pd.read_csv('../results/guney_gr.csv', index_col=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "characteristic models"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ks_c = pd.read_csv('../results/kongur_shan_char.csv', index_col=0)\n",
      "gm_c = pd.read_csv('../results/gurla_mandhata_char.csv', index_col=0)\n",
      "slr_c = pd.read_csv('../results/s_lunggar_char.csv', index_col=0)\n",
      "nlr_c = pd.read_csv('../results/n_lunggar_char.csv', index_col=0)\n",
      "pqxn_c = pd.read_csv('../results/pqx_n_char.csv', index_col=0)\n",
      "pqxq_c = pd.read_csv('../results/pqx_qingdu_char.csv', index_col=0)\n",
      "lp_c = pd.read_csv('../results/leo_pargil_char.csv', index_col=0)\n",
      "nqtl_c = pd.read_csv('../results/nqtl_char.csv', index_col=0)\n",
      "at_c = pd.read_csv('../results/alto_tiberina_char.csv', index_col=0)\n",
      "dv_c = pd.read_csv('../results/death_valley_char.csv', index_col=0)\n",
      "cd_c = pd.read_csv('../results/canada_david_char.csv', index_col=0)\n",
      "sd_c = pd.read_csv('../results/sevier_desert_char.csv', index_col=0)\n",
      "dxv_c = pd.read_csv('../results/dixie_valley_char.csv', index_col=0)\n",
      "dd_c = pd.read_csv('../results/dayman_dome_char.csv', index_col=0)\n",
      "pp_c = pd.read_csv('../results/papangeo_char.csv', index_col=0)\n",
      "tk_c = pd.read_csv('../results/tokorondo_char.csv', index_col=0)\n",
      "cb_c = pd.read_csv('../results/cordillera_blanca_char.csv', index_col=0)\n",
      "kz_c = pd.read_csv('../results/kuzey_char.csv', index_col=0)\n",
      "gn_c = pd.read_csv('../results/guney_char.csv', index_col=0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "make list of faults"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f_list = ['ks', 'lp', 'gm', 'slr', 'nlr', 'pqxn', 'pqxq', 'nqtl', 'at',\n",
      "          'dv', 'pv', 'cd', 'sd', 'dxv', 'dd', 'pp', 'tk', 'cb', 'kz', 'gn']\n",
      "\n",
      "fault_list = [ks, lp, gm, slr, nlr, pqxn, pqxq, nqtl, at,\n",
      "              dv, pv, cd, sd,  dxv, dd, pp, tk, cb, kz, gn]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 37
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fault_c_list = [ks_c, lp_c, gm_c, slr_c, nlr_c, pqxn_c, pqxq_c, nqtl_c, at_c,\n",
      "              dv_c, pv_c, cd_c, sd_c, dxv_c, dd_c, pp_c, tk_c, cb_c, kz_c, gn_c]\n",
      "\n",
      "fc_list = ['ks_c', 'lp_c', 'gm_c', 'slr_c', 'nlr_c', 'pqxn_c', \n",
      "                'pqxq_c', 'nqtl_c', 'at_c', 'dv_c', 'pv_c', 'cd_c', \n",
      "                'sd_c', 'dxv_c', 'dd_c', 'pp_c', 'tk_c', 'cb_c', 'kz_c', 'gn_c']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Calculate $q$, probability of *not* observing an earthquake, for each fault"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "q = {}\n",
      "\n",
      "for i, f_ in enumerate(f_list):\n",
      "    q[f_] = fault_list[i].copy()\n",
      "    q[f_][tw_cols] = 1 - q[f_][tw_cols]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "qc = {}\n",
      "\n",
      "for i, f_ in enumerate(fc_list):\n",
      "    qc[f_] = fault_c_list[i].copy()\n",
      "    qc[f_][tw_cols] = 1 - qc[f_][tw_cols]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 43
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### Now estimate joint probabilities\n",
      "\n",
      "Make list of columns to retain in final dataframe"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prob_cols = list( chain( *['M', list(tw_cols)] ) )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 39
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "calculate $p_{joint}$ as $1-(q \\cdot q \\cdot q...)$ and make final dataframe"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "all_prob_gr = 1 - np.product([qq[tw_cols] for qq in  q.values()])\n",
      "\n",
      "all_prob_gr['M'] = ks['M']\n",
      "\n",
      "all_prob_gr = all_prob_gr.reindex_axis(prob_cols, axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "all_c_prob = 1- np.product([qq[tw_cols] for qq in qc.values()])\n",
      "\n",
      "all_c_prob['M'] = ks['M']\n",
      "\n",
      "all_c_prob = all_c_prob.reindex_axis(prob_cols, axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tw=35\n",
      "nbins=501\n",
      "x_sec_lims = [0,0.4,0,400]\n",
      "\n",
      "figure()\n",
      "\n",
      "subplot2grid((3,2), (0,0), colspan=2)\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(all_prob_gr, lw=0.1, alpha=0.1)\n",
      "\n",
      "subplot2grid((3,2), (1,0), colspan=2)\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(all_c_prob, lw=0.1, alpha=0.1)\n",
      "\n",
      "subplot2grid((3,2), (2,0))\n",
      "pmf_x_secs(all_prob_gr, tw, n_bins=nbins)\n",
      "#axis(x_sec_lims)\n",
      "\n",
      "subplot2grid((3,2), (2,1))\n",
      "pmf_x_secs(all_c_prob, tw, n_bins=nbins)\n",
      "#axis([0,0.1,0,400])\n",
      "xlim([0,0.8])\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tw=35\n",
      "nbins=501\n",
      "x_sec_lims = [0,0.4,0,400]\n",
      "\n",
      "figure()\n",
      "\n",
      "subplot2grid((2,2), (0,0))\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(all_prob_gr, lw=0.1, alpha=0.1)\n",
      "\n",
      "subplot2grid((2,2), (0,1))\n",
      "axvline(tw, lw=0.5, color='grey')\n",
      "plot_probs(all_c_prob, lw=0.1, alpha=0.1)\n",
      "\n",
      "subplot2grid((2,2), (1,0))\n",
      "pmf_x_secs(all_prob_gr, tw, n_bins=nbins)\n",
      "#axis(x_sec_lims)\n",
      "\n",
      "subplot2grid((2,2), (1,1))\n",
      "pmf_x_secs(all_c_prob, tw, n_bins=nbins)\n",
      "#axis([0,0.1,0,400])\n",
      "xlim([0,0.8])\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}