massive transfer of risk-reliability book

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