sexta-feira, 13 de maio de 2016

Espaços de hipótese e algoritmos de aprendizagem

É comum durante o aprendizado de Machine Learning confundir os papéis do espaço espaço de hipóteses e do algoritmo de aprendizagem no determinação da hipótese encontrada. Para tirar esta dúvida resolvi escrever um pequeno exemplo que mostra dois algoritmos simples (Perceptron e Regressão Linear) explorando espaços de hipóteses completamente distintos.

Primeiramente, o espaço de hipóteses $\mathcal{H}$ contém a classe de funções consideradas por um algoritmo de aprendizagem, que irá escolher a função $g \in \mathcal{H}$ baseado em uma função de erro calculada em um conjunto de treinamento. A junção destas três partes é chamada de modelo de aprendizagem. Nesta aula sobre VC dimension, o prof. Mostafa argumenta que a capacidade de generalização de um modelo depende somente do espaço de hipóteses. Isto pode parecer contra-intuitivo, pois algoritmos diferentes podem resultar em soluções com desempenho muito diferente. Porém, como iremos ver nos exemplos abaixo, é o espaço de hipóteses que define a complexidade de um modelo (e, consequentemente sua capacidade de generalização) e o resultado apresentado pelo prof. Mostafa é um limitante inferior para o erro. Como não conhecemos a função objetivo $f$, que pode ser tão complexa quanto possível, existirá sempre uma parte de $f$ que não poderá ser capturada usando funções no espaço $\mathcal{H}$. Na aula sobre Bias-Variance, esta "margem" é definida como o bias.

Caso linearmente separável

Neste exemplo mostramos como diferentes algoritmos de aprendizado encontram diferentes hipóteses $g$ mesmo que estejam considerando o mesmo espaço de hipóteses $\mathcal{H} = \{f_w: f_w(x) = sign(w^T x) \}$.

Primeiramente, vamos amostrar duas populações de dois retângulos no plano.

In [9]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import scipy.stats
import sklearn
import sklearn.linear_model

dados1 = sp.stats.uniform.rvs(size=100 * 2).reshape((100, 2))
dados1[:,0] = dados1[:,0] * 3 + 5
dados1[:,1] = dados1[:,1] * 0.5 - 2

dados2 = sp.stats.uniform.rvs(size=100 * 2).reshape((100, 2))
dados2[:,0] = dados2[:,0] * 2 
dados2[:,1] = dados2[:,1] * 5 + 2

plt.xlim((-2, 10))
plt.ylim((-10, 15))
plt.plot(dados1[:,0], dados1[:,1], 'ro')
plt.plot(dados2[:,0], dados2[:,1], 'bo')
Out[9]:
[<matplotlib.lines.Line2D at 0x2179590>]
In [11]:
def plot_decision(model, color, caption, min_=-2, max_=10):
    xx = np.linspace(min_, max_, 100)
    # Graças ao scikit learn e ao python, podemos passar  
    # model como parâmetro de plot_decision
    coef = model.coef_.reshape(-1)
    a = -coef[0] / coef[1]
    yy = xx * a - model.intercept_ / coef[1]
    plt.plot(xx, yy, color, label=caption)

Para selecionar a hipótese, usaremos dois algoritmos diferentes: perceptron e regressão linear. Em seguida, chamamos a função plot_decision para desenhar a superfície de separação. Escolhemos estes dois algoritmos pois ambos são algoritmos considerados "simples" por estimarem somente decisões lineares em $w$.

In [12]:
def treina_e_plota(X1, X2, min_=-2, max_=10, intercept=True):
    X = np.r_[X1, X2]
    y = np.r_[np.ones(X1.shape[0]), -np.ones(X2.shape[0])]

    perc = sklearn.linear_model.Perceptron(fit_intercept=intercept, n_iter=500)
    perc.fit(X, y)

    linreg = sklearn.linear_model.LinearRegression(fit_intercept=intercept)
    linreg.fit(X, y)

    plt.figure(figsize=(10, 10))
    plt.plot(X1[:,0], X1[:,1], 'ro')
    plt.plot(X2[:,0], X2[:,1], 'bo')
    plot_decision(perc, 'g', 'Perceptron', min_, max_)
    plot_decision(linreg, 'y', 'Regressão linear', min_, max_)
    plt.legend()
    plt.title('Caso linearmente separável')
    print('Score Perceptron', perc.score(X, y))
    print('Score LinRegression', np.sum(np.sign(linreg.predict(X)) == y) / X.shape[0])
    return perc, linreg

