Python course University of Heidelberg Summer 2024

Solutions Day 03 (April 10th)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

09.01

In [2]:
np.logspace(-20,-10,11)
Out[2]:
array([1.e-20, 1.e-19, 1.e-18, 1.e-17, 1.e-16, 1.e-15, 1.e-14, 1.e-13,
       1.e-12, 1.e-11, 1.e-10])
In [3]:
np.array(10*[2.])
#
2*np.ones(10)
Out[3]:
array([2., 2., 2., 2., 2., 2., 2., 2., 2., 2.])
In [4]:
np.zeros(5, dtype=np.float32)
Out[4]:
array([0., 0., 0., 0., 0.], dtype=float32)

09.02

In [5]:
x=np.random.poisson(100,10)
x[1:] - x[:-1]
Out[5]:
array([  3,   3,   9, -29,  17, -18,  14,   6,   9])

09.03

In [6]:
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
                        unpack=True)
temp_red=temp[(temp>-90) & (temp<90)]
print (len(temp)-len(temp_red),'elements were removed')

60 elements were removed

10.01

In [7]:
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
                        unpack=True)
keep=(temp>-90) & (temp<90)
In [8]:
plt.plot(year[keep],temp[keep], marker='.',linewidth=0)
Out[8]:
[<matplotlib.lines.Line2D at 0x7fe0dc2a0850>]
In [9]:
plt.plot(year[keep]%1,temp[keep], marker='.',linewidth=0)
Out[9]:
[<matplotlib.lines.Line2D at 0x7fe0dc1c9f40>]

10.02

In [10]:
# define normalized Gaussian with mean 'mu' and sigma 'sig'
def gaussian(x,mu,sig):
    return (1/(np.sqrt(2*np.pi)*sig))*np.exp(-(x-mu)**2/(2*sig**2))
In [11]:
#  Gaussian with mean 10 and sigma sqrt(10)
x=np.linspace(-5,30,100)
d=np.random.normal(7,np.sqrt(7),10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))

plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
Out[11]:
[<matplotlib.lines.Line2D at 0x7fe0dc12a7f0>]
In [12]:
#  I plot here the sample from the Poisson distribution together with the Gaussian
x=np.linspace(-5,30,100)
d=np.random.poisson(7,10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))

plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
Out[12]:
[<matplotlib.lines.Line2D at 0x7fe0dc061880>]

10.03

In [13]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [14]:
x=np.arange(-10,10,0.1)
y=np.arange(-10,10,0.1)

# Note we calculate the meshgrid first here
#    this makes the computations below
#    more readable
xx,yy = np.meshgrid(x,y)


r1sq=xx**2+(yy-5)**2          # work out radii
r2sq=xx**2+(yy+5)**2

b1x,b1y   = - (yy-5)/r1sq, xx/r1sq # calculate magnetic field vectors
b2x,b2y   = + (yy+5)/r2sq, -xx/r2sq

plt.streamplot(x,y,b1x+b2x,b1y+b2y,color='k')
Out[14]:
<matplotlib.streamplot.StreamplotSet at 0x7fe0cff90a60>

Solution for the coil

In [15]:
a=5                        #radius of the coil
n=90                       #elements on the coil

x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n

xx,yy=np.meshgrid(x,y)     #generate coordinate grid
bx=0*xx
by=0*yy

for phi in np.arange(0,2*np.pi,dphi):
    cos=np.cos(phi)
                           #separation coil element-observer
    rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+yy**2) 

                           #magnetic field components (cross product)
    bx+=-a*cos*yy*dphi/rp**3       
    by+=-a*(a-cos*xx)*dphi/rp**3

plt.streamplot(x,y,bx,by,color='k')
Out[15]:
<matplotlib.streamplot.StreamplotSet at 0x7fe0ce6277f0>
In [16]:
# Example "Helmholtz coils"

a=5                        #radius of the coil
n=90                       #elements on the coil

x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n

xx,yy=np.meshgrid(x,y)     #generate coordinate grid
bx=0*xx
by=0*yy

