Exercise solution¶

In [1]:
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.1112

alternative solution

In [2]:
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.1308

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 [3]:
import numpy as np
a = np.array([1, 2, 3, 4], dtype=float)
a
Out[3]:
array([1., 2., 3., 4.])
In [4]:
a = range(1000)
%timeit [ i**2 for i in a]

25 µs ± 136 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [5]:
b = np.arange(1000)
%timeit b**2

533 ns ± 1.84 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)
In [6]:
c = np.array([[1,2],[3,4]])
c
Out[6]:
array([[1, 2],
       [3, 4]])
In [7]:
c.ndim
Out[7]:
2
In [8]:
c.shape
Out[8]:
(2, 2)
In [9]:
c = np.arange(27)
c.reshape((3,3,3))
Out[9]:
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 [10]:
np.zeros((2,2))
Out[10]:
array([[0., 0.],
       [0., 0.]])
In [11]:
np.ones((2,1))
Out[11]:
array([[1.],
       [1.]])
In [12]:
a = np.arange(27).reshape((3,3,3))
np.ones_like(a)
Out[12]:
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 [13]:
np.eye(3)
Out[13]:
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
In [14]:
a = np.arange(10)
a
Out[14]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [15]:
a[0]
Out[15]:
0
In [16]:
a[-1]
Out[16]:
9
In [17]:
a[0:3]
Out[17]:
array([0, 1, 2])
In [18]:
a[::2]
Out[18]:
array([0, 2, 4, 6, 8])
In [19]:
a = a.reshape(5,2)
a
Out[19]:
array([[0, 1],
       [2, 3],
       [4, 5],
       [6, 7],
       [8, 9]])
In [20]:
a[3,1]
Out[20]:
7
In [21]:
a[2,:]
Out[21]:
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 [22]:
a = np.arange(9)
In [23]:
b = a.reshape((3,3))
In [24]:
np.shares_memory(a,b)
Out[24]:
True
In [25]:
a[3] = -1
In [26]:
b
Out[26]:
array([[ 0,  1,  2],
       [-1,  4,  5],
       [ 6,  7,  8]])
In [27]:
np.shares_memory(a,b)
Out[27]:
True
In [28]:
b = a.copy()
np.shares_memory(a,b)
Out[28]:
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 [29]:
ene = np.random.exponential(size=10, scale=10.)
ene
Out[29]:
array([12.54888052,  0.96077347, 14.13114011,  5.59390883,  4.23607211,
       14.23889928,  2.02597793, 20.57168657,  5.11088435, 25.08232079])
In [30]:
mask = ene > 3
mask
Out[30]:
array([ True, False,  True,  True,  True,  True, False,  True,  True,
        True])
In [31]:
ene[mask]
Out[31]:
array([12.54888052, 14.13114011,  5.59390883,  4.23607211, 14.23889928,
       20.57168657,  5.11088435, 25.08232079])
In [32]:
ene[ene<3]
Out[32]:
array([0.96077347, 2.02597793])
In [33]:
ene[ene<3] = 0
ene
Out[33]:
array([12.54888052,  0.        , 14.13114011,  5.59390883,  4.23607211,
       14.23889928,  0.        , 20.57168657,  5.11088435, 25.08232079])

Simple Operations¶

In [34]:
a = np.arange(4)
a
Out[34]:
array([0, 1, 2, 3])
In [35]:
a+1
Out[35]:
array([1, 2, 3, 4])
In [36]:
10**a
Out[36]:
array([   1,   10,  100, 1000])
In [37]:
np.sin(a)
Out[37]:
array([0.        , 0.84147098, 0.90929743, 0.14112001])

matrix multiplication

In [38]:
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 [39]:
a = np.random.randint(low=0,high=10,size=4)
a
Out[39]:
array([9, 5, 6, 5])
In [40]:
np.sum(a)
Out[40]:
25
In [41]:
np.max(a), np.min(a)
Out[41]:
(9, 5)
In [42]:
np.argmax(a), np.argmin(a)
Out[42]:
(0, 1)
In [43]:
np.mean(a), np.median(a), np.std(a)
Out[43]:
(6.25, 5.5, 1.6393596310755)
In [44]:
a = a.reshape(2,2)
a
Out[44]:
array([[9, 5],
       [6, 5]])
In [45]:
np.sum(a,axis=1)
Out[45]:
array([14, 11])
In [46]:
np.all(ene<3)
Out[46]:
False
In [47]:
np.any(ene<3)
Out[47]:
True

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

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

y = 1
f(y)
print(y)

1

this is NOT true for numpy arrays!

As in the function is not

In [49]:
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 [50]:
rng = np.random.default_rng(1234)

example of usage:

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