.. _tdnote20182rst: ================================== 1A.e - Enoncé 12 décembre 2017 (2) ================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/exams/td_note_2018_2.ipynb|*` Correction du premier énoncé de l’examen du 12 décembre 2017. Celui-ci mène à l’implémentation d’un algorithme qui permet d’approximer une fonction :math:`f` par une fonction en escalier à partir d’un ensemble de points :math:`(X_i, f(X_i))`. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Q1 - échantillon aléatoire -------------------------- Générer un ensemble aléatoire de 1000 nombres :math:`(X_i,Y_i)` qui vérifie : - :math:`X_i` suit une loi uniforme sur :math:`[0,16]` - :math:`Y_i = \sqrt{X_i}[\sqrt{X_i}]` où :math:`[A]` est la partie entière de :math:`A`. On pourra se servir de la fonction ``random`` du module ``random``. .. code:: ipython3 import random X = [random.random() * 16 for i in range(0,1000)] Y = [ x**0.5 * int(x**0.5) for x in X] Q1 - dessiner le nuage de points - donnée ----------------------------------------- Le code suivant vous est donné afin de vérifier vos réponses. .. code:: ipython3 %matplotlib inline .. code:: ipython3 import matplotlib.pyplot as plt plt.plot(X, Y, '.') .. parsed-literal:: [] .. image:: td_note_2018_2_6_1.png Q2 - tri -------- Trier les points selon les :math:`X`. .. code:: ipython3 nuage = [(x,y) for x,y in zip(X,Y)] nuage.sort() Q3 - moyenne ------------ On suppose que les :math:`Y` sont triés selon les :math:`X` croissants. Calculer la moyenne des différences entre :math:`Y` et la moyenne :math:`m` des :math:`Y` (au carré) sur un intervalle :math:`[i,j]`, :math:`j` exclu. Ecrire une fonction ``def somme_diff(nuage, i, j)`` qui exécute ce calcul qui correspond à :math:`\sum_{k=i}^j (Y_k - m)^2` avec :math:`m = (\sum_{k=i}^j Y_k) / (j-i)`. .. code:: ipython3 def somme_diff(xy, i, j): m = sum(e[1] for e in xy[i:j]) / (j-i) return sum((e[1]-m)**2 for e in xy[i:j]) somme_diff(nuage, 0, 5), somme_diff(nuage, 0, len(nuage)) .. parsed-literal:: (0.0, 15754.105018618644) Q4 - distance ------------- Soit :math:`i,j` deux entiers, on coupe l’intervalle en deux : :math:`i,k` et :math:`k,j`. On calcule ``somme_diff`` sur ces deux intervalles, on fait la somme des différences (en valeurs absolues) de ces moyennes par rapport à la valeur sur le plus grand intervalle. On écrit la fonction ``def difference(nuage, i, j, k):``. .. code:: ipython3 def difference(nuage, i, j, k): m1 = somme_diff(nuage, i, k) m2 = somme_diff(nuage, k, j) m = somme_diff(nuage, i, j) return abs(m-m1) + abs(m-m2) difference(nuage, 0, len(nuage), 100) .. parsed-literal:: 19898.600443365925 Q5 - fonction comme paramètre ----------------------------- Le langage Python permet de passer une fonction à une autre fonction en tant qu’argument. Un exemple : .. code:: ipython3 def fct(x, y): return abs(x-y) def distance_list(list_x, list_y, f): return sum(f(x,y) for x,y in zip(list_x, list_y)) distance_list([0, 1], [0, 2], fct) .. parsed-literal:: 1 Ecrire la fonction précédente en utilisant la fonction ``fct``. .. code:: ipython3 def somme_diff(xy, i, j, f): m = sum(e[1] for e in xy[i:j]) / (j-i) return sum(f(e[1], m) for e in xy[i:j]) def difference(nuage, i, j, k, f): m1 = somme_diff(nuage, i, k, f) m2 = somme_diff(nuage, k, j, f) m = somme_diff(nuage, i, j, f) return abs(m1 + m2 - m) difference(nuage, 0, len(nuage), 100, fct) .. parsed-literal:: 552.6383487080093 Q6 - optimiser -------------- On veut déterminer le :math:`i` optimal, celui qui maximise la différence dans l’intervalle :math:`[i,j]`. On souhaite garder la fonction ``fct`` comme argument. Pour cela, implémenter la fonction ``def optimise(nuage, i, j, f):``. .. code:: ipython3 def optimise(nuage, i, j, f): mx = -1 ib = None for k in range(i+1,j-1): d = difference(nuage, i,j,k, f) if ib is None or d > mx: mx = d ib = k return ib, mx optimise(nuage, 0, len(nuage), fct) .. parsed-literal:: (553, 2184.8079894060775) .. code:: ipython3 import matplotlib.pyplot as plt x = nuage[553][0] plt.plot(X,Y,'.') plt.plot([x,x], [0,10]) .. parsed-literal:: [] .. image:: td_note_2018_2_19_1.png Q7 - optimisation encore ------------------------ Recommencer sur les deux intervalles trouvés .. code:: ipython3 optimise(nuage, 0, 570, fct), optimise(nuage, 570, len(nuage), fct) .. parsed-literal:: ((253, 787.5154656398129), (789, 156.7106930739271)) .. code:: ipython3 import matplotlib.pyplot as plt x = nuage[253][0] x2 = nuage[789][0] plt.plot(X,Y,'.') plt.plot([x,x], [0,10]) plt.plot([x2,x2], [0,10]) .. parsed-literal:: [] .. image:: td_note_2018_2_22_1.png Q8 - fonction récursive ----------------------- Pouvez-vous imaginer une fonction récursive qui produit toutes les séparations. Ecrire la fonction ``def recursive(nuage, i, j, f, th=0.1):``. La fonction ``recursive`` tente de placer chaque point dans un intervalle distinct et elle échoue car elle dépasse le nombre d’appels récursifs autorisés par le langage Python. .. code:: ipython3 def recursive(nuage, i, j, f): if i+1 == j or i == j: return [i] k, mx = optimise(nuage, i, j, f) if k is None: return None else: r1 = recursive(nuage, i, k, f) r2 = recursive(nuage, k, j, f) if r1 is None and r2 is None: return [k] elif r1 is None: return [k] + r2 elif r2 is None: return r1 + [k] else: return r1 + [k] + r2 # déclenche une exception r = recursive(nuage, 0, len(nuage), fct) len(r), r[:5], r[-5:] .. parsed-literal:: (946, [0, 1, 1, 2, 2], [996, 996, 997, 997, 998]) La fonction n’est pas encore parfait et elle retourne plusieurs fois le même point de coupure mais ce n’est pas le plus grave car la fonction retourne quasiment tous les points. Q9 - seuil ---------- L’algorithme produit beaucoup de points de coupures. On souhaite arrêter la récursion plus tôt en mettant un seuil sur la quantité obtenue :math:`|\Delta_{i}^k + \Delta_{k}^j - \Delta_{i}^j|` qui doit être supérieur à 50. .. code:: ipython3 def recursive(nuage, i, j, f, th=50): k, mx = optimise(nuage, i, j, f) if mx <= th: return None r1 = recursive(nuage, i, k, f, th=th) r2 = recursive(nuage, k, j, f, th=th) if r1 is None and r2 is None: return [k] elif r1 is None: return [k] + r2 elif r2 is None: return r1 + [k] else: return r1 + [k] + r2 r = recursive(nuage, 0, len(nuage), fct) len(r), r[:5] .. parsed-literal:: (5, [61, 253, 384, 553, 782]) .. code:: ipython3 import matplotlib.pyplot as plt plt.plot(X, Y, '.') for i in r: x = nuage[i][0] plt.plot([x,x], [0,10]) .. image:: td_note_2018_2_28_0.png Q10 - coût ---------- Quel est le coût de la fonction ``optimize`` en fonction de la taille de l’intervalle ? Peut-on mieux faire (ce qu’on n’implémentera pas). Tel qu’il est implémenté, le coût est en :math:`O(n^2)`, le coût peut être linéaire en triant les éléments dans l’ordre croissant, ce qui a été fait, ou :math:`n\ln n` si on inclut le coût du tri bien qu’on ne le fasse qu’une fois. Voyons plus en détail comment se débarrasser du coût en :math:`O(n^2)`. Tout d’abord la version actuelle. .. code:: ipython3 def somme_diff2(xy, i, j): m = sum(e[1] for e in xy[i:j]) / (j-i) return sum((e[1]-m)**2 for e in xy[i:j]) def difference2(nuage, i, j, k): m1 = somme_diff2(nuage, i, k) m2 = somme_diff2(nuage, k, j) m = somme_diff2(nuage, i, j) return abs(m1+m2-m) def optimise2(nuage, i, j): mx = -1 ib = None for k in range(i+1,j-1): d = difference2(nuage, i,j,k) if ib is None or d > mx: mx = d ib = k if ib is None: ib = i mx = 0 return ib, mx %timeit optimise2(nuage, 0, len(nuage)) .. parsed-literal:: 628 ms ± 6.15 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) L’instruction suivante permet de voir où le programme passe la majeure partie de son temps. .. code:: ipython3 # %prun optimise_abs(nuage, 0, len(nuage)) La fonction `sum `__ cache une boucle, avec la boucle ``for`` dans la fonction ``optimise``, cela explique le coût en :math:`O(n^2)`. Le fait qu’à chaque itération, on passe une observation d’un côté à l’autre de la coupure puis on recalcule les moyennes… Nous allons optimiser ce calcul en tenant compte du fait que la fonction de coût est :math:`f(x,y) = (x-y)^2`. Il faut également se souvenir de la formule :math:`\mathbb{V}X = \mathbb{E}(X^2) - (\mathbb{E}X)^2` ce qu’on transforme en :math:`\sum_{i=1}^n (X_i-M)^2 = \sum_{i=1}^n X_i^2 - n M^2` avec :math:`M = \frac{1}{n} \sum_{i=1}^n X_i`. .. code:: ipython3 def optimise_rapide(nuage, i, j): # On calcule les histogrammes. Y = sum(y for x,y in nuage) Y_1 = sum(y for x,y in nuage[i:i+1]) Y_2 = sum(y for x,y in nuage[i+1:j]) Y2 = sum(y**2 for x,y in nuage) Y2_1 = sum(y**2 for x,y in nuage[i:i+1]) Y2_2 = sum(y**2 for x,y in nuage[i+1:j]) m = Y2 - Y**2 / len(nuage) m1 = Y2_1 - Y_1**2 m2 = Y2_2 - Y_2**2/(len(nuage)-1) mx = -1 ib = None for k in range(i+1,j-1): d = abs(m1+m2-m) if ib is None or d > mx: mx = d ib = k # On met à jour les sommes Y?_? y = nuage[k][1] Y_1 += y Y_2 -= y Y2_1 += y**2 Y2_2 -= y**2 m1 = Y2_1 - Y_1**2 / (k-i+1) m2 = Y2_2 - Y_2**2 / (j-k-1) if ib is None: ib = i mx = 0 return ib, mx # On vérifie qu'on obtient les mêmes résultats. optimise_rapide(nuage, 0, len(nuage)), optimise2(nuage, 0, len(nuage)) .. parsed-literal:: ((553, 13082.574312018376), (553, 13082.574312018447)) .. code:: ipython3 %timeit optimise_rapide(nuage, 0, len(nuage)) .. parsed-literal:: 1.77 ms ± 43.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) Carrément plus rapide.