NumPy

Functionalitatea Python este asigurata de module. Un modul este un fisier Python ce contine definitii si instructiuni destinate implementarii unei anumite problematici. Numele fisierului este NumeModul.py.

Exista module pe care si le creaza un programator pentru a organiza codul. De exemplu pentru a implementa calculul formei scara redusa definim modulul Formascara.py, in care includem toate definitiile si functiile ce realizeaza aceasta reducere.

Pentru a folosi acest modul intr-un alt program, SistemLin.py, importam modulul, inserand in fisierul SistemLin.py linia:

import Formascara

O functie Funct din acest modul este apelata in SistemLin.py sub forma:

Formascara.Funct(argumente)

Pentru a evita scrierea unor nume lungi de module se practica redenumirea lor in linia de import:

import Formascara as Fs

In acest caz apelul functiei va fi:

Fs.Funct(argumente)

Biblioteca standard a limbajului Python contine numeroase module. De interes pentru matematica sunt modulele math si cmath.

Modulul math contine functiile matematice uzuale, de argument real (numere in virgula mobila). De exemplu daca intr-un fisier sursa Test.py sau in IPython Notebook importam modulul math, atunci o functie din acest modul se apeleaza in forma math.NumeFunctie(argument).

In [1]:
import math

x=math.pi
y=math.sin(x/4)
z=math.sqrt(2)/2
print y, z
0.707106781187 0.707106781187

Lista functiilor din modulul math se gaseste aici: http://docs.python.org/2/library/math.html

Modulul cmath contine functii de argument complex. Functiile din acest modul le vom folosi in semestrul 2.

O colectie de module se numeste pachet.

Pachetele de baza pentru calcul stiintific in Python sunt numpy si scipy.

Numele numpy este derivat de la numerical Python, iar scipy de la scientific Python.

Pachetul numpy contine printre altele:

  • definitia array-ului n-dimensional

  • functii care actioneaza asupra array-urilor

  • modulul linalg ce contine functii specifice algebrei liniare, dar si alte module

  • instrumente pentru a integra cod C/C++

Cursul de Algebra liniara opereaza cu vectori si matrici.

  • Un vector din $\mathbb{R}^n$ este reprezentat in numpy printr-un array de dimensiune 1 (il vom referi in continuare ca array 1D).

  • O matrice $A\in\mathbb{R}^{m\times n}$ este un array 2-dimensional (2D)

Modalitati de a crea array-uri

Pentru a crea si manipula array-uri intr-un program Python sau in IPython Notebook, se importa pachetul numpy:

In [2]:
import numpy as np
In [3]:
v=np.array([-2, 1, 4], dtype=int)
w=np.array([3.14, -1, 0, 5.2], dtype=float)
print v, w
[-2  1  4] [ 3.14 -1.    0.    5.2 ]

v si w sunt vectori carora li s-au precizat coordonatele si tipul elementelor, respectiv int si float

In [4]:
A=np.array([[-1,2,4], [3,5,7],[2,0,-7]], dtype=int)
print A
[[-1  2  4]
 [ 3  5  7]
 [ 2  0 -7]]

In declararea si initializarea array-urilor in acest fel se poate omite dtype (data type) scriind doar tipul:

In [5]:
u=np.array([-3,0,2], int)
B=np.array([[1,5],[-3,4]], float)
print u
print B
[-3  0  2]
[[ 1.  5.]
 [-3.  4.]]

Spre deosebire de Python standard, in numpy exista un numar mare de tipuri de date.

Elementele unui array 1D, $u\in\mathbb{R}^n$, se acceseaza u[i], iar indicele i ia valorile 0,1,..., n-1,

Elementele unui array 2D, B, de dimensiuni (m,n) se acceseaza B[i][j] sau B[i,j], i=0,1, ..., m-1, j=0,1, ..., n-1.

In [6]:
print u[0]
print B[0][1], B[1,1]
-3
5.0 4.0

Crearea unui array vid, adica a unui array de dimensiuni prescrise de un tuple, dar neinitializat, se realizeaza astfel:

In [7]:
x=np.empty(4)
B=np.empty((2,3))

x este un vector din $\mathbb{R}^{4}$, care eventual se va initializa cu elemente de tip float, iar B este un array de 2 linii si 3 coloane (vezi tuple-ul (2,3) in definitie).

