--- author: - Boolean algebra, conditional statements, loops. date: "by [Eugeniy E. Mikhailov](http://physics.wm.edu/~evmik/) and [Greg Bentsen](https://physics.wm.edu/~gbentsen/)" title: Functions and scripts --- # Loops examples ## Interest rate related example using loops Suppose a bank gives you a 50% interest rate (let's call it 'x'), and you put one dollar in. How much would you get at the end of the year? - one payment at the end of the year $$M_1=1 \times (1+x)=1 \times (1+0.5)=1.5$$ - interest payment every half a year $$M_2=1 \times (1+x/2) \times (1+x/2)=1 \times (1+0.5/2)^2=1.5625$$ - interest payment every month $$M_{12}=1 \times (1+x/12)^{12}=1.6321$$ ## Interest rate related example (execution) Now let's find how your return on investment ($M_N$) depends on the number of payments per year ::: columns ::: {.column width="45%"} ``` {.python .numberLines} import numpy as np Nmax = 1000 M = np.zeros(Nmax) # we need to preallocate space x=0.5 allN = range(1, Nmax+1) for N in allN: M[N-1] = np.pow(1+x/N, N) # position in array starts from 0! import matplotlib.pyplot as plt plt.plot(allN, M, '.-') plt.grid(True) plt.xlabel("N, number of payments per year") plt.ylabel("$M_N$, return on investment") # Note it understands LaTeX! plt.title("Return on investment vs number of payments") plt.savefig("return_on_investment.png") ``` ::: ::: {.column width="45%"} ![Resulting plot of return on investment](plots_python/return_on_investment.png) ::: ::: . . . Of course we do not need a computer to show that $M_\infty =e^{1/2}=1.6487$ but we need it to calculate something like $M_{1001}-M_{1000}=2.0572 \times 10^{-7}$ > Question: can we do it without use of loops? # Scripts and functions ## Scripts Script is the sequence of the Python expressions written in the file. So far we typed them directly to the python prompt. ``` {.python .numberLines} import numpy as np Nmax = 1000 M = np.zeros(Nmax) # we need to preallocate space x=0.5 allN = range(1, Nmax+1) for N in allN: M[N-1] = np.pow(1+x/N, N) # position in array starts from 0! import matplotlib.pyplot as plt plt.plot(allN, M, '.-') plt.grid(True) plt.xlabel("N, number of payments per year") plt.ylabel("$M_N$, return on investment") # Note it understands LaTeX! plt.title("Return on investment vs number of payments") plt.savefig("return_on_investment.png") ``` Let's save it to the file [plot_return_on_investment.py](plots_python/plot_return_on_investment.py). Now we can edit the file and assign any `x` and `Nmax`. ::: columns ::: {.column width="45%"} To execute the script directly in the shell ``` bash > python plot_return_on_investment.py ``` There is no visible output, but the figure regenerated. ::: ::: {.column width="50%"} Alternatively, we can run our script within `ipython` itself! ``` python In [1]: %run plot_return_on_investment.py ``` `\alert{will not work in python}`{=tex} ::: ::: Now we are free to massage our script to satisfy our needs. ## Scripts variable space Unlike functions `\alert{scripts modify Workspace variables}`{=tex} ``` {.python .numberLines} import numpy as np Nmax = 1000 M = np.zeros(Nmax) # we need to preallocate space x=0.5 allN = range(1, Nmax+1) for N in allN: M[N-1] = np.pow(1+x/N, N) # position in array starts from 0! import matplotlib.pyplot as plt plt.plot(allN, M, '.-') plt.grid(True) plt.xlabel("N, number of payments per year") plt.ylabel("$M_N$, return on investment") # Note it understands LaTeX! plt.title("Return on investment vs number of payments") plt.savefig("return_on_investment.png") ``` ::: columns ::: {.column width="50%"} ``` {.python .numberLines} In [1]: Nmax=22 In [2]: %run plot_return_on_investment.py In [3]: Nmax Out[3]: 1000 ``` See that `Nmax` value is changed ::: ::: {.column width="45%"} Think about the script as it is a keyboard macro. Calling a script is equivalent to typing the scripts statements from the keyboard. ::: ::: # Functions ## Python functions Used for separation of a meaningful chunk of code ``` python def function_name(arg1, arg2, arg3, ..., argN, karg1=val1, karg2=val2, ..., kargN=valN) """Optional but super handy doc string""" # body of the function # notice the indentation! . . return (ret1, ret2, ret3, ... , retN) ``` - All `args`, `kargs`, and `rets` are optional and defined as needed. - There is also a way to make a function with arbitrary number of options or arguments - but this is advance topic - if you wonder about example, `print` function is a good example - Function can be typed at the python prompt directly - Or collections of functions could be saved in file or files - in this case we call such collections modules, packages, or libraries ## Example of a function definition ``` python In [5]: def minmax(n1, n2=0): ...: """Take two numbers and return them in ascending order. ...: The second argument is optional, if not provided 0 will be used. ...: """ ...: if n1 < n2: ...: return (n1, n2) # this is actually one value, i.e. *tuple* of 2 numbers ...: else: ...: return (n2, n1) ...: # note that python will never reach this point ...: In [6]: minmax(23) Out[6]: (0, 23) In [7]: minmax(4, -2) Out[7]: (-2, 4) In [8]: minmax(1, 5) Out[8]: (1, 5) ``` We can also ask documentation about our function ``` python In [15]: minmax? Signature: minmax(n1, n2=0) Docstring: Take two numbers and return them in ascending order. The second argument is optional, if not provided 0 will be used. File: ~/jobs/wm/2026.fall_practical_computing_for_scientists_256/lecture03/plots_python/ Type: function ``` ## Local space of variables in functions Unlike scripts, functions leave Workspace intact. `\alert{There are exceptions, in particular in numpy world!}`{=tex} ``` python # Recall that Signature: minmax(n1, n2=0) >>> n2=33 >>> n1=10 >>> mininmax(10) ``` What we will get back? . . . ``` python (0, 10) ``` ## Reusing variables which exist in workspace considered bad idea ``` python In [6]: def f(x): ...: print(n) # notice that we did not provide n as an argument ...: return x ``` . . . ``` python In [7]: f(3) --------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[7], line 1 ----> 1 f(3) Cell In[6], line 2, in f(x) 1 def f(x): ----> 2 print(n) 3 return x NameError: name 'n' is not defined ``` . . . ``` python In [8]: n=111 In [9]: f(3) 111 Out[9]: 3 ``` In big projects this "feature" is very prone to `\alert{run out of control!}`{=tex} ## Recursion: function calls itself Canonical example: factorial $$N!=N\times(N-1)\times(N-2)\cdots3\times2\times1$$ We can rewrite it as $$N!=N\times(N-1)!$$ Notice that $0!=1$ ``` python In [10]: def fact(n): ...: if n <= 1: ...: return 1 ...: return n * fact(n-1) ...: In [11]: fact(3) Out[11]: 6 In [12]: fact(5) Out[12]: 120 ``` While certain code will look more elegant with recursion, there is a `\alert{hard limit}`{=tex} for number of the recursive calls. Also generally such code is slower than if you can implement it with loops. ## Example: calculating $\sin$ via Taylor expansion You might ask how computer calculates $\sin$? One way is to do it by approximating value of $\sin(x)$ via Taylor expansion $$ \sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} \cdots = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!} $$ ``` {.python .numberLines} import math import numpy as np def mysin(x, N=10): s=0 for n in range(0,N+1): s = s + np.pow(-1.0, n) * np.pow(x, 2*n+1) / math.factorial( 2*n +1) # we could have use our own factorial implementation as well # notice that direct use of factorial is *bad* idea # since as n grows, it has to recaluclate same things (n-1)! again and again # it would be better to rethink it # There is incoming problem for large n! factorial overflow the float representation return s ``` ``` {.python .numberLines} >>> mysin(math.pi) np.float64(1.0348185903053497e-11) >>> mysin(math.pi, N=50) np.float64(3.3280566969799443e-16) >>> mysin(math.pi/2) np.float64(1.0000000000000002) >>> mysin(math.pi/2, N=50) np.float64(1.0000000000000002) ``` . . . Spoiler, this is not how trig functions are calculated for arbitrary $x$. Usually, via trig identities $x$ brought to the smallest usable value and then it is calculated. Recall that expansion requires the least members of the sum the close $x$ to 0 (since our equation was expansion of $\sin$ around 0. ## Comparing our `mysin` implementation to the `numpy.sin` one ``` {.python .numberLines} import numpy as np import matplotlib.pyplot as plt x = np.linspace(0, 2*np.pi, 100) # we create 100 equidistant points y_numpy_sin = np.sin(x) # calculate sin for all points of array in one go # unfortunately mysin cannot work with array so we have to loop over x values y_mysin = x*0 # preallocating array for idx, xp in enumerate(x): y_mysin[idx] = mysin(xp) plt.figure(1) plt.clf() plt.plot(x, y_numpy_sin, label='numpy sin') plt.plot(x, y_mysin, label='my sin') plt.legend() plt.grid(True) plt.xlabel("Angle (rad)") plt.ylabel("$\sin$ value") plt.title("Comparison of numpy and our own $\sin$ calculations") plt.figure(2) plt.clf() plt.plot(x, y_numpy_sin-y_mysin, label='difference') plt.grid(True) plt.legend() plt.xlabel("Angle (rad)") plt.ylabel("np.sin(x) - mysin(x)") plt.title("Comparison of numpy and our own $\sin$ calculations") ``` ## Now we know enough of python to ~~be dangerous~~ do useful work By now we completed equivalent of CSCI 141 class :) The rest is mastering what we know. ### What did we skip? Quite a lot. The major omitted pieces - [data structures](https://docs.python.org/3/tutorial/datastructures.html) - list, set, tuple, dictionary, etc - [classes](https://docs.python.org/3/tutorial/classes.html) - [lamda expressions](https://docs.python.org/3/tutorial/controlflow.html#lambda-expressions) We will talk about them as needed or avoid them completely. After all numerical computations usually do not require fancy things. The sliver of python iceberg can be seen in [python tutorial](https://docs.python.org/3/tutorial/).