NumPy extends lists in Python when we need to work with arrays of numbers, providing a higher performance (throught vectorization) and functionality (basically, Linear Algebra).
In [2]:
try:
import numpy as np
except:
!pip3 install numpy --user
import numpy as np
In [3]:
np.lookfor("invert")
In [7]:
np.*?
NumPy can works (create, operate and I/O) with arrays containig:
In [73]:
A = np.empty(12)
print(A) # Uninitialized garbage values, possiblely at random.
In [76]:
print(np.empty_like(A))
In [ ]:
l = [1, 2, 3]
type(l)
In [30]:
A = np.array(l)
A
In [38]:
A = np.array([[1,1.0],(1+1j,.3)])
print(A)
In [3]:
np.zeros(10)
Out[3]:
In [4]:
np.zeros((5,5))
Out[4]:
In [32]:
np.ones(10)
In [24]:
np.full((5,5), 2)
In [34]:
np.arange(10)
In [35]:
np.linspace(1., 4., 6)
In [36]:
np.random.rand(10)
In [51]:
# The random number generator is not reset in each call
print(np.random.random((5,5)))
print(np.random.random((5,5)))
In [25]:
print(np.eye(5)) # Identity matrix
In [29]:
A = np.array([i for i in range(5)]) # 1-D list comprehension
print(A, A[1], A.shape)
In [30]:
A = np.array([[j+i*4 for j in range(4)] for i in range(5)]) # 2-D list comprehension
print(A, A.shape)
In [12]:
A = np.array([1, 2], dtype=np.uint8)
A
Out[12]:
In [13]:
type(A[0])
Out[13]:
In [119]:
A = np.array([[1,2,3],[4,5,6]])
print(A)
print(A.shape)
In [120]:
np.reshape(A, (3, 2))
Out[120]:
In [108]:
A = np.empty((5,4), dtype=np.int16, order='C') # C order is the default
for row in range(A.shape[0]):
for col in range(A.shape[1]):
A[row, col] = col+row*A.shape[1]
print(A)
In [103]:
print(np.isfortran(A))
In [116]:
B = np.reshape(A, (4,5), order='F')
print(B)
In [117]:
print(np.isfortran(B))
In [118]:
print(A.flags)
print(B.flags)
In [9]:
A = np.array([], dtype=np.uint8)
print(A, A.dtype)
A
Out[9]:
In [11]:
A = np.concatenate((A, np.array([1, 2], dtype=np.uint8)))
A
Out[11]:
(see NumPy's data types)
In [252]:
print(type(A))
In [28]:
print(f"number of dimensions={A.ndim}, shape={A.shape}")
In [121]:
# https://stackoverflow.com/questions/56090021/list-comprehension-python-prime-numbers
Primes_less_than_100 = np.array([x for x in range(2,100) if not any([x % y == 0 for y in range(2, int(x/2)+1)])])
Primes_less_than_100
Out[121]:
In [122]:
A = Primes_less_than_100 # This is a copy of pointers
A
Out[122]:
In [123]:
A[0]=1
A
Out[123]:
In [124]:
Primes_less_than_100
Out[124]:
In [125]:
id(A)
Out[125]:
In [126]:
id(Primes_less_than_100)
Out[126]:
In [128]:
A is Primes_less_than_100
Out[128]:
In [129]:
print(Primes_less_than_100.shape)
In [130]:
A = Primes_less_than_100.reshape(5,5)
print(A.shape)
In [131]:
A is Primes_less_than_100
Out[131]:
In [132]:
A.base is Primes_less_than_100
Out[132]:
In [133]:
A.flags.owndata
Out[133]:
In [135]:
print(Primes_less_than_100, A)
In [138]:
A[0,0] = 2
print(Primes_less_than_100, A)
In [55]:
Primes_less_than_100 = np.array([x for x in range(2,100) if not any([x % y == 0 for y in range(2, int(x/2)+1)])])
A = np.copy(Primes_less_than_100)
In [139]:
A is Primes_less_than_100
Out[139]:
In [58]:
Primes_less_than_100[0]
Out[58]:
In [59]:
A[0] = 1
In [60]:
print(Primes_less_than_100[0], A[0])
In [61]:
timeit A=Primes_less_than_100 # This is much faster, depending on the size of the array
In [62]:
%timeit A = np.copy(Primes_less_than_100)
In [82]:
A = np.arange(20).reshape(5, 4)
print(A)
print(A[1, 2]) # [row, column]
In [83]:
print(A[1]) # Get the second row
In [32]:
A[1][2] # [row][column]. Access first to the row, and then, to the column
Out[32]:
In [34]:
timeit A[1][2]
In [35]:
timeit A[1,2] # Access directly to the element
In [156]:
A = np.arange(20).reshape(5, 4)
print(A)
dim_0_coordinates = [0, 1, 2]
dim_1_coordinates = [3, 2, 1]
print(A[dim_0_coordinates, dim_1_coordinates])
In [91]:
print(A[[i for i in range(3)], [i for i in range(3,0,-1)]])
In [38]:
print(np.arange(3))
print(np.arange(3,0,-1))
print(A[np.arange(3), np.arange(3,0,-1)])
In [178]:
A = np.arange(100).reshape(10,10)
print(A)
lst_of_rows = [1, 2, 4]
lst_of_columns = [1, 2, 5]
sub_matrix = A[lst_of_rows][:, lst_of_columns]
print(sub_matrix)
In [13]:
B = A[[0, 1, 2], [0, 1, 2]]
print(B)
In [8]:
print(A)
In [18]:
B[...] = -1
In [19]:
print(B)
In [20]:
print(A)
In [308]:
A = np.arange(20)
print(A, A.shape)
bool_idx = (A>12)
print(bool_idx, bool_idx.shape)
In [309]:
print(A[bool_idx])
In [344]:
A = (100*(0.5-np.random.rand(25))).astype(np.int16).reshape(5,5)
print(A)
print(A[A<0]) # Notice that len(A[A<0]) <= len(A)
In [345]:
A[A<0] = 0
print(A)
In [157]:
for row in A:
print(row)
In [159]:
for element in A.flat:
print(element, end=' ')
In [141]:
A = np.arange(3)
print(A)
A = np.append(A, 4)
print(A)
In [142]:
A = np.append(-1, A)
print(A)
In [145]:
B = np.concatenate((A, A), axis=None)
print(B)
In [146]:
B = np.concatenate((A, A), axis=0)
print(B)
In [148]:
A = np.arange(30).reshape(5,6)
print(A)
In [150]:
B = np.ones(5, dtype=np.int16).reshape(5,1)
print(B)
In [153]:
C = np.hstack((A, B))
print(C)
In [154]:
B = np.ones(7, dtype=np.int16)
print(B)
In [155]:
C = np.vstack((C, B))
print(C)
In [178]:
A = np.arange(20).reshape(5,4)
print(A, A.shape)
In [180]:
print(np.transpose(A), np.transpose(A).shape)
In [179]:
print(A.T, A.T.shape)
In [183]:
A = np.arange(60).reshape(3,5,4)
print(A, A.shape)
In [186]:
print(A.T, A.T.shape)
In [192]:
A = np.arange(5)
print(A, A.shape, A.ndim)
In [195]:
B = A[:, None]
print(B, B.shape, B.ndim)
In [194]:
B = A[None, :]
print(B, B.shape, B.ndim)
In [399]:
B = A[np.newaxis, :]
print(B, B.shape, B.ndim)
In [160]:
A = np.arange(10)
print(A)
In [162]:
print(A[1:3]) # [start:end] (end not included)
In [163]:
print(A[1:5:2]) # [start:end:step]
In [164]:
print(A[:3]) # By default, the first one
In [165]:
print(A[3:])
In [168]:
print(A[::3])
In [170]:
print(A, A[:], A[::])
In [171]:
print(A[::-1])
In [286]:
A = np.arange(50).reshape(5,10)
print(A)
print(A[:]) # <- works, although it is preferable ...
print(A[:,:]) # <- ... this
In [299]:
print(A[0:2,:])
In [298]:
print(A[:,1::2])
In [300]:
print(A[A.shape[0]-2:,A.shape[1]-2:]) # bottom-right 2x2 array:
In [304]:
print(A[2]) # Row extraction
print(A[2:3,:]) # Sub-matrix extraction
In [299]:
A = np.arange(10)
print(A)
B = A[1:3] # B is simply a new view of A
B[:] = 1000
print(B)
print(A)
In [300]:
A = np.arange(10)
print(A)
B = A[::-1]
B[1] = 1000
print(B)
print(A)
In [301]:
A = np.arange(10)
print(A)
B = A[::-1].copy()
B[1] = 1000
print(B)
print(A)
In [396]:
A = np.arange(27).reshape(3,3,3)
print(A)
In [397]:
print(A[1,:,:])
In [398]:
print(A[1,...])
In [187]:
A = np.random.randint(low=-1, high=+2, size=10)
print(A)
print(-A)
In [194]:
B = np.random.randint(low=-1, high=+2, size=10)
print(A)
print(B)
print("-"*31)
print(A+B)
In [198]:
print(A*B)
print(A/B)
print(A//2)
print(A>>1)
print(A*2)
print(A**2)
print(1/A)
print(np.absolute(A))
In [200]:
A = np.arange(10)
print(A)
In [232]:
In [203]:
print("Sum =", np.sum(A))
print("Max =", np.max(A))
print("Min =", np.min(A))
print("L0 norm =", np.linalg.norm(A, ord=0)) # np.max(A)
print("L1 norm =", np.linalg.norm(A, ord=1)) # np.sum(A)
print("L2 norm =", np.linalg.norm(A)) # math.sqrt(sum(A_i**2 for A_i in A))
print("L4 norm =", np.linalg.norm(A, ord=4))
In [248]:
from IPython.display import display, Math
A = np.arange(10)
print(A)
B = A[::-1]
display(Math(r"A \cdot B={}".format(np.dot(A, B))))
print("A dot product B =", np.dot(A, B))
print("sum(A_i*B_i for A_i, B_i in zip(A, B))=", sum(A_i*B_i for A_i, B_i in zip(A, B)))
print("sum(A[:]*B[:]) =", sum(A[:]*B[:]))
In [252]:
A = np.array([1, 2, 3])
B = np.array([-1, -2, -3])
display(Math(r"A \times B={}".format(np.cross(A, B))))
# https://stackoverflow.com/questions/1984799/cross-product-of-two-vectors-in-python
print(f"[A[1]*B[2] - A[2]*B[1], A[2]*B[0] - A[0]*B[2], A[0]*B[1] - A[1]*B[0]] = [{A[1]*B[2] - A[2]*B[1]} {A[2]*B[0] - A[0]*B[2]} {A[0]*B[1] - A[1]*B[0]}]")
In [255]:
A = np.arange(4).reshape(2,2)
print(A)
B = A[::-1]
print(B)
print("A dot product B = ", np.dot(A, B)) # dot[i,j] = sum(A[i,:] * B[:,j])
# https://stackoverflow.com/questions/11033573/difference-between-numpy-dot-and-inner
for i in range(A.shape[0]):
for j in range(A.shape[1]):
print(np.sum(A[i,:]*B[:,j]), end=" ")
print("")
print("A inner product B =", np.inner(A, B)) # inner[i,j] = (A[i,:] * B[j,:])
for i in range(A.shape[0]):
for j in range(A.shape[1]):
print(np.sum(A[i,:]*B[j,:]), end=" ")
print("")
In [213]:
A = np.array([[(i+j)%2 for j in range(10)] for i in range(10)])
print(A, A.shape)
In [101]:
B = np.array([[1] for i in range(10)])
print(B, B.shape)
In [102]:
C = A @ B # Product matrix-matrix
print(C)
In [106]:
print(C.T, C.T.shape, C.shape) # Transpose
In [214]:
print("Determinant =", np.linalg.det(A))
In [221]:
R = np.random.rand(5,5)
iR = np.linalg.inv(R)
print("Inverse =", np.linalg.inv(iR)) # Matrix inverse
In [222]:
print(np.round(R @ iR))
In [223]:
print(R @ iR)
In [224]:
print(np.round(iR @ R))
In [225]:
print(iR @ R)
In [226]:
R = np.random.rand(5,4)
iR = np.linalg.pinv(R) # Pseudo-inverse
print(iR)
In [227]:
print(np.round(R @ iR))
In [228]:
print(R @ iR)
In [231]:
print(np.round(iR @ R))
In [232]:
print(iR @ R)
In vectorized operations, NumPy "extends" scalars and arrays with one of its dimensions equal to the size of the other(s) array(s).
In [256]:
A = np.array([1, 2, 3])
b = 1
print(A+b) # Scalars are broadcasted
In [257]:
B = [1]
print(A+B) # When possible, arrays are broadcasted in all possible dimensions
In [262]:
A = np.ones((5,3), dtype=np.int16)
print(A)
B = np.arange(3)
print(B)
print(A+B) # Broadcasting in the axis 0
In [264]:
print(A)
B = np.arange(5).reshape((5, 1))
print(B)
print(A+B) # Broadcasting in the axis 1
ValueError: frames are not aligned is thrown.
In [127]:
B = np.arange(4)[:, None]
print(B)
In [128]:
print(A.shape)
In [129]:
print(B.shape)
In [130]:
try:
A + B
except ValueError as e:
print("ValueError exception: ", end='')
if hasattr(e, 'message'):
print(e.message)
else:
print(e)
In [266]:
A = np.array([[(i*10+j) for j in range(10)] for i in range(10)])
print(A, A.shape)
B[] to all the rows of A[][] using scalar arithmetic:
In [137]:
C = np.empty_like(A)
def add():
for i in range(A.shape[1]):
for j in range(A.shape[0]):
C[i, j] = A[i, j] + B[j]
%timeit add()
print(C)
B[] to all the rows of B[][] using vectorial computation:
In [138]:
C = np.empty_like(A)
def add():
for i in range(A.shape[1]):
C[i, :] = A[i, :] + B
%timeit add()
print(C)
B[] to all the rows of A[][] using fully vectorial computation:
In [139]:
%timeit C = A + B # <- broadcasting is faster
print(C)
In [140]:
A = np.array(['hello', 'world!'])
print(A)
print(A.shape)
print(np.char.upper(A))
In [321]:
A = np.array([("Spain", 100), ("France", 200), ("Italy", 300)])
print(A) # Notice that all the elements are srings
print(A.shape)
print(A[:,0])
print(A[A[:,0] == "France"])
print(A[A[:,0] == "France"][:,1])
print("The value associated to the key France is", A[A[:,0] == "France"][:,1][0])
In [327]:
%timeit A[A[:,0] == "France"][:,1][0]
dictionary = {"Spain":100, "France":200, "Italy":300}
print(dictionary["France"])
%timeit dictionary["France"]
In [337]:
others = [value for key, value in dictionary.items() if key != "France"]
print(others)
%timeit [value for key, value in dictionary.items() if key != "France"]
In [338]:
print(A[:,0] != "France")
print(A[A[:,0] != "France"])
print(A[A[:,0] != "France"][:,1])
print(A[A[:,0] != "France"][:,1].astype(np.int16))
%timeit A[A[:,0] != "France"][:,1].astype(np.int16)
In [142]:
X = np.array([(1, 2., "Hello"), (3, 4., "World")],
dtype=[("first", "i4"),("second", "f4"), ("third", "S10")])
# See https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
print(X)
In [143]:
print(X["first"])
In [144]:
print(X[0])
In [145]:
print(X["second"])
In [146]:
print(X["third"])
In [147]:
Data = np.array([[1., 200.], [2., 150.], [3., 250.]])
np.savetxt("data.txt", Data)
!cat Data.txt
In [148]:
np.genfromtxt('data.txt')
Out[148]:
In [149]:
ofile = open("data.float64", mode="wb")
Data.tofile(ofile)
In [150]:
print(np.fromfile("data.float64", dtype=np.float64))
In [151]:
!cat create_float64.c
!gcc create_float64.c -o create_float64
!./create_float64
In [152]:
print(np.fromfile("data.float64", dtype=np.float64))
In [153]:
print(np.fromfile("data.float64", dtype=">d"))
# (> = bit-endian, d = double, see https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html)
In [154]:
ofile = open("data.npy", mode="wb")
A = (100*np.random.rand(2,3)).astype(np.uint16)
print(A)
In [155]:
np.save(ofile, A)
In [156]:
!ls save*
In [157]:
print(np.load("data.npy"))
In [ ]:
print(A.dtype.byteorder)
In [2]:
l = list(range(0,100000)); print(type(l), l[:10])
%timeit sum(l)
In [3]:
A = np.arange(0, 100000); print(type(A), A[:10])
%timeit np.sum(A)
In [4]:
!cat sum_array.c
!gcc -O3 sum_array.c -o sum_array
%timeit !./sum_array
In [5]:
# Example extracted from https://github.com/pyHPC/pyhpc-tutorial
lst = range(1000000)
for i in lst[:10]:
print(i, end=' ')
print()
%timeit [i + 1 for i in lst] # A Python list comprehension (iteration happens in C but with PyObjects)
x = [i + 1 for i in lst]
print(x[:10])
arr = np.arange(1000000) # A NumPy list of integers
%timeit arr + 1 # Use operator overloading for nice syntax, now iteration is in C with ints
y = arr + 1
print(y[:10])