1A.e - Enoncé 12 décembre 2017 (1)

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 \(f\) en escalier à partir d’un ensemble de points \((X_i, f(X_i))\).

Q1 - échantillon aléatoire

Générer un ensemble aléatoire de 1000 nombres \((X_i,Y_i)\) qui vérifie :

• \(X_i\) suit une loi uniforme sur \([0,16]\)

• \(Y_i = \mathbb{1}_{[\sqrt{X_i}] \mod 2}\) où \([A]\) est la partie entière de \(A\).

On pourra se servir de la fonction random du module random.

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.

%matplotlib inline

import matplotlib.pyplot as plt

plt.plot(X, Y, ".")

[<matplotlib.lines.Line2D at 0x1f893f9e518>]

Q2 - tri

Trier les points selon les \(X\).

nuage = [(x, y) for x, y in zip(X, Y)]
nuage.sort()
nuage[:5]

[(0.014962888038782651, 0),
 (0.020462778257442693, 0),
 (0.022310859639508962, 0),
 (0.03078728731371605, 0),
 (0.03153252863972433, 0)]

Q3 - moyenne

On suppose que les \(Y\) sont triés selon les \(X\) croissants. Calculer la moyenne des différences entre \(Y\) et la moyenne \(m\) des \(Y\) (en valeur absolue) sur un intervalle \([i,j]\), \(j\) exclu. Ecrire une fonction def somme_diff(nuage, i, j) qui exécute ce calcul qui correspond à \(\sum_{k=i}^{j-1} |Y_k - m|\) avec \(m = (\sum_{k=i}^{j-1} Y_k) / (j-i)\).

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

(0.0, 476.2380000000092)

Q4 - distance

Soit \(i,j\) deux entiers, on coupe l’intervalle en deux : \(i,k\) et \(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):.

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)

18.56022222223197

Q5 - fonction comme paramètre

Le langage Python permet de passer une fonction à une autre fonction en tant qu’argument. Un exemple :

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)

1

Ecrire la fonction précédente en utilisant la fonction fct.

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)

494.7982222222412

Q6 - optimiser

On veut déterminer le \(i\) optimal, celui qui maximise la différence dans l’intervalle \([i,j]\). On souhaite garder la fonction fct comme argument. Pour cela, implémenter la fonction def optimise(nuage, i, j, f):.

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)

(565, 711.6476814159435)

import matplotlib.pyplot as plt

x = nuage[552][0]
plt.plot(X, Y, ".")
plt.plot([x, x], [0, 1])

[<matplotlib.lines.Line2D at 0x1f893ea57b8>]

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.

optimise(nuage, 0, 68, fct), optimise(nuage, 68, len(nuage), fct)

((1, 0.0), (565, 618.0710615624871))

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

[<matplotlib.lines.Line2D at 0x1f8943ece80>]

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):.

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

[68, 242, 565]

import matplotlib.pyplot as plt

plt.plot(X, Y, ".")
for i in r:
    x = nuage[i][0]
    plt.plot([x, x], [0, 1])

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 \(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 \(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 \(O(n^2)\). Tout d’abord la version actuelle.

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

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.

# %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 \(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 \(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 \(k\) varie, il suffit de mettre à jour les histogrammes en déduisant et en ajoutant le \(Y_k\) aux bons endroits.

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

((565, 235.4096814159292), (565, 235.40968141593424))

C’est carrément plus rapide et cela marche pour toute fonction fct.

%timeit optimise_rapide(nuage, 0, len(nuage))

1.63 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Si on ne suppose pas que les \(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 \(Y_i\), de part et d’autres de la coupure. Lorsqu’on bouge la coupure \(k\), cela revient à déplacer \(Y_k\) d’un tableau à l’autre ce qui se fera par recherche dichotomique donc en \(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 ?

import random

X2 = list(range(10))
Y2 = X2

import matplotlib.pyplot as plt

plt.plot(X2, Y2, ".")

[<matplotlib.lines.Line2D at 0x1f8944d3ba8>]

nuage2 = [(x, y) for x, y in zip(X2, Y2)]
nuage2.sort()

r = recursive(nuage2, 0, len(nuage2), fct)
len(r), r

(5, [2, 3, 5, 7, 8])

import matplotlib.pyplot as plt

plt.plot(X2, Y2, ".")
for i in r:
    x = nuage2[i][0]
    plt.plot([x, x], [0, 10])

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.

---

Notebook on github