In [8]:
y=np.zeros(3)
C=np.zeros((4,4))
print y
print C
[ 0.  0.  0.]
[[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]]

np.zeros(n) si respectiv np.zeros((m,n)) creaza array-uri si initializeaza elementele lor cu zero, implicit de tip float.

In [9]:
u=np.ones(4, int)
M=np.ones((3,3), int)
print u
print M
[1 1 1 1]
[[1 1 1]
 [1 1 1]
 [1 1 1]]

np.ones creaza un array cu toate elementele egale cu 1.

In [10]:
v=np.array([1,4,6,-3], int)
z=np.zeros_like(v, float)
print v, z
[ 1  4  6 -3] [ 0.  0.  0.  0.]

np.zeros_like(A) creaza un array de aceleasi dimensiuni ca A, dar de elemente 0.

np.eye(n) genereaza matricea unitate $I_n$, avand elementele de tip implicit float:

In [11]:
I=np.eye(3)
print I
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

Daca dorim sa generam o matrice unitate cu elemente de tip int atunci se apeleeaza np.eye astfel:

In [12]:
J=np.eye(2, dtype=int)
print J
[[1 0]
 [0 1]]

Conversia listelor si tuple-urilor in array-uri

In [13]:
t=(3,4,7)# acesta este un tuple
ta=np.asarray(t)

print t
print ta
(3, 4, 7)
[3 4 7]
In [14]:
L=[5,6,8,11]

print type(L), L
La=np.asarray(L, float)
print type(La), La
<type 'list'> [5, 6, 8, 11]
<type 'numpy.ndarray'> [  5.   6.   8.  11.]

np.asarray(Ob) converteste obiectul Ob (in cazul de mai sus un tuple si respectiv o lista de date numerice) la un array.

In [15]:
A=np.array([[1,2,3],[4,5,6]], float)
B=np.copy(A)
print A
print B

A[1][2]=-7
print A
print B
[[ 1.  2.  3.]
 [ 4.  5.  6.]]
[[ 1.  2.  3.]
 [ 4.  5.  6.]]
[[ 1.  2.  3.]
 [ 4.  5. -7.]]
[[ 1.  2.  3.]
 [ 4.  5.  6.]]

np.copy(A) creaza un nou array cu aceleasi elemente ca si A. Observam ca modificand un element in A si afisand din nou pe A si B, B nu s-a schimbat, deci este un array pentru care s-a alocat o zona separata de memorie.

De interes deosebit sunt functiile np.arange si np.linspace, care creaza array-uri speciale:

In [16]:
n=6
x=np.arange(6)
print x
xx=np.arange(2,9)
print xx
[0 1 2 3 4 5]
[2 3 4 5 6 7 8]
In [17]:
y=np.arange(3,12,2)
print y
[ 3  5  7  9 11]
In [18]:
z=np.arange(0,1, 0.3)
print z
[ 0.   0.3  0.6  0.9]
  • Daca n este un numar intreg pozitiv, atunci x=np.arange(n) este vectorul de elemente 0,1, ..., n-1.

  • xx=np,arange(m,n), cu m<n, este vectorul de elemente m, m+1, ..., n-1

  • v=np.arrange(a, b, h), cu a<b si h>0, este vectorul de elemente v[i]= a+i*h, unde i ia valori de la 0 pana la valoarea maxima, imax, pentru care a+imax*h<b.

  • a este valoarea de start, b valoarea de stop, iar h este pasul cu care se divide intervalul [a,b] prin puncte echidistante.

    • In definitia lui y avem a=3, b=12, h=2

    • In definitia lui z, a=0, b=1, h=0.3

Pasul h poate fi si negativ. In acest caz, x=np.arange(b,a,h), cu b>a, genereaza vectorul de elemente x[i]=b+i*h, unde i ia valori de la 0 pana la valoarea maxima imax pt care b+imax*h>a.

In [19]:
x=np.arange(10,4,-1)
print x
[10  9  8  7  6  5]
In [20]:
x=np.linspace(np.pi/6, 2*np.pi/3, 10)
print x
[ 0.52359878  0.6981317   0.87266463  1.04719755  1.22173048  1.3962634
  1.57079633  1.74532925  1.91986218  2.0943951 ]

x=np.linspace(a, b, n) este vectorul ale carui elemente se obtin divizand intervalul [a,b] prin n puncte echidistante. In exemplul precedent, intervalul $[\pi/6, 2\pi/3]$ a fost divizat prin 10 puncte echidistante.

