---
title: "SciPy"
subtitle: "Lecture 07"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2023"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute:
  echo: true
---

```{r config}
#| include: false
options(
  width=80
)

local({
  hook_old <- knitr::knit_hooks$get("error")  # save the old hook
  knitr::knit_hooks$set(error = function(x, options) {
    x = sub("## \n## Detailed traceback:\n.*$", "", x)
    x = sub("Error in py_call_impl\\(.*?\\)\\: ", "", x)
    hook_old(x, options)
  })
})
```

```{python setup}
#| include: false
import scipy
import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(edgeitems=3, linewidth=180, precision=4)
```

## What is SciPy

> Fundamental algorithms for scientific computing in Python

::: {.small}
| Subpackage    | Description                                           |   | Subpackage    | Description            
|:--------------|:------------------------------------------------------|---|:--------------|:-------------------------------------------
| `cluster`     | Clustering algorithms                                 |   | `odr`         | Orthogonal distance regression       
| `constants`   | Physical and mathematical constants                   |   | `optimize`    | Optimization and root-finding routines
| `fftpack`     | Fast Fourier Transform routines                       |   | `signal`      | Signal processing    
| `integrate`   | Integration and ordinary differential equation solvers|   | `sparse`      | Sparse matrices and associated routines
| `interpolate` | Interpolation and smoothing splines                   |   | `spatial`     | Spatial data structures and algorithms
| `io`          | Input and Output                                      |   | `special`     | Special functions
| `linalg`      | Linear algebra                                        |   | `stats`       | Statistical distributions and functions
| `ndimage`     | N-dimensional image processing                        |   | &nbsp;        | &nbsp;

:::


## SciPy vs NumPy

> In an ideal world, NumPy would contain nothing but the array data type and the most basic operations: indexing, sorting, reshaping, basic elementwise functions, etc. All numerical code would reside in SciPy. However, one of NumPy's important goals is compatibility, so NumPy tries to retain all features supported by either of its predecessors. Thus, NumPy contains some linear algebra functions and Fourier transforms, even though these more properly belong in SciPy. In any case, SciPy contains more fully-featured versions of the linear algebra modules, as well as many other numerical algorithms. If you are doing scientific computing with Python, you should probably install both NumPy and SciPy. Most new features belong in SciPy rather than NumPy.

::: {.aside}
From [scipy.org/faq/](https://scipy.org/faq/)
:::


# Example 1 <br/> k-means clustering

## Data

::: {.small}
```{python}
rng = np.random.default_rng(seed = 1234)

cl1 = rng.multivariate_normal([-2,-2], [[1,-0.5],[-0.5,1]], size=100)
cl2 = rng.multivariate_normal([1,0], [[1,0],[0,1]], size=150)
cl3 = rng.multivariate_normal([3,2], [[1,-0.7],[-0.7,1]], size=200)

pts = np.concatenate((cl1,cl2,cl3))
```
:::

. . .

```{python}
#| echo: false
#| out-width: 75%
#| fig-align: center
plt.cla()
plt.scatter(cl1[:,0], cl1[:,1], c="r", marker = ".")
plt.scatter(cl2[:,0], cl2[:,1], c="b", marker = "*")
plt.scatter(cl3[:,0], cl3[:,1], c="c", marker = "D")
plt.show()
```

## k-means clustering

:::: {.columns .small}
::: {.column width='50%'}
```{python}
from scipy.cluster.vq import kmeans
ctr, dist = kmeans(pts, 3)
```
```{python}
ctr
dist
```
:::

::: {.column width='50%' .fragment}
```{python}
cl1.mean(axis=0)
cl2.mean(axis=0)
cl3.mean(axis=0)
```
:::
::::

. . .

```{python}
#| echo: false
#| out-width: 66%
#| fig-align: center
plt.cla()
plt.scatter(cl1[:,0], cl1[:,1], c="r", marker = ".")
plt.scatter(cl2[:,0], cl2[:,1], c="b", marker = "*")
plt.scatter(cl3[:,0], cl3[:,1], c="c", marker = "D")
plt.scatter(ctr[:,0], ctr[:,1], c="k", marker = "x", s = 500, linewidths=3.5)
plt.show()
```



## k-means distortion plot

> The mean (non-squared) Euclidean distance between the observations passed and the centroids generated.

::: {.small}
```{python}
ks = range(1,8)
dists = [kmeans(pts, k)[1] for k in ks]
```
```{python}
np.array(dists).reshape(-1)
```

```{python}
#| echo: false
#| out-width: 66%
#| fig-align: center
plt.cla()
p = plt.plot(ks, dists, "-ok")
plt.show()
```
:::


# Example 2 <br/> Numerical integration

## Basic functions

For general numeric integration in 1D we use `scipy.integrate.quad()`, which takes as arguments the function to be integrated and the lower and upper bounds of the integral.

```{python}
from scipy.integrate import quad
```

```{python}
quad(lambda x: x, 0, 1)
```

. . .

```{python}
quad(np.sin, 0, np.pi)
```

. . .

```{python}
quad(np.sin, 0, 2*np.pi)
```

. . .

```{python}
quad(np.exp, 0, 1)
```


## Normal PDF

The PDF for a normal distribution is given by,

$$ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} \exp\left(-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2  \right) $$


```{python}
def norm_pdf(x, μ, σ):
  return (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)
```

. . .

```{python}
norm_pdf(0, 0, 1)
norm_pdf(np.Inf, 0, 1)
norm_pdf(-np.Inf, 0, 1)
```

## Checking the PDF

We can check that we've implemented a valid pdf by integrating the pdf from $-\inf$ to $\inf$,

. . .

```{python error=TRUE}
quad(norm_pdf, -np.inf, np.inf)
```

. . .

```{python error=TRUE}
quad(lambda x: norm_pdf(x, 0, 1), -np.inf, np.inf)
```

. . .

```{python error=TRUE}
quad(lambda x: norm_pdf(x, 17, 12), -np.inf, np.inf)
```


## Truncated normal PDF

$$
f(x) = \begin{cases}
\frac{c}{\sigma \sqrt{2 \pi}} \exp\left(-\frac{1}{2} \left(\frac{x-\mu}{\sigma}\right)^2  \right), & \text{for } a \leq x \leq b \\
0,                                            & \text{otherwise.} \\
\end{cases}
$$

```{python}
def trunc_norm_pdf(x, μ=0, σ=1, a=-np.inf, b=np.inf):
  if (b < a):
      raise ValueError("b must be greater than a")
  
  x = np.asarray(x).reshape(-1)
  full_pdf = (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)
  full_pdf[(x < a) | (x > b)] = 0
  return full_pdf
```


## Testing trunc_norm_pdf


```{python}
trunc_norm_pdf(0, a=-1, b=1)
trunc_norm_pdf(2, a=-1, b=1)
trunc_norm_pdf(-2, a=-1, b=1)
trunc_norm_pdf([-2,1,0,1,2], a=-1, b=1)
```

. . .

```{python}
quad(lambda x: trunc_norm_pdf(x, a=-1, b=1), -np.inf, np.inf)

quad(lambda x: trunc_norm_pdf(x, a=-3, b=3), -np.inf, np.inf)
```


## Fixing trunc_norm_pdf

::: {.small}
```{python}
def trunc_norm_pdf(x, μ=0, σ=1, a=-np.inf, b=np.inf):
  if (b < a):
      raise ValueError("b must be greater than a")
  x = np.asarray(x).reshape(-1)
  
  nc = 1 / quad(lambda x: norm_pdf(x, μ, σ), a, b)[0]
  
  full_pdf = nc * (1/(σ * np.sqrt(2*np.pi))) * np.exp(-0.5 * ((x - μ)/σ)**2)
  full_pdf[(x < a) | (x > b)] = 0
  
  return full_pdf
```
:::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
trunc_norm_pdf(0, a=-1, b=1)
trunc_norm_pdf(2, a=-1, b=1)
trunc_norm_pdf(-2, a=-1, b=1)
trunc_norm_pdf([-2,1,0,1,2], a=-1, b=1)
```
:::

::: {.column width='50%' .fragment}
```{python}
quad(lambda x: trunc_norm_pdf(x, a=-1, b=1), -np.inf, np.inf)

quad(lambda x: trunc_norm_pdf(x, a=-3, b=3), -np.inf, np.inf)
```
:::
::::


## Multivariate normal

$$
f(\bf{x}) = \det{(2\pi\Sigma)}^{-1/2} \exp{\left(-\frac{1}{2} (\bf{x}-\mu)^T \Sigma^{-1}(\bf{x}-\mu) \right)}
$$

::: {.small}
```{python}
def mv_norm(x, μ, Σ):
  x = np.asarray(x)
  μ = np.asarray(μ)
  Σ = np.asarray(Σ)
  
  return ( np.linalg.det(2*np.pi*Σ)**(-0.5) * 
           np.exp(-0.5 * (x - μ).T @ np.linalg.solve(Σ, (x-μ)) ) )
```
:::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
norm_pdf(0,0,1)
mv_norm([0], [0], [[1]])
```
:::

::: {.column width='50%'}
```{python}
mv_norm([0,0], [0,0], [[1,0],[0,1]])
mv_norm([0,0,0], [0,0,0], 
        [[1,0,0],[0,1,0],[0,0,1]])
```
:::
::::

## 2d & 3d numerical integration

are supported by `dblquad()` and `tplquad()` respectively (see `nquad()` for higher dimensions)

```{python cache=TRUE}
from scipy.integrate import dblquad, tplquad
```

```{python cache=TRUE}
dblquad(lambda y, x: mv_norm([x,y], [0,0], np.identity(2)), 
        a=-np.inf, b=np.inf, 
        gfun=lambda x: -np.inf,   hfun=lambda x: np.inf)
```

. . .

```{python cache=TRUE}
tplquad(lambda z, y, x: mv_norm([x,y,z], [0,0,0], np.identity(3)),
        a=0, b=np.inf, 
        gfun=lambda x:   0, hfun=lambda x:   np.inf,
        qfun=lambda x,y: 0, rfun=lambda x,y: np.inf)
```



# Example 3 <br/> (Very) Basic optimization

## Scalar function minimization

:::: {.columns .small}
::: {.column width='50%'}
```{python}
def f(x):
    return x**4 + 3*(x-2)**3 - 15*(x)**2 + 1
```

```{python}
#| echo: false
#| out-width: 100%
x = np.linspace(-8, 5, 100)
plt.plot(x, f(x))
```
:::

::: {.column width='50%' .fragment}
```{python}
from scipy.optimize import minimize_scalar
minimize_scalar(f, method="Brent")
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
minimize_scalar(f, method="bounded", bounds=[0,6])
```
:::

::: {.column width='50%'}
```{python}
minimize_scalar(f, method="bounded", bounds=[-8,6])
```
:::
::::


## Results

```{python}
res = minimize_scalar(f)
type(res)
dir(res)
res.success
res.x
res.fun
```


## More details

::: {.small}
```{python}
from scipy.optimize import show_options
show_options(solver="minimize_scalar")
```
:::


## Local minima

:::: {.columns .small}
::: {.column width='50%'}
```{python}
def f(x):
  return -np.sinc(x-5)
```


```{python echo=FALSE, out.width="90%"}
x = np.linspace(-20, 20, 500)
plt.cla()
p = plt.plot(x, f(x));
plt.show()
```
:::

::: {.column width='50%' .fragment}
```{python}
res = minimize_scalar(f); res
```


```{python echo=FALSE, out.width="90%"}
x = np.linspace(-20, 20, 500)
plt.cla()
p = plt.plot(x, f(x));
plt.axvline(res.x, c='red')
plt.show()
```
:::
::::


## Random starts

::: {.small}
```{python}
#| output-location: column
rng = np.random.default_rng(seed=1234)

lower = rng.uniform(-20, 20, 100)
upper = lower + 1

sols = [minimize_scalar(f, bracket=(l,u)) 
        for l,u in zip(lower, upper)]
funs = [sol.fun for sol in sols]

best = sols[np.argmin(funs)]
best
```
:::

```{python}
#| echo: false
#| out-width: 75%
#| fig-align: center
plt.cla()
p = plt.plot(x, f(x));
plt.axvline(best.x, c='red')
plt.show()
```


## Back to Rosenbrock's function

$$
f(x,y) = (1-x)^2 + 100(y-x^2)^2
$$ 

::: {.small}
```{python}
def f(x):
  return (1-x[0])**2 + 100*(x[1]-x[0]**2)**2
```
:::

```{python}
#| echo: false
from scipy.optimize import minimize
```

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
minimize(f, [0,0])
```

:::

::: {.column width='50%'}
```{python}
minimize(f, [-1,-1])
```

:::
::::


# Example 4 <br/> Spatial Tools

## Nearest Neighbors

:::: {.small}
```{python}
#| output-location: column
rng = np.random.default_rng(seed=12345)
pts = rng.multivariate_normal(
  [0,0], [[1,.8],[.8,1]], 
  size=10
)
pts
```
:::

```{python}
#| echo: false
#| out-width: 75%
#| fig-align: center
plt.cla()
plt.scatter(pts[:,0], pts[:,1], c='w')

