La librería científica Scipy
============================

Scipy es el paquete científico (es decir, un módulo que tiene otros módulos) más completo, que incluye interfases a librerías científicas muy conocidas como LAPACK, BLAS u ODR entre muchas otras . Si importamos ``scipy`` para tener su espacio de nombres y consultamos la ayuda, vemos todos los módulos que posee:

.. code:: ipython3

    In [1]: import scipy
    In [2]: help(scipy)

::

     cluster                      --- Vector Quantization / Kmeans
     fftpack                      --- Discrete Fourier Transform algorithms
     integrate                    --- Integration routines
     interpolate                  --- Interpolation Tools
     io                           --- Data input and output
     lib                          --- Python wrappers to external libraries
     lib.lapack                   --- Wrappers to LAPACK library
     linalg                       --- Linear algebra routines
     misc                         --- Various utilities that don't have
                                      another home.
     ndimage                      --- n-dimensional image package
     odr                          --- Orthogonal Distance Regression
     optimize                     --- Optimization Tools
     signal                       --- Signal Processing Tools
     sparse                       --- Sparse Matrices
     sparse.linalg                --- Sparse Linear Algebra
     sparse.linalg.dsolve         --- Linear Solvers
     sparse.linalg.dsolve.umfpack --- :Interface to the UMFPACK library:
                                      Conjugate Gradient Method (LOBPCG)
     sparse.linalg.eigen.lobpcg   --- Locally Optimal Block Preconditioned
                                      Conjugate Gradient Method (LOBPCG) [*]
     special                      --- Airy Functions [*]
     lib.blas                     --- Wrappers to BLAS library [*]
     sparse.linalg.eigen          --- Sparse Eigenvalue Solvers [*]
     stats                        --- Statistical Functions [*]
     lib                          --- Python wrappers to external libraries
                                      [*]
     lib.lapack                   --- Wrappers to LAPACK library [*]
     integrate                    --- Integration routines [*]
     ndimage                      --- n-dimensional image package [*]
     linalg                       --- Linear algebra routines [*]
     spatial                      --- Spatial data structures and algorithms
     special                      --- Airy Functions
     stats                        --- Statistical Functions


Ajustes lineales y no lineales
------------------------------

Si simplemente necesitamos hacer ajustes básicos de polinomios, lo podemos hacer fácilmente sólo con ``numpy``:

.. code:: ipython3

    # Librería numpy
    import numpy as np

    # Datos experimentales
    x = np.array([ 0.,  1.,  2.,  3.,  4.])
    y = np.array([ 10.2 ,  12.1,  15.5 ,  18.3,  20.6 ])

    # Ajuste a una recta (polinomio de grado 1)
    p = np.polyfit(x, y, 1)

    print(p)
    # imprime [ 2.7   9.94]

en este ejemplo ``np.polyfit()`` devuelve la lista de parámetros p de la recta, por lo que el modelo lineal :math:`f(x) = ax + b` de nuestros datos será:

:math:`y(x) = p_0 x + p_1 = 2.7x + 9.94`

Ahora podemos dibujar los datos experimentales y la recta ajustada:

.. code:: ipython3

    from matplotlib import pyplot as plt

    # Valores de y calculados del ajuste
    y_ajuste = p[0]*x + p[1]

    # Dibujamos los datos experimentales
    p_datos, = plt.plot(x, y, 'b.')
    # Dibujamos la recta de ajuste
    p_ajuste, = plt.plot(x, y_ajuste, 'r-')

    plt.title('Ajuste lineal por minimos cuadrados')

    plt.xlabel('Eje X')
    plt.ylabel('Eje Y')

    plt.legend(('Datos experimentales', 'Ajuste lineal'), loc="upper left")
    plt.show()


.. image:: img/polyfit.png
  :align: center
  :width: 14cm

Como se ve en este ejemplo, la salida por defecto de ``np.polyfit()`` es un array con los parámetros del ajuste. Sin embargo, si se pide una salida detalla con el parámetro full=True (por defecto ``full=False``), el resultado es una tupla con el array de parámetros, el residuo, el rango, los valores singulares y la condición relativa. Nos interesa especialmente el residuo del ajuste, que es la suma cuadrática de todos los resíduos :math:`\sum^n_{i=1} \lvert y_i - f(x_i)\rvert^2`. Para el ejemplo anterior tendríamos lo siguiente:

.. code:: ipython3

    # Ajuste a una recta, con salida completa
    resultado = np.polyfit(x, y, 1, full=True)

    print(resultado)
    """ Imprime tupla
    (array([ 2.7 ,  9.94]),                 # Parámetros del ajuste
     array([ 0.472]),                       # Suma de residuos
     2,                                     # Rango de la matriz del sistema
     array([ 2.52697826,  0.69955764]),     # Valores singulares
     1.1102230246251565e-15)                # rcond
     """


Ajuste de funciones generales
-----------------------------

Si queremos hacer un ajuste general, no necesariamente polinómico, debemos usar alguno de los métodos del paquete :mod:`optimize`, que contiene varios optimizadores locales y globales. El más común es ``leastsq`` que, al ser un optimizador, hay que definir previamente una **función residuo** que es la que realmente se va minimizar.


.. code:: ipython3

    import numpy as np
    from matplotlib import pyplot as plt
    from scipy.optimize import leastsq

    # Datos de laboratorio
    datos_y = np.array([ 2.9, 6.1, 10.9, 12.8, 19.2])
    datos_x = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

    # Función para calcular los residuos, donde
    # se calcula (datos - modelo)
    def residuos(p, y, x):
        error  = y - (p[0]*x + p[1])
        return error

    # Parámetros iniciales estimados
    # y = p0[0]*x  + p0[0]

    p0 = [2.0, 0.0]

    # Hacemos  el ajuste por minimos cuadrados con leastsq(). El primer parámetro
    # es la funcion de residuos, luego los parámetro iniciales y una tupla con los
    # argumentos de la funcion de residuos, en este caso, datos_y y datos_x en
    # ese orden, porque así se definió la función de error
    ajuste = leastsq(residuos, p0, args=(datos_y, datos_x))

    # El resultado es una lista, cuyo primer elemento es otra
    # lista con los parámetros del ajuste
    print(ajuste[0])
    # array([ 3.93, -1.41])

Veamos otro ejemplo para ajustar una función seno:

.. code:: ipython3

    import numpy as np
    from matplotlib import pyplot as plt
    from scipy.optimize import leastsq
    from scipy import random

    # Generamos unos datos artificiales para hacer el ejemplo
    # A datos_y se le añade "ruido" que simula error de
    # medida, añadiendole un valor aleatorio
    datos_x = np.arange(0, 0.1, 0.005)
    A, k, theta = 10.0, 33.3, np.pi/5.0
    y_real = A*np.sin(2*np.pi*k*datos_x + theta)
    datos_y = y_real + 2*random.randn(len(datos_x))

    # Ahora se trata de ajustar estos datos una función
    # modelo tipo senoidal A*sin(2*pi*k*x + theta)

    # Defino la funcion de residuos
    def residuos(p, y, x):
        A, k, theta = p
        error = np.abs(y - A*np.sin(2*np.pi*k*x + theta))
        return error

    # Parámetros iniciales
    p0 = [8.0, 40.0, np.pi/3]

    # hacemos  el ajuste por minimos cuadrados
    ajuste = leastsq(residuos, p0, args=(datos_y, datos_x))

    # El resultado es una lista, cuyo primer elemento es otra
    # lista con los parámetros del ajuste.
    print(ajuste[0])

    # Ahora muestro los datos y el ajuste gráficamente
    plt.plot(datos_x, datos_y, 'o')  # datos

    # Defino la funcion modelo, para representarla gráficamente
    def funcion(x, p):
        return p[0]*np.sin(2*np.pi*p[1]*x + p[2])

    # genero datos a partir del modelo para representarlo
    x1 = np.arange(0, datos_x.max(), 0.001)  # array con muchos puntos de x
    y1 = funcion(x1, ajuste[0])           # valor de la funcion modelo en los x

    plt.plot(x1, y1, 'r-')
    plt.xlabel('Eje X')
    plt.ylabel('Eje Y')
    plt.title('Ajuste de funcion seno con leastsq')
    plt.legend(('Datos', 'Ajuste lineal'))
    plt.show()

.. image:: img/leastsq.png
  :align: center
  :width: 14cm

Este ejemplo es bastante alaborado porque hemos usado un optimizador general para hacer un ajuste, pero podemos usar :func:`curve_fit` para ahorrarnos la función residuo.

La anterior es una manera artesanal de hacer el ajuste, al construir la función de error. Para un ajuste de datos

.. code:: ipython3

    from scipy.optimize import curve_fit

    curve_fit(funcion, datos_x, datos_y, p0)

    (array([ -9.787095  ,  32.91201348,  -2.3390355 ]),
     array([[ 0.20148401, -0.00715614,  0.00215931],
           [-0.00715614,  0.07184634, -0.02241144],
           [ 0.00215931, -0.02241144,  0.00925902]]))


