Boucles versus vectorisation
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]
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
Réponses
-
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 3import 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 -
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 -
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 gmshR_ = 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 -
Personnellement je note jusqu'à des facteurs 100 en pratique sur mes applications
-
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 -
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.
-
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
-
#!/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
Algebraic symbols are used when you do not know what you are talking about.
-- Schnoebelen, Philippe -
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 -
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.
Paulimport 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)
Connectez-vous ou Inscrivez-vous pour répondre.
Bonjour!
Catégories
- 163.1K Toutes les catégories
- 8 Collège/Lycée
- 21.9K Algèbre
- 37.1K Analyse
- 6.2K Arithmétique
- 53 Catégories et structures
- 1K Combinatoire et Graphes
- 11 Sciences des données
- 5K Concours et Examens
- 11 CultureMath
- 47 Enseignement à distance
- 2.9K Fondements et Logique
- 10.3K Géométrie
- 62 Géométrie différentielle
- 1.1K Histoire des Mathématiques
- 68 Informatique théorique
- 3.8K LaTeX
- 39K Les-mathématiques
- 3.5K Livres, articles, revues, (...)
- 2.7K Logiciels pour les mathématiques
- 24 Mathématiques et finance
- 312 Mathématiques et Physique
- 4.9K Mathématiques et Société
- 3.3K Pédagogie, enseignement, orientation
- 10K Probabilités, théorie de la mesure
- 772 Shtam
- 4.2K Statistiques
- 3.7K Topologie
- 1.4K Vie du Forum et de ses membres