matplotlib y SciPy

class: center, middle, inverse, title-slide

# matplotlib y SciPy
### Licenciatura en Ciencias Genómicas, UNAM
### First version: 2021-08-22; Last update: 2021-11-16

---

<table class="table">
<colgroup>
<col style="width: 25%" />
<col style="width: 75%" />
</colgroup>
<thead>
<tr class="row-odd"><th class="head">Subpackage</th>
<th class="head">Description</th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><a class="reference internal" href="../reference/cluster.html#module-scipy.cluster" title="scipy.cluster"><code class="xref py py-obj docutils literal notranslate">cluster</code></a></td>
<td>Clustering algorithms</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/constants.html#module-scipy.constants" title="scipy.constants"><code class="xref py py-obj docutils literal notranslate">constants</code></a></td>
<td>Physical and mathematical constants</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/fftpack.html#module-scipy.fftpack" title="scipy.fftpack"><code class="xref py py-obj docutils literal notranslate">fftpack</code></a></td>
<td>Fast Fourier Transform routines</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/integrate.html#module-scipy.integrate" title="scipy.integrate"><code class="xref py py-obj docutils literal notranslate">integrate</code></a></td>
<td>Integration and ordinary differential equation solvers</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/interpolate.html#module-scipy.interpolate" title="scipy.interpolate"><code class="xref py py-obj docutils literal notranslate">interpolate</code></a></td>
<td>Interpolation and smoothing splines</td>
</tr>
</tbody>
</table>

---
# High-performance Python for crystallographic computing
#### A. Boullea and J. Kiefferb
#### August 2019

https://journals.iucr.org/j/issues/2019/04/00/gj5229/

The Python programming language, combined with the numerical computing library **NumPy** and the scientific computing library **SciPy**, has become the de facto standard for scientific computing in a variety of fields.[...] In the present article it is shown how these approaches can be efficiently used to tackle different problems, with increasing complexity, that are relevant to crystallography:

+ the 2D Laue function
+ scattering from a strained 2D crystal
+ scattering from 3D nanocrystals 
+ diffraction from films and multilayers.
---
<table class="table">
<colgroup>
<col style="width: 25%" />
<col style="width: 75%" />
</colgroup>
<thead>
<tr class="row-odd"><th class="head">Subpackage</th>
<th class="head">Description</th>
</tr>
</thead>
<tbody>
<tr class="row-odd"><td><a class="reference internal" href="../reference/io.html#module-scipy.io" title="scipy.io"><code class="xref py py-obj docutils literal notranslate">io</code></a></td>
<td>Input and Output</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/linalg.html#module-scipy.linalg" title="scipy.linalg"><code class="xref py py-obj docutils literal notranslate">linalg</code></a></td>
<td>Linear algebra</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/ndimage.html#module-scipy.ndimage" title="scipy.ndimage"><code class="xref py py-obj docutils literal notranslate">ndimage</code></a></td>
<td>N-dimensional image processing</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/odr.html#module-scipy.odr" title="scipy.odr"><code class="xref py py-obj docutils literal notranslate">odr</code></a></td>
<td>Orthogonal distance regression</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/optimize.html#module-scipy.optimize" title="scipy.optimize"><code class="xref py py-obj docutils literal notranslate">optimize</code></a></td>
<td>Optimization and root-finding routines</td>
</tr>
</tbody>
</table>

---

