import numpy as np # con eso voy a poder acceder a las funciones de numpy a través de np.función()

# en ejemplo
print('El numero e = {:.3}'.format(np.e))
print('O el numero Pi = ', np.pi)

# y podemos calcular senos y cosenos entre otras cosas
print(np.sin(np.pi)) # casi cero! guarda con los floats!

"""
Explanation: Taller de Python - Estadística en Física Experimental - 2do día

NumPy: vectores, matrices y tablas de datos
NumPy es LA biblioteca para cálculo vectorial. Además de contener un nuevo tipo de dato que nos va a ser muy útil para representar vectores y matrices, nos provee de un arsenal de funciones de todo tipo.
Vamos a empezar por importar la bliblioteca numpy. La sintaxis típica de eso era import biblio as nombre:
End of explanation
"""


a = np.array([1, 2, 3, 4]) # array toma como argumento un vector_like (lista, tupla, otro array...vector_like)
b = np.array([5, 6, 7, 8])

print(type(a)) #Tipos de dato de a

print(a + b) # vector suma
print(a * b) # acá multiplicó

# en el caso de las listas esto era muy molesto
l1 = [1, 2, 3, 4]
l2 = [5, 6, 7, 8]

print(l1 + l2) # sumar concatena
# print(l1 * l2) # esto ni siquiera se puede hacer!

"""
Explanation: Todo eso está muy bien, pero lo importante de numpy son los arrays numéricos. Los arrays numéricos nos van a servir para representar vectores (el objeto matemático) o columnas/tablas de datos (el objeto oriyinezco o de laboratorio).
La idea es que es parecido a una lista: son muchos números juntos en la misma variable y están indexados (los puedo llamar de a uno dando la posición dentro de la variable). La gran diferencia con las listas de Python es que los arrays de numpy operan de la forma que todos queremos:
1. Si sumamos o restamos dos arrays, se suman componente a componente.
2. Si multiplicamos o dividimos dos arrays, se multiplican o dividen componente a componente.
Veamos ejemplos usando la función array para crear arrays básicos.
End of explanation
"""


print(a[0], a[1], a[2], a[3]) # son 4 elementos, los indices van del 0 al 3
# y más o menos vale todo lo que valía con listas
print(b[-1]) # agarro al último elemento de b
print(b[1:]) # desde el primero hasta el 3 (no incluido el final, nunca se incluye)

"""
Explanation: Y al igual que con las listas, uno puede acceder a elementos específicos de un array:
End of explanation
"""


# equiespaciados
equilin = np.linspace(0, 1, 10) # 10 número equiespaciados linealmente entre 0 y 1
print('Equiespaciado lineal:', equilin)

arange = np.arange(0, 1, 1./10) # como el range de Python pero puede venir con un paso en coma flotante
print('Como el range de las listas:', arange)

identidad = np.identity(3)
print('Identidad de 3x3:', identidad)
print()

#otros para que prueben ustedes
ceros = np.zeros((4, 4)) # todos ceros, matriz de 4x4
unos = np.ones((2,3)) # todos unos, matriz de 2x3
ojos = np.eye(5, k=0) # unos en la diagonal, como identidad
ojos2 = np.eye(5, k=2) # qué pasó acá?

print(ceros)
print() 
print(unos)
print()
print(ojos)
print()
print(ojos2)

"""
Explanation: Para facilitar la vida del usuario numpy viene con un montón de rutinas de creación de arrays típicos. En particular, matrices típicas como las identidades o de todos elementos iguales a 1 o 0 y arrays con cierta cantidad de elementos entre dos números (muy útil para crear dominios para gráficos).
Veamos ejemplos de esos:
End of explanation
"""


x = np.linspace(0, 10, 1000) # ese array tiene 1000 elementos, andar printeando es poco práctico!
print(x.dtype) # array.dtype nos dice que tipo de elementos tiene el array
ceros = np.zeros((100, 100)) # matriz de 100x100
print(ceros.shape) # array.shape nos dice cuántas filas y columnas tiene el array


# prueben que pasa cuando le piden el shape a un array con una sola fila o columna como el x

"""
Explanation: Y antes de seguir, algo que siempre puede ser útil: los arrays tienen ciertas propiedades como su shape (de cuánto por cuánto) y el dtype (qué tipo de cosas tiene adentro). Podemos acceder a estos datos de la siguiente manera:
End of explanation
"""


print("Vector concatenado: ", np.r_[np.arange(-1, 1, 0.3), 2, 3, 4, np.linspace(5, 6, 10)])

#Para concatenar tres vectores podemos usarlo
a = np.linspace(0, 5, 5)
b = np.linspace(6, 11, 5)
c = np.linspace(12, 17, 5)
print("Matriz concatenada: ", np.r_['0,2', a, b, c].T)

"""
Explanation: Para construir un vector (o array) más complejo podemos usar las funciones r_ y c_ (busquen la documentación de c_, pero r_ mostramos un ejemplo)
End of explanation
"""


# Realicen el ejercicio 1

"""
Explanation: Ejercicio 1
Creen un vector "dominio" con 20 valores entre -5 y 5, y un vector "imagen" donde la imagen sea la de la función $x^2$
End of explanation
"""


# Realicen el ejercicio 2

"""
Explanation: Ejercicio 2
Creen una matriz de 15x15 con los números del 1 al 255.
Ayuda: siempre que puedan utilicen las funciones de creación de arrays de numpy. 
Ayuda bis: la función reshape hace lo que promete (reshapear) y puede tomar tuplas como argumento. (Si, hay que googlear. ¿Nunca les pasó?)
End of explanation
"""


from matplotlib import pyplot as plt
# muestra los gráficos en el mismo notebook
%matplotlib inline

"""
Explanation: numpy es además la librería encargada de las operaciones de álgebra lineal. Ya de por si la suma y producto por escalar de arrays de dos dimensiones son operaciones matriciales. Leanse la documentación de numpy.linalg (o revisen nuestro apartado de álgebra lineal del taller numérico); la idea es que está hecho para ser MATLAB compatible.
Gráficos, datos y ajustes
Hacer gráficos es lo primero que aprendimos a hacer en Origin, así que es lo primero que vamos a aprender para reemplazarlo. Van a ver que no es nada complicado.
Primero, debemos importar las bibliotecas necesarias para graficar, numpy por si no la teníamos, y de la biblioteca matplotlib (que tiene infinitas funciones y posibilidades), solamente pyplot de donde sacaremos las funciones que necesitaremos para graficar.
End of explanation
"""


# Con esto cargamos los datos
# Si vemos los archivos son dos columnas
# que vamos a guardar como x, y
data = np.loadtxt("doble_exp.dat") 

x = data[:,0]
y = data[:,1]
error_y = data[:,2]

"""
Explanation: Importemos los datos que tenemos en el archivo datos.csv. Este archivo tiene tres columnas, x, y, errores en y. Separemoslos en tres variables
End of explanation
"""


# Ploteamos. Los datos como puntos de color rojo
plt.plot(x, y, 'ro', label = 'Datos')

# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best') 

plt.show() # si no usaron %matplotlib inline, esto abre una ventanita con el gráfico


"""
Explanation: Para tener un poco de intuición de los datos, sin el error en y, primeros grafiquemoslos con plt.plot
End of explanation
"""


# Ploteamos
plt.plot(x, y, 'go', label = 'Modelo')

# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.yscale('log')
plt.xscale('log')
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best') 

plt.show()


"""
Explanation: Notemos que dentro de la función plot pusimos como parámetros ro que significa que el color que queremos para la curva sea azul, y que el trazo sea una línea continua. Esto es customizable, pueden probar letras de otros colores (g, r, y, k) o bien otros trazos. 
Alteremos los ejes para que sean logaritmicos (porque podemos, aunque es un proceso standard cuando tenés datos "exponenciales") y los puntos de color verde
End of explanation
"""


plt.errorbar(x, y, yerr=error_y, fmt='go', label="Datos")

# Detalles del gráfico
plt.grid(True) # Para que quede en hoja cuadriculada
plt.yscale('log')
plt.xscale('log')
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best') 

plt.show()


"""
Explanation: Así aparece una recta a ajustar, con el proceso de linealizado. Esto ya lo deberían tener más o menos claro de cursar laboratorios, pero es algo muy importante; los algoritmos lineales de cuadrados mínimo (sea por la formula normal o por descendiente del gradiente) siempre convergen. No les voy a contar toda la historia porque ya lo van a ver, pero por ahora creame
Volviendo sobre la cosas para customizar los gráficos, pasense por el repositorio nuestro pero especialmente por la documentation de matplotlib
Finalmente, para graficar los errores además de los datos tenemos la función plt.errorbar (vean la documentación y entiendan por qué se usa como la usamos acá)
End of explanation
"""