# Chama os algoritmos de treinamento.

treina_e_plota(dados1, dados2)
Score perceptorn 1.0
Score LinRegression 1.0
Out[12]:
(Perceptron(alpha=0.0001, class_weight=None, eta0=1.0, fit_intercept=True,
       n_iter=500, n_jobs=1, penalty=None, random_state=0, shuffle=True,
       verbose=0, warm_start=False),
 LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False))

Como podemos ver, cada algoritmo de aprendizagem selecionou uma hipótese diferente, mas ambas separam perfeitamente as duas classes.

Caso não linearmente separável

Como visto na aula sobre modelos lineares, os modelos estimados pela Regressão Linear e pelo Perceptron são lineares em $w$, mas não necessariamente em $x$. Para estimar modelos não-lineares em $x$ criamos uma versão transformada $z$ de todas as entradas e estimamos $f(z) = w^T z$. É importante notar que ao fazer essa transformação o espaço de hipóteses muda. Iremos exemplificar esta transformação estimando uma superfície de decisão para os dados abaixo.

In [15]:
def generate_circle(loc, scale, A, B, inside, size):
    X = np.zeros((size, 2))
    center = loc + scale/2
    #print(center)
    while size > 0:
        point = scale * sp.stats.uniform.rvs(size=2) + loc
        #print(point, (point[0] - center[0])**2/A**2 + (point[1] - center[1])**2/B**2)
        if inside and (point[0] - center[0])**2/A**2 + (point[1] - center[1])**2/B**2 <= 1:
            X[size-1] = point
            size -= 1
        elif not inside and (point[0] - center[0])**2/A**2 + (point[1] - center[1])**2/B**2 > 1:
            X[size-1] = point
            size -= 1
        #size -= 1
    return X

plt.figure(figsize=(8, 8))

dados1 = generate_circle(np.array([-5, -5]), np.array([10, 10]), 5, 3, True, 500)
dados2 = generate_circle(np.array([-5, -5]), np.array([10, 10]), 5, 3, False, 500)
plt.plot(dados1[:,0], dados1[:,1], 'ro')
plt.plot(dados2[:,0], dados2[:,1], 'bo')
Out[15]:
[<matplotlib.lines.Line2D at 0x5203a90>]

Vamos explorar o espaço de hipóteses das funções de decisão elípticas $\mathcal{H} = \{ f : f(x; A, B) = sign(\frac{(x-c_x)^2}{A^2} + \frac{y - c_y)^2}{B^2} - 1) \}$ e, portanto, estimaremos o seguinte modelo linear:

$$ f(x; \alpha, \beta, \gamma) = sign( \alpha (x_1 - \overline{x_1})^2 + \beta (x_2 - \overline x_2)^2 + \gamma) = sign(\alpha z_1 + \beta z_2 + \gamma), $$

com

$$ z_i = (x_i - \overline x_i)^2 $$

Note que, dados $\alpha, \beta$ e $\gamma$, conseguimos recuperar os parâmetros $A$ e $B$ da elipse.

$$ A = \sqrt{\frac{-\gamma}{\alpha}}, B = \sqrt{\frac{-\gamma}{\beta}}. $$
In [16]:
media = dados1.mean(axis=0)

dados1_X = (dados1 - media)**2
dados2_X = (dados2 - media)**2

perc, linreg = treina_e_plota(dados1_X, dados2_X, -3, 30)
Score perceptorn 0.972
Score LinRegression 0.925

Podemos ver no gráfico acima uma diferença entre ambos algoritmos: como para a Regressão Linear a distância entre um ponto classificado incorretamente e o hiperplano faz diferença, a reta amarela é puxada para cima. Isto ocorre pois os pontos azuis estão mais espalhados verticalmente.

Note que a decisão linear estimada acima é feita no espaço transformado, não no espaço original. Veja abaixo as duas decisões acima plotadas no espaço original de características.

In [17]:
print('Perceptron', perc.coef_, perc.intercept_)
print('Lin Regression', linreg.coef_, linreg.intercept_)
A_perc, B_perc = np.sqrt(-perc.intercept_/perc.coef_)[0]
A_lin, B_lin = np.sqrt(-linreg.intercept_/linreg.coef_)

from matplotlib.patches import Ellipse

plt.figure(figsize=(8,8))
ax = plt.gca()
plt.xlim((-10, 10))
plt.ylim((-10, 10))


plt.plot(dados1[:,0], dados1[:,1], 'ro')
plt.plot(dados2[:,0], dados2[:,1], 'bo')
e_perc = Ellipse(xy=media, width=2*A_perc, height=2*B_perc)
e_perc.set_facecolor('b')
e_perc.set_alpha(0.5)
ax.add_patch(e_perc)
e_lin = Ellipse(xy=media, width=2*A_lin, height=2*B_lin)
e_lin.set_facecolor('y')
e_lin.set_alpha(0.5)
ax.add_patch(e_lin)
Perceptron [[ -63.16696553 -167.94716625]] [ 1476.]
Lin Regression [-0.03890952 -0.10278301] 1.13912219107
Out[17]:
<matplotlib.patches.Ellipse at 0x12fdf9d0>

Espero que este exemplo tenha sido significativo e que a diferença entre algoritmo de aprendizagem e espaço de hipóteses tenha sido esclarescida. Como sempre, quaisquer sugestões ou críticas são benvindas nos comentários.

quarta-feira, 27 de janeiro de 2016

Programando em GPUs usando Numba

A biblioteca Numba pode ser usada para compilar código Python que opera em arrays para código de máquina com excelente performance. Podemos utilizar Numba também para compilar funções nopython para rodar em GPUs NVIDIA. Dependendo do tipo de processamento a ser executado, a utilização de GPUs pode ser muito vantajosa e resultar em ganhos de desempenho de 10-100 vezes em relação à codigo otimizado rodando em CPUs.

Uma CPU moderna possui, em geral, entre 2 e 8 cores. Ao executar um programa multi-thread (ou usando multiprocessing, como já fizemos antes), cada core da CPU executa uma thread (ou processo) de maneira completamente independente dos outros processadores. GPUs, por sua vez, possuem centenas ou milhares de cores que estão organizados de maneira hierárquica e não são independentes uns dos outros. Os cores de uma GPU estão agrupados em Streaming Processes (SMs), que funcionam em um modo batizado de SIMT (Single Instruction, Multiple Threads). Basicamente, cada SM funciona de maneira independente dos outros SMs e contém um certo número de cores que trabalham de maneira sincronizada. A mesma instrução é executada em todos os cores de um SM, mas utilizando dados diferentes. Por esta razão, a NVIDIA batizou este tipo de processamento de data-parallel processing.

Para executar uma função (kernel) na GPU é necessário dividir o problema a ser tradato em um grid de blocos de threads, como na figura abaixo. Cada bloco executará em um SM seu conjunto de threads, fazendo com que seja possível executar vários blocos ao mesmo tempo. Dentro de um bloco, todas as threads executam em paralelo, mas de modo sincronizado. Um máximo de 1024 threads podem existir em cada bloco e cada thread tem acesso tanto à sua posição no bloco como à posição de seu bloco no grid.

Para, por exemplo, processar uma imagem de tamanho $512\times 512$, poderíamos criar um grid de tamanho $16\times16$ em que cada bloco contém 32 threads. Cada thread determina o valor de um pixel. O grid faz o papel de um duplo for que percorre todos os pontos da imagem.

Como podemos ver na imagem (1), uma outra diferença desta nova arquitetura é que a GPU possui sua própria memória, que pode estar inclusive fisicamente separada da memória da CPU. Antes de lançar um kernel é necessário copiar os dados da memória da CPU para a GPU e após o processamento é necessário copiar o resultado de volta. Logo, processar um volume pequeno de dados pode não ser vantajoso, pois a maior parte do tempo será gasta copiando os dados.

Infelizmente, não podemos executar código Python arbitrário em GPUs. Ao contrário da utilização de Numba em CPUs, que gera código em modo objeto se não for possível otimizar todas as operações, a compilação para GPUs só pode ser feita em modo nopython (veja algumas das restrições no post anterior sobre Numba). Usaremos como exemplo uma implementação do gradiente morfológico feita em GPU.

Primeiramente, é necessário instalar os drivers mais atuais da GPU usada (disponíveis aqui) e o CUDA SDK. Um pacote binário com o CUDA SDK chamado cudatoolkit pode ser instalado via conda.

In [1]:
import numba.cuda
import numpy as np
import math

@numba.cuda.jit
def grad_morpho(img, out):
    pi, pj = numba.cuda.grid(2)
    
    if pi > 2 and pi < img.shape[0]-4 and pj > 2 and pj < img.shape[1]-4:
        minv = 255; maxv = 0
        for i in range(-3, 4):
           for j in range(-3, 4):
                pix = img[pi+i, pj+j]
                if pix < minv:
                    minv = pix
                if pix > maxv:
                    maxv = pix
        
        out[pi, pj] = maxv - minv

