Estoy haciendo un experimento para calcular la entropía aproximada de una señal. Los detalles (y el código real) se pueden encontrar en su página de Wikipedia. Desafortunadamente, aunque el algoritmo en sí funciona, es muy lento para el gran conjunto de datos (por ejemplo, tarda aproximadamente 25 segundos en una señal de 2000 largos). Dado que debo hacer este cálculo con muchas señales más largas, a esta velocidad, esperaría que mi experimento dure al menos 1 mes. Me preguntaba si hay alguna forma de acelerar el algoritmo.

import numpy as np

def ApEn(U, m, r):

    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(U)

    return abs(_phi(m + 1) - _phi(m))
2
Dillion Ecmark 29 oct. 2017 a las 07:22

3 respuestas

La mejor respuesta

Si está dispuesto a mover esa función a cython y agregar algunas anotaciones de tipo, se pueden lograr importantes mejoras de rendimiento. Aquí está mi versión de este algoritmo:

Apen.pyx:

cimport cython
from libc.math cimport fabs, log
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
@cython.cdivision(True)
cdef double max_dist(double[:] x_i, double[:] x_j, int m) nogil:
    #Performs the max function described in step 4 of ApEn algorithm
    cdef double out
    cdef double dist
    out = fabs(x_i[0] - x_j[0])
    for k in range(1, m - 1):
        dist = fabs(x_i[k] - x_j[k])
        if dist > out:
            out = dist
    return out

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.initializedcheck(False)
@cython.cdivision(True)
cdef double phi(double[:] Sn, int m, int r):
    cdef int N = len(Sn)
    cdef int i
    cdef int j
    cdef int k

    cdef int c_val
    cdef int counter
    cdef double phi_sum = 0
    cdef double phi
    cdef double m_dist

    #Performs step 3 of the ApEn algorithm
    cdef double[:, :] x = np.empty((N - m  + 1, m), dtype=np.float64)
    with nogil:
        for i in range(N - m + 1):
            for j in range(0, m):
                x[i, j] = Sn[j + i]

        #Performs a combined steps 4 & 5 of the ApEn algorithm
        for i in range(N - m + 1):
            counter = 0
            for j in range(N - m + 1):
                m_dist = max_dist(x[i], x[j], m)
                c_val = 1 if m_dist <= r else 0
                counter += c_val
            phi_sum += log(counter / (N - m + 1.0))
        phi = phi_sum / (N - m + 1.0)
    return phi

cpdef double approx_entropy(double[:] Sn, int m, int r):#Passing in steps 1 & 2 of the ApEn algorithm
    cdef double ApEn = abs(phi(Sn, m, r) - phi(Sn, m + 1, r))#Performs step 6 of the ApEn algorithm
    return ApEn

Apen.pxd:

cdef double max_dist(double[:] x_i, double[:] x_j, int m) nogil
cdef double phi(double[:] Sn, int m, int r)
cpdef double approx_entropy(double[:] Sn, int m, int r)

Setup.pxd:

from distutils.core import setup
from Cython.Build import cythonize
from distutils.core import Extension
import numpy as np

extensions = [
    Extension("apen", sources=["apen.pyx"], include_dirs=[np.get_include()], extra_compile_args=["-w"]),
]

setup(
    ext_modules = cythonize(extensions)
)

Main.py:

import time
import apen
import numpy as np

start = time.time()
data = np.random.rand(2000)
#data = np.array([85, 80, 89] * 17, dtype=np.float64)
answer = apen.approx_entropy(Sn=data, m=2, r=3)
print(answer)
end = time.time()
print(end - start)

Usando este código para 2000 puntos de datos aleatorios en mi computadora portátil, el código cython calcula ApEn en 0.36 s. Por el contrario, el código de Wikipedia tarda 14.75 s. Esto equivale a un aumento de velocidad de 40x . ¡Espero que esto sea útil!

1
CodeSurgeon 29 oct. 2017 a las 10:51

No analicé todo, pero para darte un ejemplo de cómo puedes optimizar la función usando el cálculo vectorial:

def maxdist_opti(x_i,x_j):
    return max(abs(x_i-x_j))

Cuando sus datos se almacenan en matrices numpy, puede usar operadores numpy en ellos (y hay muchos de ellos, puede echar un vistazo aquí: https://docs.scipy.org/doc/numpy-1.13.0/user/index.html) y lo hará sea mucho más rápido, en el caso anterior utilicé soustraction y la función np.max en matrices numpy.

Aquí, utilizando datos aleatorios:

x_i = np.random.rand(10000)
x_j = np.random.rand(10000)

Los datos utilizados aquí no son tan largos, pero puede ver una ganancia muy buena de rendimiento:

%timeit _maxdist(x_i,x_j)
100 loops, best of 3: 3.01 ms per loop

%timeit maxdist_opti(x_i,x_j)
10000 loops, best of 3: 28 µs per loop

Puede usar la siguiente lógica para hacer solo cálculos vectoriales en toda la fórmula, y la ganancia en rendimiento será tremenda.

Tenga en cuenta que cuanto más largos sean sus datos, más optimizado será usar el cálculo vectorial.

1
Erlinska 29 oct. 2017 a las 05:06

Normalmente, al optimizar, debe comenzar con optimizaciones algorítmicas que reducen la complejidad del algoritmo, no simplemente las constantes.

Una regla general es mirar en el bucle más interno: contiene las operaciones que se ejecutan la mayoría de las veces.

No estoy seguro de haber leído el código correctamente, pero parece que U es una matriz y _maxdist hace cálculos en sus columnas. En ese caso, tiene sentido asegurarse de que el cálculo se realice solo una vez por columna.

Por ejemplo, calcule su valor para cada columna, almacénelo en una matriz y úselo en _phi.

1
maxim1000 29 oct. 2017 a las 07:42