python - combinatorial - ¿Es `scipy.misc.comb` más rápido que un cálculo binomial ad-hoc?

combinatorial in python

Question

¿Es concluyente que ahora el scipy.misc.comb es de hecho más rápido que la implementación ad-hoc?

Según una respuesta anterior, Estadísticas: combinaciones en Python , esta función homebrew es más rápida que scipy.misc.comb cuando se calculan las combinaciones nCr :

def choose(n, k):
    """
    A fast way to calculate binomial coefficients by Andrew Dalke (contrib).
    """
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in xrange(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0

Pero después de ejecutar algunas pruebas en mi propia máquina, este no parece ser el caso, usando esta secuencia de comandos:

from scipy.misc import comb
import random, time

def choose(n, k):
    """
    A fast way to calculate binomial coefficients by Andrew Dalke (contrib).
    """
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in xrange(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok // ktok
    else:
        return 0

def timing(f):
    def wrap(*args):
        time1 = time.time()
        ret = f(*args)
        time2 = time.time()
        print ''%s function took %0.3f ms'' % (f.__name__, (time2-time1)*1000.0)
        return ret
    return wrap

@timing
def test_func(combination_func, nk):
    for n,k in nk:
        combination_func(n, k)

nk = []
for _ in range(1000):
    n = int(random.random() * 10000)
    k = random.randint(0,n)
    nk.append((n,k))

test_func(comb, nk)
test_func(choose, nk)

Obtengo el siguiente resultado:

$ python test.py
/usr/lib/python2.7/dist-packages/scipy/misc/common.py:295: RuntimeWarning: overflow encountered in exp
  vals = exp(lgam(N+1) - lgam(N-k+1) - lgam(k+1))
999
test_func function took 32.869 ms
999
test_func function took 1859.125 ms

$ python test.py
/usr/lib/python2.7/dist-packages/scipy/misc/common.py:295: RuntimeWarning: overflow encountered in exp
  vals = exp(lgam(N+1) - lgam(N-k+1) - lgam(k+1))
999
test_func function took 32.265 ms
999
test_func function took 1878.550 ms

¿La prueba de perfil de tiempo mostró que el nuevo scipy.misc.comb es más rápido que la función ad-hoc choose() ? ¿Hay algún error en mi script de prueba que hace que el tiempo sea inexacto?

¿Por qué es que scipy.misc.comb es más rápido ahora? ¿Es por algunos trucos de envoltura cython / c ?

EDITADO

Después del comentario de @WarrenWeckesser:

Al usar la aproximación de coma flotante predeterminada al usar scipy.misc.comb() , el cómputo se rompe debido al desbordamiento de coma flotante.

(Consulte http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.misc.comb.html para obtener documentación)

Cuando se prueba con exact=True que computa con enteros largos en lugar de coma flotante usando la función siguiente, es mucho más lento cuando se calculan 1000 combinaciones:

@timing
def test_func(combination_func, nk):
    for i, (n,k) in enumerate(nk):
        combination_func(n, k, exact=True)

[fuera]:

$ python test.py
test_func function took 3312.211 ms
test_func function took 1764.523 ms

$ python test.py
test_func function took 3320.198 ms
test_func function took 1782.280 ms

Answer 1

En referencia al código fuente de scipy.misc.comb , la rutina de actualización del resultado es:

    val = 1
    for j in xrange(min(k, N-k)):
        val = (val*(N-j))//(j+1)
    return val

mientras que la rutina de actualización que sugirió es:

    ntok = 1
    ktok = 1
    for t in xrange(1, min(k, n - k) + 1):
        ntok *= n
        ktok *= t
        n -= 1
    return ntok // ktok

Supongo que la razón por la que la implementación de SciPy es más lenta se debe al hecho de que la subrutina implica una división entera en cada iteración, mientras que la suya solo llama a la división una vez en la declaración de devolución.