for phi in np.arange(0,2*np.pi,dphi):
    cos=np.cos(phi)
                           #separation coil element-observer
    rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy-a/2)**2) 
    rp2=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy+a/2)**2) 

                           #magnetic field components (cross product)
    bx+=-a*cos*(yy-a/2)*dphi/rp**3       
    by+=-a*(a-cos*xx)*dphi/rp**3
    bx+=-a*cos*(yy+a/2)*dphi/rp2**3       
    by+=-a*(a-cos*xx)*dphi/rp2**3

plt.streamplot(x,y,bx,by,color='k')
Out[16]:
<matplotlib.streamplot.StreamplotSet at 0x7fe0ce3a5a00>

11.1

In [17]:
import glob
import os
for file in glob.glob('data/*'):
    print (str(int(os.path.getsize(file)/1024+0.5))+'kB', file)

122kB data/munich_temperatures_average.txt
0kB data/data_new.txt
42kB data/germany.txt
4kB data/off
73kB data/autofahrt_2018.txt
2149kB data/ice_data_2019_area.npy
4kB data/ice_data
9kB data/autofahrt.txt
2kB data/problem_1_question_6.fasta
3kB data/UID_0031246_RVC_001.tbl
1kB data/problem_1_question_5.fasta
2149kB data/ice_data_area.npy
9kB data/autofahrt_nohd.txt
0kB data/problem_1_codons.txt
2kB data/problem_1_question_4_new.fasta
123kB data/munich_temperatures_average_with_bad_data.txt
5kB data/new_decay_data.txt
0kB data/columns.txt
129kB data/autofahrt2019.txt
9kB data/autofahrt_data.txt
1kB data/data.txt
3kB data/decay_data.txt

12.1, 12.2

In [18]:
def date2str(li):
    [year,month,day,hour,minute,second] = li
    date="{0:04d}-{1:02d}-{2:02d}".format(year,month,day)
    time="{0:02d}:{1:02d}:{2:02d}".format(hour,minute,second)
    s="{0} {1}".format(date,time)
    return s

def str2date(s):
    date, time =s.split()
    year,month,day = date.split('-')
    hour,minute,second = time.split(':')
    return [int(year),int(month),int(day),int(hour),int(minute),int(second)]

# foward and backward
print (str2date(date2str([2017,4,5,17,43,12])))
print (date2str(str2date('2017-04-05 17:43:12')))

[2017, 4, 5, 17, 43, 12]
2017-04-05 17:43:12

Monte Carlo Error propagation

In [19]:
# Defined constant(s)
G = 6.67384e-11  # SI units

# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 0.05e4
mean_m_2 = 30.e4
sigma_m_2 = 0.1e4
mean_r = 3.2
sigma_r = 0.01

Standard error propagation

In [20]:
# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_1) ** 2
                           + (sigma_m_2 / mean_m_2) ** 2
                           + 4. * (sigma_r / mean_r) ** 2)
print(mean_f, sigma_f)

0.7820906249999999 0.005691093987209795

Generate random samples for m_1, m_2 and r and calculate the force

In [21]:
N = 1000000
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)

# calculate force
F = G * m_1 * m_2 / r ** 2

print(np.mean(F), np.std(F))

0.7821140402913825 0.00562528538900065

Plot the Monte-Carlo histogram using the option density=True and the Gaussian from the error propagation

In [22]:
# Define range of output values for plotting
xmin = 0.75
xmax = 0.82

# use Gaussian definition from above!

x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)

plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
Out[22]:
(0.75, 0.82)

Repeat this with larger error bars:

In [23]:
# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 2.e4
mean_m_2 = 30.e4
sigma_m_2 = 10.e4
mean_r = 3.2
sigma_r = 1.0

# Define range of output values for plotting
xmin = -1.
xmax = 5.

# Define number of samples
N = 1000000

## STANDARD ERROR PROPAGATION

# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_2) ** 2
                           + (sigma_m_2 / mean_m_2) ** 2
                           + 2. * (sigma_r / mean_r) ** 2)

# Define Gaussian function
x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)

## MONTE-CARLO ERROR PROPAGATION

# Sample from initial values
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)

# Compute final values
F = G * m_1 * m_2 / r ** 2

## PLOTTING

plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
Out[23]:
(-1.0, 5.0)