Adding Chapter 6 & 7

Chapter_6/csvhisto.py

+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()
+        

Chapter_6/edosolver.py

+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[:])
+            
+            
+            
+        
+

Chapter_6/histogram.py

+
+
+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()

Chapter_6/mplhisto.py

+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()
+
+
+

Chapter_6/playing_with_edo.py

+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)

Chapter_7/Vecteur.py

+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)
+
+

Chapter_7/complexe.py

+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")

Chapter_7/not_implemented_expl.py

+
+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
+
+
+