Scientific Programming in Python

Topic 5: Accelerating Python with Cython: Writting C in Python

Notebook created by Martín Villanueva - martin.villanueva@usm.cl - DI UTFSM - May2017.



In [2]:

    
import numba
import numpy as np
from math import sqrt

%load_ext Cython









    



The Cython extension is already loaded. To reload it, use:
  %reload_ext Cython

La distancia de Hausdorff nuevamente...

En esta actividad volveremos a implementar la distancia/métrica de Hausdorff, pero ahora utilizando Cython.

La métrica de Hausdorff corresponde a un métrica o distancia ocupada para medir cuán disímiles son dos subconjuntos dados.

Esta tiene muchas aplicaciones, en particular para comparar el parecido entre imágenes. En el caso en donde los conjuntos son arreglos bidimensionales, la definición es la siguiente:

Sean $X \in \mathbb{R}^{m \times 3}$ e $Y \in \mathbb{R}^{n \times 3}$ dos matrices, la métrica/distancia de Hausdorff sobre sobre estas como:

$$ d_H(X,Y) = \max \left(\ \max_{i\leq m} \min_{j \leq n} d(X[i],Y[j]), \ \max_{j\leq n} \min_{i \leq m} d(Y[j],X[i]) \ \right) $$

donde $d$ es la distancia Euclideana clásica. ($X[i]$ indíca la i-ésima fila de X).

Ilustración unidimensional: Distancia entre funciones.

A continuación se le proveen 3 funciones que implementan tal métrica, usando Numba.



In [ ]:

    
@numba.jit('float64 (float64[:], float64[:])')
def metric_numba(x, y):
    """
    standard Euclidean distance
    """
    ret = x-y
    ret *= ret
    return np.sqrt(ret).sum()


@numba.jit('float64 (float64[:], float64[:,:])', nopython=True)
def inf_dist_numba(x, Y):
    """
    inf distance between row x and array Y
    """
    m = Y.shape[0]
    inf = np.inf
    
    for i in range(m):
        dist = metric_numba(x, Y[i])
        if dist < inf:
            inf = dist
    return inf

@numba.jit('float64 (float64[:,:], float64[:,:])', nopython=True)
def hausdorff_numba(X, Y):
    """
    Hausdorff distance between arrays X and Y
    """
    m = X.shape[0]
    n = Y.shape[0]
    sup1 = -1.
    sup2 = -1.
    
    for i in range(m):
        inf1 = inf_dist_numba(X[i], Y)
        if inf1 > sup1:
            sup1 = inf1
    for i in range(n):
        inf2 = inf_dist_numba(Y[i], X)
        if inf2 > sup2:
            sup2 = inf2
            
    return max(sup1, sup2)

Se solicita que realice lo siguiente:

Escribir el equivalente Cython de las tres funciones anteriores, ocupando todas las optimizaciones posibles: Compiler directives, Memory Views, Inline Functions, Pure C functions o cualquier otra optimización que usted considere conveniente.
Cree 10 arreglos $X,Y$ aleatorios, con cantidad creciente de filas, y realice análsis de tiempos de ejecuciones de las versiones Numba y Cython de las funciontes anteriores sobre estos arreglos.
Concluya.



In [ ]: