text:
# Tema 2: Sistemas de ecuaciones numéricas lineales

```python
!pip install -r https://raw.githubusercontent.com/alexmascension/ANMI/main/requirements.txt
```

```python
from sympy import *
from sympy.matrices import Matrix as mat
from sympy.matrices import randMatrix
from sympy import symbols
import sympy

import matplotlib.pyplot as plt

import numpy as np

from scipy.linalg import orth
```

```python
from anmi.genericas import norma, print_verbose, norma_1, norma_inf, norma_2

from numpy.linalg import cond as numero_condicion
from anmi.T2 import descomposicion_LU, descomposicion_LDU, cholesky, gram_schmidt, householder, factorizacion_QR
from anmi.T2 import metodo_iterativo, criterio_m_matriz, criterio_SOR, criterio_simetrica_definida_positiva, criterio_diagonal_dominante, criterio_radio_espectral
```

### Número de condicion

El número de condición se define como $\vert\vert A \vert\vert \cdot \vert\vert A^{-1}\vert\vert$. Un número de condición cercano a 1 implica una matriz más estable frente a métodos con elementos diferenciales. Se cumple que las matrices ortogonales tienen número de condición 1.

```python
help(numero_condicion)
```

```python
M = mat(((1, 2, 3), (2, 3, 1), (3, 2, 4)))
M
```

```python
numero_condicion(np.array(M).astype(int))
```

```python
M_ort = (orth(np.array(M).astype(int)))
M_ort
```

```python
numero_condicion(M_ort)
```

### Factorización LU y LDU*

Recordemos que para una matriz, la factorización LU es el proceso de aplicación de la simplificación de Gauss, de modo que la matriz $L$ es una matriz triangular inferior con los coeficientes de transformación, y la matriz $U$ es la matriz superior con los elementos tras las transformaciones lineales.

Además, se puede hacer que $D$ sea una matriz diagonal con los valores de la diagonal de $U$, de modo que $LU$ = $LDD^{-1}U$, y si hacemos $U^* = D^{-1}U$ entonces tenemos $LDU^*$, donde $U^*$ sigue siendo una matriz diagonal superior, pero con la diagonal igual a 1.

A la hora de aplicar la factorización LU y LDU* se puede hacer una permutación de filas, de modo que en cada iteración se coge la fila con mayor valor (de entre las que no se han *procesado*) y se permuta, garantizando una solución siempre. También, es importante tener en mente que la factorización falla si algún elemento de la diagonal (desde el principio o durante la factorización) es 0, de modo que para solucionar ese caso se aplica la permutación.

Todas las permutaciones quedan recogidas en una matriz $P$, de modo que $$LU = LDU^* = PA$$

```python
help(descomposicion_LU)
```

```python
M = mat(((1, 4, 4), (3, 2, 1), (2, 4, 1)))
```

```python
descomposicion_LU(M, permutar_max=False)
```

```python
descomposicion_LDU(M, permutar_max=False)
```

```python
descomposicion_LU(M, rhs=ones(M.shape[0], 1), permutar_max=True)
```

```python
descomposicion_LDU(M, permutar_max=True)
```

### Factorización de Cholesky

La factorización de Cholesky es una factorización que genera una matriz triangular inferior $L$ tal que $A = LL^T$. Para que una matriz sea factorizable, tiene que cumplir que sus menores principales sean positivos, y que sea simétrica.

```python
help(cholesky)
```

```python
M = mat(((2, -1, 0), (-1, 2, -2), (0, 2, 1)))
M
```

```python
cholesky(M)
```

```python
cholesky(M + M.T)
```

```python
M + M.T
```

```python
cholesky(M + M.T) * cholesky(M + M.T).T
```

```python
# Podemos hacer también Cholesky a una matriz con símbolos!
x = symbols('x')

Mx = mat(((1, x, 1), (x, 2, 1), (1, 1, 3)))
Mx
```

```python
cholesky(Mx)
```

### Ortogonalización de Gram-Schmidt

```python
help(gram_schmidt)
```

```python
M
```

```python
GS = gram_schmidt(M)
GS['P']
```

```python
GS['Pn']
```

```python
GS['c']
```

```python
GS['P'][:, 1].T * GS['P'][:, 2]
```

```python
Mx = mat(((1, 2, 3), (1, x, 1), (0, 0, 3)))
```

```python
GSx = gram_schmidt(Mx)
GSx['P']
```

### Transformación de Householder

La transformación de Householder es una transformación para pasar de un vertor $x$ a un vector $y$. Para ello se toma el vector $e$, que sería el eje de transformación, y la matriz aplicación es $H = I - 2ee^t$.

En este caso, $$e = \pm \frac{x - y}{\vert\vert x - y \vert\vert}$$

Para la transformación de Householder $x$ e $y$ **tienen que tener la misma norma**. El vector resultante $e$ tiene norma 1.

```python
help(householder)
```

```python
v0 = mat(((1, 0)))
vf = mat((1, 1))
vf = vf/norma(vf)
```

```python
H, e = householder(v0, vf, normalizar=False)
H
```

```python
e
```

```python
# Ejercicio 10
t = symbols('t')

v0 = mat(((1, 0)))
vf = mat((cos(t), sin(t)))

H, e = householder(v0, vf, normalizar=False)
```

```python
H
```

```python
e
```

### Factorización QR

La factorización QR consiste en transformar $A = QR$ donde $Q$ es ortogonal y $R$ es triangular superior.

Si $P$ es la matriz ortogonalizada, $C$ es la matriz con los factores de ortonormalización (para Gram-Schmidt, por ejemplo, es $m_{ij} = \frac{a^j\cdot p^i}{\vert\vert p^{i}\vert\vert^2}$ y $D$ es la matriz de las normas de los vectores ortogonales ($\vert\vert p^i\vert\vert$), entonces se tiene que:
$$Q = PD^{-1}$$
$$R = D(I + C)$$

```python
help(factorizacion_QR)
```

```python
M = mat(((2, -1, 0), (0, 0, -2), (0, 2, -1)))
```

```python
factorizacion_QR(M, metodo='householder')
```

```python
A1 = mat(((2, -1, 0), (0, 0, -2), (0, 2, -1)))
```

```python
A1 = mat(((1, 2, 3), (4, 5, 6), (7, 8, 9)))
A1
```

```python
dqr = factorizacion_QR(A1, metodo='householder')
```

```python
dqr['Q']
```

```python
dqr['R']
```

```python
simplify(dqr['Q'].T * dqr['Q'])
```

```python
simplify(dqr['Q'] * dqr['R'])
```

### Métodos iterativos

Los métodos iterativos son métodos para resolver el sistema $Ax=b$. Todos los métodos se basan en el mismo criterio, que es emplear una matriz $H$ tal que, para una iteración:
$$x^{(k+1)} = Hx^{(k)} + b'$$

Si tomamos la matriz $A = D - L - U$ (donde $D$ es la diagonal de $A$, y $L$ y $U$ son matrices triangulares de elementos de $A$), la ecuación anterior queda:
$$Mx^{(k+1)} = (M-A)x^{(k)} + b$$

Y simplemente hay que despejar para $x^{(k+1)}$ (pues $x^{(k)}$ lo conocemos) y repetir el procedimiento hasta convergencia:
$$x^{(k+1)} = M^{-1}(M-A)x^{(k)} + M^{-1}b$$

Según el método, la matriz difiere: para Jacobi es $M = D$, para Gauss-Seidel es $M = D-L$ y para sobre-relajación (ros) es $M = \frac{1}{\omega}D - L$.

```python
help(metodo_iterativo)
```

```python
A = mat([[2, 1, 1], [1, -2, 1], [1, 1, 2]])
b = mat([[1, 1, 1]]).T
```

```python
diter = metodo_iterativo(A, b, n_iter=200, verbose=True)
np.array(diter['x'], dtype=float)
```

```python
# EJERCICIO 12
a = symbols('a')

A = mat([[4, 1, a], [1, 4, 1], [a, 1, 4]])
metodo_iterativo(A, metodo='gs', n_iter=3)
```

```python
# EJERCICIO 14
A = mat([[-11, 20, 8], [20, 16, -8], [8, -8, 5]]) / 27
norma_1(A)
norma_inf(A)
d = [simplify(i) for i in list((A.T * A).eigenvals().keys())]
# norma_2(A)
```

```python
A
```

```python
A.T * A
```

```python
simplify((A.T * A - a * eye(3)))
```

# Ejercicios

#### Ejercicio 14
Calcular $\vert\vert A\vert\vert_1, \vert\vert A\vert\vert_2, \vert\vert A\vert\vert_\infty$ de la matriz

```python
A = Matrix([[-11, 20, 8], [20, 16, -8], [8, -8, 5]]) / 27
A
```

```python
np.sum(abs(A), 0), norma_1(A)
```

```python
np.sum(abs(A), 1), norma_inf(A)
```

```python
# Para la norma_2 podemos usar norma_2(), pero vamos a hacerlo a mano. Para ello:
# 1) Hallamos A.T * A
AA = simplify(A.T * A)
AA
```

```python
# En este punto los resultados difieren, así que usaremos sus resultados intermedios
AA = Matrix([[65, 4, -32], [4, 80, 8], [-32, 8, 17]]) / 81
```

```python
# 2) Hallamos los autovalores
factor(det(AA - symbols('lambda') * eye(3)))
```

```python
# Los autovalores son 0 y 1. Luego nos quedamos con 1
rho = max(solve(det(AA - symbols('lambda') * eye(3))))
norma_2 = sqrt(rho)
norma_2
```

#### Ejercicio 15
La matriz de Hilbert $H_n$ tiene como coeficientes $h_{ij} = \frac{1}{i + j -1}$
Determinar el límite para $\vert\vert \; \vert\vert_1$ y $\vert\vert \; \vert\vert_\infty$: $\lim_{n\to \infty}\vert\vert H_n \vert\vert$.

```python
def Hn(n):
    hn = zeros(n, n)
    for row in range(n):
        for col in range(n):
            hn[row, col] = S(1) / S(row + col + 1) # porque i y j les sumamos 1

    return hn
```

```python
# Vamos a crear unas matrices de prueba
```

```python
Hn(7)
```

Vemos que la matriz es simétrica, luego norma 1 y norma inf serán iguales.
La columna / fila con mayor suma es la primera, pues:
$A[1:, 0] = A[:-1, 1]$ y $A[0, 0] > A[-1, 1]$

Y por inducción en el resto de columnas se ve que $\sum A[:, i] > \sum A[:, i+1]$.

El valor para $n$ de la primera fila es $\sum A[:, 0] = \sum_i^n 1/i$, que es la serie harmónica.

Para $n \to \infty$, sabemos que $\sum_i^n 1/i \to \infty$, luego $\lim_{n\to \infty}\vert\vert H_n \vert\vert = \infty$

#### Ejercicio 16
Se considera la matriz

$$ A = \begin{pmatrix}a & 1+a\\0 & a\end{pmatrix}$$

con $a > 0 \in \mathbb{R}$.

1) Calcular el número de condición de $A$ en la norma $\vert \vert A \vert \vert_\infty$$

```python
# El núemro de condición es ||A||·||A^-1||
a = symbols('a')

A = Matrix([[a, 1+a], [0, a]])
A
```

```python
Ainv = A.inv()
Ainv
```

```python
# Norma_inf A es a + a + 1 = 1 + 2a
ninfA = 1 + 2*a

# Norma_inf Ainv es 1/a + 1/a + 1/a2 =
ninfAinv = 1/a + 1/a + 1/a**2

ninfA, ninfAinv
```

```python
condicion = ninfA * ninfAinv
simplify(condicion)
```

```python
simplify(1 - factor(condicion))
```

2) Estimar el error relativo de la solución del sistema lineal perturbado
$$(A + \delta A)(x + \delta x) = b + \delta b$$

Por el lema de Banach tenemos que

$$\frac{\vert\vert \delta x\vert\vert}{\vert\vert x\vert\vert} = \frac{cond(A)}{1-cond(A)\frac{\vert\vert \delta A\vert\vert}{\vert\vert A\vert\vert}}\left(
    \frac{\vert\vert \delta b\vert\vert}{\vert\vert b\vert\vert} +
    \frac{\vert\vert \delta A\vert\vert}{\vert\vert A\vert\vert}\right)$$

Para este ejercicio es
$$\frac{\vert\vert \delta x\vert\vert}{\vert\vert x\vert\vert} = \frac{
    (2a+1)^2
}{(a^2 - (2a+1)^2)\frac{\vert\vert \delta A\vert\vert}{\vert\vert A\vert\vert}}\left(
    \frac{\vert\vert \delta b\vert\vert}{\vert\vert b\vert\vert} +
    \frac{\vert\vert \delta A\vert\vert}{\vert\vert A\vert\vert}\right)$$

De lo cual se cumple que
$$\frac{\vert\vert \delta A\vert\vert}{\vert\vert A\vert\vert} < \frac{1}{cond(A)}$$

#### Ejercicio 17
Calcular la factorización LU de la matriz $A$ con pivote parcial con $0 < \alpha < \frac{\pi}{4}$.

```python
a = symbols('alpha')
A = Matrix([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, cos(a), -sin(a)], [0, 0, sin(a), cos(a)]])
A
```

En este caso tenemos que la función LU nos falla porque el máximo de la matriz "no existe". Vamos a programarlo.

Vamos a emplear la factorización de Doolitle del estilo $PA = LU$. En este caso, $P = P^{(1)}$

```python
P = Matrix([[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
P
```

```python
dict_LU = descomposicion_LU(P * A, permutar_max=False)
```

```python
dict_LU['L']
```

```python
dict_LU['U']
```

```python
simplify(dict_LU['L'] * dict_LU['U'])
```

Vemos que tras la primera permutación, la función funciona sin problemas.

### Ejercicio 18

Sea $L$ la matriz con coeficientes $l_{ij} = \max\{i-j+1, 0\}$ para todo $i, j = 1, \cdots, n$. Se sabe que $L$ es triangular inferior que aparece en la factorización de Cholesky de una matriz $A$ de $n$ filas y $n$ columnas, de la que únicamente se sabe que su coeficiente $a_{n2} = 14$. Determinar la dimensión de $A$

Es decir, $a_{n2} = 14$, y $a_{ij} = \sum_{k \le \min\{i,j\}}^n l_{ik}l_{jk}$
Luego $a_{n2} = l_{n1}l_{21} + l_{n2}l_{22} = (n) \cdot (2) + (n - 1) \cdot (1) = 2n + n - 1 = 3n - 1$

Luego $3n - 1 = 14$ y $n = 15/3 = 5$.

Por tanto, la matriz $L$ es:

```python
L = zeros(5, 5)
for i in range(5):
    for j in range(5):
        ij = max(i - j + 1, 0)
        L[i, j] = ij
L
```

Y $A$ es

```python
A = L * L.T
A
```

Y se ve que $a_{52} = a_{n2} = 14$

### Ejercicio 19

Determinar en que rango los valores de $\alpha$ de la matriz $A$ es invertible y admite una factorización de Cholesky.

```python
a = Symbol('alpha')

A = Matrix([[2, 1, a], [1, 2, 1], [a, 1, 2]])
A
```

La matriz es invertible en aquellos casos donde su determinante no sea nulo.

```python
factor(det(A))
```

Vemos que las soluciones son $\alpha = 2$ y $\alpha = -1$. En esos casos la matriz no es invertible. Además, existe una factorización de Cholesky para los casos en los que $\det(A) > 0$, que se cuple cuando $1 < \alpha < 2$. En esos casos la factorización es

```python
cholesky(A)
```

### Ejercicio 20

Determinar la factorización de Cholesky de $A_{n\times n}$ con coeficientes.
$$a_{11} = 1$$
$$a_{ii} = 1 + i^2$$
$$a_{i, i+1} = a_{i+1, i} = -i$$

Calcular el determinante de $A$

Vamos a sacar algunos ejemplos ilustrativos para hacernos una idea.

```python
n = 7
A = zeros(n, n)
for i in range(n):
    for j in range(n):
        if i == j:
            A[i,j] = (i+1) ** 2 + 1
        elif i == (1 + j):
            A[i,j] = -i
        elif j == (1 + i):
            A[i,j] = -j
        else:
            pass
A[0, 0] = 1
A
```

```python
cholesky(A)
```

Así a primera vista vemos el patrón evidente. Vamos a demostrarlo.

Sabemos que $l_{11} = 1$ y $l_{21} = -1$

Vamos a probar un par de casos.

Para $i = 2$ tenemos que
$$l_{22} = \sqrt{a_{22} - l_{21}^2} =
\sqrt{1 + 2^2 - (-1)^2} = 2$$
$$l_{32} = \frac{a_{32} - \sum_k^1l_{3k}l_{2k}}{l_{22}} =
 \frac{-2 - (-1)(0)}{2}  =  -1$$
Y $l_{31} = 0$

Para $i=3$ tenemos

$$l_{33} = \sqrt{a_{33} - \sum_{k=1}^{2}l_{3k}^2} =
\sqrt{1 + 3^2 - l_{32}^2} = \sqrt{1 + 3^2 - (-1)^2} = 3$$
$$l_{43} = \frac{-3 - l_{42}l_{32}}{l_{33}} =
 \frac{-3 - (-1)(0)}{3}  =  -1$$

Y, por inducción:
$$l_{ii} = \sqrt{a_{ii} - l_{i, i-1}^2} = \sqrt{1 + i^2 - 1} = i$$
$$l_{i+1, i} = \frac{a_{i+1,i}}{l_ii} = \frac{-i}{i} = -1$$

Luego la matriz es de la forma:
$\det{A} = \det({L})\det({L^T}) = \det({L})^2 = (n!)^2$

### Ejercicio 21

Construir una trnaformación de Householder que transforme el vector $x = (\cos\alpha, -\sin\alpha, 1, 1)^t$ en el vector $y = (1, 1, \sin\alpha, \cos\alpha)^t.

Encontrar la matriz asociada a esta transformación para $\alpha = \pi/2$

```python
vx = Matrix([cos(a), -sin(a), 1, 1])
vy = Matrix([1, 1, sin(a), cos(a)])
vx
```

```python
H, e = householder(vx, vy, normalizar=True)
simplify(H)
```

```python
simplify(e)
```

```python
H.subs(a, pi/2)
```

### Ejercicio 22

Sean $A$ y $b$. Aplicar el método de Gram-Schmidt para ortonormalizar el conjunto de vectores $\{b, Ab, A^2b\}$.

```python
A = Matrix([[1, 2, 0], [2, 1, 2], [0, 2, 1]])
b = Matrix([1, -1, 1])
```

```python
T = zeros(3, 3)
T[:, 0] = b
T[:, 1] = A * b
T[:, 2] = A ** 2 * b
T
```

```python
GS = gram_schmidt(T)
GS['P']
```

```python
GS['Pn']
```

### Ejercicio 23

Se considera el sistem lineal $Ax=b$. Hallar el valor de $\beta$ que minimice el radio espectral de $I - \beta A$

```python
A = Matrix([[1, 1, 2], [0, 2, -1], [0, -1, 2]])
b = Matrix([1, 1, 1])
beta = Symbol('beta')
```

Primero vamos a hallar el radio espectral de la matriz $I - \beta A$

```python
iba = eye(3) - beta * A
iba
```

```python
l = Symbol('lambda')
eq = factor(det(iba - eye(3) * l))
eq
```

Vemos que los autovalores de $I - \beta A$ son $\lambda = 1 - \beta$ y
$\lambda = 1 - 3\beta$

```python
x = np.linspace(-1, 1, 100)
plt.plot(x, [max(abs(1-i), abs(1-3*i)) for i in x])
```

Vemos que el menor valor está en $\beta = 0.5$ donde $\rho(I - \beta A) = 0.5$

### Ejercicio 24

Se considera el sistema lineal
$$x_{i+1} = x_i + x_{i-1} \;\;\; i = 1, 2, 3$$
$$x_0 = 1, x_4 = 5$$

Es decir:
$$\begin{align}
x_2 & =  x_1 + x_0\\
x_3 & =  x_2 + x_1\\
x_4 & = x_3 + x_2\\
\end{align}
$$

Resolviendo parámetros
$$\begin{align}
x_2 & =  x_1 + 1\\
x_3 & =  x_2 + x_1\\
5 & = x_3 + x_2\\
\end{align}
$$

Analizar la convergencia del método de jacobi cuando se aplica a la resolución de este sistema de ecuaciones.

El sistema de ecuaciones se puede reordenar de la siguiente manera:
$$\begin{align}
-x_1 + x_2 + 0 \cdot x_3 & =  1\\
-x_1 - x_2 + x_3 & =  0\\
0\cdot x_1 - x_2 - x_3 & = 5\\
\end{align}
$$

Que, simplificando signos, queda:
$$
\begin{align}
-x_1 + x_2 + 0  & =  1\\
x_1 + x_2 - x_3 & =  0\\
0 + x_2 + x_3 & = -5\\
\end{align}
$$

En sistema matricial queda como
$$
\begin{pmatrix}
-1 & 1 & 0 \\
1 & 1 & -1 \\
0 & 1 & 1 \\
\end{pmatrix}
\begin{pmatrix}
x_1\\x_2\\x_3
\end{pmatrix} =
\begin{pmatrix}
1 \\ 0 \\ -5
\end{pmatrix}
$$

```python
A = Matrix([[-1, 1, 0], [1, 1, -1], [0, 1, 1]])
D = Matrix([[-1, 0, 0], [0, 1, 0], [0, 0, 1]])
b = Matrix([1, 0, -5])
```

```python
help(metodo_iterativo)
```

```python
metodo_iterativo(A=A, b=b, metodo='jacobi')
```

Por el criterio del radio espectral, tenemos que los autovalores son $\pm\sqrt{2}i$, y el radio espectral es $\sqrt{2} > 1$, luego el método de Jacobi no converge.
De hecho vemos que la matriz de diferencias crece con cada iteración.

### Ejercicio 25

Analizar la convergencia de los métodos de Jacobi, Gauss-Seidel y SOR cuando se aplica el sistema lineal

$$
\begin{pmatrix}
4 & 1 & 0 & \cdots & 0 & 0 & 0 \\
1 & 4 & 1 & \cdots & 0 & 0 & 0 \\
0 & 1 & 4 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 4 & 1 & 0 \\
0 & 0 & 0 & \cdots & 1 & 4 & 1 \\
0 & 0 & 0 & \cdots & 0 & 1 & 4 \\
\end{pmatrix}
\begin{pmatrix}
x_1\\x_2\\x_3\\\vdots\\x_{n-2}\\x_{n-1}\\x_n
\end{pmatrix}
=
\begin{pmatrix}
1\\2\\3\\\vdots\\n-2\\n-1\\n
\end{pmatrix}
$$

```python
n = 8
A = eye(n) * 4
b = ones(n, 1)
for i in range(n-1):
    A[i+1, i], A[i, i+1], b[i+1, 0] = 1, 1, i+2
A
```

```python
b
```

```python
dict_jac = metodo_iterativo(A=A, b=b, metodo='jacobi')
N(dict_jac['x']), N(Matrix(dict_jac['diff']))
```

```python
dict_gs = metodo_iterativo(A=A, b=b, metodo='gs')
N(dict_gs['x']), N(Matrix(dict_gs['diff']))
```

```python
dict_sor = metodo_iterativo(A=A, b=b, metodo='sor')
N(dict_sor['x']), N(Matrix(dict_sor['diff']))
```

Vemos que las matrices son convergentes, y los criterios inmediatos de convergencia se justifican. Vamos a estudiar cada uno de ellos en detalle.

**Criterio de radio espectral**: vemos que hay autovalores mayores que 1, así que no podemos asegurar la convergecia del sistema por este criterio.

```python
solve(factor(det(A - eye(n) * l)), l)
```

**Criterio de diagonal dominante**: 4 > 2, en base al criterio, los métodos de Jacoby y Gs son convergentes.

**Criterio SOR**: Si la matriz es simétrica definida positiva y 0 < $\omega$ < 2, el método SOR converge.

Para probar que este criterio se cumple tenemos que comprobar que el determinante de la matriz es > 0 para todo $n$. Para $n = 1, 2, 3$ los determinantes son $4, 14, 56$.

Asumamos que para una matriz de $n-1$ y $n-2$ su determinante es positivo. Veamos para el caso de $n$. Desarrollando el determinante por filas tenemos que

$$det(A_n) = 4 \cdot det(A_{n-1}) - A_{n-2}$$

Lo de $A_{n-2}$ lo podemos hacer desarrollando un poco los términos:
$$
det\begin{pmatrix}
\color{gray}{4} & \color{gray}{1} & \color{gray}{0} & \cdots & \color{gray}{0} & \color{gray}{0} & \color{gray}{0} \\
1 & \color{gray}{4} & 1 & \cdots & 0 & 0 & 0 \\
0 & \color{gray}{1} & 4 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & \color{gray}{0} & 0 & \cdots & 4 & 1 & 0 \\
0 & \color{gray}{0} & 0 & \cdots & 1 & 4 & 1 \\
0 & \color{gray}{0} & 0 & \cdots & 0 & 1 & 4 \\
\end{pmatrix}_{n\times n}
\rightarrow
det\begin{pmatrix}
1 & 1 & 0 & \cdots & 0 & 0 & 0 \\
0 & 4 & 1 & \cdots & 0 & 0 & 0 \\
0 & 1 & 4 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 4 & 1 & 0 \\
0 & 0 & 0 & \cdots & 1 & 4 & 1 \\
0 & 0 & 0 & \cdots & 0 & 1 & 4 \\
\end{pmatrix}_{n-1\times n-1}
=
det\begin{pmatrix}
4 & 1 & 0 & \cdots & 0 & 0 & 0 \\
1 & 4 & 1 & \cdots & 0 & 0 & 0 \\
0 & 1 & 4 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 4 & 1 & 0 \\
0 & 0 & 0 & \cdots & 1 & 4 & 1 \\
0 & 0 & 0 & \cdots & 0 & 1 & 4 \\
\end{pmatrix}_{n-2\times n-2}$$

Si asumimos que $\det(A_{n-2}) < \det(A_{n-1})$ entonces se cumple que $A_n > A_{n-1} > A_{n-2} > 0$ y el criterio SOR se cumple.

max_stars_repo_name: alexmascension/ANMI

max_issues_repo_name: alexmascension/ANMI

max_forks_repo_name: alexmascension/ANMI


text:
# 足先（エンドエフェクタ）のジョイント角度ベクトルによるヤコビ行列の導出

制御したいもの -> 足先の2次元位置 + 足先の地面となす角度

足先の状態変数を次のようにおく．

```python
import sympy as sy
from sympy import pi, cos, sin, tan
from IPython.display import display
```

```python
# 式を定義
q1 = sy.Symbol("q1")
q2 = sy.Symbol("q2")
q3 = sy.Symbol("q3")

qvec = sy.Matrix([[q1, q2, q3,]]).T

l1 = sy.Symbol("l1")
l2 = sy.Symbol("l2")
l3 = sy.Symbol("l2")
l4 = sy.Symbol("l2")

D1 = sy.Matrix([
        [cos(q1), -sin(q1), 0],
        [sin(q1), cos(q1), 0],
        [0, 0, 1],
    ])

D2 = sy.Matrix([
        [1, 0, l1],
        [0, 1, 0],
        [0, 0, 1],
    ])

D3 = sy.Matrix([
        [cos(q2), -sin(q2), 0],
        [sin(q2), cos(q2), 0],
        [0, 0, 1],
    ])

D4 = sy.Matrix([
        [1, 0, l2],
        [0, 1, 0],
        [0, 0, 1],
    ])

theta = 270/180*pi
D5 = sy.Matrix([
        [cos(theta), -sin(theta), 0],
        [sin(theta), cos(theta), 0],
        [0, 0, 1],
    ])

D6 = sy.Matrix([
        [1, 0, l3],
        [0, 1, 0],
        [0, 0, 1],
    ])

D7 = sy.Matrix([
        [cos(q3), -sin(q3), 0],
        [sin(q3), cos(q3), 0],
        [0, 0, 1],
    ])

D8 = sy.Matrix([
        [1, 0, l4],
        [0, 1, 0],
        [0, 0, 1],
    ])
```

```python
# 足先（エンドエフェクター）の位置
D_end = D1 * D2 * D3 * D4 * D5 * D6 * D7 * D8
D_end = sy.simplify(D_end)

phi = sy.Symbol('phi')
phi = q1 + q2 + q3 -2*pi

ee = sy.Matrix([
    [D_end[0, 2]],
    [D_end[1, 2]],
    [phi]
])
display(ee)
```

$\displaystyle \left[\begin{matrix}l_{1} \cos{\left(q_{1} \right)} + \sqrt{2} l_{2} \sin{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} + l_{2} \sin{\left(q_{1} + q_{2} + q_{3} \right)}\\l_{1} \sin{\left(q_{1} \right)} - \sqrt{2} l_{2} \cos{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} - l_{2} \cos{\left(q_{1} + q_{2} + q_{3} \right)}\\q_{1} + q_{2} + q_{3} - 2 \pi\end{matrix}\right]$

```python
jacobi = ee.jacobian(qvec)
jacobi = sy.simplify(jacobi)
```

結果を表示

```python
display(jacobi)
```

$\displaystyle \left[\begin{matrix}- l_{1} \sin{\left(q_{1} \right)} + \sqrt{2} l_{2} \cos{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} + l_{2} \cos{\left(q_{1} + q_{2} + q_{3} \right)} & l_{2} \left(\sqrt{2} \cos{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} + \cos{\left(q_{1} + q_{2} + q_{3} \right)}\right) & l_{2} \cos{\left(q_{1} + q_{2} + q_{3} \right)}\\l_{1} \cos{\left(q_{1} \right)} + \sqrt{2} l_{2} \sin{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} + l_{2} \sin{\left(q_{1} + q_{2} + q_{3} \right)} & l_{2} \left(\sqrt{2} \sin{\left(q_{1} + q_{2} + \frac{\pi}{4} \right)} + \sin{\left(q_{1} + q_{2} + q_{3} \right)}\right) & l_{2} \sin{\left(q_{1} + q_{2} + q_{3} \right)}\\1 & 1 & 1\end{matrix}\right]$

ext: ipynb

max_stars_repo_name: YoshimitsuMatsutaIe/quadruped

max_issues_repo_name: YoshimitsuMatsutaIe/quadruped

max_forks_repo_name: YoshimitsuMatsutaIe/quadruped


text:
$$\begin{align}
&\text{Eric Nguyen} &\qquad \text{January 16, 2019} \\
&\textbf{Calculus Ch. 1 Study Guide} \\
&\textbf{Directions: } \text{Please complete each of the following problems to the best of your ability.} \\
&\text{If you have a question, I will do my best to answer them.} \\
&\text{Work honestly and with integrity, and always } \textbf{DO YOUR BEST WORK.} \\\\
\end{align}$$

```python
%matplotlib inline
%config InlineBackend.figure_format='retina'
from numpy import e, power, log, exp, linspace
from matplotlib.pyplot import plot
```

$\text{Differentiate.}$

$1)$
$$\begin{align} \\
y &= xe^{-2x} + e^{-x} + x^3 \\
y' &= -2xe^{-2x} + e^{-2x} - e^{-x} + 3x^2 \\
&= -xe^{-2x} - e^{-x} + 3x^2
\end{align}$$

```python
y = lambda x: x * exp(-2 * x) + exp(-x) + power(x, 3)
dy = lambda x: -x * exp(-2 * x) - exp(-x) + 3 * power(x, 2)
x = linspace(0, 10, 100)
plot(x, y(x))
_ = plot(x, dy(x))
```

$2)$
$$\begin{align}
f(t) &= \ln{\frac{1 - t}{1 + t}} \\
&= \ln{(1 - t)} - \ln{(1 + t)} \\
f'(t) &= -\frac{1}{1 - t} - \frac{1}{1 + t}
\end{align}$$

```python
f = lambda t: log((1 - t) / (1 + t))
my_df = lambda t: -(1 / (1 - t)) - (1 / (1 + t))
x = linspace(-0.9, 0.9, 1000)
plot(x, f(x))
_ = plot(x, my_df(x))
```

$3)$
$$\begin{align}
y &= 7^{x^4+2} \\
y' &= \ln{7} \cdot 7^{x^4+2} \cdot 4x^3
\end{align}$$

```python
y = lambda x: power(7, power(x, 4) + 2)
dy = lambda x: log(7) * y(x) * 4 * power(x, 3)
x1 = linspace(1.8, 2, 100)
x2 = linspace(1.8, 1.93, 100)
plot(x1, y(x1))
_ = plot(x2, dy(x2))
```

$4)$
$$\begin{align}
y &= 6^x \cdot \log_7{x} \\
y' &= 6^x \cdot \frac{1}{\ln{7}} \cdot \frac{1}{x} + \log_7{x} \cdot \ln{6} \cdot 6^x \\
&= \frac{6^x}{x \cdot \ln{7}} + \log_7{x} \cdot \ln{6} \cdot 6^x
\end{align}$$

```python
log7 = lambda x: log(x) / log(7)
y = lambda x: power(6, x) * log7(x)
dy = lambda x: power(6, x) * (1 / log(7)) * (1 / x) + log7(x) * log(6) * power(6, x)
x = linspace(7, 10, 100)
plot(x, y(x))
_ = plot(x, dy(x))
```

$5)$
$$\begin{align}
y &= \ln{\left(\left(-\frac{4^{x^4}}{x^3-3}\right)^5\right)} \\
y' &= \frac{1}{u} \cdot u' \\
u &= \left(-\frac{4^{x^4}}{x^3-3}\right)^5 \\
v &= -\frac{4^{x^4}}{x^3-3} \\
v' &= -\frac{(x^3 - 3) \cdot (\ln{4} \cdot 4^{x^4} \cdot 4x^3) - (x^3 - 3) \cdot (3x^2)}{(x^3 - 3)^2} \\
&= -\frac{(\ln{4} \cdot 4^{x^4} \cdot 4x^3) - (3x^2)}{(x^3 - 3)} \\
u' &= 5v^4 \cdot v' \\
&= 5\left(-\frac{4^{x^4}}{x^3-3}\right)^4 \cdot -\frac{(\ln{4} \cdot 4^{x^4} \cdot 4x^3) - (3x^2)}{(x^3 - 3)} \\
y' &= \frac{1}{\left(-\frac{4^{x^4}}{x^3-3}\right)^5} \cdot 5\left(-\frac{4^{x^4}}{x^3-3}\right)^4 \cdot -\frac{(\ln{4} \cdot 4^{x^4} \cdot 4x^3) - (3x^2)}{(x^3 - 3)} \\
&=\frac{-5 \left((\ln{4} \cdot 4^{x^4} \cdot 4x^3) - (3x^2)\right)}{\left(-\frac{4^{x^4}}{x^3-3}\right)^5 \cdot (x^3 - 3)} \cdot \left(-\frac{4^{x^4}}{x^3-3}\right)^4 \\
\end{align}$$

```python
def v(x):
    return -((power(4, power(x, 4))) / (power(x, 3) - 3))
def u(x):
    return power(v(x), 5)
def y(x):
    return log(u(x))
def dy(x):
    return ((-5 * ((log(4) * power(4, power(x, 4)) * 4 * power(x, 3)) - (3 * power(x, 2)))) / (u(x) * (power(x, 3) - 3))) * power(v(x), 4)
x = linspace(0, 1.4, 100)
plot(x, y(x))
_ = plot(x, dy(x))
```

$\text{6) Solve for } t.$
$$\begin{align} \\
P &= P_0e^{kt} \\
\frac{P}{P_0} &= e^{kt} \\
\ln{\frac{P}{P_0}} &= kt \\
t &= \frac{\ln{\frac{P}{P_0}}}{k}
\end{align}$$

$\text{7) Find the equation of the line tangent to the graph of } f(x) = e^{-x} \text{ at the point (0, 1)}$
$$\begin{align} \\
&& f'(x) &= -e^{-x} \\
&& f'(0) &= -e^{0} \\
&& &= -1 \\
&& y &= mx + b \\
&& y &= f'(0)(x - 0) + f(0) \\
&& y &= -x + 1
\end{align}$$

```python
f = lambda x: exp(-x)
tangent_line = lambda x: -x + 1
x = linspace(-2, 2, 100)
plot(x, f(x))
_ = plot(x, tangent_line(x))
```

$\text{8) Find the equation of the line tangent to the graph of } f(x) = e^{3x} \cdot \ln{(4x)} \text{ at } x = 1$
$$\begin{align} \\
f'(x) &= e^{3x} \cdot \frac{4}{4x} + \ln{4x} \cdot 3e^{3x} \\
&= e^{3x} \left(\frac{1}{x} + 3\ln{4x}\right) \\
f'(1) &= e^{3(1)} \left(\frac{1}{(1)} + 3\ln{4(1)}\right) \\
&= e^3 \left(1 + 3 \ln{4}\right) \\
y &= f'(1)(x - 1) + f(1) \\
&= e^3 \left(1 + 3 \ln{4}\right)(x - 1) + e^3 \cdot \ln{4}
\end{align}$$

```python
f = lambda x: exp(3 * x) * log(4 * x)
def df(x):
    return exp(3 * x) * ((1 / x) + 3 * log(4 * x))
tangent_line = lambda x: df(1) * (x - 1) + f(1)
x = linspace(0.5, 1.5, 100)
plot(x, f(x))
_ = plot(x, tangent_line(x))
```

$\text{9) At a price of x dollars, the demand, in thousands of units, for a certain MP3 player is given by the demand function } q = 240e^{-0.003x}.$
$\text{a) How many MP3 players will be bought at a price of \$250? Round to the nearest thousand.}$
$\text{b) Find the marginal demand } q'(x).$
$\text{c) Interpret the meaning of the derivative.}$

```python
q = lambda x: 240 * exp(-0.003 * x)
```

```python
q(250)
```

    113.36797265784352

$\text{9a) Approximately 113 MP3 players will be bought at a price of \$250.}$

$\text{9b)}$

$$\begin{align}
q'(x) &= -0.003 \cdot 240e^{-0.003x} \\
&= -0.72e^{-0.003x}
\end{align}$$

```python
dq = lambda x: -0.72 * exp(-0.003 * x)
x = linspace(-1000, 750, 1000)
plot(x, q(x))
_ = plot(x, dq(x))
```

$\text{9c) The derivative shows how the number of MP3 players sold changes slowly as the price increases.}$

$\text{10) Since 2000, sales of bottled water have increased at the rate of approximately 9.3% per year.}$
$\text{That is, the volume of bottled water sold, } G \text{, in billions of gallons, } t \text{ years after 2000 is growing at a rate given by}$
$$\frac{dG}{dt} = 0.093G.$$
$\text{a) Find the function that satisfies the equation, given that approximately 4.7 billion gallons of bottled water were sold in 2000.}$
$\text{b) Predict the number of gallons of water that will be sold in 2010.}$
$\text{c) What is the doubling time for } G(t)?$

$\text{a) } G(t) = 4.7 * 10^9e^{0.093t}$

$\text{b) } G(10) = 4.7 * 10^9e^{0.093(10)}$

```python
4.7 * power(10, 9) * exp(0.093 * 10)
```

    11912193134.803917

$\approx 12 \text{ billion}$

$\text{c) } \dfrac{\ln{2}}{0.093}$

```python
log(2) / 0.093
```

    7.453195489891885

$\approx 7.45 \text{ years}$

$\text{11) Iodine-131 has a decay rate of 9.6% per day. The rate of change of an amount N of iodine-131 is given by}$
$$\frac{dN}{dt} = -0.096N.$$
$\text{a) Let } N_0 \text{ represent the amount of iodine-131 present at } t = 0. \text{ Find the exponential function that models the situation.}$
$\text{b) Suppose that 500g of iodine-131 is present at } t = 0. \text{ How much will remain after 4 days?}$
$\text{c) After how many days will half of the 500g iodine-131 will remain?}$

$\text{a) } N(t) = N_0e^{-0.096t}$

$\text{b) } N(4) = 500e^{-0.096(4)}$

```python
500 * exp(-0.096 * 4)
```

    340.56571358977357

$\approx 340.566 \text{g}$

$\text{c) } \dfrac{\ln{2}}{0.096}$

```python
log(2) / 0.096
```

    7.2202831308327635

$\approx 7.2 \text{ days}$

max_stars_repo_name: airicbear/calculus-homework

max_issues_repo_name: airicbear/calculus-homework

max_forks_repo_name: airicbear/calculus-homework


text:
## 6. Nonlinear ensemble filtering for the Lorenz-63 problem

In this notebook, we are interested in the sequential inference

References:

[1] Evensen, G., 1994. Sequential data assimilation with a nonlinear quasi‐geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans, 99(C5), pp.10143-10162.

[2] Asch, M., Bocquet, M. and Nodet, M., 2016. Data assimilation: methods, algorithms, and applications. Society for Industrial and Applied Mathematics.

[3] Bishop, C.H., Etherton, B.J. and Majumdar, S.J., 2001. Adaptive sampling with the ensemble transform Kalman filter. Part I: Theoretical aspects. Monthly weather review, 129(3), pp.420-436.

[4] Lorenz, E.N., 1963. Deterministic nonperiodic flow. Journal of atmospheric sciences, 20(2), pp.130-141.

[5] Spantini, A., Baptista, R. and Marzouk, Y., 2019. Coupling techniques for nonlinear ensemble filtering. arXiv preprint arXiv:1907.00389.

