tpmni/3/TP3_voisin.py

#!/usr/bin/env python3
import sys, os
import matplotlib.pyplot as plt
import numpy as np
from time import sleep
import argparse
from math import * # pour utliser les fonctions mathématiques dans l'option -f

"""
    MNI TP3 - 27 janvier 2020
    Dylan Voisin - M1 Informatique Luminy
    Franck Palacios - M1 Informatique Luminy
    Cyril Colin - M1 Informatique Luminy

    Notes:
        - u1 et u2 ne semblent pas orthogonaux, ceci est dû à la non utilisation
          d'une échelle 1:1 pour l'affichage des points
        - u1 et u2 ne sont pas proportionnels par rapport à lambda mais à la
          distance entre le barycentre et la limite de la fenêtre
        - u1 est orienté «vers la droite» et u2 «vers le haut» dans ce programme
        - u1 et u2 apparaîssent comment des traits plutôt que des flèches car
          matplotlib gère mal l'affichage de vecteurs, et plt.arrow peut dans
          certains cas donner des têtes de flèche assez étranges selon le rapport
          y:x du nuage de point
        - le code n'a pas besoin d'être touché pour changer le comportement lors
          du lancement, exemples:
            
            * afficher un nuage de point selon cos(x) entre 0 et 1 (pas de 0.1 par
              défaut et inchangeable sauf dans le code):
              > python3 TP3_voisin.py -f "cos(x)" -m 1 -s

            * idem mais avec un bruit de 0.5 (-g applique le bruit -n BRUIT dès la génération)
              > python3 TP3_voisin.py -f "cos(x)" -m 1 -g -n 0.5 -s
            
            * afficher 10 boucles (avec bruit, défaut=10) avec la fonction sqrt(x)
              > python3 TP3_voisin.py -f "sqrt(x)" -m 10 -n 2 -l 10
              le délai entre chaque itération se fait avec le paramètre -d SECONDES (ex -d 2.5)
            
            * afficher les points contenus dans le fichier exemple.txt (format plus bas)
              > python3 TP3_voisin.py -i exemple.txt -s

            * les boucles, génération de bruit à l'initialisation sont également possible pour
              la lecture d'un fichier
            
            * on peut afficher l'aide via
              > python3 TP3_voisin.py -h
        
        Format du fichier d'entrée:
        x0 y0
        x1 y1
        …
        xn yn
"""

def calc_barycentre(l_pts):
    x = sum(l_pts[:,0]) / len(l_pts)
    y = sum(l_pts[:,1]) / len(l_pts)
    return x, y

def centrer_points(l_pts, barycentre):
    return l_pts - barycentre

def calc_matrice_corrélation(l_pts):
    c_pts = centrer_points(l_pts, calc_barycentre(l_pts))
    A = np.array([
        [sum(pt[0] ** 2 for pt in c_pts), sum(pt[0] * pt[1] for pt in c_pts)],
        [sum(pt[0] * pt[1] for pt in c_pts), sum(pt[1] ** 2 for pt in c_pts)]
    ])
    return A

def diagonaliser_matrice22(m):
    a, b, c, d = m[0][0], m[0][1], m[1][0], m[1][1]
    discr = (-d - a) ** 2 - 4 * (a * d - b**2)
    lambda1 = ((d + a) + discr ** .5) / 2
    lambda2 = ((d + a) - discr ** .5) / 2
    if lambda1 != 0:
        if a - lambda1 != 0 and b != 0:
            u1 = np.array([-b, (a-lambda1)])
            u1 *= -1 if b > 0 else 1
        else:
            u1 = np.array([lambda1, 0])
        u2 = np.array([-u1[1], u1[0]])
    else:
        if a - lambda2 != 0 and b != 0:
            u2 = np.array([-b, (a-lambda2)])
            u2 *= -1 if b < 0 else 1
        else:
            u2 = np.array([lambda2, 0])
        u1 = np.array([u2[1], -u2[0]])
    return u1, u2, lambda1, lambda2

