From 31f5a04bbc6a0fd03f267e9d41a2466f962152b5 Mon Sep 17 00:00:00 2001
From: Martin Delcourt <martin.delcourt@uclouvain.be>
Date: Tue, 29 Oct 2019 16:43:53 +0100
Subject: [PATCH] Adding Chapter 6 & 7
---
Chapter_6/csvhisto.py | 53 ++++++++++++
Chapter_6/edosolver.py | 86 +++++++++++++++++++
Chapter_6/histogram.py | 76 +++++++++++++++++
Chapter_6/mplhisto.py | 20 +++++
Chapter_6/playing_with_edo.py | 44 ++++++++++
Chapter_7/Vecteur.py | 133 ++++++++++++++++++++++++++++++
Chapter_7/complexe.py | 102 +++++++++++++++++++++++
Chapter_7/not_implemented_expl.py | 76 +++++++++++++++++
8 files changed, 590 insertions(+)
create mode 100644 Chapter_6/csvhisto.py
create mode 100644 Chapter_6/edosolver.py
create mode 100644 Chapter_6/histogram.py
create mode 100644 Chapter_6/mplhisto.py
create mode 100644 Chapter_6/playing_with_edo.py
create mode 100644 Chapter_7/Vecteur.py
create mode 100644 Chapter_7/complexe.py
create mode 100644 Chapter_7/not_implemented_expl.py
diff --git a/Chapter_6/csvhisto.py b/Chapter_6/csvhisto.py
new file mode 100644
index 0000000..8ba596b
--- /dev/null
+++ b/Chapter_6/csvhisto.py
@@ -0,0 +1,53 @@
+from histogram import Histogram
+import csv
+
+class CSVHisto(Histogram):
+
+ def save(self, file_name):
+ """Save histogram as CSV file"""
+ with open(file_name,'w') as file:
+ writer = csv.writer(file)
+ writer.writerow(["bin id","bin center", "n entries"])
+ bin_size = (self._b - self._a)*1./self._n_bins
+ for bin_id, bin_content in enumerate(self._data):
+ writer.writerow([bin_id, self._a+(bin_id+0.5)*bin_size, bin_content])
+
+
+ def load(self, file_name):
+ """Read histogram from CSV file"""
+ self._data = []
+ line_id = 0
+ first_bin_center = 0
+ second_bin_center = 0
+ with open(file_name,'r') as file:
+ reader = csv.reader(file)
+ for row in reader:
+ line_id += 1
+ if line_id == 1:
+ #Skipping line containing column titles
+ continue
+ else:
+ if line_id == 2:
+ first_bin_center = float(row[1])
+ elif line_id == 3:
+ second_bin_center = float(row[1])
+
+ self._data.append(float(row[2]))
+ bin_width = second_bin_center - first_bin_center
+ self._a = first_bin_center-0.5*bin_width
+ self._b = self._a + len(self._data) * bin_width
+ self._n_bins = len(self._data)
+
+
+
+if __name__ == "__main__":
+ h2 = CSVHisto(100,-5,5)
+ from random import gauss
+ for i in range(int(100000)):
+ h2.fill(gauss(0,1))
+ h2.save("test.csv")
+
+ h3 = CSVHisto()
+ h3.load("test.csv")
+ h3.draw()
+
diff --git a/Chapter_6/edosolver.py b/Chapter_6/edosolver.py
new file mode 100644
index 0000000..bf3e484
--- /dev/null
+++ b/Chapter_6/edosolver.py
@@ -0,0 +1,86 @@
+import matplotlib.pyplot as plt
+
+
+class EDOSolver:
+ def __init__(self, f, y_0, t_0, t_max, dt, store_dt = 1):
+ self.f = f
+ self.y_stored = [y_0[:]]
+ self.time = [t_0]
+ self.t_max = t_max
+ self.h = dt
+ self.store_dt = store_dt
+
+ def solve(self):
+ raise NotImplemented
+
+ def draw(self):
+ data = []
+ for i in range(len(self.y_stored[0])):
+ data.append([entry[i] for entry in self.y_stored])
+ plt.plot(self.time,data[i])
+
+ plt.show()
+
+ def draw_diff(self,other_solver):
+ data = []
+ for i in range(len(self.y_stored[0])):
+ data.append([self.y_stored[j][i] - other_solver.y_stored[j][i] for j in range(len(self.y_stored))])
+ plt.plot(self.time,data[i])
+
+ plt.show()
+
+
+class EulerSolver(EDOSolver):
+ def solve(self):
+ t = self.time[0]
+ y = self.y_stored[0][:]
+
+ while abs(t - self.t_max) >= 0.5*self.h and t < self.t_max:
+ if abs(t - (self.time[-1] + self.store_dt)) <= 0.5*self.h:
+ self.time.append(t)
+ self.y_stored.append(y[:])
+
+ dy = self.f(t,y)
+ for i,x in enumerate(dy):
+ y[i] += self.h*x
+ t += self.h
+ self.time.append(t)
+ self.y_stored.append(y[:])
+
+
+
+class RK4Solver(EDOSolver):
+ def solve(self):
+ t = self.time[0]
+ y = self.y_stored[0][:]
+ while abs(t - self.t_max) >= 0.5*self.h and t < self.t_max:
+ if abs(t - (self.time[-1] + self.store_dt)) <= 0.5*self.h:
+ self.time.append(t)
+ self.y_stored.append(y[:])
+
+ dy = self.f(t,y)
+ k1 = [self.h*x for x in dy]
+
+
+ y_k2 = [y[i] + 0.5*k1[i] for i in range(len(k1))]
+ dy_2 = self.f(t+0.5*self.h,y_k2)
+ k2 = [self.h*x for x in dy_2]
+
+
+ y_k3 = [y[i] + 0.5*k2[i] for i in range(len(k2))]
+ dy_3 = self.f(t+0.5*self.h, y_k3)
+ k3 = [self.h*x for x in dy_3]
+
+ y_k4 = [y[i] + k3[i] for i in range(len(k3))]
+ dy_4 = self.f(t+self.h, y_k4)
+ k4 = [self.h*x for x in dy_4]
+
+ y = [y[i] + 1./6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]) for i in range(len(y))]
+ t += self.h
+ self.time.append(t)
+ self.y_stored.append(y[:])
+
+
+
+
+
diff --git a/Chapter_6/histogram.py b/Chapter_6/histogram.py
new file mode 100644
index 0000000..4765178
--- /dev/null
+++ b/Chapter_6/histogram.py
@@ -0,0 +1,76 @@
+
+
+class Histogram:
+ def __init__(self, n_bins = 0, a = 0, b = 0):
+ """Crée un histograme de n_bins réparties uniformément entre a et b"""
+ self._data = [0 for i in range(n_bins)]
+ self._n_bins = n_bins
+ self._a = a
+ self._b = b
+
+ def fill(self, x):
+ """Remplis la bin correspondant à la valeur 'x'"""
+ bin_number = int((x-self._a)*self._n_bins/(self._b-self._a))
+ if bin_number >= 0 and bin_number < self._n_bins:
+ self._data[bin_number] += 1
+ return bin_number
+ else:
+ print("Warning, {} out of bounds.".format(x))
+ return -1
+
+ def get_maximum(self):
+ """Retourne le maximum de l'histogramme"""
+ maximum = self._data[0]
+ for x in self._data[1:]:
+ if x > maximum:
+ maximum = x
+ return maximum
+
+
+ def draw(self):
+ """Affiche l'histogramme en console"""
+ y_max = self.get_maximum()
+ scale = 1.
+ if y_max > 50:
+ scale = 50*1./y_max
+ y_max = 50
+ # vertical axis : 50 pxl
+ y = y_max + 1
+ while y > 0:
+ to_print = ""
+ if (y == y_max+1):
+ to_print = "y ^"
+ else:
+ to_print = " |"
+
+ for x in range(self._n_bins):
+ if self._data[x] >= y/scale:
+ to_print += "#"
+ else:
+ to_print += " "
+ print (to_print)
+ y-=1
+
+ print (" -"+"+"+(self._n_bins-1)*"-"+">")
+ print (" |"+(self._n_bins-1)*" "+"x")
+ if scale != 1:
+ print ("y-axis scale : 1 # = {}".format(1./scale))
+
+
+
+
+if __name__=="__main__":
+ h = Histogram(10,0,1)
+ h.fill(0.5)
+ h.fill(0.5)
+ h.fill(0.5)
+ h.fill(0.1)
+ h.fill(2)
+ h.fill(-1)
+ h.draw()
+
+ h2 = Histogram(100,-5,5)
+ from random import gauss
+ for i in range(int(100000)):
+ h2.fill(gauss(0,1))
+ h2.draw()
diff --git a/Chapter_6/mplhisto.py b/Chapter_6/mplhisto.py
new file mode 100644
index 0000000..98fa1a1
--- /dev/null
+++ b/Chapter_6/mplhisto.py
@@ -0,0 +1,20 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from histogram import Histogram
+
+class MPLHisto(Histogram):
+ def draw(self):
+ bin_size = (self._b - self._a)*1./self._n_bins
+ bins = [self._a + i*bin_size for i in range(len(self._data))]
+ plt.hist(bins[:],bins,weights = self._data)
+ plt.show()
+
+if __name__ == "__main__":
+ h2 = MPLHisto(100,-5,5)
+ from random import gauss
+ for i in range(int(100000)):
+ h2.fill(gauss(0,1))
+ h2.draw()
+
+
+
diff --git a/Chapter_6/playing_with_edo.py b/Chapter_6/playing_with_edo.py
new file mode 100644
index 0000000..8a647a7
--- /dev/null
+++ b/Chapter_6/playing_with_edo.py
@@ -0,0 +1,44 @@
+from edosolver import RK4Solver, EulerSolver
+from math import pi
+def free_fall_acceleration(t,y):
+ # y = [x,v]
+ dv = - 9.81 - 0.00315 * abs(y[1]) * y[1]
+ dx = y[1]
+
+ return [dx,dv]
+
+def oscillation(t,y):
+ return[y[1],-y[0]]
+
+
+
+start = 0
+stop = 30
+dt = 0.01
+dt_store = 0.1
+
+
+
+
+solution_oscillation_euler = EulerSolver(oscillation,[1,0],start,stop,dt,dt_store)
+solution_oscillation_euler.solve()
+solution_oscillation_euler.draw()
+
+solution_oscillation_rk4 = RK4Solver(oscillation,[1,0],start,stop,dt,dt_store)
+solution_oscillation_rk4.solve()
+solution_oscillation_rk4.draw()
+
+solution_oscillation_rk4.draw_diff(solution_oscillation_euler)
+
+
+
+solution_free_fall_euler = EulerSolver(free_fall_acceleration,[0,0], start, stop, dt,dt_store)
+solution_free_fall_euler.solve()
+solution_free_fall_euler.draw()
+
+
+solution_free_fall_rk4 = RK4Solver(free_fall_acceleration,[0,0], start, stop, dt,dt_store)
+solution_free_fall_rk4.solve()
+solution_free_fall_rk4.draw()
+
+solution_free_fall_rk4.draw_diff(solution_free_fall_euler)
diff --git a/Chapter_7/Vecteur.py b/Chapter_7/Vecteur.py
new file mode 100644
index 0000000..995b327
--- /dev/null
+++ b/Chapter_7/Vecteur.py
@@ -0,0 +1,133 @@
+import math
+
+class Vecteur:
+ """Classe vecteur de nombres de dimension arbitraire."""
+
+ def __init__(self, *coord):
+ """Initialise un vecteur de coordonnees "coord" """
+ self._coord = [c for c in coord]
+
+ def get_vec(self):
+ """Retourne un tuple contenant les coordonnees du vecteur"""
+ return tuple(self._coord)
+
+ def set(self,*coord):
+ """Donne les coordonnees "coord" au vecteur"""
+ self._coord = [c for c in coord]
+
+ def set_from_list(self,coord_list):
+ """Donne les coordonnées venant de la liste 'coord_list'."""
+ self._coord = [c for c in coord_list ]
+
+
+ def set_component(self, index, value):
+ """Donne la valeur 'value' a la composante 'index' du vecteur"""
+ self._coord[index] = value
+
+ def norme(self):
+ """Retourne la norme du vecteur"""
+ nn = 0
+ for x in self._coord:
+ nn += x*x
+ return math.sqrt(nn)
+
+ def prod_scal(self, vec):
+ """Retourne le produit scalaire du vecteur avec vec"""
+ if not len(vec._coord) == len(self._coord):
+ print("Error, vectors of different dimensions...")
+ return 0
+ ps = 0
+ for i in range(len(self._coord)):
+ ps+=self._coord[i] * vec._coord[i]
+ return ps
+
+ def prod_vec(self,vec):
+ """Retourne le produit vectoriel du vecteur avec vec"""
+ if not len(vec._coord) == 3 or not len(self._coord)==3:
+ raise TypeError("vector product not possible on vectors of dimension != 3 ({} and {} provided)".format(len(self._coord),len(other._coord)))
+ print("Error, only possible for dim 3 vectors.")
+ return Vecteur(0,0,0)
+ return Vecteur(
+ self._coord[1]*vec._coord[2] - self._coord[2]*vec._coord[1],
+ -self._coord[0]*vec._coord[2] + self._coord[2]*vec._coord[0],
+ self._coord[0]*vec._coord[1] - self._coord[1]*vec._coord[0])
+
+ def __repr__(self):
+ return str(self._coord)
+
+ def __abs__(self):
+ return self.norme()
+
+ def __add__(self,other):
+ try:
+ addition = Vecteur()
+ addition.set_from_list([ x+y for (x,y) in zip(self._coord, other._coord) ])
+ return addition
+ except:
+ return NotImplemented
+
+ def __radd__(self,other):
+ return self + other
+
+ def __mul__(self,other):
+ if isinstance(other,Vecteur):
+ if len(other._coord) != len(self._coord):
+ raise TypeError("Vector size inconsistent ({} and {})".format(len(self._coord),len(other._coord)))
+ else:
+ ps = 0
+ for (x,y) in zip(self._coord, other._coord):
+ ps+= x*y
+ return ps
+ else:
+ try:
+ multiplied = Vecteur()
+ multiplied.set_from_list([other*x for x in self._coord])
+ return multiplied
+ except:
+ return NotImplemented
+
+ def __rmul__(self,other):
+ return self*other
+
+ def __truediv__(self,other):
+ try:
+ divided = Vecteur()
+ divided.set_from_list([x/other for x in self._coord])
+ return divided
+ except:
+ return NotImplemented
+
+
+ def polar_coordinates(self):
+ """Returns the polar coordinates (r,theta) of the vector if the dimension is 2"""
+ if len(self._coord) != 2:
+ raise TypeError("Unable to return polar coordinates. Dimension should be 2 (was {})".format(len(self._coord)))
+ if self._coord[0] == 0:
+ return (abs(self), -math.pi/2 + 2*(self._coord[1] > 0))
+ return (abs(self), math.atan(self._coord[1]/self._coord[0]))
+
+
+
+if __name__ == '__main__':
+ print("Debugging...")
+ u = Vecteur()
+ u.set(1,2,3)
+ print("vecteur u = {} et |u| = {}".format(u,u.norme()))
+
+ v = Vecteur(4,5,6)
+ print("vecteur v = {} et |v| = {}".format(v,v.norme()))
+
+ print("u.v = {}".format(u.prod_scal(v)))
+ print("v.u = {}".format(v.prod_scal(u)))
+ w = u.prod_vec(v)
+ print("w = u x v = {}".format(w))
+ print(" v x u = {}".format(v.prod_vec(u)))
+ print("w.u = {} ; u.w = {}".format(w.prod_scal(u),u.prod_scal(w)))
+ print("w.v = {} ; v.w = {}".format(w.prod_scal(v),v.prod_scal(w)))
+
+ print(u*3)
+ print(3*u)
+ print(u*v)
+ print(u/3)
+
+
diff --git a/Chapter_7/complexe.py b/Chapter_7/complexe.py
new file mode 100644
index 0000000..c1ccf43
--- /dev/null
+++ b/Chapter_7/complexe.py
@@ -0,0 +1,102 @@
+import math
+
+class complexe:
+ def __init__(self,real = 0, imag = 0):
+ """ Initialise complex number with real and imaginary parts respectively set to real and imag"""
+ self.real = real
+ self.imag = imag
+
+ def complexe_conjugue(self):
+ """ Returns the complex conjugate"""
+ return complexe(self.real,-self.imag)
+
+ def __repr__(self):
+ if self.imag >= 0:
+ return "{}+{}i".format(self.real,self.imag)
+ else:
+ return "{}{}i".format(self.real,self.imag)
+
+ def __add__(self,other):
+ if isinstance(other,complexe):
+ return complexe(self.real + other.real, self.imag+other.imag)
+ else:
+ try:
+ return complexe(self.real+other, self.imag)
+ except:
+ return NotImplemented
+
+ def __radd__(self,other):
+ return self + other
+
+ def __mul__(self,other):
+ if isinstance(other,complexe):
+ return complexe(self.real*other.real - self.imag * other.imag,
+ self.real*other.imag + self.imag * other.real)
+ else:
+ try:
+ return complexe(other*self.real,other*self.imag)
+ except:
+ return NotImplemented
+
+ def __rmul__(self,other):
+ return self*other
+
+ def __abs__(self):
+ """Returns the module of the complex number"""
+ return math.sqrt(self.real**2 + self.imag**2)
+
+ def __truediv__(self,other):
+ if isinstance(other,complexe):
+ return self*self.complexe_conjugue()*1./(abs(self)**2)
+ else:
+ try:
+ return complexe(self.real/other,self.imag/other)
+ except:
+ return NotImplemented
+
+ def __rtruediv__(self,other):
+ return self/other
+
+ def __eq__(self,other):
+ if isinstance(other,complexe):
+ return self.real == other.real and self.imag == other.imag
+ else:
+ try:
+ return self.real == other and self.imag == 0
+ except:
+ return NotImplemented
+ def __ne__(self,other):
+ return not self == other
+
+ def __lt__(self,other):
+ if isinstance(other,complexe):
+ return abs(self) < abs(other)
+ else:
+ try:
+ return abs(self) < other
+ except:
+ return NotImplemented
+ def __le__(self,other):
+ if isinstance(other,complexe):
+ return abs(self) <= abs(other)
+ else:
+ try:
+ return abs(self) <= other
+ except:
+ return NotImplemented
+
+ def __gt__(self,other):
+ return other < self
+
+ def __ge__(self,other):
+ return other <= self
+
+ def polar(self):
+ if self.real == 0:
+ return (abs(self),math.pi/2)
+ else:
+ return (abs(self),math.atan(self.imag/self.real))
+
+
+if __name__ == "__main__":
+ print("For explanation on all the 'NotImplemented in the code, please check-out not_implemented_expl.py")
diff --git a/Chapter_7/not_implemented_expl.py b/Chapter_7/not_implemented_expl.py
new file mode 100644
index 0000000..e44b830
--- /dev/null
+++ b/Chapter_7/not_implemented_expl.py
@@ -0,0 +1,76 @@
+
+print("Imaginons deux programmeurs python différents.")
+print("Le premier tente de construire une classe 'vecteur à deux dimension' rudimentaire")
+print("Il crée la classe 'first_class' et définit l'addition entre deux vecteurs")
+print("Et il a pensé à rajouter le 'NotImplemented' dans ses opérations.")
+
+
+
+class first_class:
+ def __init__(self, x, y):
+ self.x = x
+ self.y = y
+
+ def __add__(self,other):
+ try:
+ return(first_class(self.x+other.x,self.y+other.y))
+ except:
+ print("Cannot use first_class.__add__ ---> returning NotImplemented")
+ return NotImplemented
+
+ def __radd__(self,other):
+ return self+other
+
+
+print("Un second programmeur crée une autre classe, 'second_class' qui a pour seul but de compter le nombre de fois qu'un tel objet aura été dans une addition")
+print("On voudrait que cet objet puisse être additionné avec n'importe quel autre objet.")
+
+class second_class:
+ def __init__(self):
+ self.addition_counter = 0
+
+ def __repr__(self):
+ return "Was added {} time(s)".format(self.addition_counter)
+
+ def __add__(self,other):
+ new_addition = second_class()
+ new_addition.addition_counter = self.addition_counter + 1
+ return new_addition
+
+ def __radd__(self,other):
+ print("second_class.__radd__ called")
+ return self + other
+
+
+print("En particulier, on voit que ça devrait marcher avec la classe de l'autre programmeur\n\n")
+v1 = first_class(0,0)
+counter = second_class()
+print("v1 + counter va appeler v1.__add__, mais v1 ne sait pas que counter existe !")
+print("Il renvoie donc 'NotImplemented' et on appelle counter.__radd__, qui lui fonctionne bien.")
+counter = v1 + counter
+print(counter)
+print("\n\n\n")
+print("Maintenant, regardons ce même exemple si le premier programmeur n'avait pas renvoyé 'NotImplemented'")
+
+
+class first_class_without_protection:
+ def __init__(self, x, y):
+ self.x = x
+ self.y = y
+
+ def __add__(self,other):
+ print("Calling first_class_without_protection.__add__")
+ return(first_class_without_protection(self.x+other.x,self.y+other.y))
+
+ def __radd__(self,other):
+ return self+other
+
+
+
+v2 = first_class_without_protection(0,0)
+counter2 = second_class()
+print("Dans ce cas-ci, la méthode __add__ de la première classe va planter lors de l'addition, sans laisser la possibilité à second_class d'effectuer __radd__")
+counter = v2 + counter
+
+
+
--
GitLab