np.arange(a, b, h) si np.linspace(a,b,n) se vor folosi pentru a genera grafice de functii, pentru a trasa curbe.

Crearea unui array dintr-un fisier

Exemplele de mai sus au fost ilustrate prin array-uri cu numar redus de elemente. In problemele practice se lucreaza cu array-uri 1D sau 2D cu numar mare de elemente. In acest caz acestea se citesc din fisiere.

Functia numpy.loadtxt citeste un fisier text ce contine date din care se pot construi array-uri.

In [21]:
v=np.loadtxt("Myfiles/vector.txt")
print v
print v.size
[ 2.  3.  5.  7.  1. -2.  3.  0.  1.  6.  3.  2.  6.  3.  5.  2.  4.]
17

In fisierul vector.txt fiecare numar este scris pe o alta linie. Evident ca ele se pot scrie pe aceeasi linie separate de spatiu.

Citind un astfel de fisier se genereaza un vector, care apoi se poate redimensiona la un array 2D, folosind metoda reshape ( a se vedea mai jos).

In [22]:
A=np.loadtxt("Myfiles/matrice.txt")
print A.shape
print  A[:,1]
(5L, 5L)
[  7.   4.   5.   2.  21.]

In fisierul matrice.txt datele sunt inregistrate matricial, adica pe 5 linii si in fiecare linie sunt 5 valori.

Functia numpy.fromfile construieste un array dintr-un fisier de date text sau binar, continand date de tip cunoscut.

Exemple de utilizare a acestei functii vom da mai tarziu.

Metode asociate unui array

O metoda asociata unui array se apeleaza astfel: NumeArray.metoda.

NumeArray.shape returneaza un tuple cu un element pentru un array 1D, respectiv cu doua elemente pentru un array 2D. Din punctul de vedere al Algebrei liniare, elementul/elementele tuple-ului indica spatiul $\mathbb{R}^n$, respectiv $\mathbb{R}^{m\times n}$, caruia apartine array-ul respectiv.

In [23]:
v=np.arange(2,21,5)
t=v.shape
print v, ',t=', t

A=np.array([[-2,3,5], [1,4,0]], int)
print 'A are dimensiunile', A.shape
[ 2  7 12 17] ,t= (4L,)
A are dimensiunile (2L, 3L)

NumeArray.size returneaza numarul de elemente ale array-ului NumeArray:

In [24]:
print v.size
print A.size
4
6

Redimensionarea unui array

In [25]:
a=np.array([2,4,6,3,5,7,8,10,12,9,11,13], int)
b=a.reshape((3,4))
print a
print b
[ 2  4  6  3  5  7  8 10 12  9 11 13]
[[ 2  4  6  3]
 [ 5  7  8 10]
 [12  9 11 13]]

Remarcam ca NumeArray.reshape(t) returneaza un array de dimensiunile setate de tuple-ul t, iar array-ul initial ramane neschimbat.

Metoda reshape este foarte utila atunci cand un array este construit din date citite dintr-un fisier.

Transpusa unui array 2D se poate genera apeland metoda NumeArray.transpose()

In [26]:
A=np.arange(10).reshape((2,5))
print A

print A.transpose()
[[0 1 2 3 4]
 [5 6 7 8 9]]
[[0 5]
 [1 6]
 [2 7]
 [3 8]
 [4 9]]

Daca v este un array 1D, deci un vector linie, apeland v.transpose() nu se genereaza un vector coloana:

In [27]:
v=np.arange(1,5)
print v, v.transpose()
[1 2 3 4] [1 2 3 4]

deoarece in Python un vector este unic afisat, pe linie.

Pentru a a afisa/manipula coordonatele lui v pe coloana se defineste v ca un array 2D, cu o singura linie, astfel:

In [28]:
v=np.array([[1,2,4]], float)
print v
print v.transpose()
print v.shape, v.transpose().shape
[[ 1.  2.  4.]]
[[ 1.]
 [ 2.]
 [ 4.]]
(1L, 3L) (3L, 1L)

Dintr-un array 2D se obtine un array 1D, concatenand linie dupa linie, apeland metoda NumeArray.flatten()