def get_u1_u2_affichage(u1, u2, bary, dist_x, dist_y, pb, pg):
    """renvoie u1 et u2 prêt pour l'affichage en fonction du barycentre, 
       de la distance x et y sur l'ensemble des points, le point le plus bas et
       le point le plus à gauche"""
    if np.linalg.norm(u1) == np.linalg.norm(u2) == 0:
        # un seul point
        u1 = np.array([1, 0]) / 500
        u2 = np.array([0, 1]) / 500
        return u1, u2
    bary = np.array((bary[0] - pg[0], bary[1] - pb[1]))
    if u1[0] != 0 and u1[1] != 0:
        if u1[0] > 0: # vers la droite
            t1x = (dist_x - bary[0]) / u1[0]
        else:
            t1x = -bary[0] / u1[0]
        if u1[1] > 0: # vers le haut
            t1y = (dist_y - bary[1]) / u1[1]
        else:
            t1y = -bary[1] / u1[1]
        t1 = min(t1x, t1y)
    elif u1[0] != 0:
        if u1[0] > 0: # vers la droite
            t1 = (dist_x - bary[0]) / u1[0]
        else:
            t1 = -bary[0] / u1[0]
    else:
        if u1[1] > 0: # vers le haut
            t1 = (dist_y - bary[1]) / u1[1]
        else:
            t1 = -bary[1] / u1[1]
    if u2[0] != 0 and u2[1] != 0:
        if u2[0] > 0: # vers la droite
            t2x = (dist_x - bary[0]) / u2[0]
        else:
            t2x = -bary[0] / u2[0]
        if u2[1] > 0: # vers le haut
            t2y = (dist_y - bary[1]) / u2[1]
        else:
            t2y = -bary[1] / u2[1]
        t2 = min(t2x, t2y)
    elif u2[0] != 0:
        if u2[0] > 0: # vers la droite
            t2 = (dist_x - bary[0]) / u2[0]
        else:
            t2 = -bary[0] / u2[0]
    else:
        if u2[1] > 0: # vers le haut
            t2 = (dist_y - bary[1]) / u2[1]
        else:
            t2 = -bary[1] / u2[1]
    u1 = u1 * t1 * 2/3
    u2 = u2 * t2 * 1/3
    return u1, u2

def get_info_affichage(l_pts):
    """renvoie les vecteurs u1 et u2 pour l'affichage ainsi que xmin, xmax, ymin, ymax"""
    bary = calc_barycentre(l_pts)
    u1, u2 = diagonaliser_matrice22(calc_matrice_corrélation(l_pts))[:2]
    pb = min(l_pts, key = lambda x: x[1])
    pg = min(l_pts, key = lambda x: x[0])
    dist = min(max(l_pts[:,0]) - bary[0], max(l_pts[:,1]) - bary[1])
    dist_x = max(l_pts[:,0]) - min(l_pts[:,0])
    dist_y = max(l_pts[:,1]) - min(l_pts[:,1])
    x_range = min(l_pts[:,0]) - 0.1 * dist_x, max(l_pts[:,0]) + 0.1 * dist_x
    y_range = min(l_pts[:,1]) - 0.1 * dist_y, max(l_pts[:,1]) + 0.1 * dist_y
    u1, u2 = get_u1_u2_affichage(u1, u2, bary, dist_x, dist_y, pb, pg)
    return u1, u2, x_range, y_range


def visualiser_nuage(l_pts):
    bary = calc_barycentre(l_pts)
    u1, u2 = get_info_affichage(l_pts)[:2]
    ax = plt.axes()
    # plot du nuage
    plt.scatter(l_pts[:,0], l_pts[:,1])
    ax.arrow(*bary, *u1, ec="red", fc="white")
    ax.arrow(*bary, *u2, ec="orange", fc="white")
    plt.show()