def versao_numba_gpu(img):
    out = img.copy()
    img2 = img.copy()
    threads_per_block = (16, 16)
    nblocks = (math.ceil(img.shape[0] / threads_per_block[0]),
               math.ceil(img.shape[1] / threads_per_block[1]))
    grad_morpho[nblocks, threads_per_block](img2, out)           
    return out

Na função call_grad_morpho definimos o número de threads por bloco e o número de blocos no grid e chamamos a função grad_morho, que é executada na GPU de maneira paralela. Isto fica explícito no código ao notarmos que não existe um duplo for que percorre os pontos da imagem. Cada que thread executa grad_morpho sabe a sua dentro de seu bloco e sua posição absoluta no grid. Uma maneira conveniente de obter esta posição é chamar a fução numba.cuda.grid(ndim), que calcula os índices para um grid de ndim=1,2,3 dimensões.

Outro ponto importante é que nem sempre o grid coincide com as dimensões dos dados, então é muito importante checar se o ponto pi, pj que recebemos não ultrapassa as dimensões das imagens tratadas.

Como comparação, colocamos abaixo o código otimizado Cython que fizemos anteriormente (colado abaixo) e executamos ambas funções em uma imagem tamanho $4096\times4096$.

In [2]:
%load_ext cython
In [3]:
%%cython
import cython

@cython.boundscheck(False)
@cython.nonecheck(False)
cdef void morpho_gradient_cython(long[:,:] img, long[:,:] out, int pi, int pj) nogil:
    cdef int i, j, minv, maxv
    minv = 255
    maxv = 0
    
    for i in range(-3, 4):
        for j in range(-3, 4):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv


@cython.boundscheck(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef void versao_cython(long[:,:] img, long[:,:] out) nogil:
    cdef int i, j, k, w, h
    h = img.shape[0]; w = img.shape[1]
    for i in range(3, h-4):
        for j in range(3, w-4):
            morpho_gradient_cython(img, out, i, j)
In [4]:
import scipy
import scipy.misc
import numpy as np

img_temp = scipy.misc.lena()
img2 = np.c_[img_temp, img_temp]
img3 = np.r_[img2, img2]
img4 = np.c_[img3, img3]
img = np.r_[img4, img4]
img = np.r_[img, img]
img = np.c_[img, img]
out = img.copy()
print('Versão Cython')
%timeit versao_cython(img, out)
print('Versão Numba GPU')
%timeit versao_numba_gpu(img)

versao_cython(img, out)
out2 = versao_numba_gpu(img)
print('Imagens iguais?', np.all(out == out2))
Versão Cython
1 loops, best of 3: 1.95 s per loop
Versão Numba GPU
1 loops, best of 3: 253 ms per loop
Imagens iguais? True

Como podemos ver, a diferença entre os tempos de execução é muito grande. Utilizar Numba para executar código em GPU pode ser muito eficiente se o problem for altamente paralelizável.

sábado, 2 de janeiro de 2016

Porque escrevo em português

Tanto as Ciências como a área de Tecnologia (C&T) são campos muito globalizados em que a grande maioria da produção técnica e inovativa é feita em inglês. Neste contexto, escrever um blog em inglês parece uma escolha natural, tanto pelo audiência maior quanto pelo possível reconhecimento obtido. Apesar disto, decidi escrever em português por diversas razões que apresentarei neste post.

A principal razão para escrever em português é o estímulo da produção de cultura nacional. Acho muito importante que um país possua uma população que seja criativa e que produza conhecimento relevante para sua realidade. A transformação de um cidadão consumidor de conhecimento para um cidadão produtor de conhecimento é importantíssima para o desenvolvimento intelectual de um país e a disponibilidade de material em nossa língua é um pré requisito para isto. Espero contribuir para um cíclo virtuoso em que as pessoas busquem e leiam textos em português, interajam com estas publicações e, eventualmente, produzam novos conhecimentos a partir disto.

É claro que já existe muito material bom em inglês, mas exigir o domínio de um idioma estrangeiro como pré-requisito para estudar C&T é colocar uma barreira de entrada muito grande para iniciar nestas áreas. Aprender um idioma pode levar anos e mesmo pessoas que se sentem confortáveis com línguas estrangeiras frequentemente preferem (e buscam) informações em sua língua nativa. Negligenciar nossa própria língua significa excluir uma porção considerável de nossa população. É verdade que o inglês é sim muito importante no mundo atual e que aqueles que conhecem esta língua podem ter diversas vantagens, mas isto não significa que o aprendizado dos fundamentos destas áreas deva ser feito em inglês. Se aprender C&T já é desafiador, fazê-lo em um idioma não familiar é mais complicado ainda.

Por último, por vezes escrevo sobre tecnologias que aprendi durante meu trabalho e que muitas vezes possuem pouco ou nenhum material em português. Outras vezes, publico resumos de artigos científicos ou de algum assunto relacionado à minha área de trabalho atual (aprendizagem computacional). Em geral, as ideias ou tecnologias apresentadas já existem, então o foco é na apresentação do conteúdo. Logo, não existe nenhum prejuízo em escrever em português, pois nenhuma ideia inovadora está sendo perdida por causa de barreiras de linguagem.

A junção dos três argumentos me fez decidir por escrever em português. Uma produção inovadora que pode ter impacto significativo deve, na minha opinião, ser escrita em inglês, mas este blog não se encaixa nesta definição. Este blog busca a divulgação do conhecimento e um possível aumento na qualidade da produção técnica nacional. Tenho plena consciência que o meu impacto é (e provavelmente continuará sendo) pequeno, mas gosto de acreditar que se outros fizerem o mesmo a junção destes esforços pode resultar em algo de grande impacto. No pior dos casos, pelo menos uma pessoa (eu) está aprendendo e se beneficiando de tudo o que está escrito aqui e isto neste momento já é o suficiente.

sexta-feira, 25 de dezembro de 2015

Acelerando programas usando Cython - parte 4

Série de posts sobre Cython:

  1. Introdução ao desenvolvimento usando Cython
  2. Distribuição de extensões em Cython
  3. Análise de eficiência de operações
  4. Paralelismo e NoGIL

Na parte 4 da série sobre Cython falaremos sobre funções que liberam completamente o interpretador e podem ser paralelizadas usando OpenMP. No texto anterior vimos que chamadas de funções cdef e cpdef são muito rápidas pois elas podem ser feitas "por fora" do interpretador. Também vimos quais operaçõe podem ser feitas sem tocar no interpretador. Podemos levar esta idéia mais além e criar funções (ou contextos) nogil, que liberam explicitamente o GIL e permitem a execução simultânea de várias threads em Python. Já vimos nos textos sobre multiprocessing (partes 1 e 2) que a utilização de vários processadores pode trazer ganhos de desempenho significativos. Neste texto exploramos algumas situações em que podemos fazer processamento paralelo em Cython.

In [7]:
%load_ext cython
The cython extension is already loaded. To reload it, use:
  %reload_ext cython
In [8]:
%%cython
# exemplo de funcoes nogil e de contextos nogil
cdef void funcao_nogil(int a, float b, char[:] c) nogil:
    cdef double d
    pass
    # toda função nogil precisa ter todos os tipos declarados e só usar operações "brancas"
    with gil:
        pass
        # esta parte requer o GIL, ou seja, é executada SEMPRE de maneira serial.
    
def funcao_gil(a, b, c):
    cdef double d
    with nogil:
        pass
        # aqui só são permitidas operações "brancas" com variáveis tipadas explicitamente.

Funções (e contextos) nogil possuem diversas limitações (parecidas com as do modo nopython que vimos no post sobre Numba). Não é possível interagir com nenhum tipo de objeto em Python, somente com variáveis que possuem tipo declarado e que o tipo está disponível diretamente em C. Isto inclui os tipos C habituais (int, float, double, char, etc) e suas versões memory view (int[:], float[:, :], etc). É importante notar que todos os tipos permitidos acessam a memória diretamente e não passam pela contagem de referências do interpretador.

A maneira mais fácil de checar se um trecho de código pode ser executado como nogil é simplesmente executando o cython -a como temos feito nos últimos posts. Toda linha completamente branca pode ser executada como nogil. Qualquer traço de amarelo indica interação com o interpretador e, portanto, é necessário obter o GIL.

Apresentamos abaixo duas versões do gradiente morfológico que usamos nos últimos posts. Primeiramente, adicionamos um tipo de retorno à função morpho_gradient_cython e a marcamos como nogil, o que significa que ela não interage nenhuma vez com o interpretador e explicitamente o libera ao ser executada. A versão versao_cython_omp usa a construção cython.parallel.prange(link), que executa cada iteração do for em uma thread diferente. Esta construção só funciona em modo nogil.

Para compilar extensões usando OpenMP no ipython notebook é necessário adicionar as flags abaixo na %%cython magic. Se você está usando um setup.py, é necessário adicioanar as opções

In [9]:
%%cython --compile-args=-fopenmp --link-args=-fopenmp --force
import numpy as np
import cython

from cython.parallel cimport prange
cimport openmp


@cython.boundscheck(False)
@cython.nonecheck(False)
cdef void morpho_gradient_cython(long[:,:] img, long[:,:] out, int pi, int pj) nogil:
    cdef int i, j, minv, maxv
    minv = 255
    maxv = 0
    
    for i in range(-3, 4):
        for j in range(-3, 4):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

@cython.boundscheck(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef void versao_cython_omp(long[:,:] img, long[:,:] out) nogil:
    cdef int i, j, k, w, h
    h = img.shape[0]; w = img.shape[1]
    n = (h-2)*(w-2)
    for k in prange(n, nogil=True):
        i = k / w
        j = k % w
        morpho_gradient_cython(img, out, i, j)

@cython.boundscheck(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cpdef void versao_cython(long[:,:] img, long[:,:] out) nogil:
    cdef int i, j, k, w, h
    h = img.shape[0]; w = img.shape[1]
    n = (h-2)*(w-2)
    for k in range(n):
        i = k / w
        j = k % w
        morpho_gradient_cython(img, out, i, j)

Para comparar os resultados, precisamos de imagens maiores para sentir a diferença entre as implementações. Processamento paralelo sempre involve uma quantidade significativa de overhead e o código em Cython já é muito rápido para imagens 512x512, como é o caso da Lena. Logo, "copiamos" a Lena em uma imagem bem maior.

In [10]:
import scipy
import scipy.misc
import numpy as np

img_temp = scipy.misc.lena()
img2 = np.c_[img_temp, img_temp]
img3 = np.r_[img2, img2]
img4 = np.c_[img3, img3]
img = np.r_[img4, img4]
img = np.r_[img, img]
img = np.c_[img, img]
print(img.shape)
(4096, 4096)
In [11]:
out = img.copy()
print('Versão Cython')
%timeit versao_cython(img, out)
print('Versão Cython OMP')
%timeit versao_cython_omp(img, out)
Versão Cython
1 loops, best of 3: 4.9 s per loop
Versão Cython OMP
1 loops, best of 3: 2.3 s per loop

Como podemos ver, existe uma diferença significativa entre as duas implementações. Este é principal caso de uso de paralelismo em Cython: execução de um loop em paralelo em que as iterações são independentes e poderiam ser executadas em qualquer ordem. Os ganhos de desempenho podem ser supreendentes, principalmente em processadores com um grande número de cores (8, 16, 24, ...).

Um dos grandes problemas da utilização de OpenMP em Cython é que é proibida qualquer interação com objetos Python. Porém, dependendo do tipo de tarefa podemos mesmo assim explorar código paralelo usando o módulo multiprocessing, que já exploramos em textos anteriores. Toda função cpdef ou def acelerada usando Cython pode ser chamada normalmente usando o multiprocessing. No exemplo abaixo, aplicamos o gradiente morfológico em um grande número de imagens e comparamos a versão OpenMP com uma versão Cython "regular" (sem OpenMP) mas utilizando multiprocessing.

Este problema possui as mesmas propriedades do anterior: é um conjunto de tarefas independentes que podem ser executadas em qualquer ordem. Porém, desta vez analisamos a situação em que não conseguimos paralelizar o código de cada tarefa individualmente, então executamos várias tarefas em paralelo. Um ponto importante é que não existem benefícios em criar um número maior de processos/threads que o número de processadores disponíveis. Isto pode, inclusive, tornar a execução mais lenta por causa da competição entre os processos. Logo, executar uma função com OpenMP usando multiprocessing provavelmente não é uma boa ideia.

In [12]:
import os
import numpy as np
import scipy as sp
import scipy.ndimage

from multiprocessing.pool import Pool

def le_e_processa_omp(t):
    imgpath = t
    img = sp.ndimage.imread(imgpath, mode='L').astype('long')
    out = img.copy()
    versao_cython_omp(img, out)

def serial_omp(path, n=-1):    
    images = [x for x in os.listdir(path) if x[-3:] == 'png'][:n]
    images = ['%s/%s'%(path, img) for img in images]
    list(map(le_e_processa_omp, images))

def le_e_processa(t):
    imgpath = t
    img = sp.ndimage.imread(imgpath, mode='L').astype('long')
    out = img.copy()
    versao_cython(img, out)

def parallel_thread(path, n=-1):
    images = [x for x in os.listdir(path) if x[-3:] == 'png'][:n]
    images = ['%s/%s'%(path, img) for img in images]
    pool = Pool()
    pool.map(le_e_processa, images)
    pool.close()
    del pool

    
path = '/media/igor/Data1/datasets/staffs/test'
%timeit -r 1 -n 1 serial_omp(path, 50)
%timeit -r 1 -n 1 parallel_thread(path, 50)
1 loops, best of 1: 50.8 s per loop
1 loops, best of 1: 45.3 s per loop

Como podemos ver, o desempenho dos dois testes é bem parecido, mesmo usando uma versão muito mais lenta do gradiente morfológico na versão multiprocessing. Logo, seja usando OpenMP, seja usando multiprocessing, processamento paralelo com Cython pode trazer ganhos significativos de performance.

segunda-feira, 7 de dezembro de 2015

Acelerando Python usando Numba

Implementações em Python puro são frequentemente muito lentas e qualquer tentativa de acelerar estes programas impacta diretamente a produtividade de um cientista que uma Python. Falamos anteriormente de Processamento Paralelo com multiprocessing e Cython como maneiras de acelerar experimentos. Ambas soluções, entretanto, exigem modificações não triviais (e potencialmente extensivas) no código, aumentando sua complexidade.

Neste texto irei apresentar uma alternativa chamada Numba, um compilador Just in Time(JIT) que compila um subconjunto da linguagem Python em código de máquina eficiente. Numba é muito mais limitado que Cython e não é capaz de traduzir todas (a maioria das?) construções Python de maneira eficiente. As funcionalidades suportadas, porém, podem ser aceleradas significativamente sem grande esforço.

Primeiramente, um compilador JIT traduz, durante a execução, código escrito em uma linguagem de alto nível como Python em instruções nativas da CPU em que está sendo executado. Otimizações são aplicadas somente nas partes mais lentas do código. A vantagem disto é que não é necessário compilar nada antes da execução do programa, porém as primeira execuções de uma função "jitted" são mais lentas.

Uma das grandes vantagens de se usar Numba é que o código compilado é Python puro. Não é necessário adicionar nenhum tipo de anotação ou código extra. Tudo é feito automaticamente pelo Numba e o código gerado é bastante rápido em alguns casos específicos. A maneira mais fácil de instalar este pacote é utilizando o Anaconda Python. Um simples conda install numba instala todo o necessário sem nenhuma complicação, mesmo em Windows.

Vamos usar como exemplo nosso já cansado gradiente morfológico. Também manteremos a versão Cython que criamos em um post anterior.

Versão Python

In [1]:
import numpy as np
import scipy as sp
import scipy.ndimage
import time
import matplotlib.pyplot as plt

def morpho_gradient_py(img, out, p):
    pi, pj = p
    minv = 255
    maxv = 0
    
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

def versao_python(img):
    out = np.zeros(img.shape, img.dtype)
    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            morpho_gradient_py(img, out, (i, j))
    return out

Versão Cython

In [2]:
%load_ext cython
In [3]:
%%cython
import numpy as np
import scipy as sp
import scipy.ndimage
import time
import matplotlib.pyplot as plt
import cython

@cython.boundscheck(False)
@cython.nonecheck(False)
cdef morpho_gradient_cython(long[:,:] img, long[:,:] out, p):
    cdef int i, j, pi, pj, minv, maxv
    
    pi, pj = p
    minv = 255
    maxv = 0
    
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

@cython.boundscheck(False)
@cython.nonecheck(False)
cpdef versao_cython(long[:,:] img):
    cdef int i, j, w, h
    cdef long[:,:] out
    h = img.shape[0]; w = img.shape[1]
    out = np.zeros((h, w), np.int)
    for i in range(1, h-1):
        for j in range(1, w-1):
            morpho_gradient_cython(img, out, (i, j))
    return out

Versão Numba

In [4]:
import numba

@numba.jit
def morpho_gradient_numba(img, out, p):
    pi, pj = p
    minv = 255
    maxv = 0
    
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

@numba.jit
def versao_numba(img):
    out = np.zeros(img.shape, img.dtype)
    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            morpho_gradient_numba(img, out, (i, j))
    return out

A versão Numba é basicamente a versão Python com um decorador adicionado às funções. Diferentemente de Cython, não adicionamos nenhum tipo manualmente nem precisamos compilar nosso código explicitamente. Vamos agora aos resultados numéricos.

In [5]:
import scipy

img = scipy.misc.lena()

print('Versão Python')
%timeit versao_python(img)
print('Versão Cython')
%timeit versao_cython(img)
print('Versão Numba')
%timeit versao_numba(img)
Versão Python
1 loops, best of 3: 3.23 s per loop
Versão Cython
10 loops, best of 3: 39.7 ms per loop
Versão Numba
The slowest run took 25.12 times longer than the fastest. This could mean that an intermediate result is being cached 
1 loops, best of 3: 26.6 ms per loop

Ambas as versões aceleradas obtiveram performance muito superior ao programa em Python puro. Surpreendentemente, o código em Numba obteve um desempenho superior ao código em Cython. Isto ocorre, principalmente, pois usamos somente funcionalidades suportadas no modo nativo do Numba, chamado de nopython. Isto inclui, basicamente, manipulação de arrays com Numpy (funções suportadas), chamadas de função com decorador @jit e controle de fluxo padrão de Python (if, for i in range(..), while). Existe suporte limitado para listas e é melhor evitar funções recursivas. Uma lista completa das funcionalidades suportadas pode ser vista na documentação do projeto. Funções nopython podem ser compiladas para código nativo bastante rápido e o próprio processo de compilação não só é invisível para o programador com também é bem rápido, diferentemente de Cython.

Porém, o modo nopython só é ativado se a função só utilizar as funcionalidades permitidas. Se uma única operação não for suportada o Numba ativa o modo object e trata todas as variáveis como objetos Python. Ou seja, apesar de algumas otimizações ainda serem possíveis, o programa será muito mais lento. Não é possível, como em Cython, misturar código que interage com o interpretador e código compilado. Vejamos abaixo um exemplo do mesmo código acima compilado no modo object.

In [6]:
@numba.jit(forceobj=True)
def morpho_gradient_numba_ruim(img, out, p):
    pi, pj = p
    minv = 255
    maxv = 0
    
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

@numba.jit(forceobj=True)
def versao_numba_ruim(img):
    out = np.zeros(img.shape, img.dtype)
    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            morpho_gradient_numba_ruim(img, out, (i, j))
    return out


print('Versão Numba Object')
%timeit versao_numba_ruim(img)
Versão Numba Object
1 loops, best of 3: 18.3 s per loop

Como podemos ver, o resultado pode ser desastroso. Apesar da aparente limitação, Numba é muito útil para o que se propõe a fazer: acelerar código em Python que contém diversos loops e acessa arrays do Numpy. Uma boa maneira de checar se um código poderá ser acelerado é usar o decorador @jit(nopython=True). Esta opção força a compilação em modo nopython e levanta uma exceção se o código não puder ser compilado deste modo.

Mesmo que a função inteira não possa ser compilada, se houver algum loop que só possui operações nopython o Numba pode executar uma otimização chamada loop lifting. Vejamos o exemplo abaixo.

In [7]:
@numba.jit
def morpho_gradient_numba_lift(img, out, p):
    pi, pj = p
    minv = 255
    maxv = 0
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

@numba.jit
def versao_numba_lift(img):
    out = np.zeros(img.shape, np.uint8)
    d = {}

    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            morpho_gradient_numba_lift(img, out, (i, j))

    d[5] = 0

    return out


print('Versão Numba Object')
%timeit versao_numba_lift(img)
Versão Numba Object
The slowest run took 19.17 times longer than the fastest. This could mean that an intermediate result is being cached 
1 loops, best of 3: 25.3 ms per loop

A função versao_numba_lift não pode ser compilada em modo nopython, mas os dois for podem. Logo, ele extrai esta parte da função, aplica loop lifting e termina com um código tão rápido quanto o anterior. Porém, isto não funcionaria bem se fizéssemos o mesmo na função morpho_gradient_numba_lift. Apesar do loop da outra função ser realmente acelerado, uma funcionalidade importante do Numba é que a chamada de funções em modo nopython é muito rápida para outras funções nopython, mas lenta para funções Python (ref). Como em Cython, é possível ganhar muita velocidade se todas as nossas funções forem compiladas no modo acelerado.


Como pudemos ver, Numba pode ser uma excelente alternativa para acelerar código escrito em Python. Em um próximo post apresentarei as funções (e limitações) de Numba para rodar código em GPUs NVidia usando CUDA.