<table class="table">
<colgroup>
<col style="width: 25%" />
<col style="width: 75%" />
</colgroup>
<thead>
<tr class="row-odd"><th class="head">Subpackage</th>
<th class="head">Description</th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><a class="reference internal" href="../reference/signal.html#module-scipy.signal" title="scipy.signal"><code class="xref py py-obj docutils literal notranslate">signal</code></a></td>
<td>Signal processing</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/sparse.html#module-scipy.sparse" title="scipy.sparse"><code class="xref py py-obj docutils literal notranslate">sparse</code></a></td>
<td>Sparse matrices and associated routines</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/spatial.html#module-scipy.spatial" title="scipy.spatial"><code class="xref py py-obj docutils literal notranslate">spatial</code></a></td>
<td>Spatial data structures and algorithms</td>
</tr>
<tr class="row-odd"><td><a class="reference internal" href="../reference/special.html#module-scipy.special" title="scipy.special"><code class="xref py py-obj docutils literal notranslate">special</code></a></td>
<td>Special functions</td>
</tr>
<tr class="row-even"><td><a class="reference internal" href="../reference/stats.html#module-scipy.stats" title="scipy.stats"><code class="xref py py-obj docutils literal notranslate">stats</code></a></td>
<td>Statistical distributions and functions</td>
</tr>
</tbody>
</table>

---
# Ejemplo
**DNA diffraction pattern**

https://scipython.com/book/chapter-8-scipy/examples/dna-diffraction-pattern/

Funcion de Bessel