def animation_nuage(l_pts, loops=10, time_sleep=2, delta_bruit=3):
    plt.ion()
    fig = plt.figure()
    ax = fig.add_subplot(111)
    fig.canvas.draw()
    # nuage courant (pas de do while en Python)
    bary = calc_barycentre(l_pts)
    u1, u2 = get_info_affichage(l_pts)[:2]
    u1l = np.linalg.norm(u1)
    sc = ax.scatter(l_pts[:,0], l_pts[:,1])
    arr1 = ax.arrow(*bary, *u1, ec="red", fc="white")
    arr2 = ax.arrow(*bary, *u2, ec="orange", fc="white")
    plt.title("itération 0")
    fig.canvas.flush_events()
    fig.canvas.draw()
    fig.canvas.flush_events()
    sleep(time_sleep) # attente avant animation suivante
    for i in range(loops):
        plt.title("itération " + str(i+1))
        # aléa
        ran = np.random.random(l_pts.shape) * 2 * delta_bruit - delta_bruit
        l_pts += ran
        # recalcul des données
        bary = calc_barycentre(l_pts)
        u1, u2, x_range, y_range = get_info_affichage(l_pts)
        u1l = np.linalg.norm(u1)
        ax.set_xlim(*x_range)
        ax.set_ylim(*y_range)
        # plot
        sc.set_offsets(l_pts) # repositionnement des points
        arr1.remove() # on enlève les vecteurs
        arr2.remove()
        arr1 = ax.arrow(*bary, *u1, ec="red", fc="white")
        arr2 = ax.arrow(*bary, *u2, ec="orange", fc="white")
        fig.canvas.flush_events()
        fig.canvas.draw() # on redessine
        fig.canvas.flush_events()
        sleep(time_sleep) # attente avant animation suivante


def main(argv):
    parser = argparse.ArgumentParser()
    gps = parser.add_mutually_exclusive_group(required=True)
    gpi = parser.add_mutually_exclusive_group(required=True)
    gpi.add_argument("-i", "--input",
                        help="input file for dots - if none, generate dots")
    gpi.add_argument("-f", "--function",
                        help="function to use to generate dots, ex -f \"3*x**2+cos(x)\"")
    parser.add_argument("-m", "--max-dot-xval", type=float, default=10,
                        help="x range of generated dots (from 0 to MAX_DOT_VAL, step=0.1) from function")
    gps.add_argument("-s", "--show", action="store_true",
                        help="show cloud of dots with u1 and u2")
    gps.add_argument("-l", "--loop", type=int,
                    help="number of loop to show/perform on randomized dots from originals ones")
    parser.add_argument("-g", "--generate-with-noise", action="store_true",
                        help="apply noise on dots generation")
    parser.add_argument("-d", "--delay", type=float, default=2,
                    help="amount of time between each loop")
    parser.add_argument("-n", "--noise", type = float, default=10,
                    help="maximum value of the noise applied between each step and on generation if option -g and -f called")
    if not argv:
        parser.print_usage()
        sys.exit(0)
    args = parser.parse_args(argv)
    if args.input:
        if not os.path.isfile(args.input):
            print("File", args.input, "not found.", file=sys.stderr)
            sys.exit(1)
        else:
            l_pts = np.loadtxt(args.input)
    else:
        n = args.max_dot_xval
        f = lambda x: eval(args.function) # fonction génératrice de la ditribution
        l_pts = np.array([[i, f(i)] for i in np.arange(0, n, 0.1)])
    if args.generate_with_noise:
            delta_bruit = args.noise # bruit généré, 0 si on en veut aucun
            bruit = np.random.random(len(l_pts)) * 2 * delta_bruit - delta_bruit            
            l_pts[:,1] += bruit # application du bruit
    if args.show:
        visualiser_nuage(l_pts)
    elif args.loop:
        animation_nuage(l_pts, args.loop, args.delay, args.noise)


if __name__ == "__main__":
    main(sys.argv[1:])