.. _tdnote20181rst: ================================== 1A.e - Enoncé 12 décembre 2017 (1) ================================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/exams/td_note_2018_1.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 de retrouver une fonction :math:`f` 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 = \mathbf{1\!\!1}_{[\sqrt{X_i}] \mod 2}` 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 = [ int(x**0.5) % 2 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_1_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() nuage[:5] .. parsed-literal:: [(0.014962888038782651, 0), (0.020462778257442693, 0), (0.022310859639508962, 0), (0.03078728731371605, 0), (0.03153252863972433, 0)] 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` (en valeur absolue) 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-1} |Y_k - m|` avec :math:`m = (\sum_{k=i}^{j-1} 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(abs(e[1]-m) for e in xy[i:j]) somme_diff(nuage, 0, 5), somme_diff(nuage, 0, len(nuage)) .. parsed-literal:: (0.0, 476.2380000000092) Q4 - distance ------------- Soit :math:`i,j` deux entiers, on coupe l’intervalle en deux : :math:`i,k` et :math:`k,j`. On calcule la ``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(m1+m2-m) difference(nuage, 0, len(nuage), 100) .. parsed-literal:: 18.56022222223197 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) # On a modifié les fonctions précédentes pour calculer # une fonction d'erreur "custom" ou définie par l'utilisateur. 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(m - m1) + abs(m - m2) difference(nuage, 0, len(nuage), 100, fct) .. parsed-literal:: 494.7982222222412 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 if ib is None: # Au cas où l'intervalle est vide, on retourne une coupure # égale à i. ib = i mx = 0 return ib, mx optimise(nuage, 0, len(nuage), fct) .. parsed-literal:: (565, 711.6476814159435) .. code:: ipython3 import matplotlib.pyplot as plt x = nuage[552][0] plt.plot(X,Y,'.') plt.plot([x,x], [0,1]) .. parsed-literal:: [] .. image:: td_note_2018_1_19_1.png Le premier point de coupure trouvé (le trait orange) correspond à un des bords d’un des escaliers. Q7 - optimisation encore ------------------------ Recommencer sur les deux intervalles trouvés. La question était juste histoire que le résultat à la question précédente est reproductible sur d’autres intervalles. .. code:: ipython3 optimise(nuage, 0, 68, fct), optimise(nuage, 68, len(nuage), fct) .. parsed-literal:: ((1, 0.0), (565, 618.0710615624871)) .. code:: ipython3 import matplotlib.pyplot as plt x = nuage[58][0] x2 = nuage[552][0] plt.plot(X,Y,'.') plt.plot([x,x], [0,1]) plt.plot([x2,x2], [0,1]) .. parsed-literal:: [] .. image:: td_note_2018_1_23_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):``. .. code:: ipython3 def recursive(nuage, i, j, f, th=0.1): 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) r .. parsed-literal:: [68, 242, 565] .. code:: ipython3 import matplotlib.pyplot as plt plt.plot(X, Y, '.') for i in r: x = nuage[i][0] plt.plot([x,x], [0,1]) .. image:: td_note_2018_1_26_0.png Q9 - 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_diff_abs(xy, i, j): m = sum(e[1] for e in xy[i:j]) / (j-i) return sum(abs(e[1]-m) for e in xy[i:j]) def difference_abs(nuage, i, j, k): m1 = somme_diff_abs(nuage, i, k) m2 = somme_diff_abs(nuage, k, j) m = somme_diff_abs(nuage, i, j) return abs(m1+m2-m) def optimise_abs(nuage, i, j): mx = -1 ib = None for k in range(i+1,j-1): d = difference_abs(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 optimise_abs(nuage, 0, len(nuage)) .. parsed-literal:: 503 ms ± 21.1 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… Il y a deux façons d’optimiser ce calcul selon qu’on tient compte du fait que les valeurs de :math:`Y` sont binaires ou non. Dans le premier cas, il suffit de compter les valeurs 0 ou 1 de part et d’autres de la coupure (histogrammes) pour calculer la moyenne. Lorsque :math:`k` varie, il suffit de mettre à jour les histogrammes en déduisant et en ajoutant le :math:`Y_k` aux bons endroits. .. code:: ipython3 def histogramme_y(xy, i, j): d = [0, 0] for x, y in xy[i:j]: d[y] += 1 return d def somme_diff_histogramme(d): m = d[1] * 1.0 / (d[0] + d[1]) return (1-m) * d[1] + m * d[0] def optimise_rapide(nuage, i, j): # On calcule les histogrammes. d1 = histogramme_y(nuage, i, i+1) d2 = histogramme_y(nuage, i+1, j) d = d1.copy() d[0] += d2[0] d[1] += d2[1] m = somme_diff_histogramme(d) m1 = somme_diff_histogramme(d1) m2 = somme_diff_histogramme(d2) 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 histogrammes. On ajoute d'un côté, on retranche de l'autre. y = nuage[k][1] d1[y] += 1 d2[y] -= 1 m1 = somme_diff_histogramme(d1) m2 = somme_diff_histogramme(d2) 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)), optimise_abs(nuage, 0, len(nuage)) .. parsed-literal:: ((565, 235.4096814159292), (565, 235.40968141593424)) C’est carrément plus rapide et cela marche pour toute fonction ``fct``. .. code:: ipython3 %timeit optimise_rapide(nuage, 0, len(nuage)) .. parsed-literal:: 1.63 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) Si on ne suppose pas que les :math:`Y_i` sont binaires et qu’ils sont quelconques, les histogrammes contiendront plus de deux éléments. Dans ce cas, il faut conserver deux tableaux triés des :math:`Y_i`, de part et d’autres de la coupure. Lorsqu’on bouge la coupure :math:`k`, cela revient à déplacer :math:`Y_k` d’un tableau à l’autre ce qui se fera par recherche dichotomique donc en :math:`O(\ln n)`. La mise à jour de la moyenne des valeurs absolues est immédiate si la fonction ``fct=abs(x-y)`` et pas forcément immédiate dans le cas général. Lorsque c’est une valeur absolue, il faut utiliser quelques résultats sur la `régression quantile `__. Q10 - autre nuage de points --------------------------- Comment l’algorithme se comporte-t-il lorsque tous les points sont distincts ? .. code:: ipython3 import random X2 = list(range(10)) Y2 = X2 .. code:: ipython3 import matplotlib.pyplot as plt plt.plot(X2,Y2,'.') .. parsed-literal:: [] .. image:: td_note_2018_1_39_1.png .. code:: ipython3 nuage2 = [(x,y) for x,y in zip(X2,Y2)] nuage2.sort() .. code:: ipython3 r = recursive(nuage2, 0, len(nuage2), fct) len(r), r .. parsed-literal:: (5, [2, 3, 5, 7, 8]) .. code:: ipython3 import matplotlib.pyplot as plt plt.plot(X2,Y2,'.') for i in r: x = nuage2[i][0] plt.plot([x,x], [0,10]) .. image:: td_note_2018_1_42_0.png La fonction va placer un point dans chaque intervalle, il y aura quasiment autant de points de coupures que de points. Presque autant car l’implémentation a quelques effets de bords comme la boucle ``for k in range(i+1,j-1):`` qui ne considère pas les extrémités d’un intervalle comme de potentiels points de coupures.