started stripping down risk-and-reliability to create template

book/notebooks/Untitled.ipynb

@@ -0,0 +1,153 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "9d7ffe9f-441f-4119-b5e4-0e872922a78d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from sympy import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "406dc83e-b114-46be-baf8-073a72b47a36",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t, a = symbols('t a')\n",
+    "theta = Function('theta')(t)\n",
+    "r = Function('r')(t)\n",
+    "U_c = Function('U_c')(t)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "df8fe246-6986-4ffd-9c39-6bd259176e42",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle - \\frac{a^{2} \\operatorname{U_{c}}{\\left(t \\right)} \\cos{\\left(\\theta{\\left(t \\right)} \\right)}}{r{\\left(t \\right)}}$"
+      ],
+      "text/plain": [
+       "-a**2*U_c(t)*cos(theta(t))/r(t)"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "phi = - (a**2*U_c)/r * cos(theta)\n",
+    "phi"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "a6e55454-52c6-493b-a753-f901cee357bc",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{a^{2} \\operatorname{U_{c}}{\\left(t \\right)} \\sin{\\left(\\theta{\\left(t \\right)} \\right)} \\frac{d}{d t} \\theta{\\left(t \\right)}}{r{\\left(t \\right)}} + \\frac{a^{2} \\operatorname{U_{c}}{\\left(t \\right)} \\cos{\\left(\\theta{\\left(t \\right)} \\right)} \\frac{d}{d t} r{\\left(t \\right)}}{r^{2}{\\left(t \\right)}} - \\frac{a^{2} \\cos{\\left(\\theta{\\left(t \\right)} \\right)} \\frac{d}{d t} \\operatorname{U_{c}}{\\left(t \\right)}}{r{\\left(t \\right)}}$"
+      ],
+      "text/plain": [
+       "a**2*U_c(t)*sin(theta(t))*Derivative(theta(t), t)/r(t) + a**2*U_c(t)*cos(theta(t))*Derivative(r(t), t)/r(t)**2 - a**2*cos(theta(t))*Derivative(U_c(t), t)/r(t)"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "phi.diff(t)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "52970004-5d18-4905-a61c-01c5d9a7f81c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{a^{2} \\operatorname{U_{c}}{\\left(t \\right)} \\cos{\\left(\\theta{\\left(t \\right)} \\right)}}{r^{2}{\\left(t \\right)}}$"
+      ],
+      "text/plain": [
+       "a**2*U_c(t)*cos(theta(t))/r(t)**2"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "phi.diff(r)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "3bc6b6e8-cdba-4eb9-b754-e07cc74c6063",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/latex": [
+       "$\\displaystyle \\frac{a^{2} \\operatorname{U_{c}}{\\left(t \\right)} \\sin{\\left(\\theta{\\left(t \\right)} \\right)}}{r^{2}{\\left(t \\right)}}$"
+      ],
+      "text/plain": [
+       "a**2*U_c(t)*sin(theta(t))/r(t)**2"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "phi.diff(theta)/r"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bcae9e54-d38d-4128-b872-e4bcb11d19c5",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

book/notebooks/contamination/exercise_setup.py

@@ -0,0 +1,247 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.special import erfc
+from scipy import stats
+import seaborn as sns
+
+def process_negative_values(q, n, a):
+    """Prevents problems with negative values in sqrt of C_C0 function."""
+    N_q = sum(q<0)
+    N_n = sum(n<0)
+    N_a = sum(a<0)
+    set_value = 0.001
+    if N_q > 0:
+        for i in range(len(q)):
+            if q[i]<0:
+                q[i] = set_value
+        print(f'q has {N_q} negative values that have been set to {set_value:.1e}.')
+    if N_n > 0:
+        for i in range(len(n)):
+            if n[i]<0:
+                n[i] = set_value
+        print(f'n has {N_n} negative values that have been set to {set_value:.1e}.')
+    if N_a > 0:
+        for i in range(len(a)):
+            if a[i]<0:
+                a[i] = set_value
+        print(f'a has {N_a} negative values that have been set to {set_value:.1e}.')
+    return q, n, a
+
+class RandomVariableDistribution():
+    """Define a N or LN probability to use in MUDE week 2.7 exercise.
+    
+    Usage:
+    -----
+    Create an instance of the class
+        by default, the random variable is ~N(0,1)
+    Specify distribution using 1st and 2nd moments for N or LN:
+        --> use method define_distribution to change it!
+             (the method prints a summary of the change)
+        for example, copy/paste the following:
+            test = RandomVariableDistribution()
+            test.print_status()
+            test.define_distribution(mean=5, stddev=1, dist_type='LN')
+    An instance will include the following attributes:
+        - distribution moment and (scipy) parameters shape, location, scale
+        - dist_object: distribution object (scipy.stats)
+        - get_sample: method for getting random samples from the distribution
+        - plot_pdf_or_cdf:
+    """
+    
+    def __init__(self):
+        self.mean = 0.
+        self.stddev = 1.
+        self.dist_type = 'N'
+        
+    def print_status(self):
+        """Returns current parameter values and distribution type."""
+        print(f'The mean is {self.mean:.3e}')
+        print(f'The standard deviation is {self.stddev:.3e}')
+        if self.dist_type == 'N':
+            distribution_name = 'Normal'
+        elif self.dist_type == 'LN':
+            distribution_name = 'Lognormal'
+        print('The distribution is ' + distribution_name)
+        
+    def define_distribution(self,
+                            mean=None,
+                            stddev=None,
+                            dist_type=None,
+                            print_new_dist=False):
+        """Compute distribution parameters and RV object from moments.
+        
+        Usage:
+        -----
+        Optional arguments will override current dist:
+            --> default distribution is standard normal, ~N(0,1)
+            mean: float
+            stddev: float
+            dist_type: 'N' or 'LN'
+            print_new_dist=True to print distribution changes.
+        """
+        if mean is None:
+            mean = self.mean
+        if stddev is None:
+            stddev = self.stddev
+        if dist_type is None:
+            dist_type = self.dist_type
+        if dist_type == 'N':
+            shape = None
+            location = mean
+            scale = stddev
+            dist_object = stats.norm(loc=location, scale=scale)
+        elif dist_type == 'LN':
+            """
+            Finds distribution parameters of LN-distributed X,
+            which are equivalent to the mean and standard deviation
+            of the log of X: Y = ln(X), then converts to the
+            parameterization used by scipy.stats.lognorm.
+            """
+            stddev_y = np.sqrt(np.log(1 + (stddev/mean)**2))
+            mean_y = np.log(mean) - 0.5*stddev_y**2
+            
+            shape = stddev_y
+            location = 0
+            scale = np.exp(mean_y)
+            dist_object = stats.lognorm(s=shape,
+                                           loc=location,
+                                           scale=scale)
+        else:
+            raise ValueError("You can only use N or LN for self.dist_type.")
+        self.mean = mean
+        self.stddev = stddev
+        self.dist_type = dist_type
+        self.shape = shape
+        self.location = location
+        self.scale = scale
+        self.dist_object = dist_object
+        self.get_sample = lambda size : self.dist_object.rvs(size=size)
+        if print_new_dist:
+            self.print_status()
+        
+    def plot_pdf_or_cdf(self, plot_type='pdf', plot_limits=None):
+        """Plots the pdf or cdf (specify choice with plot_type).
+        
+        Usage:
+        -----
+        No inputs required, but there are optional inputs:
+            plot_type='cdf' to plot cdf (pdf is default value)
+            plot_limits: list with 2 values for x-axis
+        """
+        if plot_limits is None:
+            if self.dist_type == 'N':
+                plot_limits = [self.mean - 3*self.stddev,
+                               self.mean + 3*self.stddev]
+            elif self.dist_type == 'LN':
+                'Avoid negative or zero value for lognormal dist.'
+                lower_value = self.mean - 3*self.stddev
+                if lower_value <= 0:
+                    lower_value = 0.01
+                plot_limits = [lower_value,
+                               self.mean + 4*self.stddev]
+            else:
+                raise ValueError("You can only use N or LN for self.dist_type.")
+
+        plot_x = np.linspace(plot_limits[0], plot_limits[1], 1000)
+        if plot_type == 'pdf':
+            plot_y = self.dist_object.pdf(plot_x)
+            y_label = 'PDF [-]'
+        elif plot_type == 'cdf':
+            plot_y = self.dist_object.cdf(plot_x)
+            y_label = 'CDF [-]'
+        else:
+            raise ValueError("You must specify pdf or cdf for plot_type.")
+        fontsize = 15
+        plt.figure()
+        plt.plot(plot_x, plot_y, 'k-', linewidth=3)
+        plt.xlabel('Random Variable, X', fontsize=fontsize)
+        plt.ylabel(y_label, fontsize=fontsize)
+        plt.show()
+        
+def pdf_of_function_of_RV(cases, case_C=True, case_name=None, plot_limits=None):
+    """Plot the PDF of a sampled function of random variables (C_C0 or Z).
+    
+    Usage:
+    -----
+        cases: a list of arrays
+        case_C=True (default) assumes C_C0 is given. False is for Z.
+        case_name: list of strings with the labels for each array.
+                   e.g., ['Original', 'High sigma_a', 'Low sigma_a']
+                       (default will return Case i for all arrays)
+        plot_limits: list of two values that override default limits
+                     for x-axis.
+    """
+    bw_adjust = 3 # adjusts smoothness of Seaborn kernel density contours
+    labels = []
+    
+    plt.figure()
+    
+    for i in range(len(cases)):
+        if case_name == None:
+            labels.append(f'Case {i+1}')
+        else:
+            labels.append(case_name[i])
+        sns.kdeplot(cases[i], label=labels[i], ax=plt.gca(), bw_adjust=bw_adjust)
+    plt.legend()
+    if case_C:
+        plt.xlabel('Concentration ratio, $C/C_0$ [$-$]')
+        if plot_limits == None:
+            plt.xlim(0,1)
+        else:
+            plt.xlim(plot_limits)
+    else:
+        plt.xlabel('Limit-state function, $Z=\mathrm{threshold}-C/C_0$ [$-$]')
+        if plot_limits == None:
+            plt.xlim(-2,1)
+        else:
+            plt.xlim(plot_limits)
+    plt.ylabel('PDF [$]$')
+    plt.grid()
+    plt.show()
+    
+def monte_carlo(Z):
+    """Monte Carlo analysis for sampled function of random variables.
+    
+    Returns a summary of the analysis and a plot of pf and c.o.v..
+    The only input is an array of the sampled limit state-function.
+    Assumes that failure is when the limit-state function Z<0.
+    """
+    N = len(Z)
+    Nf = np.sum(Z<0)
+    
+    if Nf == 0:
+        raise ValueError('No failure realizations in the sample!')
+    elif Nf == N:
+        raise ValueError('All realizations in the sample failed!')
+    
+    pf = Nf/N
+    
+    N_f_c = np.cumsum(Z<0)
+    N_i = np.arange(1,len(N_f_c)+1,1)
+    p_f_c = N_f_c / N_i
+    with np.errstate(divide='ignore'):
+        V = 1/(np.sqrt(N_i*p_f_c))
+    
+    print(f'MONTE CARLO RESULTS\n'
+          f'--------------------------------------------------\n'
+          f'    N = {N}, number of samples\n'
+          f'   Nf = {Nf}, number of failed realizations\n'
+          f'   pf = {pf:.3e}, estimated failure probability\n'
+          f'    V = {V[-1]:.3f}, c.o.v. of estimated pf\n')
+    
+    plt.loglog(N_i,p_f_c)
+    plt.xlim(10,N)
+    plt.ylim(10**(np.floor(np.log10(pf)) - 1),1)
+    plt.grid()
+    plt.xlabel('Realization No.')
+    plt.ylabel('Estimated failure probability, $\hat{p}_f$')
+    plt.show()
+    
+    plt.plot(N_i,V)
+    plt.xscale('log')
+    plt.grid()
+    plt.xlabel('time step')
+    plt.ylabel('c.o.v.')
+    plt.show()
+    
+    return Nf,pf

book/notebooks/flood-risk/functions.py

@@ -0,0 +1,163 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+
+def find_scenarios_with_segment_i(segments_fixed, df):
+    """Given a list of segments, return set of dependent scenarios."""
+    idx = []
+    segments_fixed = set(segments_fixed)
+    for i, row in enumerate(df['failed_segments'].values):
+        segments_row = set([int(seg) for seg in row.split(",")])
+        intersection = list(segments_fixed.intersection(segments_row))
+        if len(intersection) > 0:
+            idx.append(True)
+        else:
+            idx.append(False)
+    segment_table = df.iloc[idx]
+    return segment_table
+
+def assemble_analysis_cases(df):
+    """From scenario dataframe collect damage and fatality information."""
+    econ_damage_case_1 = df.iloc[:,3] #  Unexpected, no evacuation
+    econ_damage_case_2 = df.iloc[:,5] #  Unexpected, unorganised evacuation
+    econ_damage_case_3 = df.iloc[:,7] #  Expected, unorganised evacuation
+    econ_damage_case_4 = df.iloc[:,9] #  Expected, organised evacuation
+
+    fatalities_case_1 = df.iloc[:,4]
+    fatalities_case_2 = df.iloc[:,6]
+    fatalities_case_3 = df.iloc[:,8]
+    fatalities_case_4 = df.iloc[:,10]
+
+    scenario_prob = df.iloc[:,1] #  scenario probability
+
+    cases = {'damage': [econ_damage_case_1, econ_damage_case_2, econ_damage_case_3, econ_damage_case_4],
+             'fatality': [fatalities_case_1, fatalities_case_2, fatalities_case_3, fatalities_case_4]}
+    return cases
+
+def exceedance_prob(cases, scenario_prob):
+    """Compute exceedance probability."""
+    case_1 = cases[0]
+    case_2 = cases[1]
+    case_3 = cases[2]
+    case_4 = cases[3]
+    
+    key = case_1.index
+    
+    # Combine Evacuation Scenarios into new Dataframe
+    scenarios = [pd.DataFrame({'Nr': key, 'N/D': case_1}),
+                 pd.DataFrame({'Nr': key, 'N/D': case_2}),
+                 pd.DataFrame({'Nr': key, 'N/D': case_3}),
+                 pd.DataFrame({'Nr': key, 'N/D': case_4})]
+    
+    scen_combi = pd.concat(scenarios, ignore_index=True)
+    
+    # Sort values by key ('Nr')
+    scen_combi_sorted = scen_combi.sort_values(by=['Nr']) # useful if scen_combi_nodup.sort_values(by=['Nr']) in case of no duplicates
+#     scen_combi_sorted = scen_combi.reset_index(drop=True)  # Reset index starting from 0
+    
+    # Create new Dataframe for scenario prob with corresponding key
+    prob = pd.DataFrame({'Nr': key, 'scen_prob':scenario_prob})
+    
+    # Join both Dataframes based on common key
+    join_scen_prob = scen_combi_sorted.join(prob.set_index('Nr'), on='Nr')
+    
+    # Determine Exceedance Probability
+    sorted_N_D = scen_combi.sort_values(by=['N/D'])  # Sort N/D based on ascending order
+    complete_table = sorted_N_D.join(prob.set_index('Nr'), on='Nr')    
+    
+    N_D_prob_ordered = complete_table.iloc[:,2]
+    
+#     N_D = join_scen_prob.dropna(axis=0).iloc[:,1].reset_index(drop=True)  # N/D
+#     N_D_prob = join_scen_prob.dropna(axis=0).iloc[:,2].reset_index(drop=True)  # Corresponding scenario prob
+#     print('N_D', N_D)
+#     print('N_D_PROB', N_D_prob)
+#     order = N_D.argsort()  # Sort N/D based on ascending order
+    
+#     N_D_prob_ordered = N_D_prob[order]
+    N_D_prob_cs = np.cumsum(N_D_prob_ordered[::-1])  # Take cumulative sum of reversed order
+    
+    exceed_prob = (1/4) * N_D_prob_cs[::-1]  # Reverse result (1-FN(n)) and /4 as 4 evacuation scenarios are considered
+    complete_table = complete_table.join(exceed_prob, lsuffix='_orignal', rsuffix='_new') #  complete resulting table
+#     N_D_ordered = N_D[order].reset_index(drop=True)
+    N_D_ordered = complete_table.iloc[:,1]
+    
+    return N_D_ordered, exceed_prob, complete_table
+
+def sensitivity_analysis(cases, scenario_prob, idx_segment, N_D_ordered_original, exceed_prob_original, complete_table):
+    """Docstring incomplete."""
+    case_1 = cases[0]
+    case_2 = cases[1]
+    case_3 = cases[2]
+    case_4 = cases[3]
+    
+    # Dropping segment(s)
+    case_1_exc = case_1.drop(idx_segment, axis=0).reset_index(drop=True)
+    case_2_exc = case_2.drop(idx_segment, axis=0).reset_index(drop=True)
+    case_3_exc = case_3.drop(idx_segment, axis=0).reset_index(drop=True)
+    case_4_exc = case_4.drop(idx_segment, axis=0).reset_index(drop=True)
+    
+    scenario_prob_exc = scenario_prob.drop(idx_segment, axis=0).reset_index(drop=True)
+    
+    # Determine new exceedance probability - Excluding segment(s)
+    N_D_ordered_new, exceed_prob_new, complete_table_new = exceedance_prob([case_1_exc, case_2_exc, case_3_exc, case_4_exc], scenario_prob_exc)
+    
+    # Plotting
+    plot_curves(N_D_ordered_original, exceed_prob_original)
+    plt.plot(N_D_ordered_new.dropna(axis=0), exceed_prob_new.dropna(axis=0), marker="o", linewidth=2, markersize=5, color= 'darkblue', label='Excluding segment(s)', zorder=10)
+    
+    # Plot dropped points
+    x_point = []
+    y_point = []
+    for i in range(len(idx_segment)):
+        row_i = complete_table[complete_table['Nr']== (idx_segment[i] + 1)]
+        for j in range(len(row_i)):
+            x_point.append(row_i.iloc[j, 1])
+            y_point.append(row_i.iloc[j, 3])
+    plt.scatter(x_point, y_point, label='Dropped Points', s=45, color='darkred', zorder=15)  # Plot the points corresponding to deleted segment
+    plt.legend()
+    plt.show()
+    return N_D_ordered_new, exceed_prob_new, complete_table_new
+
+def plot_curves(N_D, exceedance_prob):
+    """Plot FN or FD risk curve."""
+    plt.figure(figsize=(10,6))
+    set_TUDstyle()
+    
+    # Plot FN/FD curve
+    plt.plot(N_D, exceedance_prob, marker="o", linewidth=2, markersize=5, label='FN/D-curve',zorder=0)
+    
+    # Add acceptable limit line
+    ND = np.linspace(1,100000,100000)  # fatalities
+    alpha = 1
+    C = 0.010
+    limit_line = (C/ND**alpha)
+    plt.plot(limit_line, linewidth=2, color='darkgreen', label='Acceptable Limit Line', zorder=5)
+    
+    # Styling
+    plt.yscale('log')
+    plt.xscale('log')
+    plt.grid(True, which="both", axis='x')
+    plt.grid(True, which="major", axis='y')
+    plt.title('FN/D-curve')
+    plt.ylabel('Probability of Exceedance, $1-F_N(n)$')
+    plt.xlabel('Economic damage [meuro] / Fatalities')
+    plt.xlim(1,100000)
+    plt.ylim(10**(-7),10**(-3))
+    plt.legend()
+
+def compute_benefit(table_old, table_new):
+    """Compute new and old consequences and net benefit."""
+    old_consequences = np.sum(table_old["N/D"].values * table_old["scen_prob_new"].values)
+    new_consequences = np.sum(table_new["N/D"].values * table_new["scen_prob_new"].values)
+    benefit = old_consequences - new_consequences
+    return old_consequences, new_consequences, benefit
+    
+def set_TUDstyle():
+    TUcolor = {"cyan": "#00A6D6", "darkgreen": "#009B77", "purple": "#6F1D77",
+               "darkred": "#A50034", "darkblue": "#0C2340",
+               "orange": "#EC6842", "green": "#6CC24A",
+               "lightcyan": "#00B8C8", "red": "#E03C31", "pink": "#EF60A3",
+               "yellow": "#FFB81C", "blue": "#0076C2"}
+    plt.rcParams.update({'axes.prop_cycle': plt.cycler(color=TUcolor.values()),
+                         'font.size': 16, "lines.linewidth": 4})
+    return TUcolor

book/notebooks/overview.md

@@ -0,0 +1,3 @@
+# Notebooks
+
+This chapter contains notebooks that can be used to put your risk and reliability knowledge to the test in real life situations. They have been rendered as individual pages to show you their contents, and each page will have a download link near the top to allow you to download all files needed for the exercise directly in a zip file.