tpmni/3/TP3_voisin.py

#!/usr/bin/env python3
import sys, os
import matplotlib.pyplot as plt
from random import random, choice
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

    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 le point le plus éloigné selon la
          composante x ou y
        - u1 est orienté «vers la droite» et u2 «vers le haut» dans ce programme
"""


def acquérir_depuis_fichier(filename):
    # l_pts = []
    # with open(filename) as f:
    #     for line in f:
    #         split = line.split()
    #         if len(split) == 2:
    #             l_pts.append(list(map(float, split)))
    # return np.array(l_pts)
    return np.loadtxt(filename)

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 = max(map(lambda x: np.linalg.norm(x-bary), l_pts))
    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)
    # print(u1, u2)
    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")
    # ax.set_aspect('equal')
    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 = acquérir_depuis_fichier(args.input)
    else:
        l_pts = []
        n = args.max_dot_xval
        f = lambda x: eval(args.function) # fonction génératrice de la ditribution
        if args.generate_with_noise:
            delta_bruit = args.noise # bruit généré, 0 si on en veut aucun
        else:
            delta_bruit = 0
        for i in np.arange(0, n, 0.1):
            y = f(i) + choice((-1, 1)) * random() * delta_bruit/2
            l_pts.append([i, y])
        l_pts = np.array(l_pts)
    if args.show:
        visualiser_nuage(l_pts)
    elif args.loop:
        animation_nuage(l_pts, args.loop, args.delay, args.noise)

    

if __name__ == "__main__":
    # if len(sys.argv) == 2:
    #     l_pts = acquérir_depuis_fichier(sys.argv[1])
    # else:
    #     l_pts = []
    #     n = 500 # nombre de points
    #     f = lambda x: x/2 + 3 # fonction génératrice de la ditribution
    #     delta_bruit = 10 # bruit généré, 0 si on en veut aucun
    #     for i in range(n):
    #         y = f(i) + choice((-1, 1)) * random() * delta_bruit/2
    #         l_pts.append([i, y])
    #     l_pts = np.array(l_pts)
    # # np.savetxt("pts.txt", l_pts, "%f")
    # u1, u2 = diagonaliser_matrice22(calc_matrice_corrélation(l_pts))[:2]
    # # print(f"{u1/np.linalg.norm(u1)=}\n{u2/np.linalg.norm(u2)=}")
    # visualiser_nuage(l_pts)
    # animation_nuage(l_pts, 20, delta_bruit=100)
    main(sys.argv[1:])