![](https://upload.wikimedia.org/wikipedia/commons/5/5d/Bessel_Functions_%281st_Kind%2C_n%3D0%2C1%2C2%29.svg)
---

```python
import numpy as np
from scipy.special import jn # Función Bessel
#...

x = np.linspace(xmin, xmax, npts)

# Diffraction pattern along each line layer: |Jn(x)|^2
# for n = 0, 1, ..., nlayers-1
*layers = np.array([jn(i, x)**2 for i in range(nlayers)])

# Obtain the indexes of the maxima in each layer
maxi = [(np.diff(np.sign(np.diff(layers[i,:]))) < 0).nonzero()[0] + 1
 for i in range(nlayers)]
```

---
Crystallographic photo of Sodium Thymonucleate, Type B. "Photo 51." May 1952. (Large Version)

Original held in the Ava Helen and Linus Pauling Papers.

Creator: Rosalind Franklin, Raymond G. Gosling

Date: May 1952

---
# Ejemplo 
**The SIR epidemic model**

https://scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/

`odeint` Integrate a system of ordinary differential equations.

```python
import numpy as np
from scipy.integrate import odeint 
import matplotlib.pyplot as plt

def deriv(y, t, N, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt
  
# Initial conditions vector
y0 = S0, I0, R0
# Integrate the SIR equations over the time grid, t.
*ret = odeint(deriv, y0, t, args=(N, beta, gamma))
S, I, R = ret.T
```
---

+ **SI**: Enfermedades víricas que causan infección vitalicia, como el VIH.
+ **SIS**: Enfermedades que no confieren inmunidad tras la infección
+ **SIR**: Enfermedades víricas en las que una vez infectado, se obtiene inmunidad vitalicia

---
### Gráfica interactiva de los modelos 
Basado en el código de modelos epidémicos de: https://blog.csdn.net/youcans 
(https://blog.csdn.net/youcans/article/details/117843875)

+ Importamos las librerías

```python
from scipy.integrate import odeint 
import numpy as np 
import matplotlib.pyplot as plt 
from matplotlib.widgets import Slider, RadioButtons  # Interactivo
```
--
+ Definimos las ecuaciones diferenciales

```python
def dySIS(y, t, lamda, mu):  # SI/SIS model
    dy_dt = lamda*y*(1-y)-mu*y 
    return dy_dt
def dySIR(y, t, lamda, mu):  # SIR model, 
    i, s = y
    di_dt = lamda*s*i-mu*i 
    ds_dt = -lamda*s*i 
    return np.array([di_dt,ds_dt])
```
---
+ Definimos los parámetros

```python
number = 1e5 # total number of people
lamda = 0.2 # Daily contact rate, the average number of susceptible persons who are effectively in contact with the sick each day
sigma = 2.5 # Number of contacts during infectious period

mu = lamda/sigma # Daily cure rate, the ratio of the number of patients cured each day to the total number of patients
tEnd = 200 # Forecast date length
t = np.arange(0.0,tEnd,1) # (start,stop,step)

i0 = 1e-4 # Initial value of the proportion of patients
s0 = 1-i0 # Initial value of the proportion of susceptible persons

Y0 = (i0, s0) # Initial value of the differential equation system
```
--
+ Resolvemos las ecuaciones diferenciales con **`SciPy`**

```python
ySI = odeint(dySIS, i0, t, args=(lamda,0)) # SI model
ySIS = odeint(dySIS, i0, t, args=(lamda,mu)) # SIS model
ySIR = odeint(dySIR, Y0, t, args=(lamda,mu)) # SIR model
```

---
+ Graficamos

```python
fig, ax = plt.subplots()
plt.subplots_adjust(left=0.25, bottom=0.4)

plt.title("Comparison among SI, SIS and SIR models")
plt.xlabel('time')
plt.axis([0, tEnd, -0.1, 1.1])

si_plt, = plt.plot(t, ySI,':g', label='i(t)-SI')

sis_plt, = plt.plot(t, ySIS,'--g', label='i(t)-SIS')

sir_i_plt, = plt.plot(t, ySIR[:,0],'-r', label='i(t)-SIR')
sir_s_plt, = plt.plot(t, ySIR[:,1],'-b', label='s(t)-SIR')
sir_r_plt, = plt.plot(t, 1-ySIR[:,0]-ySIR[:,1],'-m', label='r(t)-SIR')

plt.legend(loc='best') 
```
---

---
+ Agregamos las barras interactivas

```python
axcolor = 'lightgoldenrodyellow'

# Generamos el área de las barras
axlambda = plt.axes([0.25, 0.1, 0.65, 0.03], facecolor=axcolor)
axsigma = plt.axes([0.25, 0.18, 0.65, 0.03], facecolor=axcolor)
axi0 = plt.axes([0.25, 0.26, 0.65, 0.03], facecolor=axcolor)

# Agregamos la información
slambda = Slider(axlambda, 'Daily contact rate', 0.1, 1,
               valinit=lamda, color="green")
ssigma = Slider(axsigma, 'Contacts during\ninfectious period', 0.1, 10,
                valinit=sigma)
si0 = Slider(axi0, 'Initial proportion\nof patients', 1e-4, 5e-1,
                valinit=i0, color="orange")
```

---
+ Creamos una función que actualice los datos con base en la información modificada en las barras

```python
def update(val, ):
  
    lamda = slambda.val
    sigma = ssigma.val
    i0 = si0.val

mu = lamda / sigma
    s0 = 1 - i0
    Y0 = (i0, s0)

ySI = odeint(dySIS, i0, t, args=(lamda, 0))  # SI model
    ySIS = odeint(dySIS, i0, t, args=(lamda, mu))  # SIS model
    ySIR = odeint(dySIR, Y0, t, args=(lamda, mu))  # SIR model

si_plt.set_ydata(ySI)
    sis_plt.set_ydata(ySIS)
    sir_i_plt.set_ydata(ySIR[:, 0])
    sir_s_plt.set_ydata(ySIR[:, 1])
    sir_r_plt.set_ydata(1 - ySIR[:, 0] - ySIR[:, 1])

fig.canvas.draw_idle()
```
---
+ Aplicamos la función

```python
slambda.on_changed(update)
ssigma.on_changed(update)
si0.on_changed(update)
```

---
+ Agregamos botones para ver solo un tipo de módelo

```python
rax = plt.axes([0.025, 0.5, 0.15, 0.15], facecolor=axcolor)
radio = RadioButtons(rax, ('SI', 'SIS', 'SIR'), active=0)

lines = {'SI':[si_plt], 'SIS':[sis_plt],
             'SIR':[sir_i_plt, sir_s_plt, sir_r_plt]}

def select_model(label):

for line_m in lines[label]:
        line_m.set_alpha(1)

for others in set(lines.keys()) - set([label]):
        for line_m in lines[others]:
            line_m.set_alpha(0)
    fig.canvas.draw_idle()

radio.on_clicked(select_model)
plt.show()
```
---