Boucles versus vectorisation

paul18

Pour faire écho à un autre post, par le passé j'avais mis en place des cas tests pour comparer les temps de calculs utilisant des boucles et la vectorisation. Personnellement je note jusqu'à des facteurs 100 en pratique sur mes applications

Sur le sujet, je ne saurais trop vous conseiller le lire les documents ici et là
Paul

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]

import numpy as np
import time
 
def square_function(x):
    return x**2
 
 
n = 10_000_000; # 10 million itérations
 
 
# case 1: using loops
y1 = np.zeros(n, dtype = np.float) # initialization
 
t0 = time.time()
for i in range(n):
    y1[ i] = 0.1*i**2
 
t1 = time.time()
duration1 = t1 - t0
print("Duration using loops = {} seconds".format(duration1))
 
#case 2: using vectorization
y2 = np.zeros(n, dtype = np.float)      # initialization
 
t0 = time.time()
i = np.arange(n, dtype = np.float)    # vector of indexes is created
y2 = 0.1*i**2
t1 = time.time()
duration2 = t1 - t0
print("Duration using Vectorization = {} seconds".format(duration2))
 
## ratio calculation
print("Vectorization is {} times faster than loops use!!!".format(duration1 / duration2))
print("we check there's not difference: min(y2-y1) = {} and max(y2-y1) = {}".format(np.min(y2-y1),np.max(y2-y1)))
del y2; del y1
 
# case 3: using a function works as well
y3 = np.zeros(n, dtype = np.float)      # initialization
 
t0 = time.time()
i = np.arange(n, dtype = np.float)     # vector of indexes is created
y3 = square_function(i)
t1 = time.time()
duration3 = t1 - t0
print("Duration using Vectorization + function = {} seconds".format(duration3))
del y3

paul18

Bonjour

Dans la poursuite de ce post, je propose d'initier un fil traitant de la vectorisation en Python (la démarche reste similaire sous Scilab, et sûrement sous Matlab aussi); l'objectif se veut d'échanger sur des méthodologies et façons de gagner en efficacité et rapidité (loin de moi l'idée d'avoir les meilleurs, bien au contraire).

Quelques remarques préalables qui me viennent à l'esprit:

• bien sûr l'objectif est de n'utiliser aucune boucle
• Pour chaque exemple, le temps CPU est mesuré; bien que dépendant de la machine et du système d'exploitation, cela reste le seul indicateur de performance à mon sens
• la librarie Numba ne sera pas utilisée, ni la parallélisation ni le multithreading
• un code source est proposé, et chacun est libre de le tester et de l'améliorer; en retour toute amélioration est partagée
• la taille des matrices est une variable dès le commencement (de façon à pousser plus loin le code)
• j'espère qu'en retour d'autres proposeront leurs cas tests
• je suis sur Python 3.6 avec Spyder comme IDE

Aller je me lance:

Cas 1:
Extraire dans une nouvelle matrice "Extraction" toutes les lignes où "47" apparaît dans la colonne 3

import numpy as np
import time

t0 = time.time()
n = 10
A = np.random.randint(66, size = (600*n,n), dtype = np.int32)
indice = np.where(A[:,2] == 47)
Extraction = np.copy(A[indice[0],:])
del indice
t1 = time.time()
print("Durée: {}".format(t1-t0))

Quelques temps CPU sur mon (vieux) portable avec 6 Go de mémoire:
n = 10 $\Rightarrow$ moins de 3 ms
n = 100 $\Rightarrow$ moins de 0.2 s
n = 1000 $\Rightarrow$ moins de 20 s

paul18

Ah j'ai remis la main sur le "benchmark" que je cherchais : voir site de la Nasa

paul18

Cas 2:
Renuméroter les éléments de la première colonne de $1$ à $10 \; 000n$

n = 10
t0 = time.time()
A = np.random.random( (10_000*n,int(0.1*n)) )
i = np.arange(1,10_000*n+1)
A[:,0] = i
t1 = time.time()
print("Durée: {}".format(t1-t0))

Quelques temps CPU:
$n=10$ $\Rightarrow$ moins de 3 ms
$n=100$ $\Rightarrow$ moins de 0.2 s
$n=500$ $\Rightarrow$ moins de 4.6 s

paul18

Bonjour

Jamais deux sans trois comme il se dit

Un exemple applicatif sur le calcul du volume d'une sphère à partir de son maillage de peau. Je pensais utiliser matplotlib et la triangulation de Delaunay, mais je rencontre des problèmes dans le cadre des surfaces fermées (maillage devient volumique et non surfacique); j'ai donc fait le choix d'utiliser la librairie trimesh et de générer des fichiers STL avec GMSH

Etape 1: génération du maillage sous gmsh