In [29]:
A=np.array([[1,0,0,1],[0,1,0,1],[1,1,1,0]], np.uint8)
v=A.flatten()
print A
print v
[[1 0 0 1]
 [0 1 0 1]
 [1 1 1 0]]
[1 0 0 1 0 1 0 1 1 1 1 0]

Slicing

Modalitatile de slicing pentru siruri se extind si pentru array-uri:

In [30]:
v=np.arange(10)
print v
print v[4:]
print v[:4]
print v[3:9:2]
[0 1 2 3 4 5 6 7 8 9]
[4 5 6 7 8 9]
[0 1 2 3]
[3 5 7]
In [31]:
A=np.array([[2,3,4],[1,5,-2]], float)
print A
print 'A[0,:]=', A[0,:]
print 'A[:,2]=', A[:,2]
print A[0,:].shape, A[:,2].shape
[[ 2.  3.  4.]
 [ 1.  5. -2.]]
A[0,:]= [ 2.  3.  4.]
A[:,2]= [ 4. -2.]
(3L,) (2L,)

A[i,:] este linia i din array-ul A, iar A[:, j] coloana j din A.

Remarcam ca atat o linie cat si o coloana sunt afisate ca array-uri 1D, adica pe linie.

Pentru a genera dintr-un array 2D blocuri de sub-array-uri se folosesc regulile de slicing atat pentru primul indice, cat si pentru al doilea:

In [32]:
A=np.array([[6,-1,2, -3],[1,3, -4, 7],[-1, 0,-2,3], [2,2,-5,1]], float)

print A
print A[:,1:]
print A[:3,:]
print A[1:3, :2]
[[ 6. -1.  2. -3.]
 [ 1.  3. -4.  7.]
 [-1.  0. -2.  3.]
 [ 2.  2. -5.  1.]]
[[-1.  2. -3.]
 [ 3. -4.  7.]
 [ 0. -2.  3.]
 [ 2. -5.  1.]]
[[ 6. -1.  2. -3.]
 [ 1.  3. -4.  7.]
 [-1.  0. -2.  3.]]
[[ 1.  3.]
 [-1.  0.]]

numpy.linalg

linalg este un modul in pachetul numpy. El se apeleaza numpy.linalg sau daca numpy a fost redefinit ca np, atunci np.linalg. O functie din acest modul se va apela ca np.linalg.NumeFunctie(argumente)

Nu toate functiile ce au corespondent in Algebra liniara sunt incluse in modulul linalg. Un astfel de exemplu este functia ce returneaza produsul scalar.

Produsul scalar a doi vectori $v, w\in\mathbb{R}^n$ se calculeaza apeland functia np.dot(v,w):

In [33]:
v=np.array([-1,2,5,3], float)
w=np.array([0,-3,1,4], float)
prod=np.dot(v,w)
print prod
11.0

Daca cei doi vectori nu au celasi numar de elemente, incercarea de calcul a produsului scalar genereaza eroare:

In [34]:
u=np.arange(3,6)
p=np.dot(v,u)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-34-10918cc0956f> in <module>()
      1 u=np.arange(3,6)
----> 2 p=np.dot(v,u)

ValueError: matrices are not aligned

Daca $A\in\mathbb{R}^{m\times p}, B\in\mathbb{R}^{p\times n}$, atunci produsul lor se calculeaza apeland, la fel np.dot(A,B):

In [35]:
A=np.arange(12).reshape((4,3))
B=np.arange(1,10).reshape((3,3))
C=np.dot(A,B)
print A
print B
print C
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]
[[1 2 3]
 [4 5 6]
 [7 8 9]]
[[ 18  21  24]
 [ 54  66  78]
 [ 90 111 132]
 [126 156 186]]

Norma unui vector:

In [36]:
v=np.array([-1,2,5], float)
no=np.linalg.norm(v) # norm este o functie din modulul linalg al pachetului numpy
print no, ', Versorul lui v=', v/no
5.47722557505 , Versorul lui v= [-0.18257419  0.36514837  0.91287093]

Determinantul si inversa unui array 2D:

In [37]:
A=np.array([[2, -5], [1, -7]], float)
print A
d=np.linalg.det(A)
print 'Determinantul lui A este d=', np.linalg.det(A)
B=np.linalg.inv(A)
print  B
[[ 2. -5.]
 [ 1. -7.]]
Determinantul lui A este d= -9.0
[[ 0.77777778 -0.55555556]
 [ 0.11111111 -0.22222222]]

