Éléments de corrigé

Première partie - questions 1 à 11

# Introduction to visualization
# Question 1
import matplotlib.pyplot as plt

def first_test() :
  plt.plot([1,2,3,4])  # Single array: x coord is 0..n
  plt.show()

  # Two plots on the same figure
  plt.plot([1,2,3,4], [1,4,9,16])
  plt.plot([1,2,3,4], [1,8,27,64])
  plt.show()

  # Two plots on the same figure, alternative version
  plt.plot([1,2,3,4], [1,4,9,16], [1,2,3,4], [1,8,27,64])
  plt.show()

# first_test()

# Title and so on
import numpy as np

def example() :
  city=['Delhi','Beijing','Washington','Tokyo','Moscow']
  Gender=['Male','Female']
  pos = np.arange(len(city))
  bar_width = 0.35
  Happiness_Index_Male=[60,40,70,65,85]
  Happiness_Index_Female=[30,60,70,55,75]
  plt.bar(pos, Happiness_Index_Male, bar_width, color='blue', edgecolor='black')
  plt.bar(pos+bar_width, Happiness_Index_Female, bar_width, color='pink', edgecolor='black')
  plt.xticks(pos, city)
  plt.xlabel('City', fontsize=16)
  plt.ylabel('Happiness_Index', fontsize=16)
  plt.title('Group Barchart - Happiness index across cities By Gender', fontsize=18)
  plt.legend(Gender,loc=2)
  plt.show()

# example()

# Visualizing functions
# Question 2
import math

def center_axes(axes):
  # Move left y-axis and bottom x-axis to centre, passing through (0,0)
  axes.spines['left'].set_position('zero')
  axes.spines['bottom'].set_position('zero')
  # Eliminate upper and right axes
  axes.spines['right'].set_color('none')
  axes.spines['top'].set_color('none')
  # Show ticks in the left and lower axes only
  axes.xaxis.set_ticks_position('bottom')
  axes.yaxis.set_ticks_position('left')

# Without numpy
def plt_sine() :
  x = [xx/100 for xx in range(-200, 201, 2)]
  fy = [math.sin(math.pi * xx) for xx in x]
  gy = [math.sin(2 * math.pi * xx) for xx in x]
  #
  center_axes(plt.axes())
  plt.plot(x, fy, x, gy)
  plt.show()

# Without numpy
def plt_sine_np() :
  x = np.linspace(-2, 2, 200)
  fy = np.sin(math.pi * x)
  gy = np.sin(2 * math.pi * x)
  #
  center_axes(plt.axes())
  plt.plot(x, fy, x, gy)
  plt.show()

# plt_sine()
# plt_sine_np()

# Question 3
def two_figures() :
  x = np.linspace(-2, 2, 200)
  plt.subplot(211) # 2 rows, 1 column, figure #1
  plt.plot(x, np.sin(math.pi * x))
  plt.subplot(212) # 2 rows, 1 column, figure #2
  plt.plot(x, np.sin(2 * math.pi * x))
  plt.show()

# two_figures()

# Question 4
def two_figures_style() :
  x = np.linspace(-2, 2, 200)
  plt.subplot(211) # 2 rows, 1 column, figure #1
  plt.plot(x, np.sin(math.pi * x), 'rv')
  plt.subplot(212) # 2 rows, 1 column, figure #2
  plt.plot(x, np.sin(2 * math.pi * x), 'b+')
  plt.show()

# two_figures_style()

# Digits of pi
# Question 5

def plain_red_square() :
  n = 40
  img = [[[128, 29, 58] for j in range(n)] for i in range(n)]
  plt.imshow(img)
  plt.axis('off') # Do not display the axes
  plt.show()

# plain_red_square()

# Question 6
colors = [
  [230 ,230 ,230],  # white
  [255 ,65 ,54],    # red
  [255 ,53 ,27],    # orange
  [255 ,220 ,0],    # yellow
  [1 ,255 ,112],    # light green
  [61 ,153 ,112],   # dark green
  [57 ,204 ,204],   # light blue
  [0 ,116 ,217],    # dark blue
  [177 ,13 ,201],   # purple
  [240 ,18 ,190]    # pink
]

def color_square_digits(dig, n):
  img=[[[0,0,0] for i in range(n)] for j in range(n)]
  i = 0
  j = 0
  for k in range(len(dig)):
    img[i][j] = colors[int(dig[k])]
    j = (j + 1) % n
    if j == 0:
      i = (i + 1) % n
      if i == 0 :
        break
  plt.imshow(img)
  plt.axis('off')
  plt.show()

def frac(a, b, n) :
  """
  return the n first digits of the fractional part of a/b
  """
  digits = ""
  for d in range(n) :
    a *= 10
    digits += str((a // b) % 10)
  return digits

# color_square_digits(frac(1,7,1600), 40)

# Question 7
big_pi = """\
31415926535897932384626433832795028841971693993751058209749445923078164062862\
089986280348253421170679821480865132823066470938446095505822317253594081284811\
174502841027019385211055596446229489549303819644288109756659334461284756482337\
867831652712019091456485669234603486104543266482133936072602491412737245870066\
063155881748815209209628292540917153643678925903600113305305488204665213841469\
519415116094330572703657595919530921861173819326117931051185480744623799627495\
673518857527248912279381830119491298336733624406566430860213949463952247371907\
021798609437027705392171762931767523846748184676694051320005681271452635608277\
857713427577896091736371787214684409012249534301465495853710507922796892589235\
420199561121290219608640344181598136297747713099605187072113499999983729780499\
510597317328160963185950244594553469083026425223082533446850352619311881710100\
031378387528865875332083814206171776691473035982534904287554687311595628638823\
537875937519577818577805321712268066130019278766111959092164201989380952572010\
654858632788659361533818279682303019520353018529689957736225994138912497217752\
834791315155748572424541506959508295331168617278558890750983817546374649393192\
550604009277016711390098488240128583616035637076601047101819429555961989467678\
374494482553797747268471040475346462080466842590694912933136770289891521047521\
620569660240580381501935112533824300355876402474964732639141992726042699227967\
823547816360093417216412199245863150302861829745557067498385054945885869269956\
909272107975093029553211653449872027559602364806654991198818347977535663698074\
265425278625518184175746728909777727938000816470600161452491921732172147723501\
414419735685481613611573525521334757418494684385233239073941433345477624168625\
189835694855620992192221842725502542568876717904946016534668049886272327917860\
857843838279679766814541009538837863609506800642251252051173929848960841284886\
269456042419652850222106611863067442786220391949450471237137869609563643719172\
"""

# color_square_digits(big_pi, 40)
# color_square_digits(str(math.pi).replace('.', ''), 40)


# Koch snowflake
import cmath  # complex numbers

# Question 8 to 11
def iso_triangle() :
  a = complex(0,0)
  b = complex(1,0)
  ac = (b-a) * cmath.rect(1, math.pi/3)
  c = a + ac
  return ([a.real, c.real, b.real, a.real],
          [a.imag, c.imag, b.imag, a.imag])

# plt.axes().set_aspect('equal')
# isot = iso_triangle()
# plt.plot(isot[0], isot[1])
# plt.show()

def koch_seg(segx, segy, pos = 0) :
  # We use complex numbers to represent points
  a = complex(segx[pos], segy[pos])         # get point A
  b = complex(segx[pos + 1], segy[pos + 1]) # and point B
  ab = b - a      # Compute vector AB
  ac = ab / 3     # Cut the segment at 1/3 of its length
  c = a + ac      # to find point C
  # The next point is at 1/3 of AB, turned by 60° (pi/3), relative to C
  cd = ac * cmath.rect(1, math.pi/3)  # Can also be written ac * complex(math.cos(math.pi/3), math.sin(math.pi/3))
  d = c + cd
  # The last point is at 2/3 of AB
  ae = 2 * ab / 3
  e = a + ae
  return ([a.real, c.real, d.real, e.real, b.real],
          [a.imag, c.imag, d.imag, e.imag, b.imag])

def koch_iter(segsx, segsy) :
  finalx = []
  finaly = []
  for pos in range(len(segsx) - 1) :
    nsegs = koch_seg(segsx, segsy, pos)
    finalx = finalx + nsegs[0]
    finaly = finaly + nsegs[1]
  return (finalx, finaly)

def vonKoch(n) :
  k = iso_triangle()
  for i in range(n) :
    k = koch_iter(k[0], k[1])
  plt.axes().set_aspect('equal')
  plt.plot(k[0], k[1])
  plt.show()

# vonKoch(3)
# vonKoch(5)

Deuxième partie - questions 12 à 27

# Pi

# Question 12
import random

def approx_pi(n) :
  inside = 0
  total = 0
  while total < n :
    x = random.uniform(-1, 1)
    y = random.uniform(-1, 1)
    total += 1
    if math.hypot(x, y) < 1 :
      inside += 1
  almost_pi = 4 * (inside / total)
  rel_error = math.fabs(almost_pi - math.pi) / math.pi
  return (almost_pi, rel_error)

# print(approx_pi(1000))
# print(approx_pi(100000))

# Questions 13 and 14
import matplotlib.patches as patches

def plot_monte_carlo() :
  fig = plt.figure()     # Create a figure
  fig.patch.set_facecolor('white') # Set its background to white
  axes = plt.axes()      # Create the axes of the figure
  axes.set_aspect('equal')   # Set their aspect ratio to 1
  axes.set_axis_off()    # Do not display the axes
  axes.set_xlim(left=-1, right=1.1)  # The .1 is for avoiding to clip the rectangle
  axes.set_ylim(bottom=-1.1, top=1)
  # Add the rectangle and the circle
  axes.add_patch(patches.Circle((0, 0), 1, linewidth=1, fill=False, edgecolor='blue'))
  axes.add_patch(patches.Rectangle((-1, -1), 2, 2, linewidth=1, fill=False))
  return (fig, axes)

def approx_pi_plot(n) :
  (fig, axes) = plot_monte_carlo()
  inside = 0
  total = 0
  while total < n :
    x = random.uniform(-1, 1)
    y = random.uniform(-1, 1)
    total += 1
    if math.hypot(x, y) < 1 :
      inside += 1
      plt.plot(x, y, 'ro', markersize=2)
    else:
      plt.plot(x, y, 'go', markersize=2)
    plt.pause(0.01)  # Pause for animation
  almost_pi = 4 * (inside / total)
  rel_error = math.fabs(almost_pi - math.pi) / math.pi
  axes.set_title("n = {0}  π = {1:f}  error = {2:.2%}".format(n, almost_pi, rel_error))
  plt.pause(0.01)  # To force the update of the title
  plt.show()
  return (almost_pi, rel_error)

# print(approx_pi_plot(500))

##
# Mandelbrot set
##

# Question 15
def module_real_img(z) :
  return math.hypot(z.real, z.imag)

def module_conj(z) :
  return math.sqrt((z.conjugate()*z).real)

# Question 16
def mandelfunc(c, z) :
  return z*z + c

def mandelbrot(c, n_max) :
  z = 0
  n = 0
  while n < n_max and abs(z) < 2 :
    n += 1
    z = mandelfunc(c, z)
  return n

# Question 17
def mandelbrot_set(xmin, xmax, ymin, ymax, xnum, ynum, n_max) :
  x = [xmin + i * (xmax - xmin)/xnum for i in range(xnum+1)]
  y = [ymin + i * (ymax - ymin)/ynum for i in range(ynum+1)]
  z = [[mandelbrot(complex(x[i], y[j]), n_max) for i in range(xnum)] for j in range(ynum)]
  return z

# Version numpy (much slower than the raw Python version)
import numpy as np

def mandelbrot_set_np(xmin, xmax, ymin, ymax, xnum, ynum, n_max) :
  x = np.linspace(xmin, xmax, xnum)
  y = np.linspace(ymin, ymax, ynum)
  z = np.empty((ynum, xnum))
  for i in range(xnum) :
    for k in range(ynum) :
      z[k,i] = mandelbrot(x[i] + 1j * y[k], n_max)
  return z


# Question 18
import matplotlib.colors as colors

def plot_mandelbrot(xmin, xmax, ymin, ymax, xnum, ynum, n_max) :
  z = mandelbrot_set(xmin, xmax, ymin, ymax, xnum, ynum, n_max)
  mycolors = [(1,0,0), (0,1,0), (0,0,1), (0,0,0)]
  colmap = colors.LinearSegmentedColormap.from_list("mandelbrot", mycolors, 16)
  plt.imshow(z, cmap=colmap, interpolation='bicubic')
  plt.show()

# plot_mandelbrot(-2.25, 0.75, -1.25, 1.25, 1500, 1250, 64)

# Ranking the web

# Question 19
def subdiv_one(n):
  elems = [random.random() for _  in range(n)]
  norm = sum(elems)
  return [elem / norm for elem in elems]

# Question 20
def create_transition_matrix(size) :
  # First, we create a matrix whose lines sum up to 1
  m = np.matrix([subdiv_one(size) for i in range(size)])
  # and we transpose it to have columns that sum up to 1
  return m.transpose()

# Questions 21 and 22
# https://en.wikipedia.org/wiki/Matrix_norm
import numpy.linalg as la

def page_rank(size, epsilon) :
  m = create_transition_matrix(size)
  s0 = [[1/size] for _ in range(size)]  # 1 column matrix
  rank = 0
  norm = la.norm((m ** (rank + 1) - m ** rank), np.inf)
  while (norm/size > epsilon) :
    rank += 1
    norm = la.norm((m ** (rank + 1) - m ** rank), np.inf)
  return (m ** rank) * s0

# More efficient version, which does not recompute m**k at
# each iteration. It also compare the norm directly to epsilon,
# without dividing by the number of pages
def page_rank_lowcost(size, epsilon) :
  m = create_transition_matrix(size)
  mk = m
  s0 = [[1/size] for _ in range(size)]
  rank = 0
  norm = la.norm(m, np.inf)
  while (norm > epsilon) :
    rank += 1
    mk2 = mk * m
    norm = la.norm((mk2 - mk), np.inf)
    mk = mk2
  return mk * s0

# Prey-predator

# Question 24
# The function that gives the derivative from the current values
def f(x, y, a, b) :
  return x * (a - b * y)

# Question 25
# The function that gives the next value of z (which is either x or y)
# from the current value of z, of the other variable (t), the parameters
# of the problem (a and b) and the integration step (Euler explicit scheme).
def next_step(zn, tn, a, b, step) :
  return zn + step * f(zn, tn, a, b)

# Question 26
def populations(x0, y0, alpha, beta, delta, gamma, step, n) :
  x, y = x0, y0
  lx, ly = [x0], [y0]
  for _ in range(n) :
    # We have dx/dt =  alpha*x - beta*x*y
    # and     dy/dt = -gamma*y + delta*y*x
    # hence the parameters of next_step
    x, y = next_step(x, y, alpha, beta, step), \
           next_step(y, x, -gamma, -delta, step)
    lx.append(x)
    ly.append(y)
  return lx, ly

# Question 27
def display_populations() :
  # Initial populations
  x0, y0 = 0.9, 0.9
  # Parameters
  alpha, beta, delta, gamma = 2/3, 4/3, 1, 1
  # Temporal resolution
  resolution = 1000
  # Compute until that time
  time_max = 20
  # Time step size
  step = 1 / resolution
  # Total number of steps
  n = time_max * resolution
  ## Compute the populations
  prey, pred = populations(x0, y0, alpha, beta, delta, gamma, step, n)
  # Plot the result
  time = [step * i for i in range(n+1)]
  plt.subplot(121)
  # First subplot for the evolution of populations through time
  plt.title("Population versus time")
  plt.plot(time, prey, label="Preys")
  plt.plot(time, pred, label="Predators")
  plt.xlabel("Time")
  plt.ylabel("Population")
  plt.legend(loc='upper left')
  # Second subplot for the phase diagram
  plt.subplot(122)
  plt.title("Predators versus preys")
  plt.plot(prey, pred)
  plt.xlabel("Preys")
  plt.ylabel("Predators")
  plt.show()

# display_populations()