creando piccoli array in Cython prende una quantità colossale di tempo

Question

stavo scrivendo un nuovo generatore di numeri casuali per NumPy che produce numeri casuali in base ad una distribuzione arbitraria, quando mi sono imbattuto in questo comportamento davvero strano:

questo è prova. pisside

setup.py

from distutils.core import setup 
from Cython.Build import cythonize 
setup(name = "simple cython func",ext_modules = cythonize('test.pyx'),) 

il profiling del codice

#!/usr/bin/python 
from __future__ import division 

import subprocess 
import timeit 

#Compile the cython modules before importing them 
subprocess.call(['python', 'setup.py', 'build_ext', '--inplace']) 

sstr=""" 
import test 
import numpy 
u=numpy.random.random(10) 
a=numpy.random.random(10) 
a=numpy.cumsum(a) 
a/=a[-1] 
r=numpy.empty(10,int) 
""" 

print "binary search: creates an array[N] and performs N binary searches to fill it:\n",timeit.timeit('numpy.searchsorted(a,u)',sstr) 
print "Simple replacement for binary search:takes the same args as np.searchsorted and similarly returns a new array. this performs only one trivial operation per element:\n",timeit.timeit('test.BSReplacement(a,u)',sstr) 

print "barebones function doing nothing:",timeit.timeit('test.BareBones(a,u,r)',sstr) 
print "Untyped inputs and doing N iterations:",timeit.timeit('test.UntypedWithLoop(a,u,r)',sstr) 
print "time for just np.empty()",timeit.timeit('numpy.empty(10,int)',sstr) 

L'implementazione della ricerca binaria è nell'ordine di len(u)*Log(len(a)) tempo di esecuzione. La banale funzione cython richiede l'esecuzione di len(u). Entrambi restituiscono un array int 1D di len (u).

tuttavia, anche questa non banale implementazione computazionale richiede più tempo della ricerca binaria completa nella libreria numpy. (È stato scritto in C: https://github.com/numpy/numpy/blob/202e78d607515e0390cffb1898e11807f117b36a/numpy/core/src/multiarray/item_selection.c vedi PyArray_SearchSorted)

I risultati sono:

binary search: creates an array[N] and performs N binary searches to fill it: 
1.15157485008 
Simple replacement for binary search:takes the same args as np.searchsorted and similarly returns a new array. this performs only one trivial operation per element: 
3.69442796707 
barebones function doing nothing: 0.87496304512 
Untyped inputs and doing N iterations: 0.244267940521 
time for just np.empty() 1.0983929634 

Perché la np.empty() passo prendendo così tanto tempo? e cosa posso fare per ottenere un array vuoto che posso restituire?

La funzione C esegue questo AND esegue un intero gruppo di controlli di integrità e utilizza un algoritmo più lungo nel ciclo interno. (Ho rimosso tutta la logica tranne la stessa ciclo fro mio esempio)

---

Aggiorna

Sembra che vi sono due problemi distinti:

1. Il np.empty (10) la chiamata da sola ha un sovraccarico enorme e richiede tutto il tempo necessario per searchsorted per creare un nuovo array ed eseguire 10 ricerche binarie su di esso
2. Solo dichiarando la sintassi del bufferIlha anche un sovraccarico che richiede più tempo rispetto alla ricezione delle variabili non tipizzate E iterando 50 volte.

risultati per 50 iterazioni:

binary search: 2.45336699486 
Simple replacement:3.71126317978 
barebones function doing nothing: 0.924916028976 
Untyped inputs and doing N iterations: 0.316384077072 
time for just np.empty() 1.04949498177 

Answer 1

Creazione np.empty all'interno della funzione Cython ha un certo overhead come già visto. Qui si vedrà un esempio su come creare la matrice vuota e passarlo al modulo Cython per riempire con i valori corretti:

n=10:

numpy.searchsorted: 1.30574745517 
cython O(1): 3.28732016088 
cython no array declaration 1.54710909596 

n=100:

numpy.searchsorted: 4.15200545373 
cython O(1): 13.7273431067 
cython no array declaration 11.4186086744 

Come già sottolineato, la versione numpy si adatta meglio poiché è O(len(u)*long(len(a))) e questo algoritmo qui è O(len(u)*len(a)) ...

Ho anche provato a utilizzare Memoryview, modificando sostanzialmente np.ndarray[double, ndim=1] per double[:], ma in questo caso la prima opzione era più veloce.

Il nuovo file .pyx è:

from __future__ import division 
import numpy as np 
cimport numpy as np 
cimport cython 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def JustLoop(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] u, 
      np.ndarray[int, ndim=1] r): 
    cdef int i,j 
    for j in range(u.shape[0]): 
     if u[j] < a[0]: 
      r[j] = 0 
      continue 

     if u[j] > a[a.shape[0]-1]: 
      r[j] = a.shape[0]-1 
      continue 

     for i in range(1, a.shape[0]): 
      if u[j] >= a[i-1] and u[j] < a[i]: 
       r[j] = i 
       break 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def WithArray(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] u): 
    cdef np.ndarray[np.int_t, ndim=1] r=np.empty(u.shape[0],dtype=int) 
    cdef int i,j 
    for j in range(u.shape[0]): 
     if u[j] < a[0]: 
      r[j] = 0 
      continue 

     if u[j] > a[a.shape[0]-1]: 
      r[j] = a.shape[0]-1 
      continue 

     for i in range(1, a.shape[0]): 
      if u[j] >= a[i-1] and u[j] < a[i]: 
       r[j] = i 
       break 
    return r 

Il nuovo .py del file:

import numpy 
import subprocess 
import timeit 

#Compile the cython modules before importing them 
subprocess.call(['python', 'setup.py', 'build_ext', '--inplace']) 
from test import * 

sstr=""" 
import test 
import numpy 
u=numpy.random.random(10) 
a=numpy.random.random(10) 
a=numpy.cumsum(a) 
a/=a[-1] 
a.sort() 
r = numpy.empty(u.shape[0], dtype=int) 
""" 

print "numpy.searchsorted:",timeit.timeit('numpy.searchsorted(a,u)',sstr) 
print "cython O(1):",timeit.timeit('test.WithArray(a,u)',sstr) 
print "cython no array declaration",timeit.timeit('test.JustLoop(a,u,r)',sstr) 

Answer 2

C'è una discussione di questo sulla lista Cython che potrebbe avere alcuni suggerimenti utili: https://groups.google.com/forum/#!topic/cython-users/CwtU_jYADgM

In generale, anche se provo ad allocare piccoli array al di fuori di Cython, passali in un d riutilizzarli nelle chiamate successive al metodo. Capisco che questa non è sempre un'opzione.