### The basic steps
To carry out sequential inference in `AdaptiveTransportMap`, we need to carry out a few basic steps:
* **Specify the problem**: Define the state-space model: initial condition, dynamical and observation models (including process and observation noise)
* **Specify the inflation parameters**: Determine the levels of covariance inflation to properly balance the dynamical system and the observations from the truth system
* **Specify the filter**: Choose the ensemble filter to assimilate the observations in the state estimate
* **Perform the sequential inference**: Perform the sequential inference

We will go through all of these here.

```julia
using Revise
using LinearAlgebra
using AdaptiveTransportMap
using Statistics
using Distributions
using OrdinaryDiffEq
```

    ┌ Info: Precompiling AdaptiveTransportMap [bdf749b0-1400-4207-80d3-e689c0e3f03d]
    └ @ Base loading.jl:1278
    ┌ Warning: Type annotations on keyword arguments not currently supported in recipes. Type information has been discarded
    └ @ RecipesBase ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:116
    ┌ Warning: Type annotations on keyword arguments not currently supported in recipes. Type information has been discarded
    └ @ RecipesBase ~/.julia/packages/RecipesBase/92zOw/src/RecipesBase.jl:116

```julia
using DelimitedFiles
```

Load some packages to make nice figures

```julia
using Plots
default(tickfont = font("CMU Serif", 9),
        titlefont = font("CMU Serif", 14),
        guidefont = font("CMU Serif", 12),
        legendfont = font("CMU Serif", 10),
        grid = false)
pyplot()

using LaTeXStrings
PyPlot.rc("text", usetex = "true")
PyPlot.rc("font", family = "CMU Serif")
# gr()
using ColorSchemes
```

The Lorenz-63  model is a three dimensional system that models the atmospheric convection [4]. This system is a classical benchmark problem in data assimilation. The state $\boldsymbol{x} = (x_1, x_2, x_3)$ is governed by the following set of ordinary differential equations:

\begin{equation}
\begin{aligned}
&\frac{\mathrm{d} x_1}{\mathrm{d} t}=\sigma(x_2-x_1)\\
&\frac{\mathrm{d} x_2}{\mathrm{d} t}=x_1(\rho-x_2)-x_2\\
&\frac{\mathrm{d} x_3}{\mathrm{d} t}=x_1 x_2-\beta x_3,
\end{aligned}
\end{equation}

where $\sigma = 10, \beta = 8/3, \rho = 28$. For these values, the system is chaotic and behaves like a strange attractor. We integrate this system of ODEs with time step $\Delta t_{dyn} = 0.05$. The state is fully observed $h(t,\boldsymbol{x}) = \boldsymbol{x}$ with $\Delta t_{obs}=0.1$. The initial distribution $\pi_{\mathsf{X}_0}$ is the standard Gaussian. The process noise is Gaussian with zero mean and covariance $10^{-4}\boldsymbol{I}_3$. The measurement noise has a Gaussian distribution with zero mean and covariance $\theta^2\boldsymbol{I}_3$ where $\theta^2 = 4.0$.

### Simple twin-experiment

Define the dimension of the state and observation vectors

```julia
Nx = 3
Ny = 3
```

    3

Define the time steps $\Delta t_{dyn}, \Delta t_{obs}$  of the dynamical and observation models. Observations from the truth are assimilated every $\Delta t_{obs}$.

```julia
Δtdyn = 0.05
Δtobs = 0.2
```

    0.2

Define the time span of interest

```julia
t0 = 0.0
tf = 100.0
Tf = ceil(Int64, (tf-t0)/Δtobs)
```

    500

 Define the distribution for the initial condition $\pi_{\mathsf{X}_0}$

```julia
π0 = MvNormal(zeros(Nx), Matrix(1.0*I, Nx, Nx))
```

    FullNormal(
    dim: 3
    μ: [0.0, 0.0, 0.0]
    Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    )

We construct the state-space representation `F` of the system composed of the deterministic part of the dynamical and observation models.

The dynamical model is provided by the right hand side of the ODE to solve. For a system of ODEs, we will prefer an in-place syntax `f(du, u, p, t)`, where `p` are parameters of the model.
We rely on `OrdinaryDiffEq` to integrate the dynamical system with the Tsitouras 5/4 Runge-Kutta method adaptive time marching.

We assume that the state is fully observable, i.e. $h(x, t) = x$.

```julia
h(x, t) = x
F = StateSpace(lorenz63!, h)
```

    StateSpace(AdaptiveTransportMap.lorenz63!, h)

Define the additive inflation for the dynamical and observation models

```julia
### Process and observation noise
σx = 1e-1
σy = 2.0

ϵx = AdditiveInflation(Nx, zeros(Nx), σx)
ϵy = AdditiveInflation(Ny, zeros(Ny), σy)
```

    AdditiveInflation(3, [0.0, 0.0, 0.0], [4.0 0.0 0.0; 0.0 4.0 0.0; 0.0 0.0 4.0], [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 2.0])

```julia
model = Model(Nx, Ny, Δtdyn, Δtobs, ϵx, ϵy, π0, 0, 0, 0, F);
```

To perform the nonlinear ensemble filtering, we first need to estimate the transport map $\boldsymbol{S}^{\boldsymbol{\mathcal{X}}}$.

In this notebook, we are going to assume that the basis of features does not change over time, but solely the coefficients $c_{\boldsymbol{\alpha}}$ of the expansion.

To estimate the map, we generate joint samples $(\boldsymbol{y}^i, \boldsymbol{x}^i), \; i = 1, \ldots, N_e$ where $\{\boldsymbol{x}^i\}$ are i.i.d. samples from pushforward of the standard Gaussian distribution by the flow of the Lorenz-63 system.

```julia
# Number of runs
Ne = 500

# Time span
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz63!, zeros(Nx), tspan)
```

    ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
    timespan: (0.0, 100.0)
    u0: 3-element Array{Float64,1}:
     0.0
     0.0
     0.0

Create cache to store the different final solutions

```julia
X = zeros(Ny + Nx, Ne);
```

Solve Lorenz-96 model for the different initial conditions and store the solutions as the columns of `X`, and compute the associated observation vector.

```julia
J = ceil(Int64, (tf-t0)/Δtobs)
@inbounds for i=1:Ne
    xtmp = rand(π0)
    for j=1:J
        tspan = (t0 + (j-1)*model.Δtobs, t0 + j*model.Δtobs)
        prob = remake(prob, tspan = tspan, u0 = xtmp)
        sol = solve(prob, Tsit5(), dt = Δtdyn, dense = false, save_everystep = false)
        xtmp .= copy(sol.u[end])
        # Add process noise
        model.ϵx(xtmp)
    end

    # Store the state
    X[Ny+1:Ny+Nx,i] .= copy(xtmp)
end
# Generate the observations
observe(h, X, tf, Ny, Nx)

# Apply the observation noise
ϵy(X, 1, Ny)
```

```julia
# writedlm( "lorenzdata.csv",  X, ',')
# X = readdlm( "lorenzdata.csv",  ',')
```

```julia
# Initialize the map
m = 15
S = HermiteMap(m, X; diag = true, α = 1.0e-6, b = "CstProHermite")
```

    Hermite map of dimension 6:
    Map component of dimension 1 with Nψ = 1 active features
    Map component of dimension 2 with Nψ = 1 active features
    Map component of dimension 3 with Nψ = 1 active features
    Map component of dimension 4 with Nψ = 1 active features
    Map component of dimension 5 with Nψ = 1 active features
    Map component of dimension 6 with Nψ = 1 active features

```julia
# plt = plot(layout = grid(6,1), size = (400, 800))

# for i=1:6
#       histogram!(plt[i,1], sort(transform(X)[i,:]), bins = 100; normalize = :pdf, legend = false)
# end
# plt
# savefig("normalized_histogram.pdf")
```

The optimization works up to 22 active features, after the Cholesky decomposition crashes.

```julia
@time optimize(S, X, "kfolds"; withqr = true, verbose = false, hessprecond = true)
```

    Progress: 100%|█████████████████████████████████████████| Time: 0:00:04

      4.931657 seconds (3.22 M allocations: 6.953 GiB, 4.87% gc time, 2.47% compilation time)

    Hermite map of dimension 6:
    Map component of dimension 1 with Nψ = 7 active features
    Map component of dimension 2 with Nψ = 15 active features
    Map component of dimension 3 with Nψ = 14 active features
    Map component of dimension 4 with Nψ = 15 active features
    Map component of dimension 5 with Nψ = 15 active features
    Map component of dimension 6 with Nψ = 13 active features

```julia
# Conditioning number of the basis evaluation on the different samples
cond(evaluate_basis(S[1].I.f.f,  transform(X[1:1,:])))
```

```julia
stor = Storage(S[1].I.f, transform(X[1:1,:]))
q = QRscaling(stor)
```

    QRscaling([-22.360679774997894 -0.02271836544929029 … -0.4481692087090723 4.339431572165843; 0.0 -22.360668234108555 … -6.983473057759276 0.6699844510861149; … ; 0.0 0.0 … -4.10644447220294 -0.6093826510547249; 0.0 0.0 … 0.0 1.9420001106391365], [-0.0447213595499958 4.543675434959188e-5 … 0.011974935248436213 0.1935157652036135; 0.0 -0.04472138263178639 … -0.19161651915998346 -0.07252363472327122; … ; 0.0 0.0 … -0.24351966933173722 -0.07641434254733631; 0.0 0.0 … 0.0 0.5149330293657334], [0.5264092155084994 0.0 … 0.0 0.0; 0.0 0.38100370855816695 … 0.0 0.0; … ; 0.0 0.0 … 0.255700543071432 0.0; 0.0 0.0 … 0.0 0.2505113994360087], [1.899662791871952 0.0 … 0.0 0.0; 0.0 2.624646368362927 … 0.0 0.0; … ; 0.0 0.0 … 3.9108247013798563 0.0; 0.0 0.0 … 0.0 3.991834312735308], [-11.77086789859341 -0.008655781488559328 … -0.11459711005480372 1.0870770759000647; 0.0 -8.519497523034158 … -1.7856778533937605 0.16783874244194877; … ; 0.0 0.0 … -1.0500200816349714 -0.15265730070774408; 0.0 0.0 … 0.0 0.4864931654210938], [-0.0849555027390544 8.631451162134577e-5 … 0.022748338926530182 0.36761469879793357; 0.0 -0.11737781451268703 … -0.5029256011315958 -0.1903488944969133; … ; 0.0 0.0 … -0.9523627380944125 -0.29884309837382456; 0.0 0.0 … 0.0 2.0555273353828727], [0.00721743744564548 -7.332892728467383e-6 … -0.001932596569981773 -0.031230891550644506; -7.332892728467383e-6 0.013777558789969677 … 0.05903427143507071 0.022374467714145294; … ; -0.001932596569981773 0.05903427143507071 … 2.8723906329846187 1.5420809935725495; -0.031230891550644506 0.022374467714145294 … 1.5420809935725495 23.134650012276875])

```julia
heatmap(log10.(abs.(q.R)), yflip =true, ratio = 1, colorbartitle = L"\log_{10}(R)")
# savefig("R_matrix.pdf")
```

```julia
cond(q.R)
```

    55.27446630792376

```julia
getcoeff(S[3])
```

    14-element Vector{Float64}:
      20.53124818580839
       2.7009792548197797
     -23.727192396955544
       7.708553712290142
       3.4144837968947552
       8.690484275268181
       0.40763706739599176
       1.3189876767571778
       3.4113165048303453
      -0.29754797524571996
      -0.1196049538348405
      -7.221557230952003
      -0.17410702253430407
      -3.0920766713334977

```julia
cond(evaluate_basis(S[6].I.f, X[1:6,:]))
```

    1.1152075461128705e13

```julia
P = zeros(ncoeff(S[1]), ncoeff(S[1]))
precond!(P, getcoeff(S[1]), stor, S[1], X[1:1,:])

Pqr = zeros(ncoeff(S[1]), ncoeff(S[1]))
qrprecond!(Pqr, q.U*getcoeff(S[1]), q, stor, S[1], X[1:1,:])
```

    7×7 Matrix{Float64}:
      0.000887827   9.47872e-5    0.000427699  …  -0.00265218  -0.0241934
      9.47872e-5    0.000480127   6.41453e-5      -0.00591474  -0.00895739
      0.000427699   6.41453e-5    0.000785304     -0.00329108  -0.0231205
      0.000183272   0.000648211   0.000196036     -0.0322071   -0.036709
      0.00440998    0.00107088    0.00329822      -0.0573009   -0.320624
     -0.00265218   -0.00591474   -0.00329108   …   0.331632     0.458023
     -0.0241934    -0.00895739   -0.0231205        0.458023     2.09767

```julia
cond(P)
```

    61.15180729408676

```julia
cond(Pqr)
```

    7756.545093750413

```julia
plot(S)
```

```julia
Ωhat = zeros(Ny + Nx, Ny + Nx)
cache = zeros(1, Ny + Nx, Ny + Nx)
for i=1:Ne
    hess_x_log_pdf!(cache, S, X[:,i:i])
    Ωhat .+= copy(cache[1,:,:]).^2
end
rmul!(Ωhat, 1/Ne)
```

    6×6 Matrix{Float64}:
         1.94563e8     6.93826e6  3955.05       …  41370.9            0.0
         6.93826e6     8.88628e6  1509.89              0.0            0.0
      3955.05       1509.89          2.30398e8         1.71195e5      2.32168e8
         6.79744e7     2.06612e5     1.56237e5         2.57697e7  82169.1
     41370.9           0.0           1.71195e5         1.76886e7      1.67259e5
         0.0           0.0           2.32168e8  …      1.67259e5      4.16718e8

```julia
plt = plot()

heatmap!(plt,log10.(Ωhat), ratio = 1, yflip = true,
        colorbar = true, color = :plasma, clim = (2.0, Inf), colorbar_title = L"\log_{10}(\hat{\Omega})",
        xlabel = "Index", ylabel = "Index", background_color_inside = palette(:plasma)[1],
        yticks = (reverse(collect(0:10:Nx))))
plt
```

The rest is work in progress

Set initial condition of the true system