for i in range(10):
    plt.annotate(str(i), (pts[i,0], pts[i,1]), weight="bold", size=16, ha='center', va='center')

plt.show()
```


## KD Trees

```{python}
#| echo: false
#| out-width: 75%
#| fig-align: center
plt.cla()
plt.scatter(pts[:,0], pts[:,1], c='w')

for i in range(10):
    plt.annotate(str(i), (pts[i,0], pts[i,1]), weight="bold", size=16, ha='center', va='center')

plt.show()
```


::: {.small}
```{python}
from scipy.spatial import KDTree
kd = KDTree(pts)
```
:::

. . .

:::: {.columns .small}
::: {.column width='50%'}

```{python}
dist, i = kd.query(pts[6,:], k=3)
i
dist
```
:::

::: {.column width='50%'}
```{python}
dist, i = kd.query(pts[2,:], k=5)
i
```
:::
::::


## Convex hulls


::: {.small}
```{python}
from scipy.spatial import ConvexHull
hull = ConvexHull(pts)
```

```{python}
#| fig-align: center
#| out-width: 75%
hull.vertices
scipy.spatial.convex_hull_plot_2d(hull)
```
:::

## Delaunay triangulations

:::: {.small}
```{python}
from scipy.spatial import Delaunay
tri = Delaunay(pts)
```

```{python}
#| fig-align: center
#| out-width: 66%
tri.simplices.T
scipy.spatial.delaunay_plot_2d(tri)
```
::::


## Voronoi diagrams

:::: {.small}
```{python}
from scipy.spatial import Voronoi
vor = Voronoi(pts)
```
```{python}
#| fig-align: center
#| out-width: 66%
vor.vertices.T
scipy.spatial.voronoi_plot_2d(vor)
```
::::


# Example 5 <br/> statistics

## Distributions 

Implements classes for 104 continuous and 19 discrete distributions,

* `rvs` -  Random Variates

* `pdf` -  Probability Density Function

* `cdf` -  Cumulative Distribution Function

* `sf` -  Survival Function (1-CDF)

* `ppf` -  Percent Point Function (Inverse of CDF)

* `isf` -  Inverse Survival Function (Inverse of SF)

* `stats` -  Return mean, variance, (Fisher’s) skew, or (Fisher’s) kurtosis

* `moment` -  non-central moments of the distribution


## Basic usage

```{python}
from scipy.stats import norm, gamma, binom, uniform
```

. . .

```{python}
norm().rvs(size=5)
```

. . .

```{python}
uniform.pdf([0,0.5,1,2])
```

. . .

```{python}
binom.mean(n=10, p=0.25)
binom.median(n=10, p=0.25)
```

. . .

```{python}
gamma(a=1,scale=1).stats()
norm().stats(moments="mvsk")
```


## Freezing

Model parameters can be passed to any of the methods directory, or a distribution can be constructed using a specific set of parameters, which is known as freezing.

:::: {.columns .small}
::: {.column width='50%'}
```{python}
norm_rv = norm(loc=-1, scale=3)
norm_rv.median()
```

```{python}
unif_rv = uniform(loc=-1, scale=2)
unif_rv.cdf([-2,-1,0,1,2])
unif_rv.rvs(5)
```
:::

::: {.column width='50%' .fragment}
```{python}
#| fig-align: center
#| out-width: 66%
g = gamma(a=2, loc=0, scale=1.2)

x = np.linspace(0, 10, 100)
plt.plot(x, g.pdf(x), "k-")
plt.axvline(x=g.mean(), c="r")
plt.axvline(x=g.median(), c="b")
```
:::
::::


## MLE

Maximum likelihood estimation is possible via the `fit()` method,

```{python}
x = norm.rvs(loc=2.5, scale=2, size=1000, random_state=1234)
norm.fit(x)
norm.fit(x, loc=2.5) # provide a guess for the parameter
```

. . .

```{python}
x = gamma.rvs(a=2.5, size=1000)
gamma.fit(x) # shape, loc, scale

y = gamma.rvs(a=2.5, loc=-1, scale=2, size=1000)
gamma.fit(y) # shape, loc, scale
```