text:
# Parabolic PDEs

A classic example of a parabolic partial differential equation (PDE) is the one-dimensional unsteady heat equation:
\begin{equation}
\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial t^2}
\end{equation}
where $T(x,t)$ is the temperature varying in space and time, and $\alpha$ is the thermal diffusivity: $\alpha = k / (\rho c_p)$, which is a constant.

We can solve this using finite differences to represent the spatial derivatives and time derivatives separately.
First, let's rearrange the PDE slightly:
\begin{equation}
\frac{\partial^2 T}{\partial x^2} = \frac{1}{\alpha} \frac{\partial T}{\partial t}
\end{equation}

## Explicit scheme

Let's use a *central difference* for the spatial derivative with a spacing of $\Delta x$, and a *forward difference* for the time derivative with a time-step size of $\Delta t$. With these choices, we can obtain an approximation to the PDE that applies at time $t^k$ and spatial location $x_i$:
\begin{equation}
\frac{T_{i-1}^k - 2 T_i^k + T_{i+1}^k}{\Delta x^2} = \frac{1}{\alpha} \left( \frac{T_i^{k+1} - T_i^k}{\Delta t} \right)
\end{equation}
where $T_i^k$ is the temperature at time $t^k$ and spatial location $x_i$. The following figure shows the stencil of points involved in the PDE, for a domain with five points in the $x$-direction.

:::{figure-md} fig-stencil-explicit

Stencil for explicit solution to heat equation
:::

To solve the heat equation for a one-dimensional domain over $0 \leq x \leq L$, we will need both initial conditions at $t = 0$ and boundary conditions at $x=0$ and $x=L$ (for all time). In terms of our nodal values, this means we need $T_i^{k=1}$ for $i = 1 \ldots n$, where $n$ is the number of points, as well as information about $T_1^k$ and $T_n^k$ for all times $k$.

We can rearrange the above equation to obtain our recursion formula:
\begin{equation}
T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \frac{\alpha \Delta t}{\Delta x^2} + T_i^k  \left( 1 - 2 \frac{\alpha \Delta t}{\Delta x^2} \right) \;.
\end{equation}
This is an **explicit** scheme in time, similar to the Forward Euler method we used for ordinary differential equations, and like that method it may have stability issues. The combination of terms we see repeated is also known as the Fourier number: $\text{Fo} = \frac{\alpha \Delta t}{\Delta x^2}$, and governs the stability.
We can rewrite the recursion formula using this:
\begin{equation}
T_i^{k+1} = \left( T_{i+1}^k + T_{i-1}^k \right) \text{Fo} + T_i^k  \left( 1 - 2 \text{Fo} \right) \;.
\end{equation}

The term in parentheses there must be greater than or equal to zero for stability ($1 - 2 \text{Fo} \geq 0$); if not, the solution may become unstable and blow up. This gives us some conditions on our choice of step sizes:
\begin{align}
1 - 2 \text{Fo} &\geq 0 \\
1 & \geq 2 \text{Fo} \\
\text{Fo} &\leq \frac{1}{2} \\
\frac{\alpha \Delta t}{\Delta x^2} &\leq \frac{1}{2}
\end{align}
This is the stability criterion for the explicit method: the Fourier number must be smaller than 0.5.
For a given thermal diffusivity and chosen spatial step size, this also gives us the limit on the time-step size: $\Delta t \leq \Delta x^2 / (2 \alpha)$.

Let's look at an example where the initial temperature is 200, the temperature at the boundaries are 50, the thermal diffusivity is $\alpha = 2.3 \times 10^{-1}$ m$^2 /$ s, and $L = 1$.
In other words,
\begin{align}
T(x, t=0) &= 200 \\
T(x=0, t) &= 50 \\
T(x=L, t) &= 50
\end{align}
We'll integrate out to $t = 1$, using a Fourier number of 0.25 to be comfortably below the stability limit (`Fo = 0.25`):

```matlab
clear all

dx = 0.1;
alpha = 2.3e-1;

% for stability, set the Fourier number at 0.25 (half the stability limit of 0.5)
Fo = 0.25;
% then choose the time step based on the Fourier number
dt = Fo * dx^2 / alpha;

x = [0 : dx : 1]; n = length(x);
t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    for i = 1 : n
        if i == 1
            T(k+1, 1) = 50;
        elseif i == n
            T(k+1, n) = 50;
        else
            T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
        end
    end
    plot(x, T(k+1,:))
    axis([0 1 50 200]);
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% If you are working interactively, you can use this to make a movie in Matlab
%fig = figure;
%movie(fig, F, 2)

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i));
    [imind,cm] = rgb2ind(im,256);
    % Write to GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
    end
end
```

:::{figure-md} fig-transient-heat

Animated solution to 1D transient heat transfer PDE
:::

This shows the temperature decaying exponentially from the initial conditions, constrained by the boundary conditions.

What happens if we tried to use a Fourier number larger than 0.5, or arbitrarily chose a time-step size that was too large (and resulted in $\text{Fo} > 0.5$)?

```matlab
clear all

dx = 0.1;
alpha = 2.3e-1;

%% Purposely choose a Fourier number that is past the stability limit:
Fo = 0.75;
dt = Fo * dx^2 / alpha;

x = [0 : dx : 1]; n = length(x);
t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    for i = 1 : n
        if i == 1
            T(k+1, 1) = 50;
        elseif i == n
            T(k+1, n) = 50;
        else
            T(k+1, i) = (T(k,i+1) + T(k,i-1))*Fo + T(k,i)*(1 - 2*Fo);
        end
    end
    plot(x, T(k+1,:))
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% If you are working interactively, you can use this to make a movie in Matlab
%fig = figure;
%movie(fig, F, 2)

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_unstable_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i));
    [imind,cm] = rgb2ind(im,256);
    % Write to GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',2, 'DelayTime',1e-3);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
    end
end
```

:::{figure-md} fig-unstable-solution

Animated unstable solution
:::

In this case, the solution becomes unstable and blows up, leading to unphysical results.

For this explicit scheme, the choice of $\Delta t$ is limited by the stability criterion. This means that we may be stuck using a small time-step size.

Rather than being forced to use a very small time-step size, we can also explore *implicit schemes* that are unconditionally stable.

## Implicit scheme

```matlab
clear all

alpha = 2.3e-1;
dx = 0.1;
x = [0 : dx : 1]; n = length(x);

% choose a Fourier number that is deliberately past the explicit method stability limit
Fo = 0.75;
dt = Fo * dx^2 / alpha;

t = [0 : dt : 1]; m = length(t);

T = zeros(m, n);

% Initial conditions
T(1,:) = 200;

plot(x, T(1,:))
axis([0 1 50 200]);
xlabel('Distance'); ylabel('Temperature');
F(1) = getframe(gcf);

for k = 1 : m - 1
    A = zeros(n,n);
    b = zeros(n,1);
    for i = 1 : n
        if i == 1
            A(1,1) = 1;
            b(1) = 50;
        elseif i == n
            A(n,n) = 1;
            b(n) = 50;
        else
            A(i,i-1) = Fo;
            A(i,i) = -2*Fo - 1;
            A(i,i+1) = Fo;
            b(i) = -T(k,i);
        end
    end

    T(k+1, :) = A \ b;
    plot(x, T(k+1,:))
    axis([0 1 50 200]);
    xlabel('Distance'); ylabel('Temperature');
    F(k+1) = getframe(gcf);
end
close

%% This generates a GIF of the results (for use in Jupyter Notebook)
filename = 'parabolic_implicit_animated.gif';
for i = 1 : length(F)
    im = frame2im(F(i));
    [imind,cm] = rgb2ind(im,256);
    % Write to GIF
    if i == 1
        imwrite(imind,cm,filename,'gif', 'Loopcount',inf, 'DelayTime',1e-3);
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append', 'DelayTime',1e-3);
    end
end
```

:::{figure-md} fig-implicit-solution

Solution to 1D heat equation with implicit method. Fo = 0.75
:::


text:
# Nonlinear Equations

We want to find a root of the nonlinear function $f$ using different methods.

1. Bisection method
2. Newton method
3. Chord method
4. Secant method
5. Fixed point iterations

```python
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import sympy as sym
```

```python
t = sym.symbols('t')

f_sym = t/8. * (63.*t**4 - 70.*t**2. +15.) # Legendre polynomial of order 5

f_prime_sym = sym.diff(f_sym,t)

f = sym.lambdify(t, f_sym, 'numpy')
f_prime = sym.lambdify(t,f_prime_sym, 'numpy')

phi = lambda x : 63./70.*x**3 + 15./(70.*x)
#phi = lambda x : 70.0/15.0*x**3 - 63.0/15.0*x**5
#phi = lambda x : sqrt((63.*x**4 + 15.0)/70.)

# Let's plot
n = 1025

x = linspace(-1,1,n)
c = zeros_like(x)

_ = plot(x,f(x))
_ = plot(x,c)
_ = grid()
```

```python
# Initial data for the variuos algorithms

# interval in which we seek the solution
a = 0.7
b = 1.

# initial points
x0 = (a+b)/2.0
x00 = b
```

```python
# stopping criteria
eps = 1e-10
n_max = 1000
```

## Bisection method

$$
x^k = \frac{a^k+b^k}{2}
$$
```
                   if (f(a_k) * f(x_k)) < 0:
                      b_k1 = x_k
                      a_k1 = a_k
                   else:
                      a_k1 = x_k
                      b_k1 = b_k
```

```python
def bisect(f,a,b,eps,n_max):
    assert f(a)*f(b) < 0
    a_new = a
    b_new = b
    x = mean([a,b])
    err = eps + 1.
    errors = [err]
    it = 0
    while (err > eps and it < n_max):
        if ( f(a_new) * f(x) < 0 ):
            # root in (a_new,x)
            b_new = x
        else:
            # root in (x,b_new)
            a_new = x

        x_new = mean([a_new,b_new])

        #err = 0.5 *(b_new -a_new)
        err = abs(f(x_new))
        #err = abs(x-x_new)

        errors.append(err)
        x = x_new
        it += 1

    semilogy(errors)
    print(it)
    print(x)
    print(err)

    return errors

%time errors_bisect = bisect(f,a,b,eps,n_max)
```

```python
# is the number of iterations coherent with the theoretical estimation?
```

In order to find out other methods for solving non-linear equations, let's compute the Taylor's series of $f(x^k)$ up to the first order

$$
f(x^k) \simeq f(x^k) + (x-x^k)f^{\prime}(x^k)
$$
which suggests the following iterative scheme
$$
x^{k+1} = x^k - \frac{f(x^k)}{f^{\prime}(x^k)}
$$

The following methods are obtained applying the above scheme where

$$
f^{\prime}(x^k) \approx q^k
$$

## Newton's method
$$
q^k = f^{\prime}(x^k)
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q^k}
$$

```python
def newton(f,f_prime,x0,eps,n_max):

    it = 0

    err = eps + 1.
    errors = [err]

    x = x0
    q = f_prime(x)

    err = abs(f(x))
    errors.append(err)

    while (it < n_max and f(x) > eps):

        assert f_prime(x) != 0

        q_next = f_prime(x)
        x_next = x - f(x)/q

        err = abs(f(x_next))
        errors.append(err)

        x = x_next
        q = q_next

        it += 1

    semilogy(errors)
    print(it)
    print(x)
    print(err)

    return errors

%time errors_newton = newton(f,f_prime,1.1,eps,n_max)
```

## Chord method

$$
q^k \equiv q = \frac{f(b)-f(a)}{b-a}
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q}
$$

```python
def chord(f,a,b,x0,eps,n_max):
    pass # TODO

errors_chord = chord (f,a,b,x0,eps,n_max)
```

## Secant method

$$
q^k = \frac{f(x^k)-f(x^{k-1})}{x^k - x^{k-1}}
$$

$$
x^{k+1} = x^k - \frac{f(x^k)}{q^k}
$$

Note that this algorithm requirs **two** initial points

```python
def secant(f,x0,x00,eps,n_max):
    pass # TODO

errors_secant = secant(f,x0,x00,eps,n_max)
```

## Fixed point iterations

$$
f(x)=0 \to x-\phi(x)=0
$$

$$
x^{k+1} = \phi(x^k)
$$

```python
def fixed_point(phi,x0,eps,n_max):
    pass # TODO

errors_fixed = fixed_point(phi,0.3,eps,n_max)
```

## Comparison

```python
# plot the error convergence for the methods
loglog(errors_bisect, label='bisect')
loglog(errors_chord, label='chord')
loglog(errors_secant, label='secant')
loglog(errors_newton, label ='newton')
loglog(errors_fixed, label ='fixed')
_ = legend()
```

```python
# Let's compare the scipy implmentation of Newton's method with our..
```

