Moyens efficaces pour itérer sur plusieurs tableaux numpy et traiter les éléments courants et précédents?

Question

J'ai lu beaucoup de choses sur différentes techniques pour itérer sur les tableaux numpy récemment et il semble que le consensus soit de ne pas itérer du tout (par exemple, voir un commentaire ici). Il y a plusieurs questions similaires sur SO, mais mon cas est un peu différent car je dois combiner "l'itération" (ou ne pas itérer) et l'accès aux valeurs précédentes.

Disons qu'il y a N (N est petit, généralement 4, peut être jusqu'à 7) tableaux numpy 1-D de type float128 dans une liste X, tous les tableaux sont de la même taille. Pour vous donner un petit aperçu, il s'agit de données d'intégration EDP, chaque tableau représente une fonction, et j'aimerais appliquer une section de Poincaré. Malheureusement, l'algorithme devrait être à la fois efficace en mémoire et en temps, car ces tableaux font parfois ~1 Go chacun, et il n'y a que 4 Go de RAM à bord (je viens de découvrir l'utilisation de memmap pour les tableaux numpy et envisage maintenant de les utiliser à la place des tableaux classiques).

Un de ces tableaux est utilisé pour "filtrer" les autres, donc je commence par secaxis = X.pop(idx). Maintenant je dois localiser les paires d'indices où (secaxis[i-1] > 0 et secaxis[i] < 0) ou (secaxis[i-1] < 0 et secaxis[i] > 0) et ensuite appliquer des transformations algébriques simples aux autres tableaux, X (et sauvegarder les résultats). Il est à noter que les données ne doivent pas être gaspillées lors de cette opération.

Il y a plusieurs façons de le faire, mais aucune ne me semble efficace (et assez élégante). L'une est une approche de type C, où vous itérez simplement dans une boucle for :

import array # mieux que les listes
res = [ array.array('d') for _ in X ]
for i in xrange(1,secaxis.size):
  if condition: # voir ci-dessus
    co = -secaxis[i-1]/secaxis[i]
    for j in xrange(N):
      res[j].append( (X[j][i-1] + co*X[j][i])/(1+co) )

C'est clairement très inefficace et en plus ce n'est pas une approche pythonique.

Une autre façon est d'utiliser numpy.nditer, mais je n'ai pas encore compris comment on accède à la valeur précédente, même s'il permet d'itérer sur plusieurs tableaux en même temps :

# sans secaxis = X.pop(idx)
it = numpy.nditer(X)
for vec in it:
  # vec[idx] est la valeur actuelle, comment obtenir la précédente (ou la suivante) ?

La troisième possibilité est de d'abord trouver les indices recherchés avec des tranches numpy efficaces, puis de les utiliser pour des multiplications/additions en vrac. Je préfère celle-ci pour l'instant :

res = []
inds, = numpy.where((secaxis[:-1] < 0) * (secaxis[1:] > 0) +
                   (secaxis[:-1] > 0) * (secaxis[1:] < 0))
coefs = -secaxis[inds] / secaxis[inds+1] # tableau de coefficients
for f in X: # la boucle est faite seulement N-1 fois, c'est-à-dire, de 3 à 6
    res.append( (f[inds] + coefs*f[inds+1]) / (1+coefs) )

Mais cela semble être fait en 7 + 2*(N - 1) passes, de plus, je ne suis pas sûr de l'adressage secaxis[inds] (ce n'est pas du slicing et généralement il doit trouver tous les éléments par indices comme dans la première méthode, n'est-ce pas ?).

Enfin, j'ai également essayé d'utiliser itertools et cela a résulté en des structures monstrueuses et obscures, ce qui pourrait venir du fait que je ne suis pas très familier avec la programmation fonctionnelle :

def filt(x):
  return (x[0] < 0 and x[1] > 0) or (x[0] > 0 and x[1] < 0)
import array
from itertools import izip, tee, ifilter
res = [ array.array('d') for _ in X ] 
iters = [iter(x) for x in X]   # N-1 itérateurs dans une liste
prev, curr = tee(izip(*iters)) # 2 itérateurs similaires, dont chacun
                               # est composé de N-1 itérateurs
next(curr, None) # l'un d'eux est maintenant pour la valeur actuelle
seciter = tee(iter(secaxis))
next(seciter[1], None)
for x in ifilter(filt, izip(seciter[0], seciter[1], prev, curr)):
  co = - x[0]/x[1]
  for r, p, c in zip(res, x[2], x[3]):
    r.append( (p+co*c) / (1+co) )

Ce n'est pas seulement très laid, cela prend aussi un temps terrible pour s'exécuter.

Donc, j'ai les questions suivantes :

1. De toutes ces méthodes, est-ce que la troisième est effectivement la meilleure ? Si c'est le cas, que peut-on faire pour améliorer la dernière ?
2. Y a-t-il d'autres méthodes meilleures encore ?
3. Par simple curiosité, y a-t-il un moyen de résoudre le problème en utilisant nditer ?
4. Finalement, serais-je mieux loti en utilisant des versions de memmap des tableaux numpy, ou cela ralentirait-il probablement beaucoup les choses ? Peut-être devrais-je seulement charger le tableau secaxis en RAM, laisser les autres sur le disque et utiliser la troisième méthode ?
5. (question bonus) Une liste de N tableaux numpy 1-D de même longueur provient du chargement de N fichiers .npy dont les tailles ne sont pas connues à l'avance (mais N l'est). Serait-il plus efficace de lire un tableau, puis allouer de la mémoire pour un tableau numpy 2-D (légère surcharge mémoire ici) et lire les autres dans ce tableau 2-D ?

Answer 1

La version numpy.where() est assez rapide, vous pouvez l'accélérer un peu en utilisant method3(). Si la condition > peut être changée en >=, vous pouvez également utiliser method4().

import numpy as np

a = np.random.randn(100000)

def method1(a):
    idx = []
    for i in range(1, len(a)):
        if (a[i-1] > 0 and a[i] < 0) or (a[i-1] < 0 and a[i] > 0):
            idx.append(i)
    return idx

def method2(a):
    inds, = np.where((a[:-1] < 0) * (a[1:] > 0) +
                       (a[:-1] > 0) * (a[1:] < 0))
    return inds + 1

def method3(a):
    m = a < 0
    p = a > 0
    return np.where((m[:-1] & p[1:]) | (p[:-1] & m[1:]))[0] + 1

def method4(a):
    return np.where(np.diff(a >= 0))[0] + 1

assert np.allclose(method1(a), method2(a))
assert np.allclose(method2(a), method3(a))
assert np.allclose(method3(a), method4(a))

%timeit method1(a)
%timeit method2(a)
%timeit method3(a)
%timeit method4(a)

le résultat de %timeit:

1 boucle, la meilleure de 3: 294 ms par boucle
1000 boucles, la meilleure de 3: 1,52 ms par boucle
1000 boucles, la meilleure de 3: 1,38 ms par boucle
1000 boucles, la meilleure de 3: 1,39 ms par boucle

Answer 2

Je devrai lire votre message plus en détail, mais je commencerai par quelques observations générales (à partir des questions des itérations précédentes).

Il n'y a pas de façon efficace d'itérer sur des tableaux en Python, même s'il y a des choses qui ralentissent les choses. J'aime distinguer entre le mécanisme d'itération (nditer, pour x dans A :) et l'action (alist.append(...), x[i+1] += 1). Le grand consommateur de temps est généralement l'action, effectuée de nombreuses fois, et non pas le mécanisme d'itération lui-même.

Laisser numpy itérer dans du code compilé est le plus rapide.

 xdiff = x[1:] - x[:-1]

est bien plus rapide que

 xdiff = np.zeros(x.shape[0]-1)
 for i in range(x.shape[0]:
     xdiff[i] = x[i+1] - x[i]

Le np.nditer n'est pas plus rapide.

nditer est recommandé comme un outil d'itération générale dans du code compilé. Mais sa principale valeur réside dans la gestion du broadcasting et de la coordination de l'itération sur plusieurs tableaux (entrée/sortie). Et vous devez utiliser la mise en mémoire tampon et le code c pour obtenir la meilleure vitesse de nditer (je vais chercher une récente question sur SO).

https://stackoverflow.com/a/39058906/901925

N'utilisez pas nditer sans étudier la page de tutoriel sur l'itération pertinente (celle qui se termine par un exemple cython).

\=========================

Juste en se basant sur l'expérience, cette approche sera la plus rapide. Oui, elle va itérer sur secaxis un certain nombre de fois, mais toutes ces actions sont réalisées dans du code compilé et seront bien plus rapides que n'importe quelle itération en Python. De plus, l'itération for f dans X : n'a lieu que quelques fois.

res = []
inds, = numpy.where((secaxis[:-1] < 0) * (secaxis[1:] > 0) +
                   (secaxis[:-1] > 0) * (secaxis[1:] < 0))
coefs = -secaxis[inds] / secaxis[inds+1] # tableau des coefficients
for f in X: 
    res.append( (f[inds] + coefs*f[inds+1]) / (1+coefs) )

@HYRY a exploré des alternatives pour rendre l'étape de où plus rapide. Mais comme vous pouvez le voir, les différences ne sont pas si grandes. D'autres ajustements possibles

inds1 = inds+1
coefs = -secaxis[inds] / secaxis[inds1]
coefs1 = coefs+1
for f in X:
    res.append(( f[inds] + coefs*f[inds1]) / coefs1)

Si X était un tableau, res pourrait être un tableau aussi.

res = (X[:,inds] + coefs*X[:,inds1])/coefs1

Mais pour un petit N, je suppose que la liste res est tout aussi bonne. Pas besoin de rendre les tableaux plus grands que nécessaire. Les ajustements sont mineurs, juste essayer d'éviter de recalculer les choses.

\=================

Cet utilisation de np.where est juste np.nonzero. Cela fait en réalité deux passages du tableau, une fois avec np.count_nonzero pour déterminer combien de valeurs il retournera et créer la structure de retour (liste de tableaux de longueur désormais connue). Et une deuxième boucle pour remplir ces indices. Donc de multiples itérations sont bonnes si elles rendent l'action simple.