Vectorul constituit din diagonala principala a unui array 2D:

In [38]:
A=np.array([[-2,3,6],[1,5,2],[0,-7,4]], float)
d=np.diag(A)
print A
print d
[[-2.  3.  6.]
 [ 1.  5.  2.]
 [ 0. -7.  4.]]
[-2.  5.  4.]

Matricea diagonala asociata unui vector:

In [39]:
v=np.array([-1,4,3], float)
D=np.diag(v)
print D
[[-1.  0.  0.]
 [ 0.  4.  0.]
 [ 0.  0.  3.]]

Alte functii din numpy si numpy.linalg vor fi prezentate pe masura ce la cursul de Algebra liniara este discutat suportul lor teoretic.

Caracteristici ale calculului cu array-uri in numpy

In Python standard daca se dau doua liste de n elemente numerice, atunci interpretand fiecare lista ca un vector din $\mathbb{R}^n$, suma lor s-ar calcula astfel:

In [40]:
def sumaL(v,w):
    n=len(v)
    s=[]
    for i in range(n):
          s.append(v[i]+w[i])
    return s
    

v=[2*i for i in range(1000)]
w=[3*j+1 for j in range(1000)]
%timeit sumaL(v,w)    
    
1000 loops, best of 3: 315 µs per loop

Pentru a estima timpul de rulare al functiei sumaL am folosit comanda magica %timeit din IPython. %timeit calculeaza media timpului de rulare a functiei. Media este calculata pentru mai multe executii ale functiei.

In numpy suma a a doua array-uri este implementata element cu element si operatia de adunare se noteaza simplu v+w:

In [41]:
n=1000
v=np.empty(n, int)
w=np.empty(n, int)
for i in range(n):
    v[i]=2*i
    w[i]=3*i+1
    
%timeit  s=v+w       
100000 loops, best of 3: 4.28 µs per loop

Comparand timpul mediu de executie a sumei listelor si respectiv array-urilor 1D in numpy realizam performanta acestei operatii in numpy.

Operatiile cu vectori si matrici in numpy sunt vectorizate.

  • Produsul Pr dintre un scalar a si un array A se calculeaza invocand Pr=a*A;

  • Daca doua array-uri A, B sunt de tip (m,n), atunci P=A*B este array-ul de tip (m,n) si elemente $p_{ij}=a_{ij}*b_{ij}$ (inmultirea se efectueaza element cu element):

In [42]:
a=2.5
v=np.arange(1,5)
print a*v

w=np.linspace(2,6,4)
print v*w
[  2.5   5.    7.5  10. ]
[  2.           6.66666667  14.          24.        ]

Functiile matematice implementate in modulul math sunt reimplementate in numpy, dar in timp ce cele din math au ca argument un scalar, cele din numpy admit ca argument array-uri.

Mai precis, daca $x=[x_0, x_1, \ldots, x_{n-1}]$, $n\geq 1$, este un array 1D, atunci $f(x)$ este array-ul $f(x)=[f(x_0), f(x_1), \ldots, f(x_{n-1})]$ si analog, daca $A$ este un array 2D de dimensiuni (m,n), atunci $B=f(A)$ este un array 2D de aceleasi dimensiuni si elementele sale sunt $b_{ij}=f(a_{ij})$.

Mai precis, sa experimentam evaluarea functiei $f(x)=\sin(x)$ apeland math.sin si np.sin, transmitand argumente scalare si respectiv vectoriale:

In [43]:
print math.sin(math.pi/3)
print np.sin(np.pi/3)
0.866025403784
0.866025403784
In [44]:
x=np.arange(0,np.pi, 0.5)
y=np.sin(x)
print x
print y.round(3)# afiseaza elementele lui y  cu 3 zecimale
[ 0.   0.5  1.   1.5  2.   2.5  3. ]
[ 0.     0.479  0.841  0.997  0.909  0.598  0.141]

Functiile matematice din modulul math nu sunt insa vectorizate:

In [45]:
X=[0, 0.5, 1,1.5, 2]
Y=math.sin(X)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-45-81498eb9146a> in <module>()
      1 X=[0, 0.5, 1,1.5, 2]
----> 2 Y=math.sin(X)

TypeError: a float is required
In [1]:
from IPython.core.display import HTML
def  css_styling():
    styles = open("./custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[1]: