<jupyter_start><jupyter_text>TKT4196 - Aspects of Structural Safety P00-Introduction to Python (version 210822)__<NAME> (2020), updated <NAME> (2021)__ This lecture aims to provide the student with a basic introduction to Python pogramming. Nevertheless, this is just a ''starter-pack'' and the student must develop confidence in the use of Python up to a certain degree of proficiency by practicing on her own. When the student faces issues, it should first resort to Python-dedicated forums. Most likely, somebody has had and shared the same questions and issues that you encounter. It is also recommended to skim through the handbook: *<NAME> (2016) A Primer on Scientific Programming with Python (Fifth Edition)* We will use Spyder to write and run code. Spyder is an integrated development environment (IDE) specifically developed for scientific programming.  Console The spyder console is used to run commands. The __global variables__ in use are collected in the *variable explorer*. Global variables, as opposed to __local variables__ are accessible from the terminal as well as from scripts. This is an important feature to bear in mind and it will be further elaborated later on.   Modules in PythonPython gets access to a number of bult-in functions. We often need to use additional modules, which can be used with the command *import*. E.g:<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>This makes all functions of numpy available under the name np. If you want to see which functions or methods are available within a module, this is stored in a variable called __all__. E.g., for the module we have just imported, you can type:<jupyter_code>np.__all__<jupyter_output><empty_output><jupyter_text>To acess the documentation of a certain function, use the function help():<jupyter_code>help(np.array)
print(np.array.__doc__)<jupyter_output>array(object, dtype=None, *, copy=True, order='K', subok=False, ndmin=0,
          like=None)

    Create an array.

    Parameters
    ----------
    object : array_like
        An array, any object exposing the array interface, an object whose
        __array__ method returns an array, or any (nested) sequence.
    dtype : data-type, optional
        The desired data-type for the array.  If not given, then the type will
        be determined as the minimum type required to hold the objects in the
        sequence.
    copy : bool, optional
        If true (default), then the object is copied.  Otherwise, a copy will
        only be made if __array__ returns a copy, if obj is a nested sequence,
        or if a copy is needed to satisfy any of the other requirements
        (`dtype`, `order`, etc.).
    order : {'K', 'A', 'C', 'F'}, optional
        Specify the memory layout of the array. If object is not an array, the
        newly created array will be in C order (row major) unless '[...]<jupyter_text>Most of the modules that we will be using can be seen in the file tktreq.txt  Value assignment and data structuresAssigning a value to a variable is done using the symbol =. Variables can be different types of data structures, e.g.: scalars, strings, lists, numerical arrays, tuples, dictionaries... * scalars:<jupyter_code>a = 1     # integer; same as a = int(1)
b = 1.2   # float<jupyter_output><empty_output><jupyter_text>* strings:<jupyter_code>a_string = "Hello, world!" 
print(a_string)<jupyter_output>Hello, world!<jupyter_text>* lists:<jupyter_code>a_list = []       # Empty list
b_list = [1,2,4]  # List with three numerical items
c_list = [1,[1,2],'Hola']
print(c_list)<jupyter_output>[1, [1, 2], 'Hola']<jupyter_text>* numerical arrays:<jupyter_code>import numpy as np
a_array1 = np.array([1,2,3,4])    # a vector (direct input of values)
print(a_array1)

a_array2 = np.linspace(1,4,4)     # a vector (using linspace function)
print(a_array2)

a_array3 = np.arange(1,5,1)       # a vector (using arange function)
print(a_array3)
b_array = np.array([[1,2,3],[2,3,1]])    # a matrix
print(b_array)
b_array.shape                  # dimension of matrix<jupyter_output><empty_output><jupyter_text>* tuples: (immutable sequences)<jupyter_code>a_tuple = (1,2,3)
a_tuple1 = a_tuple
print(a_tuple1)<jupyter_output>(1, 2, 3)<jupyter_text>* dictionaries:<jupyter_code>a_dict = {}
a_dict['emne'] = 'TKT 4196'
a_dict['name'] = "Aspects of structural safety"
print(a_dict)
print(a_dict['emne'])<jupyter_output>{'emne': 'TKT 4196', 'name': 'Aspects of structural safety'}
TKT 4196<jupyter_text>Probability functionsIn this course we will use some built-in statistical distributions. They will be explained in more detail as we move forward with the course. Statistical functions can be imported from scipy.stats module. Here, an example of the most used one is shown, i.e. the Normal distribution:<jupyter_code>import scipy.stats as sps
sps.norm.pdf(0,5,2)      # normal density function with mean 5 and standard deviation 2, evaluated at x=0.
sps.norm.cdf(0,5,2)            # normal cumulative distribution function with mean 5 and standard deviation 2, evaluated at x=0.
sps.norm.ppf(0.1,5,2)    # inverse normal distribution with mean 5 and standard deviation 2, evaluated at p=0.1
np.random.normal(0,1,10) # 10 random samples from standard normal distribution<jupyter_output><empty_output><jupyter_text>Can you guess the output of sps.norm.ppf(sps.norm.cdf(4,5,2),5,2) produce without running it? Other distributions are present in scipy.stats: pdf('name',param) function can be used, although it is not optimized. It is recommended to use the specific functions.  PlottingWe use the module matplotlib to generate plots.<jupyter_code>import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>plt.plot function * __Example 1:__ plot sinus function from 0 to 2*pi radians<jupyter_code>fig,ax = plt.subplots(1,2)       # generate figure with 2 plots, organized in 1 row with 2 columns

# Subplot 1: coarse discretization
x1 = np.linspace(0,2*np.pi,10)    # vector that goes from 0 to 2pi in steps of pi/2 rads. Same as x1 = np.arange(0,2*np.pi+np.pi/2,np.pi/2) 
y1 = np.sin(x1)                  # sinus function taken from numpy
ax[0].plot(x1,y1,'k')            # 'k' = color black
ax[0].set_xlabel('$x$')
ax[0].set_ylabel('$y$');
ax[0].set_xlim(0,2*np.pi)
ax[0].set_ylim(-1,1)

# Subplot 2: fine discretization
x2 = np.linspace(0,2*np.pi,101)  # vector that goes from 0 to 2pi in steps of pi/50 rads
y2 = np.sin(x2)
ax[1].plot(x2,y2,'k')
ax[1].set_xlabel('$x$')
ax[1].set_ylabel('$y$');
ax[1].set_xlim(0,2*np.pi)
ax[1].set_ylim(-1,1)
plt.tight_layout()               # fit size of subplots to window<jupyter_output><empty_output><jupyter_text>* __Example 2:__ Histogram of the standard normal distribution<jupyter_code>fig2,ax2 = plt.subplots(1,1)
pmf = np.random.randn(int(1e6))      # generate 1e6 normal distributed random numbers
kwargs = dict(bins=20, density=True, color='blue',edgecolor='black', stacked=True)  #define a dictionary with keyword arguments
ax2.hist(pmf,**kwargs) #apply the keyword arguments
ax2.set_xlabel('$x$')
ax2.set_ylabel('PDF');
ax2.set_xlim(-5,5)<jupyter_output><empty_output><jupyter_text>Now try to plot the *continuous* probability density function as introduced above and plot it over the bar plot<jupyter_code>from IPython import display
# Input x2,y2 to display the standard normal probability density function 
fig3,ax3 = plt.subplots(1,1)
x2 = 
y2 = 
ax3.plot(x2,y2,'r',LineWidth=1)
display.display(fig3)<jupyter_output><empty_output><jupyter_text>3D-plots (plt.contour and plt.surface) * __Surface:__ Plotting the bivariate standard normal distribution<jupyter_code>from matplotlib import cm
X = np.linspace(-5,5,100)
X1, X2 = np.meshgrid(X,X)
pos = np.dstack((X1, X2))   # pairs of [X1,X2]
MU = [0, 0]                 # Mean vector
C = [[1, 0], [0, 1]]        # Covariance matrix
fxx = sps.multivariate_normal(MU, C)
fig = plt.figure()
ax1 = fig.add_subplot(111, projection='3d')
surf = ax1.plot_surface(X1, X2, fxx.pdf(pos),cmap=cm.coolwarm)
# Enhance the plot
fig.colorbar(surf, shrink=0.5, aspect=10);
ax1.set_xlabel('$X_1$')
ax1.set_ylabel('$X_2$')
ax1.set_zlabel('PDF')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>* __Contour__: Same using contour plots<jupyter_code>fig, ax1 = plt.subplots(1, 1)
ax1.contour(X1, X2, fxx.pdf(pos))   # another option is contourf to add colour
# Enhance the plot
ax1.set_xlim(-3,3)
ax1.set_ylim(-3,3)
ax1.set_xlabel('$X_1$')
ax1.set_ylabel('$X_2$')
ax1.set_aspect('equal')
ax1.set_title('Bivariate standard normal distribution')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Conditionals (if statements) and loops  for loops for loops are used to compute repetitive task. For example,<jupyter_code>number = 5
factorial = 1
for i in range(1,number+1):
    factorial*=i
print("{n:g}! = {f:g}\n".format(n=number,f=factorial))<jupyter_output>5! = 120<jupyter_text>while loops while loops are to some extent similar to for loops and many tasks can be done using either of them. While loops are different in the sense that the repetitive task will be performed until a condition is reached.<jupyter_code>number = 5
factorial = 1
i=1
while i < number+1:
    factorial*=i
    i+=1
print("{n:g}! = {f:g}\n".format(n=number,f=factorial))<jupyter_output>5! = 120<jupyter_text>Lets try to compute the probability of a die drawing a 6 from a frequentistic approach. We use a combination of a __for__ loop with a conditional statement (__if__):<jupyter_code>import random
exact_sol = 1/6
count = 0
NoS = 100
for i in range(0,NoS):
    dice_roll = random.randint(1,6)
    if dice_roll == 6:
        count += 1;
sol = count/NoS
error = abs(sol-exact_sol)/(exact_sol)
print("a = {a:.2f}% the die rolled a 6\n".format(a=sol*100))
print("The error associated with {n:g} simulations is {e:.2f}%\n".format(n=NoS,e=error*100))<jupyter_output>a = 23.00% the die rolled a 6

The error associated with 100 simulations is 38.00%<jupyter_text>Estimate, using a while loop, how many simulations are needed to achieve ca. 0.1% accuracy<jupyter_code>NoS = 100
tol = 0.1/100
error = tol+1
while error>tol:
    count = 0
    for i in range(0,NoS):
        dice_roll = random.randint(1,6)
        if dice_roll == 6:
            count += 1;
    sol = count/NoS
    error = abs(sol-1/6)/(1/6)
    print("The error associated with {n:e} simulations is {e:.2f}%\r".format(n=NoS,e=error*100))
    NoS*=10
print("{n:g} simulations are needed to achieve the desired accuracy\n".format(n=NoS/10))<jupyter_output>The error associated with 1.000000e+02 simulations is 32.00%
The error associated with 1.000000e+03 simulations is 5.20%
The error associated with 1.000000e+04 simulations is 4.34%
The error associated with 1.000000e+05 simulations is 1.37%
The error associated with 1.000000e+06 simulations is 0.29%
The error associated with 1.000000e+07 simulations is 0.08%
1e+07 simulations are needed to achieve the desired accuracy<jupyter_text>Functions We have been using so far several built-in functions. They are used to automatize a process, which can be called when needed.   An important characteristic of functions is that they use local variables. This means that it does not have access to the *Variable explorer* in the main session. To use variables of the main session we need to input them into function.  Creating our own functions Let's exemplify how to create a function by creating one that performs the above task of generating n rolls of a dice<jupyter_code>def dicerolls(n):
    """ This function inputs an integer number n and generate n dice rolls"""
    n = int(n)   # we make sure that the number is an integer
    dice_rolls = []
    for i in range(0,n):
        dice_rolls.append(random.randint(1,6))
    return dice_rolls<jupyter_output><empty_output><jupyter_text>Note that the text at the beginning is what will be output when we call help() on the function<jupyter_code>help(dicerolls)
print(dicerolls(10))<jupyter_output>[6, 3, 5, 5, 5, 4, 3, 3, 1, 2]<jupyter_text>Let's use this function to create another one that uses M different dice<jupyter_code>def Mdice(M,n):
    """ This function rolls M dice n times. It outputs an np.array shape MxN"""
    M = int(M)
    n = int(n)
    dices_rolls = np.zeros((M,n))
    for j in range(0,M):
            dices_rolls[j] = dicerolls(n)
    return dices_rolls
Mdice(2,10)<jupyter_output><empty_output><jupyter_text>Use this function to calculate what is the probability of getting two sixs by using 2 dices<jupyter_code>NoS = 1e4
rolls = Mdice(2,NoS)
boolean = np.all(rolls ==6,axis=0)
res = sum(boolean)/NoS
exact = 1/36
error = abs(res-exact)/exact
print("The error associated with {n:e} simulations is {e:.2f}%\r".format(n=NoS,e=error*100))<jupyter_output>The error associated with 1.000000e+04 simulations is 1.36%<jupyter_text>Lambda functions When functions are simple input output relations, they can sometimes be simply defined in one line using functions handles:<jupyter_code>f = lambda x: (x+1)**2 + 2*x   # a value f is given using the provided formula by inputting x 
print(f(2))
fig3,ax3=plt.subplots()
xarray = np.linspace(x_min-0.5,x_min+0.5)
ax3.plot(xarray,f(xarray))
ax3.set_xlabel('$x$')
ax3.set_ylabel('$f(x)$')<jupyter_output><empty_output><jupyter_text>Some useful built-in functions A type of function that we will often use are solvers. Let's try to find the zero of the function f(x)=(x+1)**2 + 2*x<jupyter_code>import scipy as sp
x0 = 1.  # Initial guess (sometimes very important but always required)
f_zero = sp.optimize.fsolve(f,x0)[0]
print("f({x:.3f}) = 0\r".format(x=f_zero))<jupyter_output>f(-0.268) = 0<jupyter_text>and its minimum<jupyter_code>solmin = sp.optimize.minimize(f,x0,method='nelder-mead',options={'xatol': 1e-8, 'disp': True})
f_min = float(solmin.fun)
x_min = float(solmin.x)
print("minima is f({x:.3f}) = {fm:.3f}\r".format(x=x_min,fm=f_min))
# Input here a line of code that plots the minimum of the function f that we have just computed

display.display(fig3)<jupyter_output><empty_output><jupyter_text>Integrate the previous function from -10 to 10<jupyter_code>f_int,int_err = sp.integrate.quad(f,-10,10)
print("Integral of f from -10 to 10 is {:.2f}".format(f_int))<jupyter_output>Integral of f from -10 to 10 is 686.67<jupyter_text>Integrate the bivariate normal distribution:<jupyter_code>fx = lambda x,y: fxx.pdf([x,y])
fx_int = sp.integrate.nquad(fx,[[-5,5],[-5,5]])[0]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11207434634856786, 'fraction_numerical': 0.035508703793605656, 'mean_word_length': 4.823909531502423, 'pattern_counts': {'":': 0, '<': 108, '<?xml version=': 0, '>': 108, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '576786', 'license': 'MIT', 'n_tokens_mistral': 5231, 'n_tokens_neox': 4909, 'n_words': 1605, 'path': 'P00-Introduction to Python.ipynb', 'repo_name': 'jochenkohler/TKT4196_material', 'size': '14443'}