Exercise solution¶

In [3]:
import random
from math import sqrt

def f(x):
    return sqrt(1-x**2)

n = 0
N = 10000
for i in range(N):
    x = random.uniform(-1.,1.)
    y = random.uniform(0.,1.)
    if y < f(x):
        n += 1
res = n/N * 4
print(res)

3.1348

alternative solution

In [4]:
import random
from math import sqrt

n = 0
N = 10000
for i in range(N):
    x = random.uniform(0.,1.)
    y = random.uniform(0.,1.)
    if (x**2 + y**2) <1 :
        n += 1
res = n/N * 4
print(res)

3.14

Scientific python stack¶

Get full documentation here.
Get a tutorial here

The SciPy library contains several packages to perform specialized scientific calculations:

  • Special functions (scipy.special)
  • Integration (scipy.integrate)
  • Optimization (scipy.optimize)
  • Interpolation (scipy.interpolate)
  • Fourier Transforms (scipy.fftpack)
  • Signal Processing (scipy.signal)
  • Linear Algebra (scipy.linalg)
  • Sparse Eigenvalue Problems with ARPACK
  • Compressed Sparse Graph Routines (scipy.sparse.csgraph)
  • Spatial data structures and algorithms (scipy.spatial)
  • Statistics (scipy.stats)
  • Multidimensional image processing (scipy.ndimage)
  • File IO (scipy.io)

Numpy¶

It is the foundation of python scientific stack.
The basic building block is the numpy.array data structure. It can be used as a python list of numbers, but it is a specialized efficient way of manipulating numbers in python.

In [5]:
import numpy as np
a = np.array([1, 2, 3, 4], dtype=float)
a
Out[5]:
array([1., 2., 3., 4.])
In [59]:
a = range(1000)
%timeit [ i**2 for i in a]

313 µs ± 10.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [7]:
b = np.arange(1000)
%timeit b**2

1.16 µs ± 20.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
In [8]:
c = np.array([[1,2],[3,4]])
c
Out[8]:
array([[1, 2],
       [3, 4]])
In [9]:
c.ndim
Out[9]:
2
In [10]:
c.shape
Out[10]:
(2, 2)
In [11]:
c = np.arange(27)
c.reshape((3,3,3))
Out[11]:
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]],

       [[18, 19, 20],
        [21, 22, 23],
        [24, 25, 26]]])
In [12]:
np.zeros((2,2))
Out[12]:
array([[0., 0.],
       [0., 0.]])
In [13]:
np.ones((2,1))
Out[13]:
array([[1.],
       [1.]])
In [14]:
a = np.arange(27).reshape((3,3,3))
np.ones_like(a)
Out[14]:
array([[[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],

       [[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]],

       [[1, 1, 1],
        [1, 1, 1],
        [1, 1, 1]]])
In [15]:
np.eye(3)
Out[15]:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
In [16]:
a = np.arange(10)
a
Out[16]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [17]:
a[0]
Out[17]:
0
In [18]:
a[-1]
Out[18]:
9
In [19]:
a[0:3]
Out[19]:
array([0, 1, 2])
In [20]:
a[::2]
Out[20]:
array([0, 2, 4, 6, 8])
In [21]:
a = a.reshape(5,2)
a
Out[21]:
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7],
       [8, 9]])
In [22]:
a[3,1]
Out[22]:
7
In [23]:
a[2,:]
Out[23]:
array([4, 5])

Important: Views or copy¶

A slice or reshape is a view, simply a re-organization of the same data in memory, thus changing one element changes the same element in all views

In [24]:
a = np.arange(9)
In [25]:
b = a.reshape((3,3))
In [26]:
np.shares_memory(a,b)
Out[26]:
True
In [27]:
a[3] = -1
In [28]:
b
Out[28]:
array([[ 0,  1,  2],
       [-1,  4,  5],
       [ 6,  7,  8]])
In [29]:
np.shares_memory(a,b)
Out[29]:
True
In [30]:
b = a.copy()
np.shares_memory(a,b)
Out[30]:
False

Boolean masks for extracting values¶

A typical operation done in your daily physics data analysis is to extract from an array the values that match a condition. Consider an array of the energies of particles, and assume you want to use only the energies above a given threshold. Boolean masking comes at a rescue

In [31]:
ene = np.random.exponential(size=10, scale=10.)
ene
Out[31]:
array([1.84417793e+01, 1.01995057e-01, 1.15973376e+01, 3.91436666e-04,
       1.13769252e+01, 8.49667565e+00, 4.97730541e+00, 8.38086765e+00,
       1.65839490e+01, 5.95578226e+00])
In [32]:
mask = ene > 2
mask
Out[32]:
array([ True, False,  True, False,  True,  True,  True,  True,  True,
        True])
In [33]:
ene[mask]
Out[33]:
array([18.44177932, 11.59733762, 11.37692521,  8.49667565,  4.97730541,
        8.38086765, 16.58394899,  5.95578226])
In [34]:
ene[ene<2]
Out[34]:
array([0.10199506, 0.00039144])
In [35]:
ene[ene<2] = 0
ene
Out[35]:
array([18.44177932,  0.        , 11.59733762,  0.        , 11.37692521,
        8.49667565,  4.97730541,  8.38086765, 16.58394899,  5.95578226])

Values with list of indexes¶

Similarly to boolean masks, it is possible to access and modify values directly to an array using a list of indexes

In [36]:
status = np.random.randint(low=0,high=10,size=10)
status
Out[36]:
array([2, 2, 2, 6, 2, 0, 2, 4, 8, 5])
In [37]:
status[[0, 3, 5]]
Out[37]:
array([2, 6, 0])
In [38]:
status[[0, 3, 5]] = -1
status
Out[38]:
array([-1,  2,  2, -1,  2, -1,  2,  4,  8,  5])

Simple Operations¶

In [39]:
a = np.arange(4)
a
Out[39]:
array([0, 1, 2, 3])
In [40]:
a+1
Out[40]:
array([1, 2, 3, 4])
In [41]:
10**a
Out[41]:
array([   1,   10,  100, 1000])
In [42]:
np.sin(a)
Out[42]:
array([0.        , 0.84147098, 0.90929743, 0.14112001])

matrix multiplication

In [43]:
a = np.arange(1,5).reshape(2,2)
b = np.arange(5,9).reshape(2,2)
print(a)
print(b)
print("result: \n",a@b)

[[1 2]
 [3 4]]
[[5 6]
 [7 8]]
result: 
 [[19 22]
 [43 50]]

Reductions¶

In [44]:
a = np.random.randint(low=0,high=10,size=4)
a
Out[44]:
array([2, 4, 0, 9])
In [45]:
np.sum(a)
Out[45]:
15
In [46]:
np.max(a), np.min(a)
Out[46]:
(9, 0)
In [47]:
np.argmax(a), np.argmin(a)
Out[47]:
(3, 2)
In [48]:
np.mean(a), np.median(a), np.std(a)
Out[48]:
(3.75, 3.0, 3.344772040064913)
In [49]:
a = a.reshape(2,2)
a
Out[49]:
array([[2, 4],
       [0, 9]])
In [50]:
np.sum(a,axis=1)
Out[50]:
array([6, 9])
In [51]:
m1 = a>3
m1
Out[51]:
array([[False,  True],
       [False,  True]])
In [52]:
np.all(m1)
Out[52]:
False
In [53]:
np.any(m1)
Out[53]:
True

do you remember that in python function arguments are passed by copy?

In [54]:
def f(x):
    x+=1

y = 1
f(y)
print(y)

1

this is NOT true for numpy arrays!

In [55]:
def f(x):
    x+=1

y = np.array([1])
f(y)
print(y)

[2]

numpy arrays are passed by reference, pay attention!

Exercise¶

write another code to compute pi using numpy (avoid to use a for loop)

initialize a pseudo-random number generator with:

In [56]:
rng = np.random.default_rng(1234)

example of usage:

In [58]:
rng.uniform(low=0., high=1., size=10)
Out[58]:
array([0.97669977, 0.38019574, 0.92324623, 0.26169242, 0.31909706,
       0.11809123, 0.24176629, 0.31853393, 0.96407925, 0.2636498 ])