




```javascript
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')
```

# SciPy – A collection of scientific algorithms and data structures

Much of this lecture is based on the excellent work of  J.R. Johansson who published his work here: 
J.R. Johansson (jrjohansson at gmail.com) [http://jrjohansson.github.io](http://jrjohansson.github.io).
Some examples have been pulled from https://github.com/scipy/scipy-cookbook/

The SciPy framework builds on top of the low-level NumPy framework for multidimensional arrays, and provides a large number of higher-level scientific algorithms. Some of the topics that SciPy covers are:

* Special functions ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))
* Integration and ODE  ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))
* Optimization ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))
* Interpolation ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))
* Fourier Transforms ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))
* Signal Processing ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))
* Linear Algebra ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))
* Sparse Eigenvalue Problems ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))
* Statistics ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))
* Multi-dimensional image processing ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))
* File IO ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))

Each submodule in SciPy  provides a number of functions and classes that can be used to solve problems in their respective domain.


<h1 id="tocheading">Table of Contents</h1>
<div id="toc"></div>


```python
#some code to help display things in this notebook
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

sns.set(context='notebook')
```


```python
# plot style so figures are more readable. Hopefully
plt.rcParams['figure.figsize'] = (11, 7.5)
plt.rcParams['font.size'] = 15
plt.rcParams['lines.linewidth'] = 3
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
plt.rcParams['axes.labelsize'] = 15
plt.rcParams['axes.titlesize'] = 15
plt.rcParams['legend.fontsize'] = 14
```

# Linear algebra

The linear algebra module contains a lot of matrix related functions, including linear equation solving, eigenvalue solvers, matrix functions (for example matrix-exponentiation), a number of different decompositions (SVD, LU, cholesky), etc. 

Detailed documetation is available at: http://docs.scipy.org/doc/scipy/reference/linalg.html

Here we will look at how to use some of these functions:



## Linear equation systems

Linear equation systems on the matrix form

$A x = b$

where $A$ is a matrix and $x,b$ are vectors can be solved like:


```python
from scipy import linalg
```


```python
A = (np.arange(100) + np.random.normal(0, 10, size=(100))).reshape(10, 10)
b = np.arange(10)
A.shape, b.shape
```


```python
plt.matshow(A)
plt.grid(False) # turn off the grid for this plot
```


```python
x = linalg.solve(A, b)
x
```


```python
# check whether we found the correct x. All entries should be close to zero.
np.dot(A, x) - b
```

We can also do the same with

$A X = B$

where $A, B, X$ are matrices:


```python
A = np.random.rand(5,5)
B = np.random.rand(5,5)
f, [ax1, ax2] = plt.subplots(1, 2)
ax1.matshow(A)
ax2.matshow(B)
# plt.matshow(B)
ax1.grid()
ax2.grid()
```


```python
X = linalg.solve(A, B)
X
```


```python
# check
linalg.norm(np.dot(A, X) - B)
```

## Eigenvalues and eigenvectors

The eigenvalue problem for a matrix $A$:

$\displaystyle A v_n = \lambda_n v_n$

where $v_n$ is the $n$th eigenvector and $\lambda_n$ is the $n$th eigenvalue.

To calculate eigenvalues of a matrix, use the `eigvals` and for calculating both eigenvalues and eigenvectors, use the function `eig`. Both return complex numbers.


```python
A = np.random.rand(4, 4)

eigen_values = linalg.eigvals(A)
eigen_values
```


```python
values, vectors = linalg.eig(A)
values, vectors
```

There are also more specialized solvers, like the `eigh` for Hermitian matrices. 

## Exercise:  Simple Linear Algebra

1. Create a 10x10 matrix $A$ with random entries

2. Create a symmetric matrix from $A$ using its transpose (use A.T to transpose) and then plot the new matrix using `plt.matshow`

3. Use `linalg.eigh` to get the eigenvalues and vectors.


```python
# %load simple_linalg_solution.py
```

# Scipy Spatial

Methods and data structures to help deal with spatial data. Spatial data means that we can define some sort of metric on the data. 


```python
from scipy import spatial
N = 5000
points = np.random.uniform(0, 1, (N, 3)) # create 5000 points in 3D

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(points[:, 0], points[:, 1], points[:, 2])
```

## Exercise: SciPy Spatial

Find all points that have a distance of 0.4 to the center point (0.5, 0.5, 0.5). 

1. Read the first paragraph of https://en.wikipedia.org/wiki/K-d_tree. 
2. Create a KDtree using `scipy.spatial.KDTree(points)`
3. Check the scipy documentation of `scipy.spatial.KDTree` for a method finding all points within a given radius.


```python
# %load solution_spatial.py
```

Now we test the SciPy solution versus the NumPy solution from earlier for finding the nearest neighbour for each point in the point cloud. (See NumPy execrises from Tamas earlier).


```python
#%%timeit
#distances, indices = kdtree.query(points, k=2)
#distances[:, 1]
```


```python
#%%timeit
#distances = np.sum((points.reshape(N, 1, 3) - points)**2, axis=2)
#distances[np.diag_indices_from(distances)] = np.inf
#np.argmin(distances, axis=1)
```

# Special functions

A large number of mathematical special functions are important for many computional physics problems. SciPy provides implementations of a very extensive set of special functions. For details, see the list of functions in the reference documention at http://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special. 

- Airy functions
- Elliptic functions and integrals
- Bessel functions
- Gamma functions
- etc....

## Polynomials

The module contains a large number of functions for evaluating polynomials and calculating their roots.

The following plot shows just a few of them.



```python
from scipy import special

functions = [special.eval_chebys, special.eval_laguerre, special.eval_legendre, special.eval_chebyu]

fig, axs = plt.subplots(2, 2, sharex=True, sharey=True )

x = np.linspace(-2., 2., 200)

for ax, f in zip(np.ravel(axs), functions):
    ax.set_title(f.__name__)
    for n in range(6):
        ax.plot(x, f(n, x))
    ax.set_ylim([-3, 3])
```

## Statistical Functions

The module also contains a host of useful functions for statistical computations. for example the so called 'error function' `scipy.special.erf`

\begin{align}
\operatorname {erf}(x)&={\frac {1}{\sqrt {\pi }}}\int _{-x}^{x}e^{-t^{2}},dt \\ &={\frac {2}{\sqrt {\pi }}}\int _{0}^{x}e^{-t^{2}}\,dt.
\end{align}

Which allows you to get the gaussian probability dsitribution
$$
\phi(z) = \frac 1 2 \left(1 + erf\left(\frac{z}{sqrt(2)}\right) \right)
$$

The module `scipy.stats` contains nicer methods to handle statistical methods. More on that later.


```python
x = np.linspace(-3, 3)
plt.plot(x, special.erf(x), label='$erf$')
plt.plot(x, 0.5 * (1 + special.erf(x/np.sqrt(2))), label='$phi$')
plt.xlabel('$x$')
plt.legend()
```

## Bessel Functions
To demonstrate the typical usage of special functions we will look in more detail at the Bessel functions.
The `scipy.special` module includes a large number of Bessel-functions
Here we will use the functions $j_n$ and $y_n$, which are the Bessel functions 
of the first and second kind and real-valued order.

Bessel functions are a family of solutions to Bessel’s differential equation with real or complex order alpha:
$$
x^2 \frac{d^2 y}{dx^2} +x  \frac{dy}{dx} +(x^2−α^2)y=0
$$


```python
from scipy import special
x = np.linspace(0, 10, 100)

for n in range(4):
    plt.plot(x, special.jn(n, x), label=f'$J_{n}(x)$')
plt.legend();
```

Among other uses, these functions arise in wave propagation problems such as the vibrational modes of a thin drum head. Here is an example of a circular drum head anchored at the edge.

https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane


```python
# fix the radius of the drum
DRUM_RADIUS = 1

def drumhead_height(n, k, distance, angle, t):
    '''
    Solution to the Drumhead problem with a fixed membrane on the edge of the drum
    See https://en.wikipedia.org/wiki/Vibrations_of_a_circular_membrane
    '''
    kth_zero = special.jn_zeros(n, k)[-1] # fix the drum on the outter edge to be zero.
    return np.cos(t) * np.cos(n*angle) * special.jn(n, distance*kth_zero/DRUM_RADIUS)


# create two lists containing all combinations of radius and theta we want to plot
r = np.linspace(0, DRUM_RADIUS, 50)
t = np.linspace(0, 2*np.pi, 50)
radius, theta = np.meshgrid(r, t) # this creates all combinations we are interested in

x = radius * np.cos(theta)
y = radius * np.sin(theta)
z = drumhead_height(1, 1, radius, theta, 0.5)
```


```python
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, cmap='YlGnBu')
```


```python
# some code to create a nice animation within the notebook.
# this takes a couple of minutes on my machine so I commented it out

# from matplotlib.animation import FuncAnimation
# from functools import partial
# from tqdm import tqdm_notebook


# my_drum = partial(drumhead_height, n=2, k=1, distance=radius, angle=theta)

# frames = 200
# fps = 25

# fig = plt.figure()
# ax = Axes3D(fig)
# ax.plot_surface(x, y, my_drum(t=0), cmap='YlGnBu')

# bar = tqdm_notebook(total=frames)

# def update(i):
#    bar.update()
   
#    t = i / fps
#    ax.cla()
#    z = my_drum(t=t)
#    ax.plot_surface(x, y, z, cmap='YlGnBu', vmin=-0.5, vmax=0.5)
#    ax.set_zlim(-0.5, 0.5)
   
# ani = FuncAnimation(fig, update, frames=frames, interval=1000 / fps, blit=False, repeat=False)
# ani.save('animation_membrane.mp4', writer="ffmpeg")
```


```python
from IPython.display import Video
Video('./animation_membrane.mp4')
```

# Ordinary Differential Equations.

To demonstrate the usefulness of SciPy lets try to predict what will happen to humanity in case of a zombie outbreak.  This example comes courtesy of Christopher Campo.

We will show that, in case of a Zombie outbreak, humanity is inevitably doomed. 
As shown by [Phillip Munz et al](http://mysite.science.uottawa.ca/rsmith43/Zombies.pdf) we can model a simple outbreak scenario like so.

\begin{align}
\frac{S}{dt} &= B - \beta S Z - \delta S \\
\frac{Z}{dt} &= \beta S Z + \gamma R - \alpha S Z \\
\frac{R}{dt} &= \delta S + \alpha S Z - \gamma R
\end{align}

Where $S$ is the number of susceptible humans, $Z$ is the number of zombies and $R$ is the number of removed persons which have died either by natural causes or zombie attack. The human birthrate is assumed to be constant and modeled by $B$

Susceptibles can become zombies through an encounter with a zombie ($\beta$). Natural causes of human deaths are parameterized by $\delta$. Humans in the removed class can be resurrected and become zombies ($\gamma$).  
Zombies can enter the removed class by cutting their heads of or removing their brains ($\alpha$)

The question I'm trying to answer is 

> __Can Humanity survive the zombie Apocalypse?__

Since I don't know how to solve differential equations, I'm going to use SciPy instead.


Now we use SciPy to solve this system of differential equations. The function `solve_ivp` does the entire job for us. 

If we go to [the SciPy documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp) about this function we will read the following.


>This function numerically integrates a system of ordinary differential equations given an initial value:
>
> $y^{\prime} = f(t, y)$
>
>$y(t_0) = y_0$
>
>Here $t$ is a one-dimensional independent variable (time), $y(t)$ is an __n-dimensional vector-valued function__ (state), and an __n-dimensional vector-valued function f(t, y) determines the differential equations__. The goal is to find y(t) approximately satisfying the differential equations, given an initial value y(t0)=y0.


So we start by defining some starting values for the problem and the model function $f(t, y)$.


```python
from scipy.integrate import solve_ivp

# Type `solve_ivp?` to get inline documentation about the method in your notebook.
```


```python
# fixed model parameters
B = 0.1 # birth rate
delta = 0.0001 # rate of natural casues of death
beta = 0.00016 # transmission coefficent
gamma = 0.0001 # resurrection rate
alpha = 0.0001 # rate of destroyed zombies

# initial conditions
S0 = 5000.              # initial population
Z0 = 0                 # initial zombie population
R0 = 0                 # initial dead population
y0 = [S0, Z0, R0]     # initial condition vector


def f_model(t, y):
    Si = y[0]
    Zi = y[1]
    Ri = y[2]
    # the model equations (see Munz et al. 2009)
    f0 = B - beta*Si*Zi - delta*Si
    f1 = beta*Si*Zi + gamma*Ri - alpha*Si*Zi
    f2 = delta*Si + alpha*Si*Zi - gamma*Ri
    return [f0, f1, f2]

```

Now we solve the equation and output the result


```python
t_start, t_end = 0, 150  # 200 days of zobie acopalypse
solution = solve_ivp(f_model, (t_start, t_end), y0, t_eval=np.linspace(t_start, t_end, 200))
```


```python
s = solution.y[0, :]
z = solution.y[1, :]
plt.plot(solution.t, s, label='survivors')
plt.plot(solution.t, z, label='zombies')
plt.xlabel('t / days')
plt.ylabel('Population')
plt.legend()
```

## Exercise: Save the world.

The ODE defined above is a very pessimistic model. There is one unstable equilibrium state if $\beta = 0$.
Your task is to save the world by adding a the notion of a _cure_ to the model above. The cure 
can move individuals from the $Z$ to the $S$ group.

1. Introduce a new constant $\rho$ which models the effectivness of the cure. 
2. Add a term to the ODE to model the cure. 
3. Solve the ODE using `solve_ivp`




```python
# %load solution_zombies.py
```

## Example : Complicated  ODE

We use SciPy `solve_ivp` method to numerically solve the following system of coupled differential euqations 

\begin{aligned}{{\dot {\theta }}_{1}}&={\frac {6}{ml^{2}}}{\frac {2p_{\theta _{1}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{2}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}}\\{{\dot {\theta }}_{2}}&={\frac {6}{ml^{2}}}{\frac {8p_{\theta _{2}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{1}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}}.
\end{aligned}

\begin{aligned}{{\dot {p}}_{\theta _{1}}}&={\frac {\partial L}{\partial \theta _{1}}}=-{\tfrac {1}{2}}ml^{2}\left({{\dot {\theta }}_{1}}{{\dot {\theta }}_{2}}\sin(\theta _{1}-\theta _{2})+3{\frac {g}{l}}\sin \theta _{1}\right)\\{{\dot {p}}_{\theta _{2}}}&={\frac {\partial L}{\partial \theta _{2}}}=-{\tfrac {1}{2}}ml^{2}\left(-{{\dot {\theta }}_{1}}{{\dot {\theta }}_{2}}\sin(\theta _{1}-\theta _{2})+{\frac {g}{l}}\sin \theta _{2}\right).
\end{aligned} 

Now we plot the results for the $\theta_1$ and $\theta_2$ using 

\begin{align}
x_1 &= L \; sin(\theta_1(t))\\
y_1 &= - L\;  cos(\theta_1(t)) \\
x_2 &= x_1 + L \; sin(\theta_2(t)) \\
y_2 &= y_1 - L \; cos(\theta_2(t))
\end{align}


```python
from scipy.constants import g
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

L = 0.5
m = 0.1

t_start = 0
t_end = 20
fps = 50
frames = (t_end - t_start) * fps


def pendulum_model(t, y):
    """
    The right-hand side of the pendulum ODE
    """
    phi_1, phi_2, p_1, p_4 = y[0], y[1], y[2], y[3]

    c = (m * L**2)
    cos_delta_phi = np.cos(phi_1 - phi_2)
    sin_delta_phi = np.sin(phi_1 - phi_2)
    denominator = 16 - 9 * cos_delta_phi**2

    dphi_1 = 6 / c * (2 * p_1 - 3 * cos_delta_phi * p_4) / denominator
    dphi_2 = 6 / c * (8 * p_4 - 3 * cos_delta_phi * p_1) / denominator

    dp_1 = -c / 2 * ( dphi_1 * dphi_2 * sin_delta_phi + 3 * g / L * np.sin(phi_1))
    dp_4 = -c / 2 * (-dphi_1 * dphi_2 * sin_delta_phi + g / L * np.sin(phi_2))

    return [dphi_1, dphi_2, dp_1, dp_4]


# initial values
y0 = [np.pi / 3, -np.pi / 4, 0, 0.065]


solution = solve_ivp(
    pendulum_model,
    (t_start, t_end),
    y0,
    t_eval=np.linspace(t_start, t_end, frames),
)

x1s = + L * np.sin(solution.y[0])
y1s = - L * np.cos(solution.y[0])

x2s = x1s + L * np.sin(solution.y[1])
y2s = y1s - L * np.cos(solution.y[1])

plt.plot(x1s, y1s, '-', alpha=0.5, color='xkcd:sky blue')
plt.plot(x2s, y2s, '-', alpha=0.5, color='xkcd:dark magenta')
plt.axis('off')
plt.show()

```


```python
from matplotlib import animation
from IPython.display import HTML

fig, ax = plt.subplots(figsize=(6, 5))

ax.set_ylim([-1.1, 0.2])
ax.set_xlim([1, -1])

pendulum1, = ax.plot([], [])
pendulum2, = ax.plot([], [])


bar = tqdm_notebook(total=frames)

def init():
   pendulum1.set_data([], [])
   pendulum2.set_data([], [])
   return (pendulum1, pendulum2)

def update(n):
   bar.update(1)
   # update the line data
   pendulum1.set_data([0, x1s[n]], [0, y1s[n]])
   pendulum2.set_data([x1s[n], x2s[n]], [y1s[n], y2s[n]])
   return (pendulum1, pendulum2)


anim = animation.FuncAnimation(fig, update, init_func=init, frames=frames, blit=True, interval=1000 / fps)
anim.save('animation_pendulum.mp4', writer="ffmpeg")
```


```python
Video("./animation_pendulum.mp4")
```


```python
Video?
```

# Numerical Integration

## Numerical integration: quadrature

Numerical evaluation of a function of the type

$\displaystyle \int_a^b f(x) dx$

is called *numerical quadrature*, or simply *quadature*. SciPy provides a series of functions for different kind of quadrature, for example the `quad`, `dblquad` and `tplquad` for single, double and triple integrals, respectively.




```python
from scipy.integrate import quad, dblquad, tplquad
```

The `quad` function takes a large number of optional arguments, which can be used to fine-tune the behaviour of the function.

The basic usage is as follows:


```python
# define a simple function for the integrand
def f(x):
    return x**2
```


```python
x_lower = 0 # the lower limit of x
x_upper = 1 # the upper limit of x

val, abserr = quad(f, x_lower, x_upper)

print(f'integral value={val} , absolute error={abserr}' )
x = np.linspace(-1.2, 1.2, 200)
plt.plot(x, f(x))
plt.fill_between(x, f(x), alpha=0.3, where=(x > 0) & (x < 1))

```

If we need to pass extra arguments to integrand function we can use the `args` keyword argument:


```python
def integrand(x, n):
    return special.jn(n, x)

x_lower = 0  # the lower limit of x
x_upper = 10 # the upper limit of x

val, abserr = quad(integrand, x_lower, x_upper, args=(3,))
print(val, abserr) 
```

For simple functions we can use a lambda function (name-less function) instead of explicitly defining a function for the integrand:


```python
val, abserr = quad(lambda x: np.exp(-x**2), -np.Inf, np.Inf)
print(f'numerical= {val}, {abserr}')
print(f'analytical solution = {np.sqrt(np.pi)}')
```

As show in the example above, we can also use 'Inf' or '-Inf' as integral limits.

Higher-dimensional integration works in the same way:


```python
def integrand(x, y):
    return np.exp(-x**2 - y**2)


x, y = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(-2, 2, 100))
z = integrand(x, y)

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, cmap='YlGnBu')

# the integral boundarys can be functions when integrating in more than one dimensions
val, abserr = dblquad(integrand, -2, 2, lambda x: -2, lambda x: 2) 

print(f'numerical solution= {val}, {abserr}')
```

Note how we had to pass lambda functions for the limits for the y integration, since these in general can be functions of x.

## Exercise II:  Integrals and Special Functions

Find the circumference of an ellipse. The formula for calculating circumference $C$ is as follows.

$$
\frac{C}{4a} =  E(\epsilon)=\int _{0}^{\pi /2}{\sqrt {1-\epsilon^{2}\sin ^{2}\theta }}\ d\theta
$$

where $\epsilon$ is the eccentricity of the ellipse given by the length and width of the ellipse.

$$
\epsilon = \sqrt{1 - \frac{b^2}{a^2}}
$$

compare the solution found using `quad` with the appropriate function from the `scipy.special` package.


```python
# %load solution_2.py
```

# Signal Processing

SciPy offers a large amount of utilities to handle time series data. However the methods in the `scipy.signal` are useful for much more than just time series data.

## Convolution

The convolution of two continuos functions is defined by a simple integral transform.

\begin{aligned}
(f*g)(t)&\,{\stackrel {\mathrm {def} }{=}}\ \int _{-\infty }^{\infty }f(\tau )g(t-\tau )\,d\tau \\&=\int _{-\infty }^{\infty }f(t-\tau )g(\tau )\,d\tau .
\end{aligned}

Convolution happens whenever data gets smeared out by a detector or filter of any kind.



```python
from scipy import signal

# create two rectangular signals
f = np.repeat([0., 1., 0.], 120)
g = np.repeat([0., 1., 0.], 20)

result = signal.convolve(f, g, mode='same')  # same output shape as input, see docs
```


```python
fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True)

ax_orig.plot(f, color='gray')
ax_orig.set_title('Original pulse')
ax_orig.margins(0, 0.1)

ax_win.plot(g, color='gray')
ax_win.arrow(70, 0.5, 50.0, 0.0, width=0.08, head_length=8)
ax_win.set_title('Filter impulse response')

ax_filt.plot(result)
ax_filt.set_title('Filtered signal')
None
```

The naive algorithm for producing discrete convolutions is not very efficient. It takes $\mathcal{O}(n^2)$ steps to get a result. SciPy uses a heuristic to automagically find the fastest convolution method. It is based on the *convolution theorem*

> The convolution theorem states that
>
 $$
 \mathcal{F}\{f*g\}=k\cdot {\mathcal {F}}\{f\}\cdot {\mathcal {F}}\{g\}
 $$
>
> where 
> $\mathcal {F}\{f\}$ denotes the Fourier transform of $f$, and $\mathcal {F}\{g\}$ the Fourier transform of $g$.
>
> https://en.wikipedia.org/wiki/Convolution

## Signal Windows

SciPy conveniently provides a handful of common signal windows ready to be used for convolution or other filtering operations. All of them can be found in `scipy.signal.window` but their are also available in the `scipy.signal` name space.


```python
for g in [signal.nuttall(150), signal.hann(150), signal.cosine(150), signal.triang(150), signal.blackman(150)]:
    plt.plot(g)
plt.axis('off')
None
```

## Correlation

The correlation of two functions $f$ and $g$ looks similar to a convolution 

> For continuous functions $f$ and $g$, the cross-correlation is defined as:
> 
$$
(f\star g)(\tau )\ {\stackrel {\mathrm {def} }{=}}\int _{-\infty }^{\infty }f^{*}(t)\ g(t+\tau )\,dt,
$$
> 
> where $f^{*}$ denotes the complex conjugate of  $f$, and $\tau$ is the displacement, also known as lag, although a positive value of $\tau$ actually means that $g(t+\tau)$ leads $f(t)$.

> *https://en.wikipedia.org/wiki/Cross-correlation*

The result of the correlation is a function of the delay $\tau$. This is useful for finding the offset in shifted signals.


```python
offset = 0.25
t = np.linspace(-1, 1, 600)
f = signal.gausspulse(t, fc=3)
g = signal.gausspulse( t + offset, fc=3)

result = signal.correlate(f, g, mode='same')
offset = t[np.argmax(result)]
```


```python
grid = plt.GridSpec(2, 2, wspace=0.1, hspace=0.4, width_ratios=[1, 1])
ax1 = plt.subplot(grid[0, 0])
ax2 = plt.subplot(grid[1, 0])
ax3 = plt.subplot(grid[:, 1], )

ax1.plot(t, f, color='gray')
ax1.set_title('$f$')
ax1.set_yticks([])
ax1.set_xlabel('t')

ax2.plot(t, g, color='gray')
ax2.set_title('$g$')
ax2.set_yticks([])
ax2.set_xlabel('t')

ax3.plot(t, result)
ax3.axvline(offset, color='xkcd:tangerine')
ax3.set_title('$(f\star g)(\\tau )$')
ax3.set_yticks([])
ax3.set_xlabel('τ')
None
```

## Exercise III:  Signal Filtering

If we ever find an alien civilization we will have to communicate with them using binary encoded signals recorded by large radio telescopes.


Assuming the aliens use binary signals with 0s and 1s. 
Given the clock speed and the sample length with which the signal was created. Can you decode the following signal?

The plot below shows the shape of an example signal without any noise containing the message [0, 1 0]


```python
sample_length = 64
plt.plot(np.repeat([0, 1, 0], sample_length))
plt.text(25, 0.5, '0')
plt.text(100, 0.5, '1')
plt.text(175, 0.5, '0')
```


```python
alien_signal = np.loadtxt('alien_signal.txt')
sample_length = 128

# let's sample the signal in the middle of each period
clock = np.arange(sample_length // 2, len(alien_signal), sample_length)
plt.plot(alien_signal)
for c in clock:
    plt.axvline(c, color='gray', lw=1.5)
```


```python
# %load solution_3.py
```

# Fourier transform

Fourier transforms are one of the universal tools in computational physics, which appear over and over again in different contexts. SciPy provides functions for accessing the classic [FFTPACK](http://www.netlib.org/fftpack/) library, which is an efficient and well tested FFT library written in FORTRAN. The SciPy API has a few additional convenience functions, but overall the API is closely related to the original FORTRAN library.

To use the `fftpack` module in a python program, include it using:


```python
from scipy import fftpack 
```

To demonstrate how to do a fast Fourier transform with SciPy, we will create some artificial data first.


```python
# Number of sample points
N = 1000
L = 10
# sample spacing
T = L / N
x = np.linspace(0.0, L, N)

freqs = [8, 3, 32]
amps = [1, 0.5, 0.25]

y = np.zeros(N)
for a, f in zip(amps, freqs):
    y += a * np.sin(2 * np.pi * f * x)

y += np.random.normal(0, 0.75, len(y))
plt.plot(x, y, lw=2)
```


```python
yf = fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))

# plt.yscale('log')
# plt.xscale('log')
```

As expected, we now see a peaks in the spectrum that match the true frequencies in the signal.
Another, maybe more convenient,  way to create frequencies from a transformed signal uses the `fftfreq` function.



# Statistics

The `scipy.stats` module contains a large number of statistical distributions, statistical functions and tests. For a complete documentation of its features, see http://docs.scipy.org/doc/scipy/reference/stats.html.

Below we instantiate an object modeling a random variable following the normal distribution.


```python
from scipy import stats

mean = 0
sigma = 1
X = stats.norm(loc=mean, scale=sigma)
```

All distributions in this submodule extend either `scipy.stats.rv_continuous` or `scipy.stats.rv_discrete`.
Below is a list of functions supported by most continuous distribution objects

```
rvs(*args, **kwds) 	Random variates of given type.
pdf(x, *args, **kwds) 	Probability density function at x of the given RV.
logpdf(x, *args, **kwds) 	Log of the probability density function at x of the given RV.
cdf(x, *args, **kwds) 	Cumulative distribution function of the given RV.
logcdf(x, *args, **kwds) 	Log of the cumulative distribution function at x of the given RV.
sf(x, *args, **kwds) 	Survival function (1 - cdf) at x of the given RV.
logsf(x, *args, **kwds) 	Log of the survival function of the given RV.
ppf(q, *args, **kwds) 	Percent point function (inverse of cdf) at q of the given RV.
isf(q, *args, **kwds) 	Inverse survival function (inverse of sf) at q of the given RV.
moment(n, *args, **kwds) 	n-th order non-central moment of distribution.
stats(*args, **kwds) 	Some statistics of the given RV.
entropy(*args, **kwds) 	Differential entropy of the RV.
expect([func, args, loc, scale, lb, ub, …]) 	Calculate expected value of a function with respect to the distribution.
median(*args, **kwds) 	Median of the distribution.
mean(*args, **kwds) 	Mean of the distribution.
std(*args, **kwds) 	Standard deviation of the distribution.
var(*args, **kwds) 	Variance of the distribution.
interval(alpha, *args, **kwds) 	Confidence interval with equal areas around the median.
__call__(*args, **kwds) 	Freeze the distribution for the given arguments.
fit(data, *args, **kwds) 	Return MLEs for shape (if applicable), location, and scale parameters from data.
fit_loc_scale(data, *args) 	Estimate loc and scale parameters from data using 1st and 2nd moments.
nnlf(theta, x) 	Return negative loglikelihood function.

```


```python
N = 1200
measurments = X.rvs(size=N)
x = np.linspace(-3.5, 3.5, 200)

f, [ax1, ax2] = plt.subplots(1, 2, figsize=(16, 8))
ax1.hist(measurments, bins=40, density=True)
ax1.plot(x, X.pdf(x), color='red')


ax2.hist(measurments, bins=40, cumulative=True, density=True, )

ax2.plot(x, X.cdf(x), color='red')
None
```

Discrete functions work in a very similar manner.


```
rvs(*args, **kwargs) 	Random variates of given type.
pmf(k, *args, **kwds) 	Probability mass function at k of the given RV.
logpmf(k, *args, **kwds) 	Log of the probability mass function at k of the given RV.
cdf(k, *args, **kwds) 	Cumulative distribution function of the given RV.
logcdf(k, *args, **kwds) 	Log of the cumulative distribution function at k of the given RV.
sf(k, *args, **kwds) 	Survival function (1 - cdf) at k of the given RV.
logsf(k, *args, **kwds) 	Log of the survival function of the given RV.
ppf(q, *args, **kwds) 	Percent point function (inverse of cdf) at q of the given RV.
isf(q, *args, **kwds) 	Inverse survival function (inverse of sf) at q of the given RV.
moment(n, *args, **kwds) 	n-th order non-central moment of distribution.
stats(*args, **kwds) 	Some statistics of the given RV.
entropy(*args, **kwds) 	Differential entropy of the RV.
expect([func, args, loc, lb, ub, …]) 	Calculate expected value of a function with respect to the distribution for discrete distribution.
median(*args, **kwds) 	Median of the distribution.
mean(*args, **kwds) 	Mean of the distribution.
std(*args, **kwds) 	Standard deviation of the distribution.
var(*args, **kwds) 	Variance of the distribution.
interval(alpha, *args, **kwds) 	Confidence interval with equal areas around the median.
__call__(*args, **kwds) 	Freeze the distribution for the given arguments.
```


```python
# create a (discreet) random variable with poissionian distribution
Y = stats.poisson(3.5) 
```


```python
n = np.arange(0, 12)

fig, axes = plt.subplots(3,1, sharex=True)

# plot the probability mass function (PMF)
axes[0].plot(n, Y.pmf(n), 'o')

# plot the commulative distribution function (CDF)
axes[1].step(n, Y.cdf(n))

# plot histogram of 1000 random realizations of the stochastic variable X
axes[2].hist(Y.rvs(size=200), bins=np.arange(13) - 0.5, density=True, lw=2, edgecolor='w');

axes[2].set_xticks(n)
None
```

The objects support easy access to some of its statistical properties:


```python
print(X.mean(), X.std(), X.var()) # normal distribution
print(Y.mean(), Y.std(), Y.var()) # poission distribution
```

## Statistical tests

The module contains some usefull methods to perform statistical tests. 

Below we tests whether two sets of (independent) random data comes from the same distribution:


```python
# create two sets of data 
N = 500
a = X.rvs(size=N)
b = X.rvs(size=N)

plt.hist(a, bins=50, histtype='step', lw=4, density=True)
plt.hist(b, bins=50, histtype='step', lw=4, density=True)
t_statistic, p_value = stats.ttest_ind(a, b)

print(f't-statistic = {t_statistic}, p-value = {p_value}')
```

Since the p value is very large we cannot reject the hypothesis that the two sets of random data have *different* means.

# Optimization

Optimization (finding minima or maxima of a function) is a large field in mathematics, and optimization of complicated functions or in many variables can be rather complicated. Here we will only look at a few very simple cases. For a more detailed introduction to optimization with SciPy see: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

To use the optimization module in scipy first include the `optimize` module:


```python
from scipy import optimize
```

## Finding Minima

Let's first look at how to find the minima of a simple function of a single variable:


```python
def f(x):
    return 4.5*x**3 + (x-4)**2 + 0.75*x**4 - 20

x = np.linspace(-6, 3, 100)
plt.plot(x, f(x));
```

We can use the `minimize_scalar` function to find the minima of a scalar function. It provides wrapper around many optimization algorithms you can choose from using the `method` keyword. 


```python
result = optimize.minimize_scalar(f)

plt.plot(x, f(x));
plt.axvline(result.x, color='xkcd:orange')


# add some boundaries
result = optimize.minimize_scalar(f, bracket=[-6, -2])
plt.axvline(result.x, color='xkcd:dark red')
```

## Exercise: Likelihood minimization (might be a bit advanced)

Imagine you have a sensor with 400 pixels. The pixel measure the number of photons in in a fixed length of time. Fit the mean of the photon cunts using the negative log likelihood of the poisson distribution

$$
X(k, \lambda) \; \sim \; \frac{\lambda^k}{k!} e^{-\lambda}
$$


```python
true_lambda=2
data = np.random.poisson(true_lambda, (20,20))
plt.imshow(data)
plt.grid()
```


```python
# %load solution_7.py
```

## Minimzation in more dimensions.

Below you'll find the infamous *Rosenbrock* function in two dimensions. It's very popular for benchmarking global optimization problems. Its solutions are well known and easily parametrized. 


```python
from mpl_toolkits.mplot3d import Axes3D

def rosenbrock_2d(x, y, a=1, b=1):
    return (a - x)**2 + b*(y - x**2)**2


fig = plt.figure()
ax = Axes3D(fig)

x, y = np.meshgrid(np.linspace(-2 ,2, 100), np.linspace(-3. ,5, 100))
ax.plot_surface(x, y, rosenbrock_2d(x, y), cmap='magma_r')
```


```python
result = optimize.minimize(lambda x: rosenbrock_2d(x[0], x[1]), x0=[5, 2])
fig = plt.figure()
ax = Axes3D(fig)
ax.view_init(20, 30)


x, y = np.meshgrid(np.linspace(-2 ,2, 100), np.linspace(-3. ,5, 100))
ax.plot_surface(x, y, rosenbrock_2d(x, y), cmap='magma_r', alpha=0.6)

x, y, z = result.x[0], result.x[1], 20
ax.quiver(x, y, z, 0, 0, -1, length=20, arrow_length_ratio=0.04)


result
```

## Least Squares Minimization

Lets presume a linear weighted combination of variables:

$$
f(x)= \hat{y} =  \hat{\beta}_0 + \sum_{j=1}^p x_j \hat{\beta}_j
$$
where $f:\mathbb{R}^{p} \to \mathbb{R}$.


When we include a 1 as the first entry into our sample $x$ e.g. $x = (1, x_1, x_2, \ldots, x_p)$ we can rewrite
$f$ in matrix form

$$
f(x)= \hat{y} =  x^T \mathbf{\beta}
$$

where $\beta = (\beta_0, \beta_1, \beta_2, \ldots, \beta_p)$.



How do you find those weights? Like before we choose a loss function and try to opimize it.
In this case we choose a loss function called the residual sum of squares (RSS).
We calculate it over all samples $x_i$ in a matrix $\mathbf{X}$.

$$L(\beta) = RSS(\mathbf{\beta}) = \sum_{i=1}^N (y_i - x_i^T \beta)^2 $$

Here $x_i$ is a row in $\mathbf{X}$, hence the transpose.

We can now rewrite the loss function in matrix form:


$$
RSS(\beta) = (\mathbf{y} - \mathbf{X} \beta)^T (\mathbf{y} - \mathbf{X} \beta )
$$

Now we optimize the loss function just like we would any other function, by differentiating with respect to $\beta$ and setting the result equals to zero.

$$
 \mathbf{X}^T (\mathbf{y} - \mathbf{X} \beta ) \stackrel{!}{=} 0
$$

Solving for $\beta$ leads to

$$
\hat{\beta} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}
$$


We just performed  __Linear Least Squares__ regression. Below we fit some noisy a polynomial data.


```python
def f(x, beta=[4.5, 1, 0.75]):
    return beta[0]*x**3 + beta[1]*(x-4)**2 + beta[2]*x**4

x = np.linspace(-6, 3, 100)
data = f(x) + np.random.normal(0, 35, size=x.shape)

# build the answer matrix
A = np.vstack([x**3, (x - 4)**2, x**4]).T

# find the linear least squares solution for the parameters
result = optimize.lsq_linear(A, data)


plt.plot(x, data, '.')
plt.plot(x, f(x), color='gray', label='truth')
plt.plot(x, f(x, beta=result.x), color='crimson', label='fit')
plt.legend()
```

## Curve Fitting

For more general cases where the parameters are not linear there is the more convenient function `curve_fit`




```python
from scipy import stats
def f(x, scale=5, loc=4, ):
    return (1/(scale))**2*x**3 + (1/(2*loc))*x**2

x = np.linspace(-3,  1.5, 100)
data = f(x) + np.random.normal(0, 0.15, size=x.shape)

popt, pcov = optimize.curve_fit(f, x, data, p0=[5, 4, ])

plt.plot(x, data, '.')
plt.plot(x, f(x), color='gray', label='truth')
plt.plot(x, f(x, *popt), color='crimson', label='fit')

pcov
```

## Exercise VI:  Curve Fitting

Use `numpy.loadtxt` to load the file called `munich_temperatures_average.txt` and fit a model to it using `optimize.curve_fit`.




```python
# %load solution_6.py
```

# Interpolation

Interpolation is simple and convenient in scipy: The `interp1d` function, when given arrays describing X and Y data, returns and object that behaves like a function that can be called for an arbitrary value of x (in the range covered by X), and it returns the corresponding interpolated y value:


```python
from scipy import interpolate

def f(x):
    return np.sin(x)
```


```python
n = [0, 1, 2, 3, 3.5, 4, 6, 6.5, 7, 9]
x = np.linspace(0, 9, 100)

y_measured = f(n) 

linear_interpolation = interpolate.interp1d(n, y_measured)

cubic_interpolation = interpolate.interp1d(n, y_measured, kind='cubic')
# this returns a new function we can call for values of x
cubic_interpolation
```


```python
plt.plot(n, y_measured, 's', label='noisy data', color='black')
plt.plot(x, f(x), c='gray', lw=2, label='true function')
plt.plot(x, linear_interpolation(x), label='linear interpolation')
plt.plot(x, cubic_interpolation(x), label='cubic interpolation')
plt.legend();
```

# SciPy IO Module

## IO and the sound of hydrogen

SciPy provides some useful methods for reading file from MatLab, NetCDF or WAV data. Below is an example which uses the `scipy.io.wavfile` module to write a file containing the sound of hydrogen atoms.




```python
from scipy.io import wavfile
from scipy import signal

rate = 44100 #44.1 khz
duration = 10 # in sec

def rydberg(n, m):
    return (1/(n**2) -1/(m**2)) 

def lyman(m):
    return rydberg(1,m)

def balmer(m):
    return rydberg(2,m)

def sound(frequency, time):
    return 2**12*np.sin(2*np.pi*440*frequency*time)

t = np.linspace(0, duration, rate*duration)
hydrogen_sound = np.sum([sound(lyman(i), t) + sound(balmer(i), t) for i in range(2, 12)], axis=0)/6

plt.plot(t[:10000], hydrogen_sound[:10000])

wavfile.write('hydrogen.wav', rate, hydrogen_sound.astype(np.int16))
```


```python
%%HTML
<audio controls="controls" style="width:600px" >
  <source src="hydrogen.wav" type="audio/wav" />
</audio>
```

## Exercise IV:  Signal Processing and IO.

Load the file `synth_sound.wav` and plot a the frequency spectrum vs time.


```python
# %load solution_4.py
```

# Root Finding

To find the solution for a function of the form $f(x) = 0$ we can use the `fsolve` or the `root` function. It requires one or more initial guesses: 


```python
def f(x):
    return 4.5*x**3 + (x-4)**2 + 0.75*x**4 - 20

x = np.linspace(-6, 3, 100)

plt.plot(x, f(x))
plt.axhline(color='black')

result = optimize.root(f, [-6, -3, 0, 2])
for r in result.x:
    plt.axvline(r, color='xkcd:tangerine')

```

# Sparse matrices

Sparse matrices are often useful in numerical simulations dealing with large systems, if the problem can be described in matrix form where the matrices or vectors mostly contains zeros. Scipy has a good support for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc).

There are many possible strategies for storing sparse matrices in an efficient way. Some of the most common are the so-called coordinate form (COO), list of list (LIL) form,  and compressed-sparse column CSC (and row, CSR). Each format has some advantanges and disadvantages. Most computational algorithms (equation solving, matrix-matrix multiplication, etc) can be efficiently implemented using CSR or CSC formats, but they are not so intuitive and not so easy to initialize. So often a sparse matrix is initially created in COO or LIL format (where we can efficiently add elements to the sparse matrix data), and then converted to CSC or CSR before used in real calcalations.

For more information about these sparse formats, see e.g. http://en.wikipedia.org/wiki/Sparse_matrix

When we create a sparse matrix we have to choose which format it should be stored in. For example, 


```python
from scipy import sparse
```


```python
# dense matrix
M = np.array([[1,0,0,0], [0,3,0,0], [0,1,1,0], [1,0,0,1]]); M
```


```python
# convert from dense to sparse
A = sparse.csr_matrix(M)
A
```


```python
# convert from sparse to dense
A.todense()
```

More efficient way to create sparse matrices: create an empty matrix and populate with using matrix indexing (avoids creating a potentially large dense matrix)


```python
A = sparse.lil_matrix((4,4)) # empty 4x4 sparse matrix
A[0,0] = 1
A[1,1] = 3
A[2,2] = A[2,1] = 1
A[3,3] = A[3,0] = 1
A
```


```python
A.todense()
```

Converting between different sparse matrix formats:


```python
A
```


```python
A = sparse.csr_matrix(A); A
```


```python
A = sparse.csc_matrix(A); A
```

We can compute with sparse matrices like with dense matrices:


```python
A.todense()
```


```python
(A * A).todense()
```


```python
A.todense()
```


```python
A.dot(A).todense()
```


```python
v = np.array([1,2,3,4])[:,np.newaxis] 
v
```


```python
# sparse matrix - dense vector multiplication
A * v
```


```python
# same result with dense matrix - dense vector multiplcation
A.todense() * v
```

3205d6aee856efb0dfebf49165c0366e32ae3e2d
ipynb
Jupyter Notebook
scipy/scipy_general.ipynb
emilleishida/School2019
d5141df58d3240ceb0037e8f084f60fc2ef9b4c1
MIT
2019-02-02T08:33:10.000Z
2021-05-04T17:38:26.000Z
scipy/scipy_general.ipynb
emilleishida/School2019
d5141df58d3240ceb0037e8f084f60fc2ef9b4c1
MIT
2019-04-01T11:39:40.000Z
2019-04-09T12:53:33.000Z
scipy/scipy_general.ipynb
emilleishida/School2019
d5141df58d3240ceb0037e8f084f60fc2ef9b4c1
MIT
2019-03-25T18:37:26.000Z
2021-08-19T16:41:45.000Z
Qwen/Qwen-72B
1. YES
2. YES