# Guardamos y cargamos, a modo de ejemplo

#vean que crea un archivo Datos_taller.txt
np.savetxt('Datos_taller.txt', np.log(data[:,:-1]), delimiter = '\t')

Data = np.loadtxt('Datos_taller.txt', delimiter = '\t') 

#plt.plot(Data[:,0],Data[:,1], 'r.') # Veamos que son los mismos datos, pero linealizados!
plt.errorbar(Data[:,0],Data[:,1], yerr=1/y*error_y, fmt='ro')
plt.show()

"""
Explanation: Si necesitamos guardar cambios de los datos, por ejemplo con el logaritmo aplicado, podemos usar la función np.savetxt(que gaurdar en formato csv)
End of explanation
"""


from scipy.optimize import curve_fit

"""
Explanation: Obtenidos los datos y con un poco de inteligencia adquirida de ellos, queremos efectuar un ajuste. Para eso importamos la biblioteca con la función que usaremos, que aplica cuadrados mínimos para obtener los coeficientes.
End of explanation
"""


f = lambda x, A, B: A * x + B

# Ajustamos, pero con las funciónes logaritmicas. Usamos propagación de errores
popt, pcov = curve_fit(f, np.log(x), np.log(y), sigma = 1/y * error_y, 
                       absolute_sigma=True)

print(popt)
print(pcov)

"""
Explanation: El algoritmo de cuadrados mínimos necesita la función con la que queremos ajustar (que como ya linealizamos, es lineal), que vamos a definir como como función lambda, dominio, los datos, un vector con los valores iniciales de los parámetros desde donde debe comenzar a iterar. Los parámetros iniciales son importantes para el ajuste correcto (aún para el caso lineal, aunque acá el método hace un guess correcto).
La función nos devolverá 2 cosas. Primero, los parámetros optimizados por este algoritmo, ordenados como los pusimos en la función lambda cuando la definimos, que lo guardamos en el vector popt. Por otro lado nos dará la matriz de covarianza (recuerden, que tiene en su diagonal los $\sigma^2$ de cada parámetro).
End of explanation
"""


t = np.linspace(min(x), max(x), 1000) #Las funciones "viejas" de python siguen funcionando!
t = np.log(t) #Aplico logaritmo
plt.plot(np.log(x), np.log(y), 'ro', label = 'Datos')
plt.plot(t, f(t, *popt), 'g-', label = 'Ajuste') #grafico la función

# Detalles del gráfico
plt.grid(True)
plt.title('Grafico ejemplo')
plt.xlabel('Valores en x')
plt.ylabel('Valores en y')
plt.legend(loc = 'best')

plt.show()

"""
Explanation: Listo, ahora chequeamos con un gráfico que haya ajustado
End of explanation
"""


# Realicen el ejercicio 3

"""
Explanation: ¿Se les ocurre otra forma de obtener el ajuste que no sea por curve_fit? Revisen la librería optimize, ahí pueden encontrar la respuesta (y con lo que van a ver en las clases)
Ejercicio 3

Hagan un ajuste sobre la función $f(x) = Acos(\omega x)$ con $A = 2$ y $\omega = 3$ para 40 valores en $Dom = [-\pi, \pi]$ con valores que varían el 15% del valor dado por el modelo, y compare los parámetros obtenidos con los dados.
Bonus track: Se puede escribir en LaTeX sobre los gráficos. Averiguen qué biblioteca hace falta importar y presenten los parámetros ajustados en el título.
End of explanation
"""


mu, sigma = 100, 15 # mu es mi valor medio, sigma la desviación
x = mu + sigma * np.random.randn(10000) # Ya deben saber, que z = (x - mu)/s es una N(0,1)

n, bins, patches = plt.hist(x, bins=50, normed=1, facecolor='green', 
                            edgecolor='black', alpha=0.75)
# en la variable n se encuentran los datos del histograma
# bins es un vector con los bordes de los rangos de datos
# patches no nos interesa en general

"""
Explanation: Datos aleatorios. Histogramas
Hasta ahora hicimos un análisis de datos adquiridos en el laboratorio, pero no hablamos de datos aleatorios (lo que manejamos cuando hablamos de estadísticas!)
La computadora tiene internamente un mecanismo de creación de números pseudoaleatorios (usando un algoritmo bastante elegante) que permite obtener muestras de una variable con distribución uniforme "continua" (¿me explican por qué las comillas?). Sobre eso, Python crea un montón de librerías para muestrear variables en otras distibuciones
Y luego, vamos a crear histogramas, es decir un gráfico donde dividimos una serie de datos en rangos y contamos cuántos de nuestros datos caen en cada rango. A esos rangos se los llama bins. 
La función que vamos a usar es plt.hist toma como argumentos un array de números, en cuántos bins queremos dividir a nuestro eje x y algunas otras opciones de color como constante de normalización y color de las barras.
Hagamos un histograma simple de un set gaussiano. Para eso, creemos datos alrededor de algún valor medio usando randn de NumPy:
End of explanation
"""


n, bins, patches = plt.hist(x, bins=100, normed=1, facecolor='green', alpha=0.75)
n, bins, patches = plt.hist(x, bins=10, normed=1, facecolor='red', alpha=0.75)

"""
Explanation: Y ya que estamos, para mostrar cómo afecta la elección de bins, graficamos dos histogramas uno arriba del otro.
End of explanation
"""


bins = np.arange(0, 200, 5)
plt.hist(x, bins = bins, normed=1, facecolor='green', 
                            edgecolor='black', alpha=0.75);

"""
Explanation: Los bins los podemos construir como un array, con linspace, pero tenemos que recordar que estamos dando el borde derecho de los bins
End of explanation
"""


import scipy as sp
from scipy import stats #Tengo que importar la sublibrería

X = sp.stats.chi2(df=3)

t = np.linspace(0, 10, 1000)
plt.plot(t, X.pdf(t), 'r-') #Puedo usar la PDF de la distribución así

"""
Explanation: No solo podemos muestrear variables aleatorias con distribución normal, numpy.random tiene una gran selección de distribuciones. Aún así en la librería scipy.stats tenemos aún más distribuciones y funciones estadísticas, pero requieren otra forma de uso. Veamos como muestrar una distribución $\chi^2$ con 3 grados de libertad
End of explanation
"""


print(type(X))
#print(help(X)) #Hagan esto si están en un entorno no tan cool como Jupyter

"""
Explanation: La variable X que construimos con scipy.stats tiene propiedades útiles para la inferencia estadística. Es más veamos su tipo y además usemos la función help para ver que contiene
End of explanation
"""


print(X.mean(), X.moment(1))
print(X.var(), X.moment(2)) #Calculen si está bien la varianza y el momento no centrado
print(X.expect(lambda x: x**3), X.moment(3)) #Integral entre 0 e infinito de x**3. 

"""
Explanation: El tipo rv_frozen se llama así porque es una distribución congelada en sus parámetros. Podríamos construir una distribución para un caso puntual a partir de la función scipy.stats.chi2, variando df en cada caso. Pero para mejorar la lectura, podemos considerar a X como una variable aleatoria, donde tenemos conocimiento de su distribución, de sus momentos y podemos muestrar de ella. 
Veamos primero los primeros momentos y la capacidad de integrar funciones arbitrarias
End of explanation
"""


A = X.rvs(10000)
plt.hist(A, bins=50, normed=1, facecolor='green', 
                            edgecolor='black', alpha=0.75);

"""
Explanation: Mientras si queremos muestrar de la distribución tenemos la función rvs
End of explanation
"""


# Realicen el ejercicio 4

"""
Explanation: Ejercicio 4
Muestren de la distribución exponencial con scipy.stats (que pueden ver acá) y tomen 10000 muestras de esta distribución. ¿Se les ocurre una forma de encontrar una exponencial desde una uniforme? Verifiquen muestreando 1000 datos de la uniforme y aplicando la transformación que realmente es una exponencial
End of explanation
"""


def plot_histogram(data, bins, ticks=5, xlabel='X', ylabel='Histograma', density=True, ecolor='', ax = None):
    N = len(data)
    num_bins = len(bins)
    xmin = min(bins)
    xmax = max(bins)
    
    
    if ax is None:
        ax = plt.gca()
    
    hist, _ = np.histogram(data, bins=bins, density=density)
    
    normed = N if density else 1
    error =  np.sqrt(hist / normed)  # por qué es así el error?
    
    # Se le puede pasar un yerr y xerr, como a errorbar
    # ecolor es el color de los errores y capsize es el tamaño del fin de barra de error
    plt.bar(x=bins[:-1], height=hist, width=np.diff(bins), yerr=error, edgecolor='black', ecolor="red", capsize=2)

    plt.xticks(np.linspace(0, num_bins - 2, ticks), np.linspace(xmin, xmax, ticks), rotation='horizontal');
    plt.xlabel(xlabel);
    return ax

plt.figure(figsize=(7 * 1.78, 7))
ax = plot_histogram(A, np.linspace(0, 25, 50))
plt.grid(True, axis='y')

"""
Explanation: Con esto último ya tienen todas las herramientas para lidiar con las guías y el parcial, pero antes de terminar vamos a ver como agregarle barras de error a los histogramas.
Herramientas de ploteo estadístico
Para hacer un gráfico de histograma con errores tenemos un par de opciones. La más intuitiva es usar las herramientas pandas, que permiten crear a partir de un histograma computado con numpy.histogram un gráfico de barras. La función, que pueden usar y adaptar para esto es
End of explanation
"""


# Volvemos a cargar los datos y los guardamos en x, y
data = np.loadtxt("doble_exp.dat") 

x = data[:,0]
y = data[:,1]
error_y = data[:,2]

"""
Explanation: Yapa: Monte Carlo para estimación de parámetros
El método de Monte Carlo en estadística pertenece a una familia de métodos de remuestro, que permite obtener, de muestras pequeñas de datos, estimadores con sus valor estimado y varianza. Veamos en qué consiste este método. 
Primero volvamos a cargar los datos de double_exp.dat
End of explanation
"""


#Creo al azar datos dentro del error de y con distribución gaussiana
yp = y + np.random.randn(y.shape[0]) * error_y
plt.plot(x, yp, 'ro')
plt.yscale('log')
plt.grid()
plt.errorbar(x, y, yerr=error_y, fmt='b*');

#Creo al azar datos dentro del error de y con distribución gaussiana
M = []
logx = np.log(x)

f = lambda x, A, B: A * x + B
for i in range(10000):
    
    #Vuelvo a computar y
    yp = y + np.random.randn(y.shape[0]) * error_y
    logy = np.log(yp)
    
    #Dado que calculo el logaritmo, debo sacarme los valores NaN (que son logaritmos negativos)
    valid_idx = ~np.isnan(logy)
    
    p, cov = sp.optimize.curve_fit(f, logx[valid_idx], logy[valid_idx])
    
    M.append(p)

M = np.array(M)

#Valor medio y covarianza
p = np.mean(M, axis=0)
cov = np.cov(M.T)

#Imprimamos
print(p)
print(cov)

t = np.linspace(min(x), max(y), 1000)
plt.grid()
plt.errorbar(np.log(x), np.log(y), yerr=1/y * error_y, fmt='go')
plt.plot(np.log(t), f(np.log(t), *p), 'r-');

"""
Explanation: El método consiste en simular un nuevo conjunto de mediciones usando las propiedades del error, que en este caso asumimos gaussiano (podría ser de otra propiedad, pero por el teorema central del límite casi siempre será una distribución normal).
End of explanation
"""


counts, xedges, yedges = np.histogram2d(M[:,0], M[:,1], bins=100, normed=True);

plt.xlim((min(M[:,0]), max(M[:,0])))
plt.ylim((min(M[:,1]), max(M[:,1])))
plt.xlabel("A")
plt.ylabel("B")
plt.contour(xedges[:-1], yedges[:-1], counts);

"""
Explanation: Si revisan las cuentas de curve_fit les da un resultado sumamente parecido, lo que nos da un poco más de confianza en el método (o curve_fit hace algo parecido, fijense que scipy tiene a disposición el código fuente de forma bastante a mano).
Finalmente hagamos el histograma de ambos estimadores en 2D. Para eso usamos numpy.histogram2d que nos devuelve las cuentas, las coordenadas en x, y. Lean lo que dice histogram2d para entenderlo (y por qué usamos los bordes xedges e yedges con un elemento menos)
End of explanation
"""


#Ejercicio 5

"""
Explanation: Ahora la pregunta que deben hacerse, ¿están o no correlacionado los parámetros? ¿Tiene sentido?. 
Ejercicio 5 (de yapa)
Prueben ejecutar el ajuste no linealizado con Monte Carlo. Es decir en vez de ajustar una lineal, deben ajustar $$f(x|A,B,C,D, E) = A + B \exp(C x) + D \exp(E x)$$
Elija como parámetros iniciales $A=10$, $B=130$, $C=-0.001$, $D=960$ y $E=-0.02$
End of explanation
"""