```julia
L = LinearTransform(10)
L.L - Matrix(1.0*I, 10, 10)
```

    10×10 Matrix{Float64}:
     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  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
     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  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
     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  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
     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

```julia
x0 = rand(model.π0);
```

```julia
data = generate_lorenz63(model, x0, Tf);
```

Initialize the ensemble matrix `X` $\in \mathbb{R}^{(N_y + N_x) \times N_e}$.

```julia
Ne = 500 #ensemble size
X = zeros(model.Ny + model.Nx, Ne)

# Generate the initial conditions for the state.
viewstate(X, model.Ny, model.Nx) .= rand(model.π0, Ne)
```

    3×500 view(::Matrix{Float64}, 4:6, :) with eltype Float64:
     -0.45561   0.0652029  -1.80085   …  -0.470509   0.441291  -1.28427
     -0.987358  0.0830121  -1.23204      -1.65951    0.784583  -0.827193
      0.343795  0.0313563   0.144106      0.289434  -0.673158  -0.540165

Use the stochastic ensemble Kalman filter for the spin-up phase. There is no reason to use the stochastic map filter over the first cycles, as the performance of the inference is determined by the quality of the ensemble, not the quality of the filter.

```julia
enkf = StochEnKF(x->x, ϵy, Δtdyn, Δtobs)
tspin = 10.0
Tspin = ceil(Int64, (tspin-t0)/Δtobs)
```

    100

```julia
Xspin = seqassim(F, data, Tspin, model.ϵx, enkf, deepcopy(X), model.Ny, model.Nx, t0);
```

```julia
nb = 1
ne = Tspin
Δ = 1
plt = plot(layout = grid(3,1), xlim = (-Inf, Inf), ylim = (-Inf, Inf), xlabel = L"t",
           size = (600, 800))

for i =1:3
    plot!(plt[i,1], data.tt[nb:Δ:ne], data.xt[i,nb:Δ:ne], linewidth =  3, color = :teal,
          ylabel = latexstring("x_"*string(i)), legend = (i == 1), label = "True")
    plot!(plt[i,1], data.tt[nb:Δ:ne], mean_hist(Xspin)[i,1+nb:Δ:1+ne], linewidth = 3, grid = false,
          color = :orangered2, linestyle = :dash, label = "sEnKF")
    scatter!(plt[i,1], data.tt[nb:Δ:ne], data.yt[i,nb:Δ:ne], linewidth = 3, color = :grey,
          markersize = 5, alpha = 0.5, label  = "Observation")
end

plt
```

Define the stochastic map filter `smf`, based on the `HermiteMap` `S`.

```julia
X0 = vcat(zeros(Ny, Ne), Xspin[end])

# Add process noise
model.ϵx(X0, Ny+1, Ny+Nx)

# Generate the observations
observe(h, X0, tspin, Ny, Nx)

# Apply the observation noise
ϵy(X0, 1, Ny)
```

```julia
m = 40
S = HermiteMap(m, X0; diag = true, b = "CstLinProHermite")
```

    Hermite map of dimension 6:
    Map component of dimension 1 with Nψ = 1 active features
    Map component of dimension 2 with Nψ = 1 active features
    Map component of dimension 3 with Nψ = 1 active features
    Map component of dimension 4 with Nψ = 1 active features
    Map component of dimension 5 with Nψ = 1 active features
    Map component of dimension 6 with Nψ = 1 active features

```julia
@time optimize(S, X0, "kfolds"; withqr = true, verbose = true, hessprecond = true, start = Ny+1)
```

    1 terms - Training error: 1.4111911175683378, Validation error: 1.4449281957500124
    1 terms - Training error: 1.410623658562414, Validation error: 1.4468478274150656
    2 terms - Training error: 1.40744568459164, Validation error: 1.443157032394458
    3 terms - Training error: 1.373715606981072, Validation error: 1.413859027890204
    4 terms - Training error: 1.372491081672159, Validation error: 1.4104116965162845
    5 terms - Training error: 1.3719059100040596, Validation error: 1.406657443297734
    6 terms - Training error: 1.3716629743071087, Validation error: 1.4039588192763328
    7 terms - Training error: 1.370405716057054, Validation error: 1.4121420774474198
    8 terms - Training error: 1.3669689805412297, Validation error: 1.4259461915630753
    9 terms - Training error: 1.3664403548236859, Validation error: 1.4339933547375567
    10 terms - Training error: 1.2878974227483733, Validation error: 1.3905640241856914
    11 terms - Training error: 1.286500517183202, Validation error: 1.3955864811509817
    12 terms - Training error: 1.2768455386284392, Validation error: 1.3977791837915992
    13 terms - Training error: 1.2764965431900148, Validation error: 1.397356016388681
    14 terms - Training error: 1.2527765114256382, Validation error: 1.377323722822657
    15 terms - Training error: 1.2526121850276069, Validation error: 1.3753019819015841
    16 terms - Training error: 1.251391883080683, Validation error: 1.3870749183585374
    17 terms - Training error: 1.2497366828267065, Validation error: 1.396851816834918
    18 terms - Training error: 1.2453995054464873, Validation error: 1.4526158359259038
    19 terms - Training error: 1.2438350486097063, Validation error: 1.4401663421433117
    20 terms - Training error: 1.2372557695950597, Validation error: 1.416295910972395
    21 terms - Training error: 1.2370098661597502, Validation error: 1.4123824185354001
    22 terms - Training error: 1.2362563107487061, Validation error: 1.4204154278430259
    23 terms - Training error: 1.2355284828099617, Validation error: 1.4037288122185376
    1 terms - Training error: 1.4337153933049587, Validation error: 1.3548310928035259
    1 terms - Training error: 1.4335109770737642, Validation error: 1.3565744358645502
    2 terms - Training error: 1.403671006586731, Validation error: 1.3477352795404016
    3 terms - Training error: 1.3984021991941369, Validation error: 1.3580253628577328
    4 terms - Training error: 1.3977804165898968, Validation error: 1.3565041812051946
    5 terms - Training error: 1.3969066235764114, Validation error: 1.3564724717214534
    6 terms - Training error: 1.3960317108712443, Validation error: 1.367578537541042
    7 terms - Training error: 1.3331301423628623, Validation error: 1.2430036230340247
    8 terms - Training error: 1.3315462774938172, Validation error: 1.2360162137407686
    9 terms - Training error: 1.3299622557373636, Validation error: 1.2386671040196853
    10 terms - Training error: 1.3278365401469638, Validation error: 1.2434120608162393
    11 terms - Training error: 1.327315096574128, Validation error: 1.244767980803149
    12 terms - Training error: 1.3260509966161689, Validation error: 1.2422941902479392
    13 terms - Training error: 1.3256853822408172, Validation error: 1.24579594345804
    14 terms - Training error: 1.3055275227407606, Validation error: 1.2136562357216074
    15 terms - Training error: 1.3053107733683325, Validation error: 1.2157469118484514
    16 terms - Training error: 1.304635281488049, Validation error: 1.2135982591224506
    17 terms - Training error: 1.3026021280186466, Validation error: 1.208022365336603
    18 terms - Training error: 1.3022896127163182, Validation error: 1.210989451021638
    19 terms - Training error: 1.3019169174724354, Validation error: 1.2041344755742729
    20 terms - Training error: 1.3003321260589529, Validation error: 1.2107185549955768
    21 terms - Training error: 1.3002170479814645, Validation error: 1.2115765975036166
    22 terms - Training error: 1.2937475281040482, Validation error: 1.1823133317732546
    23 terms - Training error: 1.2925166132180632, Validation error: 1.1871536576838408
    1 terms - Training error: 1.4002861805625233, Validation error: 1.4885479437732687
    1 terms - Training error: 1.4000187976963008, Validation error: 1.4902379643858659
    2 terms - Training error: 1.3776135203373345, Validation error: 1.4193118777988651
    3 terms - Training error: 1.3736201598582958, Validation error: 1.4182240101028294
    4 terms - Training error: 1.3720768033777526, Validation error: 1.4198540255388625
    5 terms - Training error: 1.3715536667689034, Validation error: 1.4155076720195108
    6 terms - Training error: 1.3697246320278769, Validation error: 1.4161120112231627
    7 terms - Training error: 1.3690081784463546, Validation error: 1.4197911632477938
    8 terms - Training error: 1.3686942718914978, Validation error: 1.4196240164719567
    9 terms - Training error: 1.308634699158249, Validation error: 1.319593962845296
    10 terms - Training error: 1.3085412775701728, Validation error: 1.3198620523394404
    11 terms - Training error: 1.3017550146327395, Validation error: 1.3124244870089907
    12 terms - Training error: 1.3001793694707469, Validation error: 1.3243691738072327
    13 terms - Training error: 1.2992380550075422, Validation error: 1.3229082255232134
    14 terms - Training error: 1.2985710919904914, Validation error: 1.3183789271548063
    15 terms - Training error: 1.2983337459313153, Validation error: 1.3180583051833934
    16 terms - Training error: 1.2980275453853263, Validation error: 1.3211513884427277
    17 terms - Training error: 1.295379255022312, Validation error: 1.3096031160687465
    18 terms - Training error: 1.2940662302436576, Validation error: 1.3077931611056552
    19 terms - Training error: 1.2927285949447338, Validation error: 1.3092485713186712
    20 terms - Training error: 1.2903499579258255, Validation error: 1.3028961494058844
    21 terms - Training error: 1.290329307923711, Validation error: 1.3037033121005868
    22 terms - Training error: 1.2834091117809339, Validation error: 1.3144154352777926
    23 terms - Training error: 1.2825471103813446, Validation error: 1.3166274672907903
    1 terms - Training error: 1.409920346404197, Validation error: 1.4500112804065741
    1 terms - Training error: 1.409264264504158, Validation error: 1.452479418423353
    2 terms - Training error: 1.404930170068828, Validation error: 1.4632062468539169
    3 terms - Training error: 1.3779903927200077, Validation error: 1.4397795267425002
    4 terms - Training error: 1.375347502449689, Validation error: 1.4487646148249804
    5 terms - Training error: 1.3742313889958826, Validation error: 1.449662746367475
    6 terms - Training error: 1.3734340934476668, Validation error: 1.448017102416126
    7 terms - Training error: 1.3397853339538413, Validation error: 1.4269227155088842
    8 terms - Training error: 1.3393743955992006, Validation error: 1.4236824371307295
    9 terms - Training error: 1.3384952491328233, Validation error: 1.4223560729897555
    10 terms - Training error: 1.3378745688598361, Validation error: 1.4249495090205044
    11 terms - Training error: 1.33755470224358, Validation error: 1.4250959420644862
    12 terms - Training error: 1.2548120806519063, Validation error: 1.4250936929974998
    13 terms - Training error: 1.253629463896232, Validation error: 1.4364709989540796
    14 terms - Training error: 1.2461124241900838, Validation error: 1.4133574727673226
    15 terms - Training error: 1.245358784533442, Validation error: 1.4134583471776316
    16 terms - Training error: 1.244345679880001, Validation error: 1.4131964209294636
    17 terms - Training error: 1.2425838118017574, Validation error: 1.4149681558027938
    18 terms - Training error: 1.2403548630097179, Validation error: 1.409003777142902
    19 terms - Training error: 1.2399373162855165, Validation error: 1.4141544802672903
    20 terms - Training error: 1.2331413991104678, Validation error: 1.4394590049896148
    21 terms - Training error: 1.232184188945378, Validation error: 1.4414759865696467
    22 terms - Training error: 1.2270943402658785, Validation error: 1.4201376076631476
    23 terms - Training error: 1.2240103550595651, Validation error: 1.4311139009540215
    1 terms - Training error: 1.4345796281833452, Validation error: 1.3513741532899815
    1 terms - Training error: 1.434509230834942, Validation error: 1.350883426071751
    2 terms - Training error: 1.4316139702952093, Validation error: 1.3461884103785424
    3 terms - Training error: 1.3916608018170487, Validation error: 1.3438229699187603
    4 terms - Training error: 1.390961718972958, Validation error: 1.3417061792859426
    5 terms - Training error: 1.3900623037629871, Validation error: 1.3378863683881097
    6 terms - Training error: 1.388846907630216, Validation error: 1.3359216449433706
    7 terms - Training error: 1.3885377726221588, Validation error: 1.3360129220701793
    8 terms - Training error: 1.321646194474665, Validation error: 1.263762489319266
    9 terms - Training error: 1.3214848770547019, Validation error: 1.2625556643811973
    10 terms - Training error: 1.3213546592864753, Validation error: 1.2628400520915166
    11 terms - Training error: 1.319364939154128, Validation error: 1.253164079081693
    12 terms - Training error: 1.3190004867884406, Validation error: 1.252338145213979
    13 terms - Training error: 1.3152168954297874, Validation error: 1.2653297590343984
    14 terms - Training error: 1.3151920518571374, Validation error: 1.2651058089744183
    15 terms - Training error: 1.313763780556504, Validation error: 1.2690771332199
    16 terms - Training error: 1.3136739363821093, Validation error: 1.270211578966555
    17 terms - Training error: 1.3091867611296357, Validation error: 1.2725646237810424
    18 terms - Training error: 1.3088195875036084, Validation error: 1.271200466095293
    19 terms - Training error: 1.3070333971476693, Validation error: 1.2700615807626312
    20 terms - Training error: 1.306400220653484, Validation error: 1.2707225244389218
    21 terms - Training error: 1.3061568287985827, Validation error: 1.2705423798596631
    22 terms - Training error: 1.306127153293236, Validation error: 1.2710963860555122
    23 terms - Training error: 1.2944301307109531, Validation error: 1.2809939974421458
    1 terms - Training error: 1.4179385332046728
    1 terms - Training error: 1.417759194987984
    2 terms - Training error: 1.414485787236862
    3 terms - Training error: 1.3816383810977644
    4 terms - Training error: 1.3806134708569004
    5 terms - Training error: 1.3789744222463491
    6 terms - Training error: 1.3775335593644533
    7 terms - Training error: 1.377308481850296
    8 terms - Training error: 1.3095286770425674
    9 terms - Training error: 1.3091652088794257
    10 terms - Training error: 1.2887929987018334
    11 terms - Training error: 1.288554768734656
    12 terms - Training error: 1.2880399031161014

    Progress:  33%|█████████████▋                           |  ETA: 0:00:05

    13 terms - Training error: 1.2879926678990545
    14 terms - Training error: 1.2833639515438278
    15 terms - Training error: 1.2829706610879452
    1 terms - Training error: 1.433215456023802, Validation error: 1.3568308419281578
    1 terms - Training error: 1.4330367519048595, Validation error: 1.356776173868118
    2 terms - Training error: 1.4207366335295133, Validation error: 1.3547481430251966
    3 terms - Training error: 1.41809150286688, Validation error: 1.3541954847369775
    4 terms - Training error: 1.4157077575993537, Validation error: 1.3447192307541793
    5 terms - Training error: 1.2070229217218058, Validation error: 1.2114847945477996
    6 terms - Training error: 1.2017067850399292, Validation error: 1.2127661169890336
    7 terms - Training error: 1.1834104002251777, Validation error: 1.1798554538497947
    8 terms - Training error: 1.1793619770128023, Validation error: 1.1887515133571638
    9 terms - Training error: 1.173265787823672, Validation error: 1.194940769292195
    10 terms - Training error: 1.1731632311643339, Validation error: 1.195119575137552
    11 terms - Training error: 1.1729969499591932, Validation error: 1.1948995496539299
    12 terms - Training error: 1.1727713622280918, Validation error: 1.1936848934224111
    13 terms - Training error: 1.1144079756210663, Validation error: 1.1227017454104458
    14 terms - Training error: 1.1141407567010537, Validation error: 1.1224464378404397
    15 terms - Training error: 1.0956271367952468, Validation error: 1.1143930239018918
    16 terms - Training error: 1.0887621470568407, Validation error: 1.1334116666800842
    17 terms - Training error: 1.083511369593942, Validation error: 1.1147474639552668
    18 terms - Training error: 1.0834519782819878, Validation error: 1.112979329304202
    19 terms - Training error: 1.0817074627923715, Validation error: 1.119760742981184
    20 terms - Training error: 1.0808018082717237, Validation error: 1.1208687760769567
    21 terms - Training error: 1.0805095880533393, Validation error: 1.120652944419238
    22 terms - Training error: 1.048411695560202, Validation error: 1.1140402988640594
    23 terms - Training error: 1.0442195871689588, Validation error: 1.1550056960795236
    1 terms - Training error: 1.4335655085994128, Validation error: 1.3554306316257152
    1 terms - Training error: 1.4330286757218045, Validation error: 1.3569444234559873
    2 terms - Training error: 1.2148207294706128, Validation error: 1.1964680369705964
    3 terms - Training error: 1.2095161394605378, Validation error: 1.2027056451065399
    4 terms - Training error: 1.2090557462064582, Validation error: 1.2017939516691858
    5 terms - Training error: 1.1830179213666006, Validation error: 1.1889487025296244
    6 terms - Training error: 1.182225773668977, Validation error: 1.1904215576534045
    7 terms - Training error: 1.1796914301142662, Validation error: 1.184631773836422
    8 terms - Training error: 1.1767003679806085, Validation error: 1.1762810849773317
    9 terms - Training error: 1.176230741748542, Validation error: 1.1765550718093334
    10 terms - Training error: 1.1159194624909219, Validation error: 1.1156774953484576
    11 terms - Training error: 1.1150459748920225, Validation error: 1.11716267673766
    12 terms - Training error: 1.092402064725167, Validation error: 1.0966737184707107
    13 terms - Training error: 1.091369192516918, Validation error: 1.097592358099691
    14 terms - Training error: 1.0905279128349503, Validation error: 1.1037086670722243
    15 terms - Training error: 1.0886572370944256, Validation error: 1.0987963215734564
    16 terms - Training error: 1.0831502862928162, Validation error: 1.1097270068554375
    17 terms - Training error: 1.0830501741582694, Validation error: 1.111624274809518
    18 terms - Training error: 1.081862695697158, Validation error: 1.113173668756336
    19 terms - Training error: 1.0815684690409912, Validation error: 1.112234243533313
    20 terms - Training error: 1.0809356582535679, Validation error: 1.1123215125922608
    21 terms - Training error: 1.080566552257767, Validation error: 1.109576190444984
    22 terms - Training error: 1.07786568120621, Validation error: 1.104864297548806
    23 terms - Training error: 1.0776206523095078, Validation error: 1.1031755154088454
    1 terms - Training error: 1.4132201668695255, Validation error: 1.436811998545264
    1 terms - Training error: 1.4124435364288705, Validation error: 1.4403239632802831
    2 terms - Training error: 1.407498024907391, Validation error: 1.4130968559530832
    3 terms - Training error: 1.401131034920391, Validation error: 1.406398480075519
    4 terms - Training error: 1.2070499444659677, Validation error: 1.2115121486748395
    5 terms - Training error: 1.2065298423429887, Validation error: 1.2132714384783343
    6 terms - Training error: 1.1949157117916063, Validation error: 1.2056443349932127
    7 terms - Training error: 1.1938466617588939, Validation error: 1.201653063541312
    8 terms - Training error: 1.1766318416096258, Validation error: 1.158936743413208
    9 terms - Training error: 1.1722359459414562, Validation error: 1.1537862430789327
    10 terms - Training error: 1.1658689787409742, Validation error: 1.1504524381032974
    11 terms - Training error: 1.1653948108699987, Validation error: 1.1495523904054024
    12 terms - Training error: 1.097489340179023, Validation error: 1.079445874403971
    13 terms - Training error: 1.0969924087561251, Validation error: 1.0800108255241436
    14 terms - Training error: 1.0968761244869942, Validation error: 1.0790754679438943
    15 terms - Training error: 1.096382572903296, Validation error: 1.0776355893157556
    16 terms - Training error: 1.0742459577146355, Validation error: 1.030780619808358
    17 terms - Training error: 1.0741415841299258, Validation error: 1.0304773143983914
    18 terms - Training error: 1.07236875032891, Validation error: 1.03585298018868
    19 terms - Training error: 1.067434646391999, Validation error: 1.0457073048684764
    20 terms - Training error: 1.066138099964516, Validation error: 1.0507412486782366
    21 terms - Training error: 1.0646154481873458, Validation error: 1.0542082027501525
    22 terms - Training error: 1.0637325994720126, Validation error: 1.0514017976550878
    23 terms - Training error: 1.0554583919979432, Validation error: 1.0775042047856014
    1 terms - Training error: 1.3933875040410817, Validation error: 1.5161426498590396
    1 terms - Training error: 1.3925534210941761, Validation error: 1.5194848919969701
    2 terms - Training error: 1.2076866060969302, Validation error: 1.2258437761410426
    3 terms - Training error: 1.2046190840467699, Validation error: 1.2217698171255211
    4 terms - Training error: 1.2041031215662363, Validation error: 1.2209306398623547
    5 terms - Training error: 1.177633313532716, Validation error: 1.2107807970984963
    6 terms - Training error: 1.1737665949025937, Validation error: 1.2113299390088728
    7 terms - Training error: 1.1735069364746313, Validation error: 1.2099625682960078
    8 terms - Training error: 1.170344791600815, Validation error: 1.202187875829082
    9 terms - Training error: 1.1690527830768027, Validation error: 1.212220423044877
    10 terms - Training error: 1.160275775805161, Validation error: 1.19584121976917
    11 terms - Training error: 1.158270673806344, Validation error: 1.2059041790276315
    12 terms - Training error: 1.084056283300328, Validation error: 1.1710225438191337
    13 terms - Training error: 1.082499345155625, Validation error: 1.1748051366802073
    14 terms - Training error: 1.0545484177128808, Validation error: 1.1633845769063238
    15 terms - Training error: 1.0539814834270254, Validation error: 1.1655498431647842
    16 terms - Training error: 1.0529962649043272, Validation error: 1.1657393119945818
    17 terms - Training error: 1.049963131974027, Validation error: 1.181430375156636
    18 terms - Training error: 1.04971584142601, Validation error: 1.1800151353177182
    19 terms - Training error: 1.0490964773998122, Validation error: 1.1692499423416147
    20 terms - Training error: 1.0419479929354036, Validation error: 1.1414206948434074
    21 terms - Training error: 1.0388040094521132, Validation error: 1.149755462470827
    22 terms - Training error: 1.037107885306996, Validation error: 1.162963287157982
    23 terms - Training error: 1.0355156896457578, Validation error: 1.1691093860325799
    1 terms - Training error: 1.416304030489544, Validation error: 1.4244765440651899
    1 terms - Training error: 1.4160454140213155, Validation error: 1.4248042254558078
    2 terms - Training error: 1.401501956086545, Validation error: 1.4330064696571432
    3 terms - Training error: 1.3987922061495146, Validation error: 1.4326183416347629
    4 terms - Training error: 1.394854536364264, Validation error: 1.4450239467694674
    5 terms - Training error: 1.3929751031321231, Validation error: 1.434898542605503
    6 terms - Training error: 1.206407116252098, Validation error: 1.2154475480721814
    7 terms - Training error: 1.2043400542272973, Validation error: 1.2050656505505393
    8 terms - Training error: 1.1806756024786833, Validation error: 1.194696167988951
    9 terms - Training error: 1.1780244966197702, Validation error: 1.1951209417879756
    10 terms - Training error: 1.1737969300104656, Validation error: 1.1892821385735604
    11 terms - Training error: 1.1735933688635294, Validation error: 1.189737181809795
    12 terms - Training error: 1.1733867787685397, Validation error: 1.1886570198349875
    13 terms - Training error: 1.172947431672283, Validation error: 1.187888941001365
    14 terms - Training error: 1.171868167012945, Validation error: 1.192258986886634
    15 terms - Training error: 1.1717864510282823, Validation error: 1.1945513729924362
    16 terms - Training error: 1.109456764791559, Validation error: 1.140223638041706
    17 terms - Training error: 1.1091721258763898, Validation error: 1.1380867510876165
    18 terms - Training error: 1.0836608741348905, Validation error: 1.1278394501085647
    19 terms - Training error: 1.0816223126659084, Validation error: 1.1250549383862933
    20 terms - Training error: 1.0791778350003005, Validation error: 1.1199741910775882
    21 terms - Training error: 1.0776912613621668, Validation error: 1.1188313613740364
    22 terms - Training error: 1.0765349664422275, Validation error: 1.1227894220265888
    23 terms - Training error: 1.07608889542357, Validation error: 1.1205271671597496
    1 terms - Training error: 1.4179385332046734
    1 terms - Training error: 1.4177327402159592
    2 terms - Training error: 1.2109478031488603
    3 terms - Training error: 1.207693409210143
    4 terms - Training error: 1.207106271685966
    5 terms - Training error: 1.183853946637566
    6 terms - Training error: 1.1834535616005797
    7 terms - Training error: 1.180278593947676
    8 terms - Training error: 1.176073573752118
    9 terms - Training error: 1.1757458935562348
    10 terms - Training error: 1.1150219322049622
    11 terms - Training error: 1.114811889119215
    12 terms - Training error: 1.0983832288808484
    13 terms - Training error: 1.0958129525438955
    14 terms - Training error: 1.0879320986915826
    15 terms - Training error: 1.087101565196487
    16 terms - Training error: 1.086898353344848
    17 terms - Training error: 1.0867324216220517
    18 terms - Training error: 1.0858551644680814
    19 terms - Training error: 1.085356684683799
    20 terms - Training error: 1.0846135769486138

    Progress:  67%|███████████████████████████▍             |  ETA: 0:00:02

    21 terms - Training error: 1.0842214533500212
    1 terms - Training error: 1.4269957599178684, Validation error: 1.3817096263518909
    1 terms - Training error: 1.426115678927987, Validation error: 1.3882811246822826
    2 terms - Training error: 1.4094697306476276, Validation error: 1.3719322223819033
    3 terms - Training error: 1.404358477630845, Validation error: 1.3888097307285652
    4 terms - Training error: 1.3900556987317758, Validation error: 1.3542196185832056
    5 terms - Training error: 1.3884383582052091, Validation error: 1.3377142048559103
    6 terms - Training error: 1.1395874642603554, Validation error: 1.2315118615654816
    7 terms - Training error: 1.1384315557201254, Validation error: 1.2316516770437356
    8 terms - Training error: 1.137349203375671, Validation error: 1.2369755377315612
    9 terms - Training error: 1.1345829379642514, Validation error: 1.2282868373597564
    10 terms - Training error: 1.104418217815836, Validation error: 1.151673865702405
    11 terms - Training error: 1.10427376617855, Validation error: 1.1461140722776557
    12 terms - Training error: 1.1021355291941588, Validation error: 1.139329504224689
    13 terms - Training error: 1.101953555219043, Validation error: 1.1355948676198664
    14 terms - Training error: 1.1013938404638302, Validation error: 1.1355542184022611
    15 terms - Training error: 1.1012122391192962, Validation error: 1.1385609892918411
    16 terms - Training error: 1.1005393684289821, Validation error: 1.149393905407074
    17 terms - Training error: 1.1004809015699948, Validation error: 1.1480908134973777
    18 terms - Training error: 1.1003569863193916, Validation error: 1.1464947120996904
    19 terms - Training error: 1.0994066202897297, Validation error: 1.1438950098639222
    20 terms - Training error: 1.0988789719326983, Validation error: 1.1447290070946827
    21 terms - Training error: 1.0988414781725129, Validation error: 1.1446094579544832
    22 terms - Training error: 0.9752883293786028, Validation error: 1.049101145816564
    23 terms - Training error: 0.9734658311828948, Validation error: 1.063705157890561
    1 terms - Training error: 1.4093811940478076, Validation error: 1.4521678898321337
    1 terms - Training error: 1.4081885166647596, Validation error: 1.4563254484725352
    2 terms - Training error: 1.1618550718767402, Validation error: 1.1578508912446193
    3 terms - Training error: 1.1581145896891611, Validation error: 1.1530653102641548
    4 terms - Training error: 1.1574997241135956, Validation error: 1.1514482279187335
    5 terms - Training error: 1.156678542122803, Validation error: 1.1503349074185598
    6 terms - Training error: 1.1522817429049097, Validation error: 1.148842745696791
    7 terms - Training error: 1.1108909282645607, Validation error: 1.1203981628224624
    8 terms - Training error: 1.1105084507029312, Validation error: 1.1190857199877673
    9 terms - Training error: 1.1077206326435636, Validation error: 1.121389456496304
    10 terms - Training error: 1.1054245364569877, Validation error: 1.1152629886646928
    11 terms - Training error: 1.1044628293702994, Validation error: 1.1213548838690373
    12 terms - Training error: 1.0922741200579575, Validation error: 1.1011597289678636
    13 terms - Training error: 1.0916120406703658, Validation error: 1.1037607852231628
    14 terms - Training error: 0.9596001192469163, Validation error: 0.9525950165572871
    15 terms - Training error: 0.9589683574657246, Validation error: 0.9519711966198483
    16 terms - Training error: 0.9588237144122813, Validation error: 0.9530644183954915
    17 terms - Training error: 0.9577991543297883, Validation error: 0.9577327067238053
    18 terms - Training error: 0.9517436124864074, Validation error: 0.9488900905675269
    19 terms - Training error: 0.950664450090313, Validation error: 0.9443818740011674
    20 terms - Training error: 0.9501845333426282, Validation error: 0.9433935158528954
    21 terms - Training error: 0.9494697588964659, Validation error: 0.9417180523343956
    22 terms - Training error: 0.9493587518350766, Validation error: 0.9438788978487865
    23 terms - Training error: 0.8933539472937233, Validation error: 0.921946638630939
    1 terms - Training error: 1.4060875759777958, Validation error: 1.4653423621121824
    1 terms - Training error: 1.4054754048948743, Validation error: 1.4681800831360263
    2 terms - Training error: 1.3901156497364144, Validation error: 1.4265314501109392
    3 terms - Training error: 1.3792182598859573, Validation error: 1.4318219405775985
    4 terms - Training error: 1.1531693046010458, Validation error: 1.1678770575844948
    5 terms - Training error: 1.1524645610778044, Validation error: 1.1662742255160032
    6 terms - Training error: 1.149789017252491, Validation error: 1.158333377862908
    7 terms - Training error: 1.11090707810409, Validation error: 1.1165756370080835
    8 terms - Training error: 1.110563815027309, Validation error: 1.1150748989795356
    9 terms - Training error: 1.1074562155018752, Validation error: 1.1179095014781804
    10 terms - Training error: 1.1052923076702836, Validation error: 1.112850258852292
    11 terms - Training error: 1.1033510457288502, Validation error: 1.114704991266598
    12 terms - Training error: 1.102852654526372, Validation error: 1.114084432923744
    13 terms - Training error: 1.1027355744158798, Validation error: 1.1168701337189777
    14 terms - Training error: 1.1026000595127226, Validation error: 1.1179392267505577
    15 terms - Training error: 1.1025051094737777, Validation error: 1.1174920873570264
    16 terms - Training error: 1.1013174349422952, Validation error: 1.1227235897874168
    17 terms - Training error: 1.101212966976853, Validation error: 1.1225091354447037
    18 terms - Training error: 1.0988938495162017, Validation error: 1.1151776818043908
    19 terms - Training error: 1.098771304804868, Validation error: 1.116492178496776
    20 terms - Training error: 1.0985916031418632, Validation error: 1.1161772234730478
    21 terms - Training error: 1.0983848032038135, Validation error: 1.1164842920118803
    22 terms - Training error: 1.0962199995303141, Validation error: 1.1237105071013915
    23 terms - Training error: 1.0948830218319194, Validation error: 1.1244616088396011
    1 terms - Training error: 1.4133474606561875, Validation error: 1.4363028233986148
    1 terms - Training error: 1.4130421173336356, Validation error: 1.4371706010586545
    2 terms - Training error: 1.3884796154320718, Validation error: 1.4310879039597757
    3 terms - Training error: 1.3790224484579383, Validation error: 1.4310064758161618
    4 terms - Training error: 1.1606157779961397, Validation error: 1.1405656968890678
    5 terms - Training error: 1.1581311247307862, Validation error: 1.1481475729579118
    6 terms - Training error: 1.151026558616187, Validation error: 1.1604940729378483
    7 terms - Training error: 1.1130576139072115, Validation error: 1.1128337123470888
    8 terms - Training error: 1.111359552147817, Validation error: 1.118333922585575
    9 terms - Training error: 1.1081839573031296, Validation error: 1.1197134042193142
    10 terms - Training error: 1.1052821365974332, Validation error: 1.1161255951159126
    11 terms - Training error: 1.10463626252099, Validation error: 1.11719005367015
    12 terms - Training error: 1.1045869416587677, Validation error: 1.117517012065364
    13 terms - Training error: 1.1040127597381326, Validation error: 1.1197534121097956
    14 terms - Training error: 0.985544042455513, Validation error: 0.999831664861439
    15 terms - Training error: 0.9852953050749164, Validation error: 0.9999479492379981
    16 terms - Training error: 0.9516820332977646, Validation error: 0.9778945526364906
    17 terms - Training error: 0.9504960126632095, Validation error: 0.9770548262278013
    18 terms - Training error: 0.9502796320871405, Validation error: 0.9777842632986279
    19 terms - Training error: 0.9453490257556204, Validation error: 0.9750039393068132
    20 terms - Training error: 0.9429976923438735, Validation error: 0.9917323583795931
    21 terms - Training error: 0.9421587044433967, Validation error: 0.9922260598442291
    22 terms - Training error: 0.941576774007973, Validation error: 0.9925430891484164
    23 terms - Training error: 0.9385244131652936, Validation error: 0.9901376984005879
    1 terms - Training error: 1.433880675423706, Validation error: 1.3541699643285432
    1 terms - Training error: 1.4334842098067089, Validation error: 1.3555980945241495
    2 terms - Training error: 1.4010773782559391, Validation error: 1.3859656307741366
    3 terms - Training error: 1.3958093603001986, Validation error: 1.3689988664472847
    4 terms - Training error: 1.1653117082238609, Validation error: 1.1218574998893862
    5 terms - Training error: 1.1642406891024133, Validation error: 1.121989433470495
    6 terms - Training error: 1.161647734206918, Validation error: 1.1139842842639613
    7 terms - Training error: 1.1136961605922917, Validation error: 1.1088641614076957
    8 terms - Training error: 1.1132232181654442, Validation error: 1.1111694195302764
    9 terms - Training error: 1.1090266272644973, Validation error: 1.109329114798745
    10 terms - Training error: 1.108869562339202, Validation error: 1.1071858543774202
    11 terms - Training error: 1.107244924546652, Validation error: 1.1061960658321854
    12 terms - Training error: 1.1036233302773277, Validation error: 1.1077228524978182
    13 terms - Training error: 1.1024090402663567, Validation error: 1.1103682766980851
    14 terms - Training error: 1.102236853365616, Validation error: 1.110383838834608
    15 terms - Training error: 1.1022105802189572, Validation error: 1.109991844101126
    16 terms - Training error: 1.10148388201959, Validation error: 1.1141796218769138
    17 terms - Training error: 1.1013408226756087, Validation error: 1.112845649481832
    18 terms - Training error: 1.1005381644997538, Validation error: 1.112507967471846
    19 terms - Training error: 1.1002849740506988, Validation error: 1.1109913944376923
    20 terms - Training error: 1.100210997334037, Validation error: 1.1123106601067225
    21 terms - Training error: 1.0877434495957037, Validation error: 1.096957345499113
    22 terms - Training error: 1.0876468697744555, Validation error: 1.0978599411023184
    23 terms - Training error: 1.0874339322514688, Validation error: 1.0964327280889314
    1 terms - Training error: 1.4179385332046732
    1 terms - Training error: 1.4175741415636134
    2 terms - Training error: 1.1607238583563977
    3 terms - Training error: 1.1567470008998249
    4 terms - Training error: 1.1559342276204176
    5 terms - Training error: 1.155029826653794
    6 terms - Training error: 1.1511971440899542
    7 terms - Training error: 1.1119804471167507
    8 terms - Training error: 1.1113543776804857
    9 terms - Training error: 1.1093435078490443
    10 terms - Training error: 1.1066051431111779
    11 terms - Training error: 1.1062600963070934
    12 terms - Training error: 1.1051909901630292
    13 terms - Training error: 1.1045798475459598
    14 terms - Training error: 1.1044907725144792
    15 terms - Training error: 1.1013218041734065
    16 terms - Training error: 1.1012374194089298
    17 terms - Training error: 1.1011354400579394
    18 terms - Training error: 0.9763309570633508
    19 terms - Training error: 0.9757833819405233
    20 terms - Training error: 0.9757060320214088
    21 terms - Training error: 0.9453351886883915
    22 terms - Training error: 0.9451713377296139

    Progress: 100%|█████████████████████████████████████████| Time: 0:00:06

    23 terms - Training error: 0.9417235190591425
    24 terms - Training error: 0.9408816817449349
      6.098357 seconds (1.99 M allocations: 6.076 GiB, 4.47% gc time)

    Hermite map of dimension 6:
    Map component of dimension 1 with Nψ = 1 active features
    Map component of dimension 2 with Nψ = 1 active features
    Map component of dimension 3 with Nψ = 1 active features
    Map component of dimension 4 with Nψ = 15 active features
    Map component of dimension 5 with Nψ = 21 active features
    Map component of dimension 6 with Nψ = 24 active features

```julia
plot(S[6].I.f)
```

```julia
smf = StochMapFilter(x->x, ϵy, S, Ny, Nx, Δtdyn, Δtobs, 2*Δtobs, false)
```

    Stochastic Map Filter with filtered = false

```julia
mean(evaluate(S, X0; start = Ny+1); dims = 2)
```

    6×1 Matrix{Float64}:
      0.0
      0.0
      0.0
      0.0003147998678031647
     -0.0005704044283433507
     -0.00020692477228671292

```julia
cov(evaluate(S, X0; start = Ny+1); dims = 2)
```

    6×6 Matrix{Float64}:
     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   0.0           0.0           0.0
     0.0  0.0  0.0   0.997445     -0.000370127  -0.00195316
     0.0  0.0  0.0  -0.000370127   0.826134     -0.0214816
     0.0  0.0  0.0  -0.00195316   -0.0214816     0.770837

```julia
Sxpush = evaluate(S, X0; start = Ny+1)[Ny+1:end,:]
```

    3×500 Matrix{Float64}:
     -0.665326  -1.16209  -1.26467     …  0.398589  -0.465817   1.44996
     -0.829003  -1.18318  -0.605887       1.72827    0.978253   0.831009
      0.495277  -2.58625   0.00626933     1.05973   -1.63564   -0.430096

```julia
plt = plot(layout = grid(1,3), xlim = (-3, 3), ylim = (-3,3), legend = false)
    scatter!(plt[1,1], Sxpush[1,:], Sxpush[2,:], ratio = 1)
    scatter!(plt[1,2], Sxpush[2,:], Sxpush[3,:], ratio = 1)
    scatter!(plt[1,3], Sxpush[1,:], Sxpush[3,:], ratio = 1)
plt
```

```julia
plot(S)
```

```julia
plot(S.C[6].I.f)
```

```julia
S = HermiteMap()
```

```julia
smf.M[6].I.f.B
```

    MultiBasis(Basis of 40 functions: Constant, 0th -> 38th degree Probabilistic Hermite function
    , 6)

```julia
Xsmf = seqassim(F, data, 10, model.ϵx, smf, deepcopy(vcat(zeros(Ny, Ne), Xspin[end])), model.Ny, model.Nx, tspin)
```

```julia
nb = 1
ne = Tspin+10
Δ = 1
plt = plot(layout = grid(3,1), xlim = (-Inf, Inf), ylim = (-Inf, Inf), xlabel = L"t",
           size = (600, 800))

for i =1:3
    plot!(plt[i,1], data.tt[nb:Δ:ne], data.xt[i,nb:Δ:ne], linewidth =  3, color = :teal,
          ylabel = latexstring("x_"*string(i)), legend = false, label = "True")
    plot!(plt[i,1], data.tt[nb:Δ:ne], mean_hist(vcat(Xspin, Xsmf[2:end]))[i,1+nb:Δ:1+ne], linewidth = 3, grid = false,
          color = :orangered2, linestyle = :dash, label = "sEnKF")
    scatter!(plt[i,1], data.tt[nb:Δ:ne], data.yt[i,nb:Δ:ne], linewidth = 3, color = :grey,
          markersize = 5, alpha = 0.5, label  = "Observation")
end

plt

# savefig("lorenz63_kfoldoptim.pdf")
```

```julia
nb = 1
ne = Tspin+100
Δ = 1
plt = plot(layout = grid(3,1), xlim = (-Inf, Inf), ylim = (-Inf, Inf), xlabel = L"t",
           size = (600, 800))

for i =1:3
    plot!(plt[i,1], data.tt[nb:Δ:ne], data.xt[i,nb:Δ:ne], linewidth =  3, color = :teal,
          ylabel = latexstring("x_"*string(i)), legend = :bottomleft, label = "True")
    plot!(plt[i,1], data.tt[nb:Δ:ne], mean_hist(vcat(Xspin, Xsmf[2:end]))[i,1+nb:Δ:1+ne], linewidth = 3, grid = false,
          color = :orangered2, linestyle = :dash, label = "sEnKF")
    scatter!(plt[i,1], data.tt[nb:Δ:ne], data.yt[i,nb:Δ:ne], linewidth = 3, color = :grey,
          markersize = 5, alpha = 0.5, label  = "Observation")
end

plt

# savefig("lorenz63.pdf")
```

```julia
vcat(Xspin, Xsmf)
```

```julia
# Plot the first component of the state over time
nb = 1
ne = 499
Δ = 1
plt = plot(xlim = (-Inf, Inf), ylim = (-Inf, Inf), xlabel = L"t", ylabel = L"x_1")
plot!(plt, data.tt[nb:Δ:ne], data.xt[1,nb:Δ:ne], linewidth =  3, color = :teal, label = "True")
plot!(plt, data.tt[nb:Δ:ne], mean_hist(Xsmf)[1,1+nb:Δ:1+ne], linewidth = 3, grid = false,
     color = :orangered2, linestyle = :dash, label = "Stochastic Map Filter")
scatter!(plt, data.tt[nb:Δ:ne], data.yt[1,nb:Δ:ne], linewidth = 3, color = :grey, markersize = 5, alpha = 0.5, label = "Observation")
plt
```

```julia
Xsmf
```

```julia
# Plot the different component of the state over time
nb = 5500
ne = 6000
Δ = 1
plt = plot(layout = grid(3,1), xlim = (-Inf, Inf), ylim = (-Inf, Inf), xlabel = L"t",
           size = (900, 1000))

for i =1:3
    plot!(plt[i,1], data.tt[nb:Δ:ne], data.xt[i,nb:Δ:ne], linewidth =  2, color = :teal,
          ylabel = latexstring("x_"*string(i)), legend = (i == 1), label = "True")
    plot!(plt[i,1], data.tt[nb:Δ:ne], mean_hist(Xsmf)[i,1+nb:Δ:1+ne], linewidth = 2, grid = false,
          color = :orangered2, linestyle = :dash, label = "sEnKF")
    scatter!(plt[i,1], data.tt[nb:Δ:ne], data.yt[i,nb:Δ:ne], linewidth = 3, color = :grey,
          markersize = 5, alpha = 0.5, label  = "Observation")
end

plt
```

Talk about the benchmark of the different filters

### Benchmark of the EnKF filter

```julia
path = pwd()*"/lorenz63/data/"
Nx = 3
Ny = 3
Ne_array = [10, 20, 40, 60, 100, 200];
@time model, data = setup_lorenz63(path, Ne_array);
```

    Ne 10 RMSE: 3.8254663423276716
    Ne 20 RMSE: 2.0605921219961134
    Ne 40 RMSE: 0.4223537564608097
    Ne 60 RMSE: 0.40146546759545937
    Ne 100 RMSE: 0.3877197292448612
    Ne 200 RMSE: 0.4025809532902939
     15.443304 seconds (68.89 M allocations: 6.501 GiB, 6.17% gc time)

```julia
metric_list = benchmark_lorenz63(model, data, path, Ne_array, [1.0])
```

    (Ne, β) = (10, 1.0)
      0.482682 seconds (3.47 M allocations: 329.564 MiB, 11.44% gc time)
    Ne = 10
    Ne 10& β 1.0 RMSE: 0.6015590736992926

    Progress:  17%|██████▉                                  |  ETA: 0:00:05

    (Ne, β) = (20, 1.0)

    Progress:  33%|█████████████▋                           |  ETA: 0:00:04

    0.859515 seconds (6.59 M allocations: 631.994 MiB, 8.56% gc time)
    Ne = 20
    Ne 20& β 1.0 RMSE: 0.5184477645494938
    (Ne, β) = (40, 1.0)

    Progress:  50%|████████████████████▌                    |  ETA: 0:00:04

    1.662244 seconds (12.83 M allocations: 1.210 GiB, 8.65% gc time)
    Ne = 40
    Ne 40& β 1.0 RMSE: 0.4483744816193398
    (Ne, β) = (60, 1.0)

    Progress:  67%|███████████████████████████▍             |  ETA: 0:00:03

    2.463622 seconds (19.07 M allocations: 1.804 GiB, 8.51% gc time)
    Ne = 60
    Ne 60& β 1.0 RMSE: 0.4403101469791623
    (Ne, β) = (100, 1.0)

    Progress:  83%|██████████████████████████████████▏      |  ETA: 0:00:02

    4.334649 seconds (31.55 M allocations: 2.983 GiB, 8.49% gc time)
    Ne = 100
    Ne 100& β 1.0 RMSE: 0.4303647629868742
    (Ne, β) = (200, 1.0)
     13.776164 seconds (62.75 M allocations: 5.935 GiB, 7.19% gc time)
    Ne = 200
    Ne 200& β 1.0 RMSE: 0.4465248257137066

    Progress: 100%|█████████████████████████████████████████| Time: 0:00:24

    6-element Array{Any,1}:
     Metrics[Metrics(10, [0.8166106052280312, 0.8477976249831615, 0.5377435004826744, 0.2946057216215732, 0.1963372890873663, 0.1661211989653268, 0.09550122139602638, 0.10895691501209233, 0.061645037116889324, 0.1556247918798692  …  4.0942175082812975, 2.918216773645667, 2.293946056603929, 2.3368525142834664, 2.9966610036699945, 2.488711153121422, 2.060789991067101, 2.000268585410395, 2.0094430111553767, 1.633586039615903], 0.37686107778341077, 0.6015590736992926, 0.8149386709205617, [1.0125271950095414, 0.9623361052178709, 0.7792957424814254, 0.4344193482216172, 0.2873587820499424, 0.2729751953413108, 0.40382950988831384, 0.6108278943470481, 0.6592219807793599, 0.4157011617649463  …  0.6622613427311322, 0.9156031963385003, 1.1880540752197546, 1.010250360884169, 0.9055102239545725, 0.6429162631150634, 0.5592022167774022, 0.6557579770729094, 0.987851621549649, 1.1234614213819019], 0.4419068954362283, 0.5053138269227165, 0.25903240927442023, [1.0, 1.0, 0.6666666666666666, 0.6666666666666666, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.6666666666666666, 0.6666666666666666, 0.3333333333333333, 0.3333333333333333, 0.0, 0.3333333333333333, 0.0, 0.0, 0.3333333333333333, 0.6666666666666666], 1.0, 0.7261666666666667, 0.36457522740458864)]
     Metrics[Metrics(20, [0.8009234015601246, 0.8873529418152709, 0.9648186335874159, 0.2241641575750264, 0.08972122691323191, 0.40502498541772275, 0.9061432442037632, 1.2895306531554689, 0.7946468671210342, 0.3985672597736556  …  0.3346030298442835, 0.5163997563927475, 0.43517438396749325, 0.598360550796541, 0.5362700352527868, 0.4471320623829449, 0.29561670268590506, 0.18851296824592492, 0.17522643152541226, 0.07857223386207775], 0.35531061663093255, 0.5184477645494938, 0.4845044270612191, [1.3258269952984914, 1.219366525258285, 1.1067738085368013, 0.8876478449926386, 0.6257782703432062, 0.5204509248590731, 0.6103587013183903, 0.7679474941003844, 0.5675737349456604, 0.5082509495209795  …  0.5453088799482404, 0.5575510867223642, 0.5540770509755611, 0.5983031348899548, 0.50022559167807, 0.42259328666553975, 0.3198506611541602, 0.3043017231422324, 0.38395027869333503, 0.5385430054163486], 0.47192182393714555, 0.5305997869038892, 0.24119505026883653, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.3333333333333333, 0.0, 0.6666666666666666, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 1.0, 0.8525, 0.2979733949063843)]
     Metrics[Metrics(40, [0.9826817575939837, 0.8853553187784176, 0.3822879886654962, 0.1560808317176755, 0.16401163158104606, 0.2641915580627298, 0.4821786668586909, 0.7562409207457594, 0.5304897251404828, 0.2335617480836244  …  0.46302073478462163, 0.44641521503457005, 0.5198248995974374, 0.8841581293917309, 0.8944940422060007, 0.6869037355788589, 0.47199979466366954, 0.3146926290239334, 0.22052593050442004, 0.2011946269035898], 0.31988125779252113, 0.4483744816193398, 0.3888888264729722, [0.9910648265607959, 1.0389612210939443, 0.9976972004017787, 0.6975832292555886, 0.40656021680950216, 0.3543519277702079, 0.4841426130471005, 0.6802004103354949, 0.5357937586885131, 0.4567283561463507  …  0.34892484591013456, 0.4717990778913134, 0.618724675393529, 0.5828250426543808, 0.4770344985964706, 0.39083347571382326, 0.2923228793059253, 0.3054108982896446, 0.45375090190774553, 0.6625324478592381], 0.4946079848132213, 0.5558350191174204, 0.23424766693404636, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.6666666666666666, 1.0, 1.0, 0.6666666666666666, 0.6666666666666666, 0.6666666666666666, 0.6666666666666666, 0.6666666666666666, 1.0, 1.0], 1.0, 0.9485, 0.17881512390001802)]
     Metrics[Metrics(60, [0.9958246928389437, 1.1406683287468808, 0.5892665012618263, 0.06937257988643586, 0.11686273681705214, 0.2663933938825168, 0.5322621212098003, 0.9316025189352882, 0.728746365036979, 0.4082893401458479  …  0.22786860029792666, 0.31749262144454576, 0.35389057590962486, 0.5051968531575757, 0.556030690676821, 0.4877666524221944, 0.30372102936096423, 0.2222280953223254, 0.28222115198656517, 0.24913840043694432], 0.30630601607276686, 0.4403101469791623, 0.38370929842878665, [1.15339760434503, 1.111068400675063, 1.0617894821482072, 0.7246278025729431, 0.4659868724356628, 0.4125774571657407, 0.5030501098617005, 0.6339205155508028, 0.5972643721683021, 0.4689687485398599  …  0.3842024564419135, 0.48550370981036023, 0.5900843703168447, 0.6445608536580213, 0.500422315137187, 0.4156164932470485, 0.3016388798926969, 0.2916869146188145, 0.40422049044227376, 0.6320377754649462], 0.5153348159848088, 0.5717853643157081, 0.2325183875067794, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 1.0, 0.9791666666666667, 0.10627248639759496)]
     Metrics[Metrics(100, [1.1447153311331781, 1.38577826954797, 0.8011100613741654, 0.2521263581826285, 0.3020487816972422, 0.508284816604461, 0.8670734818498514, 1.1608770052370112, 0.8718084104569586, 0.5296218649707812  …  0.2398186998545124, 0.25503741685515185, 0.22177804228613865, 0.41110055435395404, 0.47496251760836805, 0.40719857492095896, 0.24764360531773363, 0.13828193095173355, 0.0826234199512165, 0.26856854675486025], 0.2893155689372151, 0.4303647629868742, 0.38191564508486936, [0.9527704170384034, 0.9321001672664082, 1.0398646440457058, 0.6935998849951813, 0.44319915626874445, 0.3912492072948739, 0.5244511435147953, 0.6086134296913539, 0.5544278894887922, 0.4791755006959619  …  0.3928515891332105, 0.5199694696873517, 0.6603156868760311, 0.7450505076075534, 0.6302553655749117, 0.5068577227060873, 0.37451309930446913, 0.35529727887821777, 0.5185904054762793, 0.7753972206416857], 0.5329209748271038, 0.5871951348369144, 0.2322523924985191, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 1.0, 0.9911666666666668, 0.07514664299502112)]
     Metrics[Metrics(200, [0.888273273483637, 1.0797033517549515, 0.3440269499593048, 0.21420101076430767, 0.16162425690189322, 0.3253767298519863, 0.7050060288401291, 1.0723841440703699, 0.7927858236953814, 0.37303796598152517  …  0.14863863331263572, 0.13098694675828024, 0.01569618827238591, 0.16380423196922922, 0.33101947558854705, 0.3183054841536702, 0.19962116984435951, 0.10929821953369317, 0.07402847742163798, 0.21688151831203697], 0.3082727000481167, 0.4465248257137066, 0.38847206086051883, [1.1025016902484062, 1.1585146386588427, 1.050280287671641, 0.7264029154997294, 0.46463567765400915, 0.4408017551630872, 0.5331142954826315, 0.6640979088648434, 0.6326136904712192, 0.5171574662091373  …  0.4291626962129586, 0.5520631256132565, 0.6739989865135746, 0.774635718849873, 0.6359763126942447, 0.5135438852089059, 0.3752157713769756, 0.3581808492199844, 0.47687145842631096, 0.6813329693353127], 0.5508320940815306, 0.6018345158581575, 0.23071327806131858, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 1.0, 0.9908333333333332, 0.07584313842648356)]

```julia
metric_list = benchmark_lorenz63(model, data, path, Ne_array, [1.0])
```

    (Ne, β) = (200, 1.0)
      7.739944 seconds (62.75 M allocations: 5.935 GiB, 6.06% gc time)
    Ne = 200
    Ne 200& β 1.0 RMSE: 0.4189442707516625

    1-element Array{Any,1}:
     Metrics[Metrics(200, [0.27418871478089074, 0.4011111920900346, 0.6224367042654201, 1.1940241372680338, 1.178706937442249, 0.629896584355798, 0.37352902606306937, 0.27598562967074297, 0.3415108073527272, 0.5385900049790774  …  0.22005926456584765, 0.11117371050261328, 0.0745833426310767, 0.06915409494524392, 0.23566185689849842, 0.4284533696688461, 0.7382482677445492, 1.2785156926466672, 0.9477650979445261, 0.8579855755647071], 0.2843509443072877, 0.4189442707516625, 0.3843942458358231, [0.2791191166770778, 0.3308824803118106, 0.5356751453956368, 0.7482658949655999, 0.6849666981780923, 0.5401293761758875, 0.3905840265207302, 0.30529026491290673, 0.3612676308031941, 0.5506587576884893  …  0.9550452598208562, 0.7275900415277436, 0.5097020874896808, 0.38068472189933705, 0.44668108906035564, 0.6711198461183089, 0.7934070557338772, 0.953289237580203, 1.0821830718436276, 1.1373384160933897], 0.529889526630428, 0.584733698922179, 0.23401428353181045, [1.0, 1.0, 1.0, 1.0, 0.6666666666666666, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], 1.0, 0.9928333333333332, 0.065876381401767)]

max_stars_repo_path: notebooks/6 Nonlinear ensemble filtering Lorenz 63.ipynb