R_ = 2.;
SetFactory("OpenCASCADE");
Mesh.CharacteristicLengthFactor = 0.05;
Sphere(1) = {0., 0., 0, R_, -Pi/2, Pi/2, 2*Pi};

En jouant sur la valeur de Mesh.CharacteristicLengthFactor appelée $n$ ici, le maillage est affiné (ce qui n'empêche pas les erreurs dues à la facettisation). En pratique:

• n=1 $\Rightarrow$ 330 triangles
• n=0.1 $\Rightarrow$ 25 293 triangles
• n=0.05 $\Rightarrow$ 100 314 triangles
• n=0.02 $\Rightarrow$ 618 972 triangles

En ligne de commande, le maillage est lancé avec:

gmsh -2 sphere.geo -o sphere.stl -format stl

Code Python:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# cet exemple demande d'installer trimesh (pyglet est optionnel)
# https://pypi.org/project/trimesh/

# pour GMSH:
# http://gmsh.info


import time, os
import numpy as np
import trimesh

PATH = str(os.getcwd())
t0 = time.time()
maillage = trimesh.load(PATH + '/sphere_01.stl') # cas sphère de rayon R = 2.
noeuds = maillage.vertices  # récupère les coordonnées de tous les noeuds du maillages
triangles = maillage.faces  # récupère les sommets de chaque triangles
print("Le maillage comporte {} éléments".format(len(triangles)))
volume_ = maillage.volume   # calcul le volume du maillage
t1 = time.time()
print("Durée de traitement par trimesh = {}\n".format(t1 - t0))

t0 = time.time()
# coordonnées des noeuds par triangle
x0 = noeuds[triangles[:,0],0]
x1 = noeuds[triangles[:,1],0]
x2 = noeuds[triangles[:,2],0]

y0 = noeuds[triangles[:,0],1]
y1 = noeuds[triangles[:,1],1]
y2 = noeuds[triangles[:,2],1]

z0 = noeuds[triangles[:,0],2]
z1 = noeuds[triangles[:,1],2]
z2 = noeuds[triangles[:,2],2]

# calcul du volume en utilisant la vectorisation avec la formule:
V_ = (x0 + x1 + x2) * ( (y1 - y0)*(z2 - z0) - (y2 - y0)*(z1 - z0))
V_ = np.sum(V_)/6.
t1 = time.time()

print("Volume avec trimesh = {}".format(volume_))
print("Volume avec vectorisation = {}".format(V_))
print("Volume théorique d'une sphère= {}".format(4*np.pi*2**3/3))
print("Durée de la vectorisation = {}".format(t1 - t0))
del x0; del x1; del x2; del y0; del y1; del y2; del z0; del z1; del z2

Les temps de calculs par vectorisation sont:

• n=1 $\Rightarrow$ moins de 1 ms
• n=0.1 $\Rightarrow$ moins de 3 ms
• n=0.05 $\Rightarrow$ moins de 0.02 s
• n=0.02 $\Rightarrow$ moins de 0.2 s

J'aurais pu continuer avec un maillage plus fin ($n$=0.01 donne 2 461 267 triangles), mais le temps de lecture par trimesh du fichier STL (plus de 460 Mo) devient rédhibitoire.

Paul

marsup

Personnellement je note jusqu'à des facteurs 100 en pratique sur mes applications

Et donc finalement, sur les exemples que tu donnes, tu vérifies encore ce gros avantage pour la vectorisation ?

paul18

Plus que jamais !

je ne me pose même plus la question car j'ai fait ce constat à maintes reprises; dans un autre post, j'ai parlé d'une démonstration que j'ai faîte à l'un de mes anciens étudiants qui avait utilisé 3 boucles imbriquées (40 minutes de calculs), contre moins d'une minute en utilisant la vectorisation.

Même constat sous Scilab.

Si comme pour Saint Thomas il faut en faire la démonstration, allons-y sur un des exemples ci-dessus

marsup

Ok, je suis d'accord. J'ai aussi remarqué la même chose sous Scilab, même si je fais des calculs beaucoup plus modestes.

paul18

sur ce petit exemple j'ai un ratio de 20.

Proposez des bouts de code à vectoriser et on verra si c'est possible et quels sont les gains, non ?

import numpy as np
import time

n = 100_000_000
t0 = time.time()
A = np.random.random( (n,1))
i = np.arange(1,n+1)
A[:,0] = i
t1 = time.time()
T1 = t1 - t0
print("Durée AVEC vectorisation: {}".format(T1))

t0 = time.time()
for i in range(n):
    A[i,0] = i+1
t1 = time.time()
T2 = t1 - t0
print("Durée SANS vectorisation: {}".format(T2))
print("ratio = {}".format(T2/T1))
# del A; del n; del i

nicolas.patrois

#!/usr/bin/python3

from time import time
from itertools import permutations

n=10
t0=time()
compte1=0
for _ in permutations(range(n)):
  pass#compte1+=1
t1=time()
T1=t1-t0
print("Durée AVEC vectorisation: {}".format(T1))

t0=time()
def liste(l):
 if len(l)==1:
  yield l
 for i in range(len(l)):
  for ll in liste(l[:i]+l[i+1:]):
   yield [l[i]]+ll

compte2=0
for _ in liste(list(range(1,n+1))):
 pass#compte2+=1

t1=time()
T2=t1-t0
print("Durée SANS vectorisation: {}".format(T2))
print("ratio = {}".format(T2/T1))
print(compte1,compte2)

Avec les variables de décompte :

Durée AVEC vectorisation: 0.8243844509124756
Durée SANS vectorisation: 18.384573221206665
ratio = 22.300970379605747

Et si les boucles tournent à vide :

Durée AVEC vectorisation: 0.4252817630767822
Durée SANS vectorisation: 17.91806960105896
ratio = 42.13223127986316

paul18

En mélangeant renumérotation et conditions, j'obtiens un ratio de plus de 160 (:P)

import numpy as np
import time

# objectifs :
# si A_i <= 0.25, la valeur est renumérotée 1
# si 0.25 < A_i <= 0.5, la valeur est renumérotée 2
# si 0.5 < A_i <= 0.75, la valeur est renumérotée 3
# si 0.75 < A_i <= 1., la valeur est renumérotée 4

n = 1_000_000
t0 = time.time()
A = np.random.random( (n,1))
B = np.copy(A)
a1 = np.where( A <= 0.25 ) # indices pour lesquels la condition est vraie
a2 = np.where( (A > 0.25) & (A <= 0.5) )
a3 = np.where( (A > 0.5) & (A <= 0.75) )
a4 = np.where( (A > 0.75) & (A <= 1.) )
A[a1] = 1.
A[a2] = 2.
A[a3] = 3.
A[a4] = 4.
t1 = time.time()
T1 = t1 - t0
print("Durée AVEC vectorisation: {}".format(T1))

t0 = time.time()
for i in range(n):
    if (B[ i] <= 0.25):
        B[ i] = 1.
    elif (B[ i] > 0.25) & (B[ i] <= 0.5):
        B[ i] = 2.
    elif (B[ i] > 0.5) & (B[ i] <= 0.75):
        B[ i] = 3.
    elif (B[ i] > 0.75) & (B[ i] <= 1.):
        B[ i] = 4.           
t1 = time.time()
T2 = t1 - t0
print("Durée SANS vectorisation: {}".format(T2))
print("ratio = {}".format(T2/T1))
# del A; del n; del i

Erratum: je n'avais pas compris le message de AD portant sur les balises avec $\left[ i \right]$, du coup le code n'apparaissait pas correctement

paul18

Un nouveau cas test en cette nouvelle d'année 2020; ici on utilise np.mod, np.unique, np.any par exemple. Je n'ai pas encore trouvé comment me passer de la boucle pour l'affichage (contrairement à Scilab par exemple qui accepte la vectorisation pour print)
Il y a sûrement d'autres façon de procéder.
Paul

import numpy as np
import time

# objectif: rechercher tous les multiples de COEF et leurs occurences
# etape 1: calcul du vecteur "Modulo"
# etape 2: extraction des multiples avec "Ou"
# etape 3: si des multiples existent, alors on recherche leur(s) occurence(s) avec "Compteur"

t0 = time.time()
n = 101 # 6 doit etre le minimum
COEF = 5

A = np.random.randint( 2, high = int(0.5*n), size = (n), dtype = np.int32)
# A = np.array([4, 12, 16]) # à tester avec COEF = 4

Modulo = np.mod(A, COEF)

Ou = np.where(Modulo == 0)

    
if not np.any(Ou): # test si Ou est vide
    print("Il n'y a pas de multiples de {}".format(COEF))
    
else:
    Multiples = np.copy(A[Ou])
    MultiplesUniques, Compteur = np.unique(Multiples, return_counts = True)
    
    if (len(Compteur) == 1):
        print("{} est l'unique multiple et il apparait {} fois".format(MultiplesUniques[0], np.shape(Multiples)[0]))
    elif ( (len(Compteur) > 1) & (np.max(Compteur) == 1) ):
        print("{} sont les multiples de {} et ils n'apparaissent qu'une seule fois chacun".format(Multiples, COEF))
    else:
        # print("{} sont les multiples de {} et ils apparaissent {} fois respectivement".format(MultiplesUniques, COEF, Compteur))
        for i in range(len(MultiplesUniques)):
            print("\t{} (multiple de {}) apparait {} fois".format(MultiplesUniques[ i], COEF, Compteur[ i]))

t1 = time.time()
print("\nDurée : ", t1-t0)