Vamo a probar esto con una ley de decaimiento exponencial:

.. code:: ipython3

    # Cargo los datos experimentales a un array
    # Los datos están delimitados por ";" y uso unpack
    # para que ponga primero las columas y pueda desenpaquetarlas
    # en  variables

    tiempo, masa, error = np.loadtxt("medidas_radio.txt", delimiter=";", unpack=True)

    # Compruebo como ser ven los datos
    plot(tiempo, masa, '.')

    def decaimiento(x, a, b):
       """Ley de decaimiento exponencial"""
       return a * exp(-b*x)

    # Ajuste de los datos
    # curve_fit devuelve dos variables, los parámetros del ajuste y
    # la matriz de covarianza
    popt, pcov = curve_fit(decaimiento, tiempo, masa)

    # Ahora creo una curva teórica a partir del modelo ajustado
    times = np.arange(0, 100, 0.1)
    model = decaimiento(times, *popt)

    plt.plot(times, model, '-r')

    plt.legend(('Medidas', 'Ajuste'))
    plt.xlabel("Tiempo (s)")
    plt.ylabel("Masa (g)")
    plt.title("Ajuste de decaimiento exponencial")

    # Guardo la grafica
    plt.savefig("ajuste.png")

    # Para mostrar la gráfica por pantalla
    plt.show()

.. image:: img/ajuste-exp.png
  :align: center
  :width: 14cm


Interpolación de datos
----------------------

El módulo :mod:`interpolate` contiene rutinas de interpolación basadas en la conocida librería FITPACK en Fortran; resulta muy útil en partes de datos donde no hay suficientes medidas. Los interpoladores funcionan creando una función interpoladora de orden predefinido, usando los datos de medidas. Luego se aplica esta función de interpolación a un array de datos más denso que la muestre. Consideremos los datos senoidales que vimos antes:

::

    from scipy.interpolate import interp1d

    interpolador_lineal = interp1d(datos_x, datos_y)
    interpolador_cubico = interp1d(datos_x, datos_y, kind='cubic')

    x_inter = linspace(0.01, 0.09, 500)
    y_inter_l = interpolador_lineal(x_inter)
    y_inter_c = interpolador_cubico(x_inter)

    plt.plot(datos_x, datos_y, 'ok', label="Datos")
    plt.plot(x_inter, y_inter_l, 'r', label="Interpolación lineal")
    plt.plot(x_inter, y_inter_c, 'y', label="Interpolación cúbico")

    # Si usamos InterpolatedUnivariateSpline podemos interpolar fuera
    # del rango de los datos
    from scipy.interpolate import InterpolatedUnivariateSpline

    # Array de valores más denso que los datos originales, pero
    # dentro del rango de los datos
    x_inter = np.linspace(-0.01, 0.11, 500)

    interpolador = InterpolatedUnivariateSpline(datos_x, datos_y, k=2)
    y_inter_u = interpolador(x_inter)

    plot(x_inter, y_inter_u, 'k', lw=0.5, label="Interpolador univariable")

    xlim(-0.003, 0.102)
    ylim(-14, 15)
    legend()

Aunque este es el método más usado para interpolar, tiene el problema que no permite interpolar puntos fuera del rango de las medidas. Esto se puede resolver con la función ``InterpolatedUnivariateSpline()``


.. image:: img/interp_univar.png
  :align: center
  :width: 14cm


Integración numérica
--------------------

Intentemos integrar numéricamente la integral:

:math:`\int_{-2}^{2} 2*\exp(\frac{-x^2}{5}) ~dx`

.. code:: ipython3

    def func1(x):
        return 2.0*exp(-x**2/5.0)

    # integración entre -2 y +2
    int1, err1 = integrate.quad(func1, -2, +2)
    print(int1,err1)
    (6.294530963693763, 6.9883332051087914e-14)

Si ahora hacemos la misma integral pero desde :math:`-\infty` hasta :math:`+\infty` obtendremos toda el área bajo la curva definida en el integrando:

.. code:: ipython3

    # integración entre -infinito y +infinito
    In [10]: int2, err2 = integrate.quad(func1, -Inf, +Inf)
    # y el resultado obtenido es:
    In [11]: print(int1, err1)
    (7.9266545952120211, 7.5246691415403668e-09)


Manipulación de arrays 2D: imágenes
-----------------------------------

El módulo :mod:`ndimage` ofrece manipulación y filtrado básicos de imágenes, dados como array de numpy.

.. code:: ipython3

    In [20]: from scipy import misc
    In [21]: from scipy import ndimage as nd

    In [22]: img = misc.face(gray=True)
    In [23]: shifted_img = nd.shift(img, (50, 50))
    In [24]: shifted_img2 = nd.shift(img, (50, 50), mode='nearest')
    In [25]: rotated_img = nd.rotate(img, 30)

    In [26]: cropped_img = img[50:-50, 50:-50]
    In [27]: zoomed_img = nd.zoom(img, 2)  # Interpola la imagen

    In [28]: zoomed_img.shape
    # (1536, 2048)


.. code:: ipython3

    In [30]: noisy_img = np.copy(img).astype(np.float)

    In [31]: noisy_img += img.std()*0.5*np.random.standard_normal(img.shape)

    In [32]: blurred_img = nd.gaussian_filter(noisy_img, sigma=3)
    In [33]: median_img = nd.median_filter(blurred_img, size=5)

    In [34]: from scipy import signal
    In [35]: wiener_img = signal.wiener(blurred_img, (5,5))

    In [40]: plt.subplot(221)
    In [41]: plt.imshow(img)

    In [42]: plt.subplot(222)
    In [43]: plt.imshow(blurred_img)

    In [44]: plt.subplot(223)
    In [45]: plt.imshow(median_img)

    In [46]: plt.subplot(224)
    In [47]: plt.imshow(wiener_img)


.. image:: img/ndimage_filters.png
  :align: center
  :width: 14cm


Para más utilidades de manipulación de imágenes más sofisticadas conviene usar
http://scikit-image.org o http://opencv.org


Módulo de constantes físicas
----------------------------

Scipy contiene un práctico módulo de constantes físicas.

.. code:: ipython3

    In [1]: from scipy import constants as C

    In [2]: C.c  # Velocidad de la luz en m/s
    Out[2]: 299792458.0

    In [4]: C.e  # Carga del electrón
    Out[4]: 1.602176565e-19

    In [6]: C.atmosphere
    Out[7]: 101325.0

    In [7]: C.mmHg
    Out[7]: 133.32236842105263

    In [8]: C.Julian_year
    Out[8]: 31557600.0

    In [9]: C.Avogadro
    Out[9]: 6.02214129e+23

    In [10]: C.parsec
    Out[10]: 3.0856775813057292e+16

    In [11]: C.Stefan_Boltzmann
    Out[11]: 5.670373e-08

    In [13]: C.convert_temperature(np.array([0, 100.0]), 'Celsius', 'Kelvin')  # Conversor temps
    Out[13]: array([ 273.15,  373.15])

    In [14]: C.day  # Dia en segundos
    Out[14]: 86400.0

    In [15]: C.pico  # Prefijos del SI
    Out[15]: 1e-12

    In [16]: C.oz
    Out[16]: 0.028349523124999998

    In [17]: # Constantes físicas de CODATA 2014
    In [18]: C.find('atm')
    Out[18]: ['standard atmosphere']
    In [19]: C.physical_constants['standard atmosphere']
    Out[19]: (101325.0, 'Pa', 0.0)


Ejercicios
----------

#. El fichero de texto medidas_PV_He.txt posee medidas de presión y volumen para 0.1 mol de helio, que se comprime sistemáticamente a temperatura constante. Este experimento se realiza a tres temperaturas distintas. Suponiendo que el gas se comporta idealmente y por tanto que se verifica que *PV=nRT*, representar los datos y realizar un ajuste lineal P(V) para cada temperatura. ¿Cuánto vale la constante de gases ideales según el experimento?

#. Para una muestra que contiene 10g de yodo 131 (semivida de 8 días), se hacen diariamente cinco medidas independientes a lo largo de 60 días. Esas medidas están en el fichero medidas_decaimiento_yodo131b.txt, donde cada fila corresponde a cada una de las 5 medidas realizadas diariamente en gramos de yodo que queda. Representar en un gráfico con puntos las cinco medidas con colores distintos para cada una y ajustar a cada una la curva teórica de decaimiento. Imprimir por pantalla los parámetros de cada uno de los cinco ajustes.

#. Hacer una función que cree una máscara circular (todos los valores son cero, excepto una zona circular de radio ``R`` y centrada en un punto arbitrario) con un array 2D para la imagen del castor usada anteriormente ``face = misc.face(gray=True)``.