```python
import scipy.optimize as opt
%time opt.newton(f, 1.0, f_prime, tol = eps)
```

ext: ipynb


text:
[The 1st IAA-CSIC Severo Ochoa School on Statistics, Data Mining and Machine Learning](https://www.granadacongresos.com/sostat)

[Zeljko Ivezic, University of Washington](http://faculty.washington.edu/ivezic/)

[This notebook](https://github.com/carmensg/IAA_School2019/tree/master/lectures/Day3-ZeljkoIvezic/notebooks/dimensionality_reduction.ipynb)
### November 6, 2019

# Lecture 4: Dimensionality Reduction

### o Introduction to multi-dimensional distributions
### o Principal Component Analysis (PCA)
### o Independent Component Analysis (ICA)
### o Non-negative matrix factorization
### o Manifold Learning (example: Local Linear Embedding)

```python
# first things first...
import numpy as np
from matplotlib import pyplot as plt
from scipy import stats
from scipy import integrate
from scipy.stats import norm
from scipy.stats import cauchy
# astroML tools
from astroML.plotting import hist
from astroML.plotting.mcmc import convert_to_stdev
from astroML.plotting import setup_text_plots
from astroML.plotting.mcmc import convert_to_stdev
setup_text_plots(fontsize=8, usetex=True)
import warnings; warnings.simplefilter('ignore')
```

## Dimensionality Reduction##

Fitting and overfitting get worse with ''curse of dimensionality'' Bellman 1961

Think about a hypersphere. Its volume is  given by

\begin{equation}
  V_D(r) = \frac{2r^D\pi^{D/2}}{D\  \Gamma(D/2)},
\end{equation}

where $\Gamma(z)$ is the complete gamma function, $D$ is the dimension, and $r$ the radius of the sphere.

If you populated a hypercube of size $2r$ how much data would be enclosed by the hypersphere
- as $D$ increases the fractional volume enclosed by the hypersphere goes to 0!

For example: the SDSS comprises a sample of 357 million sources.
- each source has 448 measured attributes
- selecting just 30 (e.g., magnitude, size..) and normalizing the data range $-1$ to $1$

probability of having one of the 357 million sources reside within a unit hypersphere 1 in 1.4$\times 10^5$.

```python
%matplotlib inline
import numpy as np
import scipy.special as sp
from matplotlib import pyplot as plt

def unitVolume(dimension, radius=1.):
    return 2*(radius**dimension *np.pi**(dimension/2.))/(dimension*sp.gamma(dimension/2.))

dim = np.linspace(1,100)

#------------------------------------------------------------
# Plot the results
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(dim,unitVolume(dim)/2.**dim)

ax.set_yscale('log')
ax.set_xlabel('$Dimension$')
ax.set_ylabel('$Volume$')

plt.show()
```

```python
D=30;unitVolume(D)/(2.**D)
```

    2.041026339664143e-14

## What can we do to reduce the dimensions?

### Principal component analysis or the Karhunen-Loeve transform

Points are  correlated along a particular direction which doesn't align with the initial choice of axes.
- we should rotate our axes to align with this correlation.
- rotation preserves the relative ordering of data

Choose  rotation to maximize the ability to discriminate between the data points
*   first axis, or <u>principal component</u>, is direction of maximal variance
*   second principal component is orthogonal to the first component and maximizes the residual variance
*   ...

## Derivation of principal component analyses##

Set of data $X$: $N$ observations by $K$ measurements

Center data by subtracting the mean **Why?**

The covariance is

>$
C_X=\frac{1}{N-1}X^TX,
$

$N-1$ as the sample covariance matrix

We want a projection, $R$,  aligned with the directions of maximal variance ($Y= X R$) with covariance

>$
C_{Y} = R^T X^T X R = R^T C_X R
$

Derive  principal component by maximizing its variance (using Lagrange multipliers and constraint)

> $
\phi(r_1,\lambda_1) = r_1^TC_X r_1 - \lambda_1(r_1^Tr_1-1).
$

derivative of $\phi(r_1,\lambda)$ with respect to $r_1$ set to 0

> $
C_Xr_1 - \lambda_1 r_1 = 0.
$

$\lambda_1$ is the root of the equation $\det(C_X -
\lambda_1 {\bf I})=0$ and the largest eigenvalue

>$
\lambda_1 =  r_1^T C_X r_1
$

Other  principal components  derived by
applying additional constraint that components are uncorrelated (e.g., $r^T_2 C_X r_1 = 0$).

## Lagrangian mulitpliers

## Computation of principal components##

Common approach is eigenvalue decomposition of the covariance or correlation matrix,
or singular value decomposition (SVD) of the data matrix

** SVD given by**

>$
U \Sigma V^T = \frac{1}{\sqrt{N - 1}} X,
$

columns of $U$ are  _left-singular vectors_

columns of $V$ are the _right-singular vectors_

The columns of $U$ and $V$ form orthonormal bases ($U^TU = V^TV = I$)

Covariance matrix is

> $
\begin{eqnarray}
  C_X &=& \left[\frac{1}{\sqrt{N - 1}}X\right]^T \left[\frac{1}{\sqrt{N - 1}}X\right]\nonumber\\
      &=& V \Sigma U^T U \Sigma V^T\nonumber\\
      &=& V \Sigma^2 V^T.
\end{eqnarray}
$

right singular vectors $V$ are the principal components. We can calculate principal components from the SVD of $X$ - we dont need $C_X$.

Singular value decomposition (SVD) can factorize an N x K matrix into $U \Sigma V^T$. There are different conventions for computing the SVD in the literature, and this figure illustrates the convention used in this text. The matrix of singular values $\Sigma$ is always a square matrix of size [R x R] where R = min(N, K). The shape of the resulting U and V matrices depends on whether N or K is larger. The columns of the matrix U are called the left-singular vectors, and the columns of the matrix V are called the right-singular vectors. The columns are orthonormal bases, and satisfy $U^T U = V^T V = I$.

## Preparing data for PCA##

- Center data by subtracting the mean of each dimension
- For heterogeneous data (e.g., galaxy shape and flux) divide by  variance (whitening). **why?**
- For spectra or images normalize each row so integrated flux of each object is one.

```python
from sklearn.decomposition import PCA
from sklearn.decomposition import RandomizedPCA

from astroML.datasets import sdss_corrected_spectra
from astroML.decorators import pickle_results

#------------------------------------------------------------
# Download data
data = sdss_corrected_spectra.fetch_sdss_corrected_spectra()
spectra = sdss_corrected_spectra.reconstruct_spectra(data)
wavelengths = sdss_corrected_spectra.compute_wavelengths(data)

#----------------------------------------------------------------------
# Compute PCA
def compute_PCA(n_components=5):
#    np.random.seed(500)
    nrows = 500
    ind = np.random.randint(spectra.shape[0], size=nrows)

    spec_mean = spectra[ind].mean(0)
#    spec_mean = spectra[:50].mean(0)

    # PCA: use randomized PCA for speed
    pca = PCA(n_components - 1)
    pca.fit(spectra[ind])
    pca_comp = np.vstack([spec_mean,
                          pca.components_])
    evals = pca.explained_variance_ratio_

    return pca_comp, evals

n_components = 5
decompositions, evals = compute_PCA(n_components)

#----------------------------------------------------------------------
# Plot the results
fig = plt.figure(figsize=(10, 8))
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05,
                    bottom=0.1, top=0.95, hspace=0.05)

titles = 'PCA components'

for j in range(n_components):
    ax = fig.add_subplot(n_components, 2, 2*j+2)

    ax.yaxis.set_major_formatter(plt.NullFormatter())
    ax.xaxis.set_major_locator(plt.MultipleLocator(1000))
    if j < n_components - 1:
        ax.xaxis.set_major_formatter(plt.NullFormatter())
    else:
        ax.set_xlabel(r'wavelength ${\rm (\AA)}$')
    ax.plot(wavelengths, decompositions[j], '-k', lw=1)

    # plot zero line
    xlim = [3000, 7999]
    ax.plot(xlim, [0, 0], '-', c='gray', lw=1)
    ax.set_xlim(xlim)

    # adjust y limits
    ylim = plt.ylim()
    dy = 0.05 * (ylim[1] - ylim[0])
    ax.set_ylim(ylim[0] - dy, ylim[1] + 4 * dy)

    ax2 = fig.add_subplot(n_components, 2, 2*j+1)
    ax2.yaxis.set_major_formatter(plt.NullFormatter())
    ax2.xaxis.set_major_locator(plt.MultipleLocator(1000))
    if j < n_components - 1:
        ax2.xaxis.set_major_formatter(plt.NullFormatter())
    else:
        ax2.set_xlabel(r'wavelength ${\rm (\AA)}$')
    ax2.plot(wavelengths, spectra[j], '-k', lw=1)

    # plot zero line
    ax2.plot(xlim, [0, 0], '-', c='gray', lw=1)
    ax2.set_xlim(xlim)

    if j == 0:
        ax.set_title(titles, fontsize='medium')

    if j == 0:
        label = 'mean'
    else:
        label = 'component %i' % j

    # adjust y limits
    ylim = plt.ylim()
    dy = 0.05 * (ylim[1] - ylim[0])
    ax2.set_ylim(ylim[0] - dy, ylim[1] + 4 * dy)

    ax.text(0.02, 0.95, label, transform=ax.transAxes,
            ha='left', va='top', bbox=dict(ec='w', fc='w'),
            fontsize='small')

plt.show()
```

## Interpreting the PCA##

Reconstruction of spectrum, ${x}(k)$, from the
eigenvectors, ${e}_i(k)$

>$ \begin{equation}
  {x}_i(k) = {\mu}(k) + \sum_j^R \theta_{ij} {e}_j(k),
\end{equation}
$

Truncating this expansion (i.e., $r<R$)

>$\begin{equation}
{x}_i(k) = {\mu}(k) + \sum_i^{r<R} \theta_i {e}_i(k),
\end{equation}
$

- eigenvectors ordered by their associated eigenvalues
- eigenvalues reflect variance  within each eigenvector (sum of the eigenvalues is total variance of the system).
- project a each spectrum onto these first few eigenspectra is a compression of the data

This is the sense in which PCA gives for dimensionality reduction.

```python
from sklearn.decomposition import PCA, KernelPCA

from astroML.datasets import sdss_corrected_spectra
from astroML.decorators import pickle_results

#------------------------------------------------------------
# Download data
data = sdss_corrected_spectra.fetch_sdss_corrected_spectra()
spectra = sdss_corrected_spectra.reconstruct_spectra(data)
wavelengths = sdss_corrected_spectra.compute_wavelengths(data)

#------------------------------------------------------------
# Compute PCA components

# Eigenvalues can be computed using PCA as in the commented code below:

#from sklearn.decomposition import PCA
#pca = PCA()
#pca.fit(spectra)
#evals = pca.explained_variance_ratio_
#evals_cs = evals.cumsum()

#  because the spectra have been reconstructed from masked values, this
#  is not exactly correct in this case: we'll use the values computed
#  in the file compute_sdss_pca.py
evals = data['evals'] ** 2
evals_cs = evals.cumsum()
evals_cs /= evals_cs[-1]
evecs = data['evecs']
spec_mean = spectra.mean(0)

#------------------------------------------------------------
# Find the coefficients of a particular spectrum
spec = spectra[1]
coeff = np.dot(evecs, spec - spec_mean)

#------------------------------------------------------------
# Plot the sequence of reconstructions
fig = plt.figure(figsize=(8, 8))
fig.subplots_adjust(hspace=0)

for i, n in enumerate([0, 4, 8, 20]):
    ax = fig.add_subplot(411 + i)
    ax.plot(wavelengths, spec, '-', c='gray')
    ax.plot(wavelengths, spec_mean + np.dot(coeff[:n], evecs[:n]), '-k')

    if i < 3:
        ax.xaxis.set_major_formatter(plt.NullFormatter())

    ax.set_ylim(-2, 21)
    ax.set_ylabel('flux')

    if n == 0:
        text = "mean"
    elif n == 1:
        text = "mean + 1 component\n"
        text += r"$(\sigma^2_{tot} = %.2f)$" % evals_cs[n - 1]
    else:
        text = "mean + %i components\n" % n
        text += r"$(\sigma^2_{tot} = %.2f)$" % evals_cs[n - 1]

    ax.text(0.01, 0.95, text, ha='left', va='top', transform=ax.transAxes)

fig.axes[-1].set_xlabel(r'${\rm wavelength\ (\AA)}$')
plt.show()
```

## Comparing dimensionality reduction (PCA vs LLE)

*  Comparisons between the eigenvectors derived from PCA and known spectral types of galaxies
have shown that these statistically orthogonal components correlate
strongly with specific physical properties
*  given the cumulative nature of the sum of
variances used in PCA, astrophysically interesting components within
the spectra (e.g., sharp spectral lines or transient features for
certain galaxy populations) may not be reflected in the largest PCA
components.
- sum of
linear components can efficiently reconstruct the features within the
data does not always hold.  An example of this is the variation in
broad emission lines (such as those from quasars). The
variation in line width is an inherently nonlinear process and can
require a large number of components to fully characterize:
for broad line quasars over 30
components are required to reproduce the underlying spectra compared
to the 10 required for quiescent and star-forming galaxies.

## PCA with missing data##

Observed spectrum, $x^o$ is the true spectrum, ${x}$ plus a wavelength-dependent weight, ${w}$.

Weight is zero where data are missing and $1/{\sigma}^2$ for rest

Minimizing the quadratic deviation between  ${x}^o$ truncated reconstruction, $\sum_i \theta_i {e}_i$
and solving for $\theta_i$ gives

>$\begin{equation}
	\sum_k \theta_i {w}(k) {e}_i(k) {e}_j(k) =
	\sum_k {w}(k) {x}^o(k) {e}_j(k),
\end{equation}
$

If $M_{ij} = \sum_k {w}(k) {e}_i(k) {e}_j(k)$ and $F_i = \sum_k {w}(k) {x}^o(k) {e}_i(k)$ then

>$\begin{equation}
	\theta_i = \sum_j M_{ij}^{-1} F_{j},
\end{equation}
$

- $F_j$ are  coefficients derived from  gappy data
- $M_{ij}^{-1}$ shows how correlated  eigenvectors are over the missing regions.

An estimate of the uncertainty on the
reconstruction coefficients is given by

>$\begin{equation}
%Cov(\theta_i,\theta_j) = \frac{1}{N}M_{ij}^{-1}
{\rm Cov}(\theta_i,\theta_j) = M_{ij}^{-1}.
\end{equation}
$

Accuracy of this reconstruction will depend on the distribution of
the gaps within the data vector.

```python
from astroML.datasets import fetch_sdss_corrected_spectra
from astroML.datasets import sdss_corrected_spectra
from matplotlib import ticker

#------------------------------------------------------------
# Get spectra and eigenvectors used to reconstruct them
data = fetch_sdss_corrected_spectra()
spec = sdss_corrected_spectra.reconstruct_spectra(data)
lam = sdss_corrected_spectra.compute_wavelengths(data)
evecs = data['evecs']
mu = data['mu']
norms = data['norms']
mask = data['mask']

#------------------------------------------------------------
# plot the results
i_plot = ((lam > 5750) & (lam < 6350))
lam = lam[i_plot]

specnums = [20, 8, 9]
subplots = [311, 312, 313]

fig = plt.figure(figsize=(12, 15))
fig.subplots_adjust(hspace=0)

for subplot, i in zip(subplots, specnums):
    ax = fig.add_subplot(subplot)

    # compute eigen-coefficients
    spec_i_centered = spec[i] / norms[i] - mu
    coeffs = np.dot(spec_i_centered, evecs.T)

    # blank out masked regions
    spec_i = spec[i]
    mask_i = mask[i]
    spec_i[mask_i] = np.nan

    # plot the raw masked spectrum
    ax.plot(lam, spec_i[i_plot], '-', color='k', lw=2,
            label='True spectrum')

    # plot two levels of reconstruction
    for nev in [10]:
        if nev == 0:
            label = 'mean'
        else:
            label = 'nev=%i' % nev
        spec_i_recons = norms[i] * (mu + np.dot(coeffs[:nev], evecs[:nev]))
        ax.plot(lam, spec_i_recons[i_plot], label=label)

    # plot shaded background in masked region
    ylim = ax.get_ylim()
    mask_shade = ylim[0] + mask[i][i_plot].astype(float) * ylim[1]
    plt.fill(np.concatenate([lam[:1], lam, lam[-1:]]),
             np.concatenate([[ylim[0]], mask_shade, [ylim[0]]]),
             lw=0, fc='k', alpha=0.2)

    ax.set_xlim(lam[0], lam[-1])
    ax.set_ylim(ylim)
    ax.yaxis.set_major_formatter(ticker.NullFormatter())

    if subplot == 311:
        ax.legend(loc=1, prop=dict(size=14))

    ax.set_xlabel('$\lambda\ (\AA)$')
    ax.set_ylabel('normalized flux')

plt.show()
```

 ## Independent Component Analysis (ICA)

Independent component analysis (ICA) solves the ``cocktail party problem'', where there are multiple microphones situated through out a room containing
$N$ people. Each microphone picks up a linear combination of the $N$
voices. The goal of ICA is to use the concept of statistical
independence to isolate (or unmix) the individual signals.

Each spectrum, $x_i(k)$, can be described by

>$x_1(k) = a_{11}s_1(k) + a_{12}s_2(k) + a_{13}s_3(k) + \cdots\\
x_2(k) = a_{21}s_1(k) + a_{22}s_2(k) + a_{23}s_3(k) + \cdots \\
x_3(k) = a_{31}s_1(k) + a_{32}s_2(k) + a_{33}s_3(k) + \cdots \\
$

where $s_i(k)$ are the individual  spectra and $a_{ij}$
the appropriate mixing amplitudes.

In matrix form,
>$ X = A S$,

 Extracting these signal spectra is equivalent to estimating the appropriate weight matrix,
$W$, such that
>$ S = W X$.

Underlying ICA is the observation that the
input signals, $s_i(k)$, should be **statistically
independent** - the joint probability distribution, $f(x,y)$, can be
fully described by  $f(x^p,y^q)=f(x^p)f(y^q)$

 The rationale for ICA is that the
sum of any two independent random variables will always be more
Gaussian than either of the individual random variables (i.e., from the
central limit theorem). If we identify an unmixing matrix, $W$, that
maximizes the non-Gaussianity of the distributions, then we would be
identifying the input signals.

```python
# Author: Jake VanderPlas
# License: BSD
#   The figure produced by this code is published in the textbook
#   "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
#   For more information, see http://astroML.github.com
#   To report a bug or issue, use the following forum:
#    https://groups.google.com/forum/#!forum/astroml-general
import numpy as np
from matplotlib import pyplot as plt

from sklearn.decomposition import FastICA

from astroML.datasets import sdss_corrected_spectra
from astroML.decorators import pickle_results

#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX.  This may
# result in an error if LaTeX is not installed on your system.  In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)

#------------------------------------------------------------
# Download data
data = sdss_corrected_spectra.fetch_sdss_corrected_spectra()
spectra = sdss_corrected_spectra.reconstruct_spectra(data)
wavelengths = sdss_corrected_spectra.compute_wavelengths(data)

#----------------------------------------------------------------------
# Compute ICA components

def compute_ICA(n_components=5):
    spec_mean = spectra.mean(0)

    # ICA treats sequential observations as related.  Because of this, we need
    # to fit with the transpose of the spectra
    ica = FastICA(n_components - 1, algorithm='deflation')
    ica.fit(spectra.T)
    ica_comp = np.vstack([spec_mean,
                          ica.transform(spectra.T).T])

    return ica_comp

n_components = 5
decompositions = compute_ICA(n_components)

#----------------------------------------------------------------------
# Plot the results
fig = plt.figure(figsize=(10, 8))
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05,
                    bottom=0.1, top=0.95, hspace=0.05)

titles = ['ICA components',]

for j in range(n_components):
    ax = fig.add_subplot(n_components, 1, j + 1)

    ax.yaxis.set_major_formatter(plt.NullFormatter())
    ax.xaxis.set_major_locator(plt.MultipleLocator(1000))
    if j < n_components - 1:
        ax.xaxis.set_major_formatter(plt.NullFormatter())
    else:
        ax.set_xlabel(r'wavelength ${\rm (\AA)}$')

    ax.plot(wavelengths, comp[j], '-k', lw=1)

    # plot zero line
    xlim = [3000, 7999]
    ax.plot(xlim, [0, 0], '-', c='gray', lw=1)
    ax.set_xlim(xlim)

    if j == 0:
        ax.set_title(titles[0])
        label = 'mean'
    else:
        label = 'component %i' % j

    ax.text(0.03, 0.94, label, transform=ax.transAxes, ha='left', va='top')

    for l in ax.get_xticklines() + ax.get_yticklines():
        l.set_markersize(2)

    # adjust y limits
    ylim = plt.ylim()
    dy = 0.05 * (ylim[1] - ylim[0])
    ax.set_ylim(ylim[0] - dy, ylim[1] + 4 * dy)

plt.show()
```

 ## Non-negative matrix factorization

Eigenvectors are defined relative to  mean data vector. Principal components that can be positive or negative but for many physical systems we know data are a linear sum of positive components (e.g. galaxy spectrum is a linear sum of stellar components

Nonnegative matrix factorization (NMF) applies positivity constraint

Assumes  any data matrix can be factored into two matrices, $W$ and $Y$, with
>$\begin{equation}
X=W Y,
\end{equation}
$

and both $W$ and $Y$ are nonnegative

$WY$ is an approximation of $X$. Minimizing the reconstruction error $|| (X - W Y)^2 ||$
nonnegative bases can be derived using a simple update rule (Lee and Seung),
>$
\begin{eqnarray}
W_{ki} &=& W_{ki}\frac{[XY^T]_{ki}}{[WYY^T]_{ki}},\\
Y_{in} &=& Y_{in}\frac{[W^TX]_{in}}{[W^TWY]_{in}},
\end{eqnarray}
$

Iterative process does not guarantee nonlocal minima (like $K$-means and EM).

Use random initialization and cross-validation

```python
# Author: Jake VanderPlas
# License: BSD
%matplotlib inline

import numpy as np
from matplotlib import pyplot as plt

from sklearn.decomposition import NMF
from sklearn.decomposition import RandomizedPCA

from astroML.datasets import sdss_corrected_spectra
from astroML.decorators import pickle_results

#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX.  This may
# result in an error if LaTeX is not installed on your system.  In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)

#------------------------------------------------------------
# Download data
data = sdss_corrected_spectra.fetch_sdss_corrected_spectra()
spectra = sdss_corrected_spectra.reconstruct_spectra(data)
wavelengths = sdss_corrected_spectra.compute_wavelengths(data)

#----------------------------------------------------------------------
# Compute PCA, and NMF components
#  we'll save the results so that they can be re-used
@pickle_results('spec_decompositions.pkl')
def compute_PCA_NMF(n_components=5):
    spec_mean = spectra.mean(0)

    # PCA: use randomized PCA for speed
    pca = RandomizedPCA(n_components - 1)
    pca.fit(spectra)
    pca_comp = np.vstack([spec_mean,
                          pca.components_])

    # NMF requires all elements of the input to be greater than zero
    spectra[spectra < 0] = 0
    nmf = NMF(n_components)
    nmf.fit(spectra)
    nmf_comp = nmf.components_

    return pca_comp, nmf_comp

n_components = 5
decompositions = compute_PCA_NMF(n_components)

#----------------------------------------------------------------------
# Plot the results
fig = plt.figure(figsize=(10, 8))
fig.subplots_adjust(left=0.05, right=0.95, wspace=0.05,
                    bottom=0.1, top=0.95, hspace=0.05)

titles = ['PCA components', 'NMF components']

for i, comp in enumerate(decompositions):
    for j in range(n_components):
        ax = fig.add_subplot(n_components, 3, 3 * j + 1 + i)

        ax.yaxis.set_major_formatter(plt.NullFormatter())
        ax.xaxis.set_major_locator(plt.MultipleLocator(1000))
        if j < n_components - 1:
            ax.xaxis.set_major_formatter(plt.NullFormatter())
        else:
            ax.set_xlabel(r'wavelength ${\rm (\AA)}$')

        ax.plot(wavelengths, comp[j], '-k', lw=1)

        # plot zero line
        xlim = [3000, 7999]
        ax.plot(xlim, [0, 0], '-', c='gray', lw=1)
        ax.set_xlim(xlim)

        if j == 0:
            ax.set_title(titles[i])

        if titles[i].startswith('PCA') or titles[i].startswith('ICA'):
            if j == 0:
                label = 'mean'
            else:
                label = 'component %i' % j
        else:
            label = 'component %i' % (j + 1)

        ax.text(0.03, 0.94, label, transform=ax.transAxes,
                ha='left', va='top')

        for l in ax.get_xticklines() + ax.get_yticklines():
            l.set_markersize(2)

        # adjust y limits
        ylim = plt.ylim()
        dy = 0.05 * (ylim[1] - ylim[0])
        ax.set_ylim(ylim[0] - dy, ylim[1] + 4 * dy)

plt.show()

```

## Manifold Learning

Real-world data sets often have very nonlinear features (e.g. QSOs with broad lines) which are hard to describe compactly using linear eigenvectors. Manifold learning techniques search for a representation of these data within a lower dimensional space

</img>

 **Local Linear Embedding:**

Step 1: define local geometry
- local neighborhoods determined from $k$ nearest neighbors.
- for each point calculate weights that reconstruct a point from its $k$ nearest
neighbors
>$
\begin{equation}
  \mathcal{E}_1(W) = \left|X - WX\right|^2
\end{equation}
$
- essentially this is finding the hyperplane that describes the local surface at each point within the data set.

Step 2: embed within a lower dimensional space
- set all $W_{ij}=0$ except when point $j$ is one of the $k$ nearest neighbors of point $i$.
- $W$ becomes very sparse for $k \ll N$.
- minimize
>$\begin{equation}
  \mathcal{E}_2(Y) = \left|Y - W Y\right|^2,
\end{equation}
$

with $W$ fixed to find an $N$ by $d$ matrix ($d$ is the new dimensionality)

Step 1 requires a nearest-neighbor search.

Step 2 requires an
eigenvalue decomposition of the matrix $C_W \equiv (I-W)^T(I-W)$,

LLE has been applied to data as diverse as galaxy spectra, stellar spectra, and photometric light curves.

```python
# Author: Jake VanderPlas
# License: BSD
import numpy as np
from matplotlib import pyplot as plt
from sklearn import manifold, neighbors

from astroML.datasets import sdss_corrected_spectra
from astroML.datasets import fetch_sdss_corrected_spectra
from astroML.plotting.tools import discretize_cmap
from astroML.decorators import pickle_results

#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX.  This may
# result in an error if LaTeX is not installed on your system.  In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)

#------------------------------------------------------------
# Set up color-map properties
clim = (1.5, 6.5)
cmap = discretize_cmap(plt.cm.jet, 5)
cdict = ['unknown', 'star', 'absorption galaxy',
         'galaxy', 'emission galaxy',
         'narrow-line QSO', 'broad-line QSO']
cticks = [2, 3, 4, 5, 6]
formatter = plt.FuncFormatter(lambda t, *args: cdict[int(np.round(t))])

#------------------------------------------------------------
# Fetch the data; PCA coefficients have been pre-computed
data = fetch_sdss_corrected_spectra()
coeffs_PCA = data['coeffs']
c_PCA = data['lineindex_cln']
spec = sdss_corrected_spectra.reconstruct_spectra(data)
color = data['lineindex_cln']

#------------------------------------------------------------
# Compute the LLE projection; save the results
@pickle_results("spec_LLE.pkl")
def compute_spec_LLE(n_neighbors=10, out_dim=3):
    # Compute the LLE projection
    LLE = manifold.LocallyLinearEmbedding(n_neighbors, out_dim,
                                          method='modified',
                                          eigen_solver='dense')
    Y_LLE = LLE.fit_transform(spec)
    print (" - finished LLE projection")

    # remove outliers for the plot
    BT = neighbors.BallTree(Y_LLE)
    dist, ind = BT.query(Y_LLE, n_neighbors)
    dist_to_n = dist[:, -1]
    dist_to_n -= dist_to_n.mean()
    std = np.std(dist_to_n)
    flag = (dist_to_n > 0.25 * std)
    print (" - removing %i outliers for plot" % flag.sum())

    return Y_LLE[~flag], color[~flag]

coeffs_LLE, c_LLE = compute_spec_LLE(10, 3)

#----------------------------------------------------------------------
# Plot the results:
for (c, coeffs, xlim) in zip([c_PCA, c_LLE],
                             [coeffs_PCA, coeffs_LLE],
                             [(-1.2, 1.0), (-0.01, 0.014)]):
    fig = plt.figure(figsize=(5, 3.75))
    fig.subplots_adjust(hspace=0.05, wspace=0.05)

    # axes for colorbar
    cax = plt.axes([0.525, 0.525, 0.02, 0.35])

    # Create scatter-plots
    scatter_kwargs = dict(s=4, lw=0, edgecolors='none', c=c, cmap=cmap)

    ax1 = plt.subplot(221)
    im1 = ax1.scatter(coeffs[:, 0], coeffs[:, 1], **scatter_kwargs)
    im1.set_clim(clim)
    ax1.set_ylabel('$c_2$')

    ax2 = plt.subplot(223)
    im2 = ax2.scatter(coeffs[:, 0], coeffs[:, 2], **scatter_kwargs)
    im2.set_clim(clim)
    ax2.set_xlabel('$c_1$')
    ax2.set_ylabel('$c_3$')

    ax3 = plt.subplot(224)
    im3 = ax3.scatter(coeffs[:, 1], coeffs[:, 2], **scatter_kwargs)
    im3.set_clim(clim)
    ax3.set_xlabel('$c_2$')

    fig.colorbar(im3, ax=ax3, cax=cax,
                 ticks=cticks,
                 format=formatter)

    ax1.xaxis.set_major_formatter(plt.NullFormatter())
    ax3.yaxis.set_major_formatter(plt.NullFormatter())

    ax1.set_xlim(xlim)
    ax2.set_xlim(xlim)

    for ax in (ax1, ax2, ax3):
        ax.xaxis.set_major_locator(plt.MaxNLocator(5))
        ax.yaxis.set_major_locator(plt.MaxNLocator(5))

plt.show()

```


text:
```
from sympy import init_session
init_session()
a , mu, k, t = symbols('a, mu, k, t')
```

    IPython console for SymPy 0.7.6 (Python 2.7.6-64-bit) (ground types: python)

    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()

    Documentation can be found at http://www.sympy.org

Given is a PDF of a Branching Erlang (or Coxian) Distribution.

```
pdf = a * mu * exp(-mu*t)+ (((1-a)*(mu**k))/(factorial(k-1))) * (t**(k-1))*exp(-mu*t)
pdf
```

```
cdf = integrate(pdf, t)
cdf
```

```
eq = diff(cdf,t) - pdf
simplify(eq)
```

```

```

max_stars_repo_name: juliusf/ipython-notebooks

max_issues_repo_name: juliusf/ipython-notebooks

max_forks_repo_name: juliusf/ipython-notebooks


text:
# Lógica proposicional

El módulo `sympy` permite hacer algunos cálculos en lógica proposicional

```python
from sympy import *
```

Vamos a decirle a sympy que renderice las salidas en LaTeX con mathjax

```python
init_printing()
```

## Satisfacibilidad

**Ejemplo** Si es cierto $a\to b$, ¿podemos inferir algo sobre $a\vee c\to b\vee c$? ¿Y sobre $a\wedge c\to b\wedge c$?

```python
a,b,c = symbols("a,b,c")
```

Vamos a usar que $\Gamma \models \alpha$ equivale a que $\Gamma\cup\{\neg \alpha\}$ es insatisfacible

```python
p = a >> b
q = (a|c)>>(b|c)
```

```python
p&~q
```

```python
satisfiable(p&~q)
```

    False

También podíamos ver si $p\to q$ es una tautología por el teorema de la deducción, o lo que es lo mismo, mostrar que siempre es verdadera

```python
simplify(p>>q)
```

Lo mismo con $a\wedge c\to b\wedge c$

```python
simplify(p>>((a&c)>>(b&c)))
```

**Ejemplo** Estudia si el conjunto $\{c\to (a\vee c), b\to(c\to a), d\wedge \neg (c\to a)\}$ es satisfacible

```python
a,b,c,d = symbols("a,b,c,d")
```

```python
satisfiable((c>>(a|b))&(b>>(c>>a))&d&~(c>>a))
```

    False

De ser satisfacible, nos daría un mundo donde lo es

```python
satisfiable(a|b)
```

    {b: True, a: True}

O bien todos los mundos donde lo es, si así se lo pedimos

```python
list(satisfiable(a|b, all_models=True))
```

    [{b: True, a: True}, {b: False, a: True}, {b: True, a: False}]

## Implicación semántica

Podemos detectar si $\Gamma\models a$ con la siguiente función

```python
def implica(gamma,a):
    """ Determina si Gamma implica semanticamente a"""
    p = True
    for x in gamma:
        p = p & x
    p=p&(~a)
    return not(satisfiable(p))
```

```python
implica([a,a>>b],b)
```

    True

```python
implica([a>>b,b>>c],a>>c)
```

    True

Volvamos al primer ejemplo

```python
implica([a>>b], (a|c)>>(b|c))
```

    True

Otro ejemplo

```python
implica([a>>b],a)
```

    False

## Tautologías

Una proposición $a$ es tautología si $\models a$

```python
implica([], a>>(b>>a))
```

    True

```python
implica([],~a>>(a>>b))
```

    True

```python
def tautologia(a):
    return not(satisfiable(~a))
```

```python
tautologia(~a>>(a>>b))
```

    True

```python
tautologia(a>>b)
```

    False

```python
(a>>b).args
```

```python
to_cnf(a>>(b&c))
```

```python
type(_)
```

    And

```python
to_cnf(a>>b)
```

```python
type(_)
```

    Or

## Davis-Putnam
### Forma clausular y operaciones con literales

Vamos a definir la forma clausular de un conjunto de proposiciones

```python
G= {c>>(a|b), b>>(c>>a), d&~(c>>a)}
```

```python
def forma_clausular(G):
    fc = set([])
    for g in G:
        q = to_cnf(g)
        if isinstance(q,And):
            fc.update(q.args)
        else:
            fc.add(q)
    return fc
```

```python
forma_clausular(G)
```

```python
def es_literal(p):
    if isinstance(p,Symbol):
        return True
    if isinstance(p,Not):
        ar = p.args
        if len(ar)==1 and isinstance(ar[0],Symbol):
            return True
    return False
```

```python
es_literal(a)
```

    True

```python
es_literal(~a)
```

    True

```python
es_literal(a&b)
```

    False

```python
v=Or(*[])
```

```python
type(v)
```

    sympy.logic.boolalg.BooleanFalse

```python
v==False
```

    True

```python
def es_clausula(p):
    if p==False:
        return True
    if es_literal(p):
        return True
    if isinstance(p,Or):
        ar = p.args
        return all(es_literal(a) for a in ar)
    return False
```

```python
es_clausula(a&b)
```

    False

```python
es_clausula(a|~b)
```

    True

```python
def complemento(p):
    if isinstance(p,Symbol):
        return ~p
    if isinstance(p,Not):
        return p.args[0]
    return None
```

```python
complemento(a)
```

```python
complemento(~a)
```

```python
l=[1,2,3]
l.remove(3)
2 in l
```

    True

```python
def quita(c,l):
    if c==False:
        return False
    if es_literal(c):
        if c==l:
            return False
        return c
    ar = list(c.args)
    if l in ar:
        ar.remove(l)
        return Or(*ar)
    return c
```

```python
quita(a|b,b)
```

```python
quita(a|b,~b)
```

```python
def primero(xs,cond):
    for x in xs:
        if cond(x):
            return x
    return None
```

```python
primero([1,2,3],lambda x:(x%2)==0)
```

```python
def clausula_a_lista(c):
    if c==False:
        return []
    if es_literal(c):
        return [c]
    return list(c.args)
```

```python
clausula_a_lista(a|b|~c)
```

### Detectando una clausula unit

```python
def clausula_unit(cs):
    l = primero(cs,es_literal)
    if l!=None:
        print("Hemos encontrado una clausula unit ", l)
        csp = [c for c in cs if not(l in clausula_a_lista(c))]
        cspq = [quita(c,complemento(l)) for c in csp]
        pprint(cspq)
        return cspq
    return cs
```

```python
clausula_unit([a,a|b, ~a|c, c|b])
```

```python
Or(a,~a,b)
```

```python
simplify(_)
```

### Quitando tautologías y cláusulas redundantes

```python
def quita_tautologias(cs):
    return [c for c in cs if simplify(c)!=True]
```

```python
quita_tautologias([a|b|~a,a|c])
```

```python
l=[1,2]
```

La siguiente función determina si todos los literales de una cláusula están contenidos en otra cláusula

```python
def contenido(c1,c2):
    return all((l in clausula_a_lista(c2)) for l in clausula_a_lista(c1))
```

```python
contenido(a|b,a|b|c)
```

    True

```python
contenido(a|b,a|c)
```

    False

```python
def quita_redundantes(cs):
    return [c for c in cs if not(any(contenido(d,c) and d!=c for d in cs))]
```

```python
quita_redundantes([a|b,a|b|c])
```

### Encontrando litereales puros

```python
def literal_puro(cs):
    literales = set([])
    for c in cs:
        literales.update(clausula_a_lista(c))
    l = primero(literales, lambda x:not(complemento(x) in literales))
    if l==None:
        return cs
    print("Hemos encontrado un literal puro", l)
    csp = [c for c in cs if not(l in clausula_a_lista(c))]
    print(csp)
    return csp
```

```python
literal_puro([a|b,~a])
```

### Dividiendo cuando no hay literales puros ni cláuslas unit

```python
def divide(cs):
    literales = set([])
    for c in cs:
        literales.update(clausula_a_lista(c))
    l = primero(literales, lambda x:(complemento(x) in literales))
    cl = complemento(l)
    csl = set([quita(c,l) for c in cs if l in clausula_a_lista(c)])
    csnl= set([quita(c,cl) for c in cs if cl in clausula_a_lista(c)])
    csc = set([c for c in cs if not(l in clausula_a_lista(c)) and not(cl in clausula_a_lista(c))])
    parte1=csl.union(csc)
    parte2=csnl.union(csc)
    print("Dividimos usando ",l)
    pprint(parte1)
    pprint(parte2)
    return [csl.union(csc),csnl.union(csc)]
```

```python
divide([a|b,~a|b,a|c, b|c])
```

### Algoritmo de Davis-Putnam

```python
def inconsistente(cs):
    """Detecta si cs es satisfacible usando Davis-Putnam"""
    css = quita_tautologias(quita_redundantes(cs))
    unitpuro=True
    while unitpuro:
        cssn = clausula_unit(css)
        if cssn== css:
            cssn = literal_puro(css)
            if cssn==css:
                unitpuro=False
        if len(cssn)==0:
            print("Hemos llegado al conjunto vacío")
            return False
        if False in cssn:
            print("Hemos encontrado la cláusula vacía")
            return True
        css=cssn

    d = divide(css)
    return inconsistente(d[0]) and inconsistente(d[1])

```

```python
inconsistente([~a,a])
```

    Hemos encontrado una clausula unit  Not(a)
    [False]
    Hemos encontrado la cláusula vacía

    True

```python
clausula_unit([a,b])
```

```python
clausula_unit([b])
```

```python
inconsistente([a,b])
```

    Hemos encontrado una clausula unit  a
    [b]
    Hemos encontrado una clausula unit  b
    []
    Hemos llegado al conjunto vacío

    False

```python
inconsistente([~a|b,a,~b])
```

    Hemos encontrado una clausula unit  a
    [b, ¬b]
    Hemos encontrado una clausula unit  b
    [False]
    Hemos encontrado la cláusula vacía

    True

Este último ejemplo es precisamente el modus ponens

```python
inconsistente(forma_clausular([a>>b,a,~b]))
```

    Hemos encontrado una clausula unit  a
    [b, ¬b]
    Hemos encontrado una clausula unit  b
    [False]
    Hemos encontrado la cláusula vacía

    True

**Ejemplo** Veamos que $(a\to (b\to c))\to(\neg (a\to \neg b)\to c)$ es una tautología. Esto equivale a probar
$$\models (a\to (b\to c))\to(\neg (a\to \neg b)\to c).$$
Por el Teorema de la deducción (dos veces), basta probar que
$$\{a\to (b\to c), \neg (a\to \neg b)\}\models c,$$
y esto es equivalente a demostrar que el conjunto
$$\{a\to (b\to c), \neg (a\to \neg b),\neg c\}$$
es insatisfacible

```python
forma_clausular([a>>(b>>c), ~(a>>~b), ~c])
```

```python
inconsistente(_)
```

    Hemos encontrado una clausula unit  Not(c)
    [¬a ∨ ¬b, b, a]
    Hemos encontrado una clausula unit  b
    [¬a, a]
    Hemos encontrado una clausula unit  Not(a)
    [False]
    Hemos encontrado la cláusula vacía

    True

**Ejemplo** Veamos que $$\{(a\to \neg b\vee d)\wedge(b\wedge \neg d\to a\vee c), (d\to (a \leftrightarrow \neg b)\vee (b\wedge \neg c)\}\models (\neg b\to (d\wedge (c\vee \neg d)))\to c\wedge d$$

```python
Gamma={(a>>(~b|d))&((b&~d)>>(a|c)), (d>>(Equivalent(a,~b))|(b|~c)),
       (~b>>(d&(c|~d))),~(c&d)}
```

```python
forma_clausular(Gamma)
```

```python
inconsistente(forma_clausular(Gamma))
```

    Dividimos usando  Not(c)
    set([¬d, b ∨ d, d ∨ ¬a ∨ ¬b])
    set([b ∨ d, b ∨ ¬d, a ∨ d ∨ ¬b, d ∨ ¬a ∨ ¬b])
    Hemos encontrado una clausula unit  Not(d)
    [b, ¬a ∨ ¬b]
    Hemos encontrado una clausula unit  b
    [¬a]
    Hemos encontrado una clausula unit  Not(a)
    []
    Hemos llegado al conjunto vacío

    False

## Resolventes

```python
Or(Or(a,b),Or(a,c))
```

```python
def resolventes(a,b):
    literales=[l for l in clausula_a_lista(a) if complemento(l) in clausula_a_lista(b)]
    return [Or(quita(a,l),quita(b,complemento(l))) for l in literales]
```

```python
resolventes(~a|b,a|~b)
```

max_stars_repo_name: lmd-ugr/LMD

max_issues_repo_name: pedritomelenas/LMD

max_forks_repo_name: pedritomelenas/LMD


text:
# Holonomic Robot

**need to fix formatting**

Robots come in a variety of types and configurations: wheeled, tracked, legs, flying, etc. Common wheeled robots typically have two wheels (directly driven) with a caster wheel to make the robot stable. There are some without the caster wheel and employ a control system to keep them upright (inverted pendulum
problem) and resemble a Segway scooter. All of these two wheeled robot are
non-holonomic systems.

A **nonholonomic system** in physics and mathematics is a system whose state depends on the path taken in order to achieve it. A car is a typical example of a nonholonomic system. It can occupy any location (x,y) in any orientation ($\phi$), but the path from one location/orientation to another position/orientation is not always linear. Thus you have to parallel park.

A **holonomic system** is not constrained by this. Thus it can move anywhere at will.

```python
%matplotlib inline
```

```python
from __future__ import division
from __future__ import print_function
import numpy as np
from numpy.linalg import norm
from math import cos, sin, pi
import matplotlib.pyplot as plt
```

## Soccer

These types of robots are great for soccer, especially the goalie. Typical configuration are 3 or 4 motors. Omni wheels allow slip perpendicular to the axis.

Due to these constraints, a holonomic robot which could travel in any direction and immediately change its position and orientation is much more desirable. There are a variety of different wheels which make this type of robot possible such as mecanum or omni wheels.

Omni wheels operate like standard wheels in that the force is produced normal to the motor's axis of rotation and as a result of friction. However, there are a series of smaller wheels which ring the main wheel and allow the wheel to slip in the direction of the motor rotational axis. Note that no force is produced parallel to the motor axis, just slippage.

## Dynamics

The dynamics for a holonomic robot with 4 omni directional wheels (can be derived using Euler-Largrange ($\mathcal{L}$) which defines a system's kinectic ($T$) and potential ($V$) energies in relation to a set of generalized coordinates ($q$) and generalized forces ($Q$):

\begin{equation}
    \newcommand{\dpar}[2]{\frac{\partial #1}{\partial #2}}
\end{equation}

\begin{equation}
	\mathcal{L}=T-V \\
	\frac{d}{dt} \left\{ \dpar{ \mathcal{L} }{\dot q} \right\} - \dpar{ \mathcal{L} }{q} = Q \\
	T = \frac{1}{2}M v_w^2+ \frac{1}{2}J \dot \psi^2 + \frac{1}{2} J_w (\dot \theta_1^2 + \dot \theta_2^2 + \dot \theta_3^2 + \dot \theta_4^2) \\
	V = 0
\end{equation}

However, the dynamics must be calculated from an inertial reference frame (${W}$) and take into account the rotating body frame dynamics (${B'}$). Now, assume the body frame is offset from the center of mass (CM) by $x_m$ and $y_m$ which compose a vector $r_m$. Thus the velocity of the robot in the rotating frame would be:

\begin{equation*}
	v_w = v_{B'} + \dot \psi \times r_m \\
	v_w = v_{B'} +
	\begin{bmatrix}
		0 & 0 & \dot \psi
	\end{bmatrix}^T
	\times
	\begin{bmatrix}
		x_m & y_m & 0
	\end{bmatrix}^T
    =
	\begin{bmatrix}
		\dot x & \dot y & 0
	\end{bmatrix}^T +
	\begin{bmatrix}
		-y_m \dot \psi & x_m \dot \psi & 0
	\end{bmatrix}^T \\
	v_{B'} = \begin{bmatrix}
		\dot x & \dot y & 0
	\end{bmatrix}^T
\end{equation*}

where $v_{B'}$ is the speed of the body frame. Now substituting that into
the above kinetic energy equation $T$, we get:

\begin{equation}
	T = \frac{1}{2}M( ( \dot x - \dot \psi y )^2 + (\dot y + \dot \psi x)^2)+ \dots \\
	T = \frac{1}{2}M( \dot x^2 - 2 \dot \psi y_m \dot x +\dot \psi^2 y_m^2 + \dot y^2 + 2 \dot \psi x_m \dot y + \dot \psi^2 x_m^2)+ \frac{1}{2}J \dot \psi^2 + \frac{1}{2} J_w (\dot \theta_1^2 + \dot \theta_2^2 + \dot \theta_3^2 + \dot \theta_4^2)  \\
	\frac{d}{dt} \left\{ \dpar{ \mathcal{L} }{\dot x} \right\} = M ( \ddot x - \ddot \psi y - \dot \psi \dot y ) \hspace{1cm} \dpar{ \mathcal{L} }{x} = M(\dot \psi \dot y + \dot \psi^2 x) \\
	\frac{d}{dt} \left\{ \dpar{ \mathcal{L} }{\dot y} \right\} = M (\ddot y + \ddot \psi x + \dot \psi \dot x) \hspace{1cm} \dpar{ \mathcal{L} }{y} = M( -\dot \psi \dot x + \dot \psi^2 y) \\
	\frac{d}{dt} \left\{ \dpar{ \mathcal{L} }{\dot \psi} \right\} = J \ddot \psi \hspace{1cm} \dpar{ \mathcal{L} }{\phi} = 0 \\
	\frac{d}{dt} \left\{ \dpar{ \mathcal{L} }{\dot \theta} \right\} = J_w \sum \limits_{i=1}^4 \ddot \theta_i \hspace{1cm} \dpar{ \mathcal{L} }{\theta} = 0
\end{equation}

Now we make the following assumptions: ${B'}$ is coincident with
${B}$, $x_m = 0$, $y_m = 0$, $\dot x = v_x$,
$\dot y = v_y$

\begin{equation*}
	F_x = M (\ddot x - 2 \dot \psi \dot y ) \\
	F_y = M (\ddot y + 2 \dot \psi \dot x) \\
	T = J \ddot \psi \\
	\tau_w = J_w \ddot \theta_1 \hspace{1cm}
	\tau_w = J_w \ddot \theta_2 \hspace{1cm}
	\tau_w = J_w \ddot \theta_3 \hspace{1cm}
	\tau_w = J_w \ddot \theta_4
\end{equation*}

\begin{equation*}
	\begin{bmatrix}
		F_x \\
		F_y \\
		T
	\end{bmatrix} =
	\begin{bmatrix}
		M & 0 & 0 \\
		0 & M & 0 \\
		0 & 0 & J
	\end{bmatrix}
	\begin{bmatrix}
		\ddot x \\
		\ddot y \\
		\ddot \psi
	\end{bmatrix} +
	\begin{bmatrix}
		0 & -2M \dot \psi & 0 \\
		2M \dot \psi & 0 & 0 \\
		0 & 0 & 0
	\end{bmatrix}
	\begin{bmatrix}
		\dot x \\
		\dot y \\
		\dot \psi
	\end{bmatrix}
	= \mathcal{M} \ddot X + \mathcal{O} \dot X = Q
\end{equation*}

## World Coordinates

Now the dynamics derived so far are all in the body frame and we could stop
here and develop a controller which performs velocity control. However, position
control is more useful and a transform needs to be performed to move the
velocities and accelerations into the world frame.

\begin{equation}
	\dot X^W = R_B^W \dot X^B \\
	R_B^W =
	\begin{bmatrix}
		\cos \psi & \sin \psi & 0 \\
		-\sin \psi & \cos \psi & 0 \\
		0 & 0 & 1
	\end{bmatrix} \\
	\ddot X^W = \dot R_B^W \dot X^B + R_B^W \ddot X^B \\
	\dot R_B^W =
	\begin{bmatrix}
		\sin \psi & -\cos \psi & 0 \\
		\cos \psi & \sin \psi & 0 \\
		0 & 0 & 1
	\end{bmatrix}
\end{equation}

Now, substituting this into the dynamics, gives dynamics in the world
coordinate system of:

\begin{equation}
	F = \mathcal{M} (\dot R \dot X + R \ddot X )  + \mathcal{O} R \dot X \\
	F = \mathcal{M} R \ddot X + (\mathcal{M}  \dot R + \mathcal{O} R) \dot X
\end{equation}

## External Forces and Torques

Now summing the forces into their body referenced $x$ and $y$ directions and the torque about the $z$ axis, gives us:

\begin{equation}
    \sum F_x=f_1 \sin(\phi) - f_2 \sin(\phi) - f_3 \sin(\phi) + f_4 \sin(\phi) \\
	\sum F_y=f_1 \cos(\phi) + f_2 \cos(\phi) - f_3 \cos(\phi) - f_4 \cos(\phi) \\
	\sum T=L(f_1+f_2+f_3+f_4)
\end{equation}

Additionally, we can simplify this by assuming all of the angles are the same
(e.g., $\phi_1 = \phi_2 = \phi_3 = \phi_4$) and can now put this into a
matrix form:

\begin{equation}
	\begin{bmatrix}
		F_x \\
		F_y \\
		T
	\end{bmatrix} =
	\begin{bmatrix}
		\sin(\phi) & 0 & 0 \\
		0 & \cos(\phi) & 0 \\
		0 & 0 & L
	\end{bmatrix}
	\begin{bmatrix}
		1 & -1 & -1 & 1\\
		1 & 1 & -1 & -1\\
		1 & 1 & 1& 1
	\end{bmatrix}
	\begin{bmatrix}
		f_1 \\
		f_2 \\
		f_3 \\
		f_4
	\end{bmatrix}
\end{equation}

where $\phi$ is again the angle of the motors, $f_i$ is the magnitude of the force produced by the motors, and $L$ is the radius of the robot.

where $pinv()$ is defined as the pseudoinverse since $A(\phi)$ is not a square matrix. Finally, substituting these into the original equation, we can calculate the torques given the desired accelerations.

\begin{equation}
	\begin{bmatrix} \tau_1 \\  \tau_2 \\  \tau_3 \\  \tau_4 \end{bmatrix} = \frac {M r_w} {4}
	\begin{bmatrix}
		-1 & 1 & 1 \\
		-1 & -1 & 1 \\
		1 & -1 & 1 \\
		1 & 1 & 1
	\end{bmatrix}
	\begin{bmatrix}
		\frac{1}{\sin(\phi)} & 0 & 0 \\
		0 & \frac{1}{\cos(\phi)} & 0 \\
		0 & 0 & \frac{1}{2}
	\end{bmatrix}
	\begin{bmatrix}
		a_x \\
		a_y \\
		R \dot \omega
	\end{bmatrix}
\end{equation}

Now looking at this equation, we notice that $\phi$ can not be equal to 0, 90, 180, 270, or 360 otherwise we get a singularity in the $A(\phi)$ matrix. This however is not an issue in the real world, since the motors would occupy the same physical space and the robot would essentially only have 2 and not 4 motors.

## Holonomic Robot Kinematics

Number of equivalent motors for any direction under linear movement only, no rotational movement allowed.

Now performing a similar exercise for what was done with the dynamics, looking
at $coordinate$, the velocity of motor 1is given by
$v_1 = -\sin(\phi) v_x + \cos(\phi) v_y + R \omega$. Performing this for
each wheel gives:

\begin{equation}
	\begin{bmatrix}
		v_1 \\
		v_2 \\
		v_3 \\
		v_4
	\end{bmatrix} =
	\begin{bmatrix}
		-\sin(\phi)  & \cos(\phi) & L \\
		-\sin(\phi) & -\cos(\phi) & L \\
		 \sin(\phi) & -\cos(\phi) & L \\
		 \sin(\phi)  & \cos(\phi) & L
	\end{bmatrix}
	\begin{bmatrix}
		v_x \\
		v_y \\
		\omega
	\end{bmatrix} =
	\begin{bmatrix}
		-1 & 1 & 1 \\
		-1 & -1 & 1 \\
		1 & -1 & 1 \\
		1 & 1 & 1
	\end{bmatrix}
	\begin{bmatrix}
		\sin(\phi) & 0 & 0 \\
		0 & \cos(\phi) & 0 \\
		0 & 0 & L
	\end{bmatrix}
	\begin{bmatrix}
		v_x \\
		v_y \\
		\omega
	\end{bmatrix}
\end{equation}

Now setting $\omega$ to zero and calculating only linear movement, we can
determine the number of equivalent motors as shown in
\figref{fig:equivalent_motors}. For example, setting
$\phi$ to 30 $^\circ$ (the red line in \figref{fig:equivalent_motors})
and traveling in the x direction only ($\begin{bmatrix} v_x & v_y & \omega \end{bmatrix}^T = \begin{bmatrix}1& 0 & 0 \end{bmatrix}^T$),
the above equation simplifies to $4 \sin(30)$ or 2 equivalent motors.
Repeating for the y direction results in $4 \cos(30)$ or 3.46 equivalent
motors.

Now it is interesting to note that when $\phi$ is set to 30 $^\circ$,
the robot has more equivalent motors when going forward or backwards, while a
$\phi$ of 60 $^\circ$ provides more equivalent motors moving left or right.
When the motors are are angled at 45 $^\circ$, movement is clearly equally
optimized for both forward/backwards and left/right ( $2 \sin(45)$ is 2.83 motors)
movement.

\figref{fig:equivalent_motors} tells us that no mater how the 4 motors are
oriented in a realistic configuration, the robot will never have the equivalent
use of all 4 motors. Movement in one direction or another can be optimized, but
then a sacrifice is made in another direction. This fact is intuitively obvious.

Another issue is these results are also ideal. This logic assumes that the wheels
will not slip and have good traction in any orientation. Unfortunately real world
results do not mimic this situation and the robot's performance will be reduced.

## Equvilent Motors

Since this robot has 4 motors with omni wheels, in certain configurations, you can get more power/speed in certain directions than others. Let's look at a configuration with the wheels oriented at the $\phi$ angles of: 30, 45, and 60 degrees.

```python
def motors(phi, angles):
    """
    in:
        phi - orientation of motors
        angles - array of angles from 0 ... 2*pi
    out:
        array of results
    """
    phi = phi*pi/180
    ans = []
    for angle in angles:
        a = np.array([
                [-sin(phi),  cos(phi), 1.0],
                [-sin(phi), -cos(phi), 1.0],
                [ sin(phi), -cos(phi), 1.0],
                [ sin(phi),  cos(phi), 1.0]
            ])
        b = np.array([
                cos(angle),
                sin(angle),
                0.0
            ])
        v = sum(abs(a.dot(b)))
        ans.append(v)
    return ans
```

```python
theta = np.arange(0,2*pi,0.1)
em30 = motors(30, theta)
em45 = motors(45, theta)
em60 = motors(60, theta)
```

```python
plt.polar(theta, em30, label='30');
plt.polar(theta, em45, label='45');
plt.polar(theta, em60, label='60');
plt.grid(True)
plt.legend(loc='upper left');
```

As seen above, 30 and 60 degrees orientations favor certain directions. However, 45 degrees is better in all directions and all around movement.

## Control

Using the equations above to transform desired movement into motor actions, the robot can move in any direction. A basic summery of robot movement to motor direction is shown above.

# References

---

* Alexander Gloye, Raul Rojas, Holonomic Control of a Robot with an Omnidirectional Drive, accepted for publication by Künstliche Intelligenz, Springer-Verlag, 2006.

* http://en.wikipedia.org/wiki/Non-holonomic_system

* http://en.wikipedia.org/wiki/Lagrangian_mechanics

* http://www.kornylak.com

* R. Balakrishna, Ashitava Ghosal, "Modeling of Slip for Wheeled Mobile Robots," lEEE TRANSACTIONS ON ROBOTICS AND AUTOMATION, VOL. I I , NO. I , FEBRUARY 1995, pp. 126-132

* J. Agullo, S. Cardona, and J. Vivancos, “Kinematics of vehicles with directional sliding wheels,” Mechanisms and Muchine Theory, vol. 22, no. 4, pp. 295-301, 1987.

* Pseudoinverse: for m > n: $A_{left}^{-1}=(A^TA)^{-1}A^T$ or m < n: $A_{right}^{-1} = A^T(AA^T)^{-1}$ such that $AA^{-1}=I$ or $A^{-1}A=I$

* Masayoshi Wada (2010). Motion Control of a Four-wheel-drive Omnidirectional Wheelchair with High Step Climbing Capability, Climbing and Walking Robots, Behnam Miripour (Ed.), ISBN: 978-953-307-030-8, [InTech](http://www.intechopen.com/books/climbing-and-walking-robots/motion-control-of-a-four-wheeldrive-omnidirectional-wheelchair-with-high-step-climbing-capability)

-----------

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.


text:
# Meta

Filename:  BasicNeuralNets.ipynb

Purpose:  To get an introductory understanding of artificial neural networks in the context of image-recognition.  To illustrate their structure, function, and effectiveness.  In practice, people generally start studying this area using the MNIST data set to construct a neural net to distinguish and correctly label handwritten digits.  This is the traditional "hello world" project in this field, so let's do that.  Let's construct basic digit-recognizing networks (in the form of FFSNNs) using gradient descent and backpropagation to improve performance from a random starting point (that is, from randomly selected initial parameters).

Abbreviation: MNIST = Modified (handwritten digits data set from the U.S.) National Institute of Standards and Technology

Abbreviation: FFSNN = Feed-Forward Sigmoid-Neuron Network

Overview: Watch this video series to get an excellent (and beautiful) overview of how neural networks (such as the one ones built here) are built and function: Neural Networks, by 3Blue1Brown ( https://www.youtube.com/watch?v=aircAruvnKk&list=PLZHQObOWTQDNU6R1_67000Dx_ZCJB-3pi )

Reference: Much of the strategy and code here is taken from Michael Nielsen's online book on "Neural Networks and Deep Learning" ( http://neuralnetworksanddeeplearning.com ).

# Network IO

The input data set originates from what's called the MNIST database of handwritten digits (found here: http://yann.lecun.com/exdb/mnist/).  This is a set of 70,000 images of hand-written digits, scanned and digitized, each image processed so the digit is size-normalized, centered, and anti-aliased, with 28x28 (=784) pixels, each pixel having a greyscale value (between 0 and 1? -- take a look at one original image example).

Originally split into a training set of 60,000 images and a test set of 10,000 images, we'll split it slightly differently into three groups: training images (for setting parameters), validation images (for setting hyperparameters), and test images (for testing the effectiveness of the resulting network).
\begin{align}
I   &=  \text{set of 70,000 input images}  \\
I1  &=  \text{set of } n1 \text{ (50,000) training images}  \\
I2  &=  \text{set of } n2 \text{ (10,000) validation images}  \\
I3  &=  \text{set of } n3 \text{ (10,000) test images}
\end{align}

# Network Specs

Each network can be visually represented as being composed of nodes and connections, which symbolize the algorithmic relationships between the input data and output data, via intermediate calculations that produce the output from the input.

The nodes (or artificial "neurons") are organized by layers, and each adjacent layer's nodes are all connected to each other.  In other words, adjacent layers are fully connected.
\begin{align}
L  &=  \text{number of layers}  \\
\mathbf{N}
   &=  \text{numbers of nodes/neurons in each layer (vector, } L\text{-dimensional, ordered by layer, from input to output)}  \\
   &=  \left< N_1, ..., N_i, ..., N_L \right>  \\
N_1
   &=  \text{number of nodes/neurons in input layer}  \\
N_i
   &=  \text{numbers of nodes/neurons in } i^{\text{th}} \text{ layer}  \\
N_L
   &=  \text{numbers of nodes/neurons in output layer}
\end{align}

Note: For any $i$ where $1<i<L$, the $i^\text{th}$ layer (or layer-$i$) is called "hidden", presumably because it is part of the eventual "black box" final-product network, where the input and output are the (more) visible parts.

Each node/neuron holds/outputs an "activation" (number).
\begin{align}
\mathbf{a}
   &=  \text{activations of nodes/neurons in each layer (vector, } L\text{-dimensional, ordered by layer, from input to output)}  \\
   &=  \left< \mathbf{a}^1, ..., \mathbf{a}^i, ..., \mathbf{a}^L \right>  \\
\mathbf{a}^i
   &=  \text{activations of nodes/neurons in } i^{\text{th}} \text{ layer (vector, } N_i\text{-dimensional, ordered by node, from top to bottom)}  \\
   &=  \left< a^i_1, ..., a^i_j, ..., a^i_{N_i} \right>
\end{align}

Each activation value in the input layer is a function of the input data (for the first node-layer) or a function of the previous layer's activation

```python

```

max_stars_repo_name: hotpocket/DigitRecognizer

max_issues_repo_name: hotpocket/DigitRecognizer

max_forks_repo_name: hotpocket/DigitRecognizer


text:
# Sensibilidad OIS

**Definición:** una curva cupón cero es una colección ordenada de plazos y tasas de interés donde cada una de las tasas es una  tasa adecuada para traer a valor presente un flujo de caja al plazo que corresponde a la tasa.

Consideremos un instrumento financiero cuyo valor depende de una curva cupón cero. Por ejemplo, un swap de tasa de interés. Matemáticamente, esta dependencia se expresa como:

$$
V = V\left(z_1,\ldots,z_n;\alpha\right)
$$

donde $z_1,\ldots,z_N$ son los valores de las tasas de la curva y $\alpha$ representa un vector de parámetros adicionales (como la tasa de cupón en un bono o un tipo de cambio en un swap de monedas).

**Definición:** la cantidad $\Delta_i$ de $V$ respecto a $z_i$ se define como:

$$
\Delta_i=\frac{\partial V}{\partial z_i}
$$

**Definición:** la sensibilidad de $V$ respecto a un cambio $\delta_i$ (positivo o negativo) en la tasa $z_i$ se define como:

$$
S\left(V,z_i,\delta_i\right)=V\left(z_1,...,z_i+\delta_i,\ldots,z_n\right)-V\left(z_1,\ldots,z_i,\ldots,z_n\right)
$$

Cuando se busca conocer la sensibilidad a un cambio $\delta_i$  en la tasa $z_i$, donde $\delta_i$ representa una magnitud sin signo, se utiliza la siguiente definición:

**Definición:** la sensibilidad de $V$ respecto a un cambio de magnitud $\delta_i$ en la tasa $z_i$ se define como:

$$
\overline{S}\left(V,z_i,\delta_i\right)=\frac{V\left(z_1,\ldots ,z_i+\delta_i,\ldots,z_n\right)-V\left(z_1,\ldots,z_i-\delta_i,...,z_n\right)}{2}
$$

En ocasiones, cuando se busca conocer la sensibilidad a un movimiento de magnitud pequeña, resulta conveniente realizar la siguiente aproximación:

$$
\overline{S}\left(V,z_i,\delta_i\right)\approx\Delta_i \cdot\delta_i
$$

este es el caso, cuando se busca conocer esta sensibilidad para un elevado número de operaciones. Es importante considerar que se puede asegurar que esta aproximación es válida cuando la función $V$ satisface los requerimientos de regularidad necesarios para aplicar el [Teorema de Taylor](https://en.wikipedia.org/wiki/Taylor%27s_theorem).

f(x0 + h) - f(x0) = f'(x0) x h + error, error --> 0, para h --> 0

## Configuración

### Librerías

```python
from finrisk import QC_Financial_3 as Qcf
import modules.auxiliary as aux
from functools import partial
from enum import Enum
import pandas as pd
```

### Variables Globales

```python
class BusCal(Enum):
    NY = 1
    SCL = 2
    TARGET = 3
```

```python
def get_cal(code: BusCal) -> Qcf.BusinessCalendar:
    """
    """
    if code == BusCal.NY:
        cal = Qcf.BusinessCalendar(Qcf.QCDate(1, 1, 2020), 20)
        for agno in range(2020, 2071):
            f = Qcf.QCDate(12, 10, agno)
            if f.week_day() == Qcf.WeekDay.SAT:
                cal.add_holiday(Qcf.QCDate(14, 10, agno))
            elif f.week_day() == Qcf.WeekDay.SUN:
                cal.add_holiday(Qcf.QCDate(13, 10, agno))
            elif f.week_day() == Qcf.WeekDay.MON:
                cal.add_holiday(Qcf.QCDate(12, 10, agno))
            elif f.week_day() == Qcf.WeekDay.TUE:
                cal.add_holiday(Qcf.QCDate(11, 10, agno))
            elif f.week_day() == Qcf.WeekDay.WED:
                cal.add_holiday(Qcf.QCDate(10, 10, agno))
            elif f.week_day() == Qcf.WeekDay.THU:
                cal.add_holiday(Qcf.QCDate(9, 10, agno))
            else:
                cal.add_holiday(Qcf.QCDate(8, 10, agno))

        cal.add_holiday(Qcf.QCDate(15, 2, 2021))

    return cal
```

```python
get_cal(BusCal.NY)
```

    <finrisk.QC_Financial_3.BusinessCalendar at 0x7f3e87e9e538>

```python
frmt = {
    'tasa': '{:.6%}',
    'df': '{:.6%}',
    'valor_tasa': '{:.4%}',
    'spread': '{:.4%}',
    'nominal': '{:,.0f}',
    'interes': '{:,.0f}',
    'amortizacion': '{:,.0f}',
    'flujo': '{:,.4f}',
}
```

```python
class TypeOis(Enum):
    SOFR = 1
    ICP = 2
```

```python
type_ois_template = {
    TypeOis.SOFR: {
        'currency': Qcf.QCUSD(),
        'periodicity': Qcf.Tenor('1Y'),
        'stub_period': Qcf.StubPeriod.SHORTFRONT,
        'settlement_lag': 0,
        'calendar': BusCal.NY,
        'bus_adj_rule': Qcf.BusyAdjRules.MODFOLLOW,
        'amort_is_cashflow': True,
        'fixed_rate': Qcf.QCInterestRate(0.0, Qcf.QCAct360(), Qcf.QCLinearWf()),
    },

    TypeOis.ICP: {
        'currency': Qcf.QCCLP(),
        'periodicity': Qcf.Tenor('6M'),
        'stub_period': Qcf.StubPeriod.SHORTFRONT,
        'settlement_lag': 0,
        'calendar': BusCal.SCL,
        'bus_adj_rule': Qcf.BusyAdjRules.MODFOLLOW,
        'amort_is_cashflow': True,
        'fixed_rate': Qcf.QCInterestRate(0.0, Qcf.QCAct360(), Qcf.QCLinearWf()),
    }
}
```

## Importa Curva Cero Cupón

Se importa la data de la curva cupón cero que fue construida en el notebook 6.

```python
df_curva = pd.read_excel('data/20201012_built_sofr_zero.xlsx')
```

```python
df_curva.head(20).style.format(frmt)
```

        	        plazo	        tasa	        df

                        0
                        1
                        0.081111%
                        99.999778%

                        1
                        7
                        0.084051%
                        99.998388%

                        2
                        14
                        0.077967%
                        99.997010%

                        3
                        21
                        0.077358%
                        99.995549%

                        4
                        33
                        0.078067%
                        99.992942%

                        5
                        61
                        0.078064%
                        99.986954%

                        6
                        92
                        0.081103%
                        99.979560%

                        7
                        125
                        0.078059%
                        99.973271%

                        8
                        152
                        0.076030%
                        99.968343%

                        9
                        182
                        0.075014%
                        99.962603%

                        10
                        212
                        0.073593%
                        99.957265%

                        11
                        243
                        0.072476%
                        99.951761%

                        12
                        273
                        0.071967%
                        99.946187%

                        13
                        306
                        0.071356%
                        99.940196%

                        14
                        335
                        0.070848%
                        99.934996%

                        15
                        365
                        0.070238%
                        99.929787%

                        16
                        547
                        0.064161%
                        99.903893%

                        17
                        730
                        0.056758%
                        99.886549%

                        18
                        1097
                        0.075010%
                        99.774813%

                        19
                        1462
                        0.125752%
                        99.497569%

```python
def get_curve_from_dataframe(
        yf: Qcf.QCYearFraction,
        wf: Qcf.QCWealthFactor,
        df_curva: pd.DataFrame) -> Qcf.ZeroCouponCurve:
    """
    Retorna un objeto Qcf.ZeroCouponCurve. Esta función requiere que `df_curva` tenga una columna
    de nombre 'plazo' y una columna de nombre 'tasa'. Se usa interpolación lineal en la curva que
    se retorna.
    """
    plazos = Qcf.long_vec()
    tasas = Qcf.double_vec()
    for row in df_curva.itertuples():
        plazos.append(row.plazo)
        tasas.append(row.tasa)
    curva = Qcf.QCCurve(plazos, tasas)
    curva = Qcf.QCLinearInterpolator(curva)
    tipo_tasa = Qcf.QCInterestRate(0.0, yf, wf)
    curva = Qcf.ZeroCouponCurve(curva, tipo_tasa)
    return curva
```

```python
zcc = get_curve_from_dataframe(
    Qcf.QCAct365(),
    Qcf.QCCompoundWf(),
    df_curva
)
```

Algunos métodos del objeto`zcc`.

```python
plazo = 900
print(f"Tasa a {plazo} días es igual a {zcc.get_rate_at(plazo):.4%}")
print(f"Factor de descuento a {plazo} días es igual a {zcc.get_discount_factor_at(plazo):.6%}")
```

    Tasa a 900 días es igual a 0.0652%
    Factor de descuento a 900 días es igual a 99.839384%

## Valorización

```python
def get_ois_using_template(
    template,
    type_ois: TypeOis,
    rp: Qcf.RecPay,
    notional: float,
    start_date: Qcf.QCDate,
    tenor: Qcf.Tenor,
    fixed_rate_value: float,
    spread: float,
    gearing: float
):
    """
    """
    template_dict = template[type_ois]
    meses = tenor.get_years() * 12 + tenor.get_months()
    end_date = start_date.add_months(meses)
    template_dict['fixed_rate'].set_value(fixed_rate_value)
    es_bono = False

    # Construye la pata fija
    fixed_rate_leg = Qcf.LegFactory.build_bullet_fixed_rate_leg(
        rp,
        start_date,
        end_date,
        template_dict['bus_adj_rule'],
        template_dict['periodicity'],
        template_dict['stub_period'],
        get_cal(template_dict['calendar']),
        template_dict['settlement_lag'],
        notional,
        template_dict['amort_is_cashflow'],
        template_dict['fixed_rate'],
        template_dict['currency'],
        es_bono)

    # Construye la pata ois
    rp = Qcf.RecPay.PAY if rp == Qcf.RecPay.RECEIVE else Qcf.RecPay.RECEIVE
    icp_clp_leg = Qcf.LegFactory.build_bullet_icp_clp2_leg(
        rp,
        start_date,
        end_date,
        template_dict['bus_adj_rule'],
        template_dict['periodicity'],
        template_dict['stub_period'],
        get_cal(template_dict['calendar']),
        template_dict['settlement_lag'],
        notional,
        template_dict['amort_is_cashflow'],
        spread,
        gearing,
        True
    )

    for i in range(icp_clp_leg.size()):
        cshflw = icp_clp_leg.get_cashflow_at(i)
        cshflw.set_start_date_icp(1.0)
        cshflw.set_end_date_icp(1.0)

    return fixed_rate_leg, icp_clp_leg
```

### Operación Ejemplo

```python
op = get_ois_using_template(
    type_ois_template,
    TypeOis.SOFR,
    Qcf.RecPay.RECEIVE,
    10000000,
    Qcf.QCDate(14, 10, 2020),
    Qcf.Tenor('2Y'),
    .01,
    0.0,
    1.0
)
op
```

    (<finrisk.QC_Financial_3.Leg at 0x7f3e85682100>,
     <finrisk.QC_Financial_3.Leg at 0x7f3e85682e68>)

#### Digresión: `functools.partial`

Supongamos que estoy en una situación en la que sólo quiero construir OIS de *SOFR*. Me gustaría no tener que repetir los argumentos `type_ois_template` y `TypeOis.SOFR` cada vez que llamo la función `get_ois_using_template`. Puedo definir una nueva función de la siguiente forma:

```python
get_ois_sofr = partial(get_ois_using_template, type_ois_template, TypeOis.SOFR)
```

```python
y = f(x, y, z)
y12 = f(1, 2, z)
```

Con esta nueva función, `get_ois_sofr`, ahora puedo construir la operación `op` de la siguiente forma:

```python
op = get_ois_sofr(
    Qcf.RecPay.RECEIVE,
    10000000,
    Qcf.QCDate(14, 10, 2020),
    Qcf.Tenor('2Y'),
    .01,
    0.0,
    1.0
)
```

```python
op
```

    (<finrisk.QC_Financial_3.Leg at 0x7f3e847e53d8>,
     <finrisk.QC_Financial_3.Leg at 0x7f3e847e51d0>)

#### Continuamos (fin digresión ...)

```python
aux.show_leg(op[0], 'FixedRateCashflow', '').style.format(frmt)
```

        	        fecha_inicial	        fecha_final	        fecha_pago	        nominal	        amortizacion	        interes	        amort_es_flujo	        flujo	        moneda	        valor_tasa	        tipo_tasa

                        0
                        2020-10-14
                        2021-10-14
                        2021-10-14
                        10,000,000
                        0
                        101,389
                        True
                        101,388.8889
                        USD
                        1.0000%
                        LinAct360

                        1
                        2021-10-14
                        2022-10-14
                        2022-10-14
                        10,000,000
                        10,000,000
                        101,389
                        True
                        10,101,388.8889
                        USD
                        1.0000%
                        LinAct360

```python
aux.show_leg(op[1], 'IcpClpCashflow', '').style.format(frmt)
```

        	        fecha_inicial	        fecha_final	        fecha_pago	        nominal	        amortizacion	        amort_es_flujo	        flujo	        moneda	        icp_inicial	        icp_final	        valor_tasa	        interes	        spread	        gearing	        tipo_tasa

                        0
                        2020-10-14
                        2021-10-14
                        2021-10-14
                        -10,000,000
                        0
                        True
                        0.0000
                        CLP
                        1.000000
                        1.000000
                        0.0000%
                        -0
                        0.0000%
                        1.000000
                        LinAct360

                        1
                        2021-10-14
                        2022-10-14
                        2022-10-14
                        -10,000,000
                        -10,000,000
                        True
                        -10,000,000.0000
                        CLP
                        1.000000
                        1.000000
                        0.0000%
                        -0
                        0.0000%
                        1.000000
                        LinAct360

#### Valor Presente Pata Fija

```python
vp = Qcf.PresentValue()
```

```python
fecha_val = Qcf.QCDate(14, 10, 2020)
```

```python
vp_fija = vp.pv(fecha_val, op[0], zcc)
print(f'El valor presente de la pata fija es: USD {vp_fija:,.2f}')
```

    El valor presente de la pata fija es: USD 10,008,917.25

**Ejercicio:** Replique el valor de la pata fija en Excel utilizando los flujos y los factores de descuento que obtiene de la curva `zcc`.

```python
# Se dan de alta las fechas finales de ambos cupones
fecha1 = Qcf.QCDate(14, 10, 2021)
fecha2 = Qcf.QCDate(14, 10, 2022)

# O de forma más automática:
fecha1 = op[0].get_cashflow_at(0).get_settlement_date()
fecha2 = op[0].get_cashflow_at(1).get_settlement_date()

# Se calcula el número de días entre la fecha de valorización (fecha_val) y las fechas
# finales de ambos cupones.
plazo1 = fecha_val.day_diff(fecha1)
plazo2 = fecha_val.day_diff(fecha2)

# Utilizando la curva zcc se calculan los df a esos plazos
df1 = zcc.get_discount_factor_at(plazo1)
df2 = zcc.get_discount_factor_at(plazo2)

# Se obtienen los flujos totales (interés y amortización) de ambos cupones
flujo1 = op[0].get_cashflow_at(0).amount()
flujo2 = op[0].get_cashflow_at(1).amount()

# Finalmente, se calcula el valor presente como el producto (escalar) entre los df y los flujos.
check_vp = df1 * flujo1 + df2 * flujo2

# Se muestra el resultado.
print(f'El valor presente a mano es: {check_vp:,.2f}')
```

    El valor presente a mano es: 10,008,917.25

#### Valor Presente Pata Flotante

```python
fwd = Qcf.ForwardRates()
```

Este valor presente sólo considera el flujo del nocional al vencimiento. Esto porque aún no hemos calculado el valor esperado de los índices (o las tasas equivalentes).

```python
print(f'VP: {vp.pv(fecha_val, op[1], zcc):,.2f}')
```

    VP: -9,988,658.09

Comprobamos.

```python
print(f'{-df2 * 10000000:,.2f}')
```

    -9,988,658.09

Ahora se calculan los valores esperados.

```python
fwd.set_rates_icp_clp_leg(fecha_val, 1.0, op[1], zcc)
```

```python
aux.show_leg(op[1], 'IcpClpCashflow', '').style.format(frmt)
```

        	        fecha_inicial	        fecha_final	        fecha_pago	        nominal	        amortizacion	        amort_es_flujo	        flujo	        moneda	        icp_inicial	        icp_final	        valor_tasa	        interes	        spread	        gearing	        tipo_tasa

                        0
                        2020-10-14
                        2021-10-14
                        2021-10-14
                        -10,000,000
                        0
                        True
                        -7,023.7778
                        CLP
                        1.000000
                        1.000702
                        0.0700%
                        -7,097
                        0.0000%
                        1.000000
                        LinAct360

                        1
                        2021-10-14
                        2022-10-14
                        2022-10-14
                        -10,000,000
                        -10,000,000
                        True
                        -10,004,327.9717
                        CLP
                        1.000702
                        1.001135
                        0.0400%
                        -4,056
                        0.0000%
                        1.000000
                        LinAct360

```python
vp_flot = vp.pv(fecha_val, op[1], zcc)
print(f'El valor presente de la pata flotante es: USD {vp_flot:,.2f}')
```

    El valor presente de la pata flotante es: USD -10,000,000.00

¿Porqué sabíamos que tenía que dar 10,000,000.00? Hint: la respuesta está en la construcción de la curva cero OIS.

## Sensibilidad

En `QC_Financial_3` al calcular el valor presente, también se calculan las derivadas del valor presente respecto a cada uno de los vértices de la curva.

### Pata Fija

f(g(x)) = y --> dy/dx = f'(g(x)) . g'(x)

```python
vp.pv(fecha_val, op[0], zcc)
der = vp.get_derivatives()
```

Tengo que ejecutar el método `get_derivatives` justo después de valorizar.

```python
for d in der:
    print(d)
```

    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    -10124.661219898855
    0.0
    -19986227.206320565
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0
    0.0

Con esas derivadas, se puede calcular la sensibilidad a cada vértice de la curva cupón cero para un movimiento de 1 punto básico.

```python
print(type(der))
```

    <class 'finrisk.QC_Financial_3.double_vec'>

```python
delta = .0001
total = 0
for i, d in enumerate(der):
    total += d * delta
    print(f"Sensibilidad en {i}: {d * delta: 0,.2f}")
print(f"Sensibilidad total: {total:,.2f}")
```

    Sensibilidad en 0:  0.00
    Sensibilidad en 1:  0.00
    Sensibilidad en 2:  0.00
    Sensibilidad en 3:  0.00
    Sensibilidad en 4:  0.00
    Sensibilidad en 5:  0.00
    Sensibilidad en 6:  0.00
    Sensibilidad en 7:  0.00
    Sensibilidad en 8:  0.00
    Sensibilidad en 9:  0.00
    Sensibilidad en 10:  0.00
    Sensibilidad en 11:  0.00
    Sensibilidad en 12:  0.00
    Sensibilidad en 13:  0.00
    Sensibilidad en 14:  0.00
    Sensibilidad en 15: -1.01
    Sensibilidad en 16:  0.00
    Sensibilidad en 17: -1,998.62
    Sensibilidad en 18:  0.00
    Sensibilidad en 19:  0.00
    Sensibilidad en 20:  0.00
    Sensibilidad en 21:  0.00
    Sensibilidad en 22:  0.00
    Sensibilidad en 23:  0.00
    Sensibilidad en 24:  0.00
    Sensibilidad en 25:  0.00
    Sensibilidad en 26:  0.00
    Sensibilidad en 27:  0.00
    Sensibilidad en 28:  0.00
    Sensibilidad en 29:  0.00
    Sensibilidad en 30:  0.00
    Sensibilidad en 31:  0.00
    Sensibilidad en 32:  0.00
    Sensibilidad total: -1,999.64

**Ejercicio:** Verifique por diferencias finitas centrales la derivada y la sensibilidad en el vértice 17. La aproximación de la derivada por diferencias finitas centrales es:

$$
\begin{equation}
f^{\prime}\left(x\right)\approx\frac{f\left(x+h\right)-f\left(x-h\right)}{2h}
\end{equation}
$$

### Pata Flotante

En este caso, al rescatar estas derivadas, obtengo sólo lo relacionado con la curva de descuento.

```python
vp.pv(fecha_val, op[1], zcc)
der = vp.get_derivatives()
```

La sensibilidad es muy parecida a la de la pata fija.

```python
total = 0
for i, d in enumerate(der):
    total += d * delta
    print(f"Sensibilidad en {i}: {d * delta: 0,.2f}")
print(f"Sensibilidad total: {total:,.2f}")
```

    Sensibilidad en 0:  0.00
    Sensibilidad en 1:  0.00
    Sensibilidad en 2:  0.00
    Sensibilidad en 3:  0.00
    Sensibilidad en 4:  0.00
    Sensibilidad en 5:  0.00
    Sensibilidad en 6:  0.00
    Sensibilidad en 7:  0.00
    Sensibilidad en 8:  0.00
    Sensibilidad en 9:  0.00
    Sensibilidad en 10:  0.00
    Sensibilidad en 11:  0.00
    Sensibilidad en 12:  0.00
    Sensibilidad en 13:  0.00
    Sensibilidad en 14:  0.00
    Sensibilidad en 15:  0.70
    Sensibilidad en 16:  0.00
    Sensibilidad en 17:  1,997.46
    Sensibilidad en 18:  0.00
    Sensibilidad en 19:  0.00
    Sensibilidad en 20:  0.00
    Sensibilidad en 21:  0.00
    Sensibilidad en 22:  0.00
    Sensibilidad en 23:  0.00
    Sensibilidad en 24:  0.00
    Sensibilidad en 25:  0.00
    Sensibilidad en 26:  0.00
    Sensibilidad en 27:  0.00
    Sensibilidad en 28:  0.00
    Sensibilidad en 29:  0.00
    Sensibilidad en 30:  0.00
    Sensibilidad en 31:  0.00
    Sensibilidad en 32:  0.00
    Sensibilidad total: 1,998.16

La estructura es la misma que para una pata fija, lo que indica que se debe también incluir la sensibilidad a la curva de proyección.

```python
import numpy as np
result = []
for i in range(op[1].size()):
    cshflw = op[1].get_cashflow_at(i)

    # Se calcularon las derivadas al calcular los valores
    # esperados de las tasas
    amt_der = cshflw.get_amount_derivatives()

    # Se calcula el factor de descuento al plazo del flujo.
    df = zcc.get_discount_factor_at(fecha_val.day_diff(cshflw.get_settlement_date()))

    #
    amt_der = [a * df * delta for a in amt_der]
    if len(amt_der) > 0:
        result.append(np.array(amt_der))
total = result[0] * 0

for r in result:
    total += r

for i in range(len(total)):
    print(f"Sensibilidad en {i}: {total[i]:0,.2f}")

print(f"Sensibilidad de proyección: {sum(total):,.2f} USD")
```

    Sensibilidad en 0: -0.00
    Sensibilidad en 1: -0.00
    Sensibilidad en 2: -0.00
    Sensibilidad en 3: -0.00
    Sensibilidad en 4: -0.00
    Sensibilidad en 5: -0.00
    Sensibilidad en 6: -0.00
    Sensibilidad en 7: -0.00
    Sensibilidad en 8: -0.00
    Sensibilidad en 9: -0.00
    Sensibilidad en 10: -0.00
    Sensibilidad en 11: -0.00
    Sensibilidad en 12: -0.00
    Sensibilidad en 13: -0.00
    Sensibilidad en 14: -0.00
    Sensibilidad en 15: -0.70
    Sensibilidad en 16: -0.00
    Sensibilidad en 17: -1,997.46
    Sensibilidad en 18: -0.00
    Sensibilidad en 19: -0.00
    Sensibilidad en 20: -0.00
    Sensibilidad en 21: -0.00
    Sensibilidad en 22: -0.00
    Sensibilidad en 23: -0.00
    Sensibilidad en 24: -0.00
    Sensibilidad en 25: -0.00
    Sensibilidad en 26: -0.00
    Sensibilidad en 27: -0.00
    Sensibilidad en 28: -0.00
    Sensibilidad en 29: -0.00
    Sensibilidad en 30: -0.00
    Sensibilidad en 31: -0.00
    Sensibilidad en 32: -0.00
    Sensibilidad de proyección: -1,998.16 USD

**Ejercicio:** Ambas sensibilidades se cancelan. ¿Porqué?

En fórmula:

Cashflow = F(z1, z2, z3) / Df(z1, z2, z3)

Al calcular las derivadas del valor presente, sólo derivamos respecto a z1, z2, z3 en el denominador.

Al calcular las derivadas del flujo, derivamos respecto a z1, z2 y z3 en el numerador, lo que coincide con derivar la función F respecto a z1, z2 y z3.