# importamos bibliotecas cómputo de matrices
import numpy as np

A = np.matrix([[1, 1],
               [1, 0]])

# para k=3 A se multiplica por sí misma tres veces:
A*A*A

u0 = [1,
      0]

# u3
np.dot(A*A*A,
       u0)

"""
Explanation: Sucesión de Fibonacci
La sucesión de Fibonacci se usa para modelar crecimiento de poblaciones de bichos. Por ejemplo: conejitos. Se empieza con un granero vacío, en $t_{0}$, en $t_{1}$ una parejita de conejos se metió quién sabe cómo y usan ese tiempo para instalarse, en $t_{1}$ sigue la misma parejita, pero una nueva parejita ya se está gestando. En $t_{3}$ nace una segunda pareja, ahora son dos: la nueva usa ese tiempo para instalarse y la primera otra vez está gestando una nueva.
Uno puede imaginarse el área que ocupa la población de conejitos viendo esta figura:
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/d/db/34%2A21-FibonacciBlocks.png/300px-34%2A21-FibonacciBlocks.png">
También puede usarse la serie de Fibonacci para aproximar la proporción áurea:
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/9/93/Fibonacci_spiral_34.svg/320px-Fibonacci_spiral_34.svg.png">
Más información en el artículo "Sucesión de Fibonacci" en la Wikipedia.
La ecuación de recurrencia de segundo orden es así:
$F_{k+2}=F_{k+1}+F_{k}$
<a name="eq1">.......eq1</a>
Con valores iniciales $F_{0}=1, F_{1}=1$ la secuencia es:
$0,1,1,2,3,5,8,13,...$
Hay dos formas de obtener el k-ésimo valor de esta serie. Usando la ecuación de recurrencia <a href="#eq1">eq1</a> o usando esta fórmula cerrada:
$F_{k}=\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^{k}-\left(\frac{1-\sqrt{5}}{2}\right)^{k}\right]$<a name="eq2">.......eq2</a>
Teorema
Sea $A$ una matriz cuadrada de dimensión n. Si $A$ tiene n valores propios distintos
$\lambda_{1},\lambda_{2},...,\lambda_{n},$
con vectores propios correspondientes 
$v_{1},v_{2},...,v_{n},$
entonces $A$ es diagonalizable, o sea existen las matrices $P$ y $D$ de dimensión n tales que $D$ es una matriz diagonal y
$P^{-1}AP=D$
$D=\left[\begin{array}{cccc}
 \lambda_{1}  &  0  &  \cdots\  &  0\
0  &  \lambda_{2}  &  \cdots\  &  0\
\vdots\  &  \vdots\  &  \ddots\  &  \vdots\
0  &  0  &  0  &  \lambda_{n} 
\end{array}\right]$
$P$ es la matriz cuyas columnas son los vectores propios.
$P=\left[v_{1}|v_{2}|...|v_{n}\right]$
Se generan los núneros de la sucesión de Fibonacci en las coordenadas de los vectores en $R2$ a partir de la condición inicial
$u_{0}=\left[\begin{array}{c}
F_{1}\
F_{0}
\end{array}\right]=\left[\begin{array}{c}
1\
0
\end{array}\right]$
Sea
$A=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]$
<a name="eq3">.......eq3</a>
entonces, si
$u_{k}=A\mathbf{u}_{k-1}$
<a name="eq4">.......eq4</a>
para $k=1,2,3...$
$u_{k}=\left[\begin{array}{c}
F_{k+1}\
F_{k}
\end{array}\right]$
$u_{k}=A^{k}\mathbf{u}_{0}$
<a name="eq5">.......eq5</a>

Se usa la ecuación <a href="#eq4">eq4</a> para generar algunos $u$.
$u_{0}=\left[\begin{array}{c}
F_{1}\
F_{0}
\end{array}\right]=\left[\begin{array}{c}
1\
0
\end{array}\right]$
$u_{1}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
1\
0
\end{array}\right]=\left[\begin{array}{c}
1\
1
\end{array}\right]$
$u_{2}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
1\
1
\end{array}\right]=\left[\begin{array}{c}
2\
1
\end{array}\right]$
$u_{3}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
2\
1
\end{array}\right]=\left[\begin{array}{c}
3\
2
\end{array}\right]$
$u_{4}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
3\
2
\end{array}\right]=\left[\begin{array}{c}
5\
3
\end{array}\right]$
$u_{5}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
5\
3
\end{array}\right]=\left[\begin{array}{c}
8\
5
\end{array}\right]$
$u_{6}=\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]
\left[\begin{array}{c}
8\
5
\end{array}\right]=\left[\begin{array}{c}
13\
8
\end{array}\right]$
A continuación pongamos a prueba la <a href="#eq5">eq5</a>.
End of explanation
"""


from math import sqrt

# eigenvalores
l1=(1+sqrt(5))/2.0
l2=(1-sqrt(5))/2.0

print("lambda1=%s, lambda2=%s" % (l1,l2))

# la matriz de transicion P
P = np.matrix([[l1, l2],
               [1,  1]])

# valores como 1.66533454e-16 son como ceros! esto pasa por no hacerlo analiticamente
(P**-1)*A*P

"""
Explanation: Polinomio característico
Como $\lambda$ es un valor propio de $A$
$(\lambda I-A)\overrightarrow{X}=\overrightarrow{0}$
un eigenvector que sea su solución depende de que
$\det(\lambda I-A)=\det\left(\left[\begin{array}{cc}
\lambda & 0\
0 & \lambda
\end{array}\right]-\left[\begin{array}{cc}
1 & 1\
1 & 0
\end{array}\right]\right)=\overrightarrow{0}$
Luego
$\det\left(\left[\begin{array}{cc}
\lambda-1 & -1\
-1 & \lambda
\end{array}\right]\right)=(\lambda-1)\lambda-(-1)(-1)=\lambda^{2}-\lambda-1$
Valores propios
Para obtener los valores propios usamos la fórmula general para encontrar valores de $\lambda$ que satisfagan el polinomio característico
$p(\lambda)=\lambda^{2}-\lambda-1$
Acá la fórmula general:
$x=\frac{-b\pm\sqrt{b^{2}-4ac}}{2a}$
$\lambda=\frac{1\pm\sqrt{-1^{2}-4(1)(-1)}}{2(1)}=\frac{1\pm\sqrt{1+4}}{2}$
Por eso
$\begin{array}{c}
\lambda_{1}=\frac{1+\sqrt{5}}{2}\
\lambda_{2}=\frac{1-\sqrt{5}}{2}
\end{array}$
<a name="eq6">.......eq6</a>
## Vectores propios
$\left[\begin{array}{cc}
\lambda-1 & -1\
-1 & \lambda
\end{array}\right]\overrightarrow{X}=\left[\begin{array}{cc}
\lambda-1 & -1\
-1 & \lambda
\end{array}\right]\left[\begin{array}{c}x_{1}\
x_{2}\end{array}\right]$
Equivale al sistema de ecuaciones
$(\lambda-1)x_{1}-x_{2}=0$....(1)
$-x_{1}+\lambda x_{2}=0$.....(2)
Despejando $x_{1}$ en (2)
$x_{1}=\lambda x_{2}$
que sustituyendo en (1)
$(\lambda-1)\lambda x_{2}-x_{2}=0$
factorizamos $x_{2}$ 
$x_{2}\left[(\lambda-1)\lambda-1\right]=0$
sea $x_{2}=1$ en (2)
$-x_{1}+\lambda=0$
$x_{1}=\lambda$
Por eso los vectores propios son
$v_{1}=\left[\begin{array}{c}
\lambda_{1}\
1
\end{array}\right],\: v_{2}=\left[\begin{array}{c}
\lambda_{2}\
1
\end{array}\right]$
<a name="eq7">.......eq7</a>
Matriz de transición P
Con $v_{1}$ y $v_{2}$ formamos la matriz $P$.
Para obtener la matriz invertida $P^{-1}$ Anton y Rorres muestran este teorema:
Si
$A=\left[\begin{array}{cc}
a & b\
c & d
\end{array}\right]$
entonces
$A^{-1}=\frac{1}{ad-bc}\left[\begin{array}{cc}
d & -b\
-c & a
\end{array}\right]$
Teniendo la matriz
$P=\left[\begin{array}{cc}
\lambda_{1} & \lambda_{2}\
1 & 1
\end{array}\right]$
entonces
$P^{-1}=\frac{1}{\lambda_{1}-\lambda_{2}}\left[\begin{array}{cc}
1 & -\lambda_{2}\
-1 & \lambda_{1}
\end{array}\right]$
Para comprobar
$P^{-1}AP=\left[\begin{array}{cc}
\lambda_{1} & 0\
0 & \lambda_{2}
\end{array}\right]=D$
basta con multiplicar las matrices $P^{-1}AP$. Como está laborioso hagámoslo con numpy en la siguiente celda.
End of explanation
"""


def fibonacci(k):
    if k <= 1:
        return k
    else:
        return fibonacci(k-1) + fibonacci(k-2)

#  el indice empieza en cero
for n in range(25):
    print("F%s=%s" % (n,fibonacci(n)))

"""
Explanation: En pos de la "fórmula cerrada"
Si $A$ es diagonalizable entonces existe 
$P^{-1}AP=D$
Como $PP^{-1}=I$ es la matriz identidad, y es el neutro multiplicativo.
$PP^{-1}APP^{-1}=PDP^{-1}$
luego
$A=PDP^{-1}$
y entonces
$A^{k}=\left(PDP^{-1}\right)^{k}$
ello significa que $PDP^{-1}$ se multiplica $k$ veces
$A^{k}=\left(PDP^{-1}\right)\left(PDP^{-1}\right)\left(PDP^{-1}\right)...$
cambiando la agrupación de los paréntesis se tiene
$A^{k}=PDDD...DP^{-1}$
lo cuál explica la siguiente ecuación:
$A^{k}=PD^{k}P^{-1}$
<a name="eq8">.......eq8</a>

Sustituyendo <a href="#eq8">eq8</a> en <a href="#eq4">eq4</a> tenemos que
$\mathbf{u}{k}=(PD^{k}P^{-1})\mathbf{u}{0}$
no olvidar que
$u_{0}=\left[\begin{array}{c}
F_{1}\
F_{0}
\end{array}\right]=\left[\begin{array}{c}
1\
0
\end{array}\right]$
Luego se trata de otra larga multiplicación de matrices.
$u_{k}=\left[\begin{array}{cc}
\lambda_{1} & \lambda_{2}\
1 & 1
\end{array}\right]\left[\begin{array}{cc}
\lambda_{1}^{k} & 0\
0 & \lambda_{2}^{k}
\end{array}\right]\frac{1}{\lambda_{1}-\lambda_{2}}\left[\begin{array}{cc}
1 & -\lambda_{2}\
-1 & \lambda_{1}
\end{array}\right]\left[\begin{array}{c}
1\
0
\end{array}\right]$
$u_{k}=\left[\begin{array}{cc}
\lambda_{1} & \lambda_{2}\
1 & 1
\end{array}\right]\left[\begin{array}{cc}
\lambda_{1}^{k} & 0\
0 & \lambda_{2}^{k}
\end{array}\right]\left[\begin{array}{c}
1\
-1
\end{array}\right]\frac{1}{\lambda_{1}-\lambda_{2}}$
$u_{k}=\left[\begin{array}{cc}
\lambda_{1} & \lambda_{2}\
1 & 1
\end{array}\right]\left[\begin{array}{c}
\lambda_{1}^{k}\
-\lambda_{2}^{k}
\end{array}\right]\frac{1}{\lambda_{1}-\lambda_{2}}$
y finalmente
$\mathbf{u}{k}=(PD^{k}P^{-1})\mathbf{u}{0}=\frac{1}{\lambda_{1}-\lambda_{2}}\left[\begin{array}{c}
\lambda_{1}^{k+1}-\lambda_{2}^{k+1}\
\lambda_{1}^{k}-\lambda_{2}^{k}
\end{array}\right]$
También
$u_{k}=\left[\begin{array}{c}
F_{k+1}\
F_{k}
\end{array}\right]$ entonces
$F_{k}=\frac{1}{\lambda_{1}-\lambda_{2}}\left(\lambda_{1}^{k}-\lambda_{2}^{k}\right)$
Por otro lado
$\frac{1}{\lambda_{1}-\lambda_{2}}=\frac{1}{\frac{1+\sqrt{5}}{2}-\frac{1-\sqrt{5}}{2}}=\frac{1}{\frac{(1+\sqrt{5})-(1-\sqrt{5})}{2}}=\frac{1}{\frac{1+\sqrt{5}-1+\sqrt{5})}{2}}=\frac{1}{\frac{2\sqrt{5}}{2}}=\frac{1}{\sqrt{5}}$
que si lo multiplicamos por $\left(\lambda_{1}^{k}-\lambda_{2}^{k}\right)$ equivale a la <a href="#eq2">eq2</a>:
$F_{k}=\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^{k}-\left(\frac{1-\sqrt{5}}{2}\right)^{k}\right]$
Implementación Recursiva
End of explanation
"""


def lucas(k):
    if k == 0:
        return 1
    elif k == 1:
        return 2
    else:
        return lucas(k-1) + lucas(k-2)

#  el indice empieza en cero
for n in range(25):
    print("F%s=%s" % (n,lucas(n)))

"""
Explanation: Sucesión de Lucas
Esta serie es igual a la de Fibonacci pero empieza con $F_{0}=1$ y $F_{1}=2$.
End of explanation
"""


# importamos bibliotecas para plotear
import matplotlib
import matplotlib.pyplot as plt

# para desplegar los plots en el notebook
%matplotlib inline

p = []
for n in range(11):
    p.append([fibonacci(n), lucas(n)])

figura = plt.plot( p )
p

"""
Explanation: ¿Cómo se ven las series de Lucas y Fibonacci juntas?
End of explanation
"""


def fibonacci_formula(k):
    return (1.0/sqrt(5))*(((1.0+sqrt(5))/2)**k-((1.0-sqrt(5))/2)**k)

#  el indice empieza en cero
for n in range(25):
    print("F%i=%i" % (n,fibonacci_formula(n)))


"""
Explanation: Fibonacci con fórmula
End of explanation
"""


import timeit

# cronometro para implementacion recursiva
tr = timeit.Timer("""def fibonacci(k):
    if k <= 1:
        return k
    else:
        return fibonacci(k-1) + fibonacci(k-2)

    fibonacci(25)""")

# cronómetro para implementación con formula
tf = timeit.Timer("""from math import sqrt

def fibonacci_formula(k):
    return (1.0/sqrt(5))*(((1.0+sqrt(5))/2)**k-((1.0-sqrt(5))/2)**k)

fibonacci_formula(25)""")

# ejecutarlas diez mil veces, medir su tiempo
tiempos=[tr.timeit(10000),
         tf.timeit(10000)]

# armar gráfica de barras
plt.bar([0,1], tiempos, align='center', alpha=0.5)
plt.title('10,000 ejecuciones, k=25')
plt.ylabel('milisegundos')
plt.xticks([0,1],['recursivo', 'formula'])
plt.show()

tiempos

"""
Explanation: Dificultades algorítmicas
A continuación se miden los tiempos de ejecución de ambas implementaciones de la sucesión de Fibonacci. Cada una se ejecuta diezmil veces y se reporta el tiempo promedio.
End of explanation
"""


