Nuage de points python

Gentil · January 2020

Bonjour
J'ai des bases en python mais pas suffisantes pour faire ce que je souhaite faire, ce serait super si vous pouviez me donner un coup de main.

1) J'aimerais afficher en rouge les points (i,j) par exemple (0.2,0.6) et en bleu des points (k,l) par exemple (0.1,0.3) et les lier par créer une fonction qui relie le point (i,j) au point (k,l) par un trait.
2) De plus j'aimerais afficher les points bleus et rouges sur le même graphique, pour qu'on y voit quelque chose j'aimerais séparer le nuage de points bleus et rouges. Je vous mets une image pour que vous voyiez ce que je veux dire.
Je vous remercie sincèrement pour toute aide. 95394

nicolas.patrois · January 2020

Ça devrait être faisable avec la bibliothèque matplotlib, non ?

Gentil · January 2020

Voici mes matrices (mu et nu) si vous voulez les voir.

#Paramétrisation des espaces de départ et d'arrivée
s1 = np.linspace(-5,5,N, False)
s2 = np.linspace(-5,5,N, False)
S1,S2 = np.meshgrid(s1,s2)          #points de discrétisation de l'espace de départ

t1 = np.linspace(-5,5,N, False)
t2 = np.linspace(-5,5,N, False)
T1,T2 = np.meshgrid(t1,t2)          #points de discrétisation de l'espace de d'arrivé

#Définition des mesures
def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return epsilon
    else:
        return np.exp(-(x1**2+x2**2)/2)/2*np.pi
mu = fmu(S1,S2)
mu = mu/np.sum(mu)

def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return epsilon
    else:
        return np.exp(-(x1**2+x2**2)/4)/4*np.pi
nu = fnu(T1,T2)
#print(np.sum(nu))
nu = nu/np.sum(nu)

Gentil · January 2020

Bonjour nicolas.patrois, comment feriez vous avec cette bibliothèque ?

Benoit RIVET · January 2020

Le code proposé est incomplet. Si on teste le code on obtient un message d'erreur :

[color=#FF0000]    s1 = np.linspace(-5,5,N, False)
NameError: name 'np' is not defined
[/color]

En rajoutant la ligne :

import numpy as np

on obtient un nouveau message d'erreur :

[color=#FF0000]    s1 = np.linspace(-5,5,N, False)
NameError: name 'N' is not defined
[/color]

En rajoutant la ligne :

N = 10

on obtient un nouveau message d'erreur :

[color=#FF0000]    mu = fmu(S1,S2)
NameError: name 'fmu' is not defined
[/color]

Tu définis par ailleurs deux fois la fonction n1(x1, x2), qui renvoient des résultats différents ; icelle fonction n'étant par ailleurs d'aucune utilité dans le fragment de code que tu proposes. Dans cette fonction apparaît une variable epsilon qui n'est définie nulle part. Tout cela n'est guère lisible pour quelqu'un d'autre que toi (et ce sera sans doute également illisible pour toi si tu relis ton code dans six mois).

Bref. Produis un code complet et minimal, spécifie clairement ce que tu souhaites, par exemple : j'ai les listes Ax, Bx d'abscisses des points A et B, les listes Ay et By des ordonnées, je veux relier le point A[k] au point B[k] par un trait bleu pour tout k, faire une croix bleue aux points A[k], une croix rouge aux points B[k]. Dans ce cas, on pourra se pencher sur ton problème.

brian · January 2020

Voici un exemple à peu près minimal. Je te laisse adapter avec tes données.

#! /usr/bin/env python3
# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
from matplotlib import collections as mc
import sys

def afficheOuSauve(fichierSortie):
    fig, ax = plt.gcf(), plt.gca()
    ax.set_aspect('equal', adjustable='datalim') # repère orthonormé

    if fichierSortie:
        fig.set_figwidth(10)
        plt.savefig(fichierSortie, dpi=100, bbox_inches="tight")
    else:
        fig.set_size_inches(20, 12)
        plt.show()

def main():
    fichierSortie = sys.argv[1] if len(sys.argv) >= 2 else ""

    l1 = [ (-7.5,1), (-5,-2), (-7,-4), (-9, 0) ]
    l2 = [ (5,2), (2,3), (4,5), (3, -3) ]
    abscisses1 = [ t[0] for t in l1 ]
    ordonnées1 = [ t[1] for t in l1 ]
    abscisses2 = [ t[0] for t in l2 ]
    ordonnées2 = [ t[1] for t in l2 ]

    lines = zip(l1, l2)
    lc = mc.LineCollection(lines, colors=(0.5, 0.9, 0.7, 0.5)) # RGBA

    plt.plot(abscisses1, ordonnées1, 'b+')
    plt.plot(abscisses2, ordonnées2, 'rx')
    ax = plt.gca()
    ax.add_collection(lc)
    afficheOuSauve(fichierSortie)

    sys.exit(0)

if __name__ == "__main__": main()

Si tu passes un paramètre comme /tmp/exemple.png, ça sauve la figure là-dedans. Si tu passes la chaîne vide ("" dans la plupart des shells), matpotlib affiche le graphe à l'écran et permet d'interagir avec la souris (zoom, translation, lecture de coordonnées...). 95398

Gentil · January 2020

Bonjour Benoit Rivet,

En effet j'ai eu des souci cette après midi avec pizo, il a gardé en mémoire des variables même d'autre ficher python compilé en parallèle. Donc quand je compilé il mélangeait tout et ça m'a fait perdre du temps. Bref j'imagine que c'est comme ça qu'on apprend ...

Ce que je voulais écrire est :

#Paramétrisation des espaces de départ et d'arrivée
s1 = np.linspace(-5,5,N, False)
s2 = np.linspace(-5,5,N, False)
S1,S2 = np.meshgrid(s1,s2)          #points de discrétisation de l'espace de départ

t1 = np.linspace(-5,5,N, False)
t2 = np.linspace(-5,5,N, False)
T1,T2 = np.meshgrid(t1,t2)          #points de discrétisation de l'espace de d'arrivé

#Définition des mesures
def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return 10**-5
    else:
        return np.exp(-(x1**2+x2**2)/2)/2*np.pi

for i in range(N):
    for j in range(N):
        mu[ i][j] = n1(s1[ i],s2[j])
mu = mu/np.sum(mu)

def n2(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return 10**-5
    else:
        return np.exp(-(x1**2+x2**2)/4)/4*np.pi

for i in range(N):
    for j in range(N):
        nu[ i][j] = n2(t1[ i],t2[j])
print(nu)
nu = nu/np.sum(nu)

Gentil · January 2020

Merci brian pour votre code.
Dans votre exemple vous séparez les bleus des rouges en prenant des abscisses négative et positive mais le problème c'est que si j'ai des points bleue et rouge de même coordonnée et très proche ça va être très vite illisible. D'où ma remarque il faudrait séparer artificiellement les échantillons.

Rescassol · January 2020

Bonsoir,

Gentil a écrit:

il a gardé en mémoire des variables même d'autre ficher python compilé en parallèle

Gentil, la troisième icone du menu du shell dans Pyzo permet de mettre la mémoire à 0, il est de ta responsabilité de l'utiliser quand c'est nécessaire.

Cordialement,

Rescassol

brian · January 2020

@Gentil

Mon code n'a que faire du signe des abscisses. J'ai dû inventer des coordonnées car tu n'en avais même pas fourni dans ta question. Il y a un premier nuage de points défini par la liste l1 et un second par la liste l2. Mon code trace un segment entre les éléments de l1 et de l2 de même indice et met les mêmes marques aux extrémités avec les mêmes couleurs que dans l'image que tu as postée. Il montre aussi comment choisir la couleur des segments, avoir un repère orthonormé, régler la taille de la figure interactive, de celle sauvée dans un fichier et de choisir à la volée l'une ou l'autre option depuis la ligne de commande. Bref, je crois que tu ferais mieux de lire attentivement le code.

Si tu as une seule liste de points, à toi de décider comment la partitionner (il y a sans doute plein de façons d'aborder le problème). Je traite ici la question de la visualisation. La partition d'un nuage de points en deux parties selon tel ou tel critère est un problème algorithmique. Si tu as besoin d'aide sur ce sujet, à mon avis il vaut mieux ouvrir un fil séparé avec un titre plus adapté.

Gentil · January 2020

Oui merci encore brian. J'ai deux ensembles pas besoin de partitionner.

Toute fois, c'est difficile de plot mes deux ensembles. Quand je le fais je ne vois rien, je pense que les points sont trop petit.

brian · January 2020

Quand on n'est pas sûr, il faut faire progressivement. Commence avec 3 points dans chaque ensemble et règle les problèmes d'échelle. Tu peux peut-être nous donner l1 et l2 avec 3 ou 4 points chacune ? J'espère au moins que tu obtiens une image. Car sinon, c'est tout autre chose et il faut s'intéresser aux messages d'erreurs.

brian · January 2020

Est-ce qu'au moins tu arrives à lancer mon script et obtenir l'image que j'ai postée ? Sinon, il faut te demander si matplotlib est installé et comment tu exécutes les instructions. C'est un script Python, il est fait pour être lancé par Python. Si tu exécutes les instructions dans une session interactive Pyzo, je ne suis pas étonné que ça ne marche pas. matplotlib supporte une utilisation interactive similaire à matlab, mais pas dans n'importe quel shell Python. D'après cette page, le shell standard Python est supporté ainsi que iPython à condition d'utiliser un backend interactif. Ceux-ci sont listés après la phrase :

And here are the user interfaces and renderer combinations supported; these are interactive backends, capable of displaying to the screen and of using appropriate renderers from the table above to write to a file

Parmi ces backends, TkAgg doit toujours être disponible si Python a été installé avec le toolkit Tk (case à cocher dans l'installation sous Windows). Qt5Agg et GTK3Agg doivent aussi être disponibles sur un système Linux décent. On peut aussi travailler avec Jupyter, mais je ne l'ai pas essayé en conjonction avec matplotlib.

Le plus simple pour commencer, ce serait quand même de lancer mon script en passant un chemin de fichier .png en paramètre. Ceci doit toujours fonctionner sans configuration particulière. Par exemple :

py -3 programme.py C:\sortie.png

ou

python3 programme.py /tmp/sortie.png

Si ça marche, remplace C:\sortie.png par "" ou par rien du tout. Si là ça ne marche plus, il y a sans doute un message d'erreur indiquant que le backend ne convient pas.

Gentil · January 2020

Oui Brian votre code fonctionne parfaitement sur ma machine.

Gentil · January 2020

En fait le raison je pense c'est que scatter ça marche que en 2D, il ne m'affiche pas les point (s1[l],s2[k],n2[l][k]) comme j'imaginais.

Gentil · January 2020

Bon ce n'est pas très grave si je ne peux pas plotter ça (vous m'excuserez l'expression :)o).

Ce que je peux faire en revanche, c'est utiliser le code de Brian et avec l1 et l2 mettre s1 et s2 par exemple l1 = (s1[k], s2[l] ) pour k, l et pareil pour l2. Je vais obtenir comme ça un quadrillage de -5 et 5.
Après comme j'ai dit le problème c'est que les points se confondent. Je pense qu'il va falloir plotter l2 + C une translation (dans chaque coordonnée de l2) afin que tout ne se mélange pas.

Une fois qu'on a fait ça, je vais tracer des traits entre les points qui m'intéressent avec deux couleurs différentes. Un pour les bons traits une pour le traits que j'ai trouvé par un certain algo.

Et là je serais pas mal. Après ce serait bien de pouvoir remplacer les croix du scatter par des ronds dont la taille est proportionnel à un certain coefficients. Mais ça peut vite être pas beau si les points sont trop proches.

Donc voila ma feuille de route, je vous envoie les avancées si j'ai des questions et que vous n'en avez pas marre de moi :-?. Merci pour votre soutient !

brian · January 2020

Je ne comprends rien aux $k$ et $l$ non définis. Une translation, c'est trivial avec NumPy et comme matplotlib travaille nativement sur des tableaux NumPy que tu utilises déjà dans ton code, ça ne devrait pas poser de problème. Pour les segments, le paramètre 'colors' de matplotlib.collections.LineCollection() peut être un vecteur (d'où la possibilité de donner une couleur différente à chaque segment). Pour les points de taille variable, en s'inspirant de l'exemple Plotting with keyword strings donné ici, on voit que l'on peut par exemple faire ceci dans mon code :

colors = len(l1)*[(0.5,0.5,1.0)]
plt.scatter(abscisses1, ordonnées1, s=10*np.abs(abscisses1), c=colors, marker='o')

afin de représenter chaque point de l1 à l'aide d'un disque bleu dont l'aire est proportionnelle à la valeur absolue de l'abscisse du point (c'est un juste exemple).

Edit : ceci suppose bien sûr d'avoir ajouté 'import numpy as np'.

Gentil · January 2020

En fait les $k$ et $l$ sont des indices. par s1[1],s2[1] ; s1[0],s2[3] ; s1[4],s2[2]
Voici un exemple de ce que j'obtiens.
À présent je vais passer à la représentation des flèches.
Merci pour le lien, à la fin j'essayerais de changer les croix en flèches.

Peut-être que je ne garderai pas tout les points si ce n'est pas beau. Je vous envoie ça prochainement. 95456

Gentil · January 2020

Alors là je trouve ça vraiment bizarre, ce matin j'excute le même code et voici ce que python m'affiche, ce n'est plus comme hier alors que pourtant c'est le même code.

#Paramétrisation générale
N = 10                    #pas de discrétisation de l'intervalle [0,1]. N[ i] correspond à la coordonnée i.
rho = 1

#Paramétrisation des espaces de départ et d'arrivée
s1 = np.linspace(-5,5,N, False)
s2 = np.linspace(-5,5,N, False)
S1,S2 = np.meshgrid(s1,s2)          #points de discrétisation de l'espace de départ

t1 = np.linspace(-5,5,N, False)
t2 = np.linspace(-5,5,N, False)
T1,T2 = np.meshgrid(t1,t2)          #points de discrétisation de l'espace de d'arrivé

#Définition des mesures
def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return 10**-5
    else:
        return np.exp(-(x1**2+x2**2)/2)/2*np.pi

for i in range(N):
    for j in range(N):
        mu[ i][j] = n1(s1[ i],s2[j])
mu = mu/np.sum(mu)

def n2(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) < 10**-3:
        return 10**-3
    else:
        return np.exp(-(x1**2+x2**2)/4)/4*np.pi

for i in range(N):
    for j in range(N):
        nu[ i][j] = n2(t1[ i],t2[j])

nu = nu/np.sum(nu)

#Matrice de coût
X1,Y1 = np.meshgrid(t1,s1)
X2,Y2 = np.meshgrid(t2,s2)
C1 = (X1-Y1)**2                    #Matrice contenant les c1[i,j]=(x1i-y1j)²
C2 = rho*(X2-Y2)**2               #Matrice contenant les c2[i,j]=lamb*(x2i-y2j)²

#Gamma = np.zeros([N,N,N,N])
#Gamma = np.round(Sinkhorn2D(N, C1,C2,mu,nu, 100, False), decimals = 3)

def Te(ep, i,j):
    return [s1[ i]*(2*ep+2), s1[j]*(2*ep+2)]/np.sqrt(2*ep**2 +
        2+4*ep)

list1 = []
for i in range(N):
    for j in range(N):
        list1.append( [s1[ i],s2[j]] )

list2 = []
for i in range(N):
    for j in range(N):
        list2.append( [t1[ i],t2[j]] )

list3 = []
for i in range(N):
    for j in range(N):
        a = Te(1,i,j)[0]
        b = Te(1,i,j)[1]
        list3.append( [a,b] )

abscisses1 = [ t[0] for t in list1 ]
ordonnées1 = [ t[1] for t in list1 ]

abscisses2 = [ t[0]+10 for t in list2 ]
ordonnées2 = [ t[1]+10 for t in list2 ]

abscisses3 = [ t[0]+10 for t in list3 ]
ordonnées3 = [ t[1]+10 for t in list3 ]

#plot mesure
plt.figure()
plt.plot(abscisses1, ordonnées1, 'b')
plt.plot(abscisses2, ordonnées2, 'r')
#plt.axes('off')

Gentil · January 2020

NB : Normalement c'est list2.append( [t1,t2[j]] ) mais latex comprend comme italique donc il y a un conflit... on devrait enlever le latex du l'option code formaté je pense.

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]

brian · January 2020

Le logiciel du forum (Phorum) interprète $[\text{i}]$ comme un début de balise « italique ». Quand tu postes du code Python ici, remplace les $[\text{i}]$ par

[ i]

et fais de même pour les $[\text{b}]$ et sans doute d'autres variantes (du genre $[\text{color=red}]$...).

Je n'ai pas lu ton script, mais s'il ne fait pas la même chose d'une fois à l'autre, c'est peut-être encore une mauvaise utilisation de Pyzo ? Ne vaudrait-il pas mieux exécuter le script en ligne de commande (flèche haut pour rappeler la commande précédente dans à peu près tous les shells et Enter pour l'exécuter, c'est assez rapide) ?

Edit : le problème est sans doute l'utilisation de 'b' et 'r' au lieu de 'b+' et 'rx' qui eux, précisent la forme des marques.

Gentil · January 2020

Bonjour Brian
Je pense que 'b' signifie bleu et 'b+' bleu foncé :)o. D'accord b+ et rx permettent de préciser autre chose que la couleur. Merci beaucoup !
Du coup là j'essaye de tracer les traits entre deux points mais pour l'instant ce n'est pas très fructueux.

def Te(ep, i,j):
    return [s1[ i]*(2*ep+2), s1[j]*(2*ep+2)]/np.sqrt(2*ep**2 +
        2+4*ep)

list1 = []
for i in range(N):
    for j in range(N):
        list1.append( [s1[ i],s2[j]] )

list2 = []
for i in range(N):
    for j in range(N):
        list2.append( [t1[ i],t2[j]] )

list3 = []
for i in range(N):
    for j in range(N):
        list3.append( Te(0,i,j) )

abscisses1 = [ t[0] for t in list1 ]
ordonnées1 = [ t[1] for t in list1 ]

abscisses2 = [ t[0]+12 for t in list2 ]
ordonnées2 = [ t[1]+12 for t in list2 ]

abscisses3 = [ t[0]+12 for t in list3 ]
ordonnées3 = [ t[1]+12 for t in list3 ]

#plot mesure
plt.figure()
#plt.subplot(1,2,1)
plt.plot(abscisses1, ordonnées1, 'b+')
plt.plot(abscisses2, ordonnées2, 'b+')
plt.plot(abscisses3,ordonnées3,'rx')
plt.title("Transport de Knothe entre Gaus1 et Gaus2")
plt.axis('off')

lines = zip(list1, list3)
lc = mc.LineCollection(lines, colors=(0.5, 0.9, 0.7, 0.5)) # RGB
ax = plt.gca()
ax.add_collection(lc)

[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD] 95484

brian · January 2020

Il y a bien des segments tracés sur la deuxième figure, non ? Tu as modifié mon commentaire RGBA, ce n'est pas bien ! RGBA, c'est quatre valeurs (utilise un moteur de recherche). Si tu gardes 4 valeurs mais que tu mets RGB en commentaire, c'est peu compréhensible. Si tu veux des retours concrets sur le code, fais en sorte qu'il tourne après un simple copier/coller dans un fichier vide, car ce n'est pas le cas pour l'instant (imports manquants, variables non définies).

Si tu travailles avec beaucoup de points, utilise la vectorisation comme a expliqué paul18 dans d'autres fils. Par exemple, ceci n'est pas efficace :

abscisses3 = [ t[0]+12 for t in list3 ]

Avec NumPy, on peut faire une translation de cette manière :

>>> import numpy as np
>>> l = np.array([0, 5, -2, 3], dtype=np.float64)
>>> l
array([ 0.,  5., -2.,  3.])
>>> l + 3.0
array([3., 8., 1., 6.])

Mais tant que tu en es au stade des essais avec peu de points, les méthodes « peu efficaces » ne posent pas de problème, bien sûr.

Gentil · January 2020

Oui tout va bien merci, pardon pour j'oublié d'écrire de la A dans RGBA mais je n'ai rien modifié. 95498

Gentil · January 2020

Ah je ne savais pas pour les listes, je modifierai ça effectivement si j'ai plus de données.
À présent, j'aimerais modifier un peu la couleur de trait pour le rendre plus foncé. Je pense que je peux le faire avec votre code mais je ne vois pas trop comment pour l'instant.
Deuxièmement j'aimerais faire des points dont l'aire varie en fonction d'une matrice de poids comme ci- dessous, au lieu des croix. Toutefois dans votre lien, il faut utiliser scatter. Nous nous utilisons plot. Vous pensez qu'on peut quand même le faire ? Parce que j'ai essayé de remplacer plot par scatter mais ça bug.

brian · January 2020

Gentil a écrit:

j'aimerais modifier un peu la couleur de trait pour le rendre plus foncé.

Tu es sérieux, là ? La couleur des segments est donnée par colors=(0.5, 0.9, 0.7, 0.5) et j'ai précisé que c'est au format RGBA. As-tu cherché RGBA sur Internet (en ajoutant éventuellement le mot-clé « color » ou « couleur ») ?

Quant aux points (disques) dont la taille et la couleur peuvent être différentes pour chacun d'entre eux, je t'ai déjà donné la réponse plus haut. Tu prends l'exemple complet minimal que j'ai donné au départ, tu ajoutes 'import numpy as np' en haut et les deux lignes que j'ai indiquées ici (par exemple, avant 'ax = plt.gca()') et tu vas voir que ça marche. Il faut un minimum d'initiative, tout de même !

Gentil · January 2020

Excusez-moi de vous avoir fait répéter, c'est nouveau pour moi donc je suis un peu lent à la détente. Mais là après quelques recherches sur internet, j'ai compris pour les couleurs merci.

Voilà j'ai quasiment fini, merci pour votre aide. Peut-être que je vais ajouter quelques bricoles.
Voilà le code :

#Paramétrisation générale
N = 10                    #pas de discrétisation de l'intervalle [0,1]. N[ i] correspond à la coordonnée i.
rho = 1

#Paramétrisation des espaces de départ et d'arrivée
s1 = np.linspace(-5,5,N, False)
s2 = np.linspace(-5,5,N, False)
S1,S2 = np.meshgrid(s1,s2)          #points de discrétisation de l'espace de départ

t1 = np.linspace(-5,5,N, False)
t2 = np.linspace(-5,5,N, False)
T1,T2 = np.meshgrid(t1,t2)          #points de discrétisation de l'espace de d'arrivé

#Définition des mesures
def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return 10**-5
    else:
        return np.exp(-(x1**2+x2**2)/2)/2*np.pi

mu = np.zeros([N,N])
for i in range(N):
    for j in range(N):
        mu[ i][j] = n1(s1[ i],s2[j])
mu = mu/np.sum(mu)

def n2(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) < 10**-3:
        return 10**-3
    else:
        return np.exp(-(x1**2+x2**2)/4)/4*np.pi

nu = np.zeros([N,N])
for i in range(N):
    for j in range(N):
        nu[ i][j] = n2(t1[ i],t2[j])

nu = nu/np.sum(nu)

def Te(ep, i,j):
    return [s1[ i]*(2*ep+2), s1[j]*(2*ep+2)]/np.sqrt(2*ep**2 +
        2+4*ep)

list1 = []
for i in range(N):
    for j in range(N):
        list1.append( [s1[ i],s2[j]] )

list2 = []
for i in range(N):
    for j in range(N):
        list2.append( [t1[ i],t2[j]] )

list3 = []
for i in range(N):
    for j in range(N):
        list3.append( Te(0,i,j)+12 )

abscisses1 = [ t[0] for t in list1 ]
ordonnées1 = [ t[1] for t in list1 ]

abscisses2 = [ t[0]+12 for t in list2 ]
ordonnées2 = [ t[1]+12 for t in list2 ]

abscisses3 = [ t[0] for t in list3 ]
ordonnées3 = [ t[1] for t in list3 ]

#plot mesure
plt.figure()
#plt.subplot(1,2,1)
plt.plot(abscisses1, ordonnées1, 'b+')
plt.plot(abscisses2, ordonnées2, 'b+')
plt.plot(abscisses3,ordonnées3,'rx')
plt.title("Transport de Knothe entre Gaus1 et Gaus2")
plt.axis('off')

lines = zip(list1, list3)
lc = mc.LineCollection(lines, colors= (0.2, 0.9, 0.1)) # RGBA
ax = plt.gca()
ax.add_collection(lc)

plt.show()
plt.figure()
l1 = list1
colors = len(l1)*[(0.5,0.5,1.0)]

list_weight1 = []
C = np.max(mu)
for i in range(N):
    for j in range(N):
        list_weight1.append( 250*mu[ i][j]/C+10)

list_weight2 = []
C = np.max(nu)
for i in range(N):
    for j in range(N):
        list_weight2.append( 250*nu[ i][j]/C+10)

plt.scatter(abscisses1, ordonnées1, s=list_weight1, c=colors, marker='o')
#plt.scatter(abscisses2, ordonnées2, s=list_weight2, c=colors, marker='o')
plt.scatter(abscisses3, ordonnées3, s=list_weight2, c='r', marker='o')

lines = zip(list1, list3)
lc = mc.LineCollection(lines, colors= (0.2, 0.9, 0.1)) # RGBA
ax = plt.gca()
ax.add_collection(lc)

plt.show()

[large][Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD][/large] 95504

Gentil · January 2020

J'ai une dernière question, j'ai fait une boucle for pour afficher différent graphique. J'aimerais bien en faire une animation, est-ce que vous savez comment faire ? Si c'est trop compliqué pour moi je comprends, mais je me permets de vous demander car vous avez l'air très calés et moi je suis curieux (:D.

brian · January 2020

C'est un classique pour ImageMagick (programme 'convert'). Il en est question ici, par exemple. La commande :

convert -delay 20 F_*.png -loop 0 movie.gif

a l'air correcte, de mémoire. Il me semble que sous Windows, il faut utiliser 'magick convert' au lieu de 'convert' sous peine d'appeler un exécutable de Windows (discuté ici récemment). On peut ensuite sans doute convertir le GIF animé vers un autre format (plus sophistiqué et économique en espace disque, tel qu'XviD ou H.264) avec FFmpeg, voire tout faire avec ce dernier.

Gentil · January 2020

Merci brian ! Sur ce coup j'ai été trop ambitieux je pense ! (:P)

Surtout que le code sur lequel vous avez aidé n'est pas fini. Il y a un détail auquel je n'avais pas pensé... Si le nuage cible est beaucoup plus gros (espacés) que le nuage source alors pour l'affichage on y verra rien !
Du coup j'ai travaillé (un peu d'initiative comme vous me l'avez dit ^^) et mon idée c'est simple :
=> Calculer la plus grande distance entre le points de l'image cible et normaliser par une constante afin de les rapprocher artificiellement.

Le problème de cette méthode c'est que :
1) Le temps de calcul j'ai des boucles en N**2.
2) Le resultat n'est pas top jugeant par vous même, Te(0, i,j) est une application linéaire de matrice $\begin{bmatrix} 3 & 1 \\ 1 & 4 \end{bmatrix}$ mais ça ne ressemble pas vraiment à une dilatation.

Voici le code vous comprendrez mieux.

[large][Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD][/large]

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 14 14:36:24 2020
@author: Roméo
"""
# Math
import numpy as np
import math
import scipy.stats as stat

# Graphe
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.colors as colors
from matplotlib import collections as mc
from matplotlib.ticker import LinearLocator, FormatStrFormatter


#Paramétrisation générale
N = 10                      #pas de discrétisation de l'intervalle [0,1]. N[ i] correspond à la coordonnée i.
#epsilon = .01
arret = .01
#rho = epsilon**2            #paramètre à valeur dans [0,1], permettant d'obtenir le transport de Knothe proche de 0

#Paramétrisation des espaces de départ et d'arrivée
x1min = -5
x1max = 5
x2min = -5
x2max = 5
y1min = -5
y1max = 5
y2min = -5
y2max = 5

s1 = np.linspace(x1min,x1max,N)
s2 = np.linspace(x2min,x2max,N)
t1 = np.linspace(y1min,y1max,N)
t2 = np.linspace(y2min,y2max,N)
S1,S2 = np.meshgrid(s1,s2,indexing='ij')          #points de discrétisation de l'espace de départ
T1,T2 = np.meshgrid(t1,t2,indexing ='ij')          #points de discrétisation de l'espace de départ


#Paramétrisation des espaces de départ et d'arrivée
s1 = np.linspace(-5,5,N, False)
s2 = np.linspace(-5,5,N, False)
S1,S2 = np.meshgrid(s1,s2,indexing='ij')          #points de discrétisation de l'espace de départ

t1 = np.linspace(-5,5,N, False)
t2 = np.linspace(-5,5,N, False)
T1,T2 = np.meshgrid(t1,t2,indexing='ij')          #points de discrétisation de l'espace de d'arrivé

## Définition des mesures

#Densité gaussienne centrée réduite
def n1(x1,x2):
    if  np.exp(-(x1**2+x2**2)/2) ==0:
        return 10**-5
    else:
        return np.exp(-(x1**2+x2**2)/2)/2*np.pi

# Mesure gaussienne 1
mu = np.zeros([N,N])
for i in range(N):
    for j in range(N):
        mu[ i][j] = n1(s1[ i],s2[j])
mu = mu/np.sum(mu)

#Densité gaussienne centrée variance :  3 1
#                                       1 4
def n2(x1,x2):
    if  np.exp(-(4*x1**2+3*x2**2 -2*x1*x2)/22)/(np.sqrt(11))*2*np.pi < 10**-3:
        return 10**-3
    else:
        return np.exp(-(4*x1**2+3*x2**2 -2*x1*x2)/22)/(np.sqrt(11))*2*np.pi

#Mesure gaussienne 2
nu = np.zeros([N,N])
for i in range(N):
    for j in range(N):
        nu[ i][j] = n2(t1[ i],t2[j])
nu = nu/np.sum(nu)

## Sinkhorn 2D
def SinkHorn2D(rho):
    #Matrice de coût
    X1,Y1 = np.meshgrid(t1,s1,indexing='ij')
    X2,Y2 = np.meshgrid(t2,s2,indexing='ij')
    C1 = rho*(X1-Y1)**2                    #Matrice contenant les c1[ i,j]=(x1i-y1j)²
    C2 = (X2-Y2)**2                         #Matrice contenant les c2[ i,j]=lamb*(x2i-y2j)²

    #Calcul de Gamma_Barre
    GammaB1 = np.exp(-C1/epsilon)      #Matrice contenant les GammaB1[ i,k] = Gamma_Barre1i,k
    GammaB2 = np.exp(-C2/epsilon)      #Matrice contenant les GammaB2[j,m] = Gamma_Barre2j,m

    a2 = np.ones([N,N])
    b2 = np.ones([N,N])
    erreur = math.inf
    n = 0
    while (erreur>arret and n<=10000):
        a1 = a2
        b1 = b2
        a2 = mu/np.dot(np.dot(GammaB1,b1),GammaB2.T)     #Dans le cas où mu est nulle à certain endroit, a2 peut être nulle, ce qui risque de créer des problèmes
        b2 = nu/np.dot(np.dot(GammaB1.T,a2),GammaB2)
        erreur = np.linalg.norm((a2-a1))/np.linalg.norm(a1/2+a2/2,np.inf)+np.linalg.norm((b2-b1))/np.linalg.norm(b1/2+b2/2,np.inf)    #l'arrêt doit être relatif car sinon ça diverge
        n+=1
        #print(erreur)
    #print(n)
    #print(erreur)
    return a1,b1


def transport(x,y, rho, epsilon):

    #Calcul de Sinkhorn
    a,b = SinkHorn2D(rho)

    #Matrice de coût
    X1,Y1 = np.meshgrid(t1,s1,indexing='ij')
    X2,Y2 = np.meshgrid(t2,s2,indexing='ij')
    C1 = rho*(X1-Y1)**2                    #Matrice contenant les c1[ i,j]=(x1i-y1j)²
    C2 = (X2-Y2)**2                         #Matrice contenant les c2[ i,j]=lamb*(x2i-y2j)²

    #Calcul de Gamma_Barre
    GammaB1 = np.exp(-C1/epsilon)      #Matrice contenant les GammaB1[ i,k] = Gamma_Barre1i,k
    GammaB2 = np.exp(-C2/epsilon)      #Matrice contenant les GammaB2[j,m] = Gamma_Barre2j,m

    #Calcul de Transport_optimal(x,y, rho)
    if x1min<=x<x1max and x2min<=y<x2max:
        i = int((x-x1min)/(x1max-x1min)*N)
        j = int((y-x2min)/(x2max-x2min)*N)
        ei = np.zeros(N)
        ei[ i] = 1
        GammaiB = np.dot(ei,GammaB1)
        ej = np.zeros(N)
        ej[j] = 1
        GammajB = np.dot(ej,GammaB2)
        return a[ i,j]*((b.T*GammaiB).T*GammajB) #L'erreur était là
    else:
        return print("le point ("+str(x)+","+str(y)+") n'est pas dans l'intervalle de départ ["+str(x1min)+","+str(x1max)+"]x["+str(x2min)+","+str(x2max)+"]")

def display(x,y):
    Vij = transport(x,y,rho, epsilon).T

    plt.close('all')
    im = imshow(Vij, cmap=cm.seismic, origin='lower')
    colorbar(im)
    title('transport du point ('+str(x)+', '+str(y)+'), $\epsilon$='+str(epsilon)+' , $\rho$='+str(rho))
    show()

def TKS(x1,x2,rho, epsilon):
    Vij = transport(x1,x2,rho, epsilon)
    scale = np.sum(np.sum(Vij))
    y1 = np.sum(np.sum(Vij*T1))/scale
    y2 = np.sum(np.sum(Vij*T2))/scale
    return y1,y2

## plot des transport

# Transport optimal
def Te(ep,i,j):
    return [s1[ i]*(3*ep+np.sqrt(11))+s2[j], ep*s1[ i]+s2[j]*(4+ep*np.sqrt(11))]/np.sqrt(3*ep**2 + 4 +
        2*ep*np.sqrt(11))

#Rho : C1 + rho*C2 et T_ep converge vers T_knothe lorsque ep tend vers 0
def plot(rho, epsilon):
    C = np.max(mu)
    D = np.max(nu)

    # Plot transport optimal

    #Grille -5;5 source
    list1 = []
    for i in range(N):
        for j in range(N):
            list1.append( [s1[ i],s2[j]] )

    abscisses1 = [ t[0] for t in list1 ]
    ordonnées1 = [ t[1] for t in list1 ]

    # Coordonée image qu'on décale de la source
    list3 = []
    for i in range(N):
        for j in range(N):
            list3.append( Te(rho,i,j) )

    abscisses3 = [ t[0] for t in list3 ]
    ordonnées3 = [ t[1] for t in list3 ]

    # Les point des images sont trop espacé on voit rien donc je les rapproche munuellement par dilatation
    E = 1

    for i in range(1,N):
        for j in range(1,N):
            for k in range(1,N):
                for l in range(1,N):
                    a0 = abscisses3[ i]
                    a1 = abscisses3[j]
                    b0 = ordonnées3[k]
                    b1 = ordonnées3[l]
                    a = np.array(a0,a1)
                    b = np.array(b0,b1)
                    if np.linalg.norm(a-b) > E:
                        E = np.linalg.norm(a-b)

        # Modification de l'image
    listTO = []
    for i in range(N**2):
        temp1 = list3[ i][0]
        temp2 = list3[ i][1]
        listTO.append( [(np.sqrt(200)/E)*temp1+12,(np.sqrt(200)/E)*temp2+12] )

    abscissesT0 = [ t[0] for t in listTO ]
    ordonnéesT0 = [ t[1] for t in listTO ]

    # Les poids
    list_weight1 = []
    for i in range(N):
        for j in range(N):
            list_weight1.append( 250*mu[ i][j]/C+10)

    MW2max = 1
    for i in range(N):
        for j in range(N):
            temp1 = Te(rho,i,j)[0]
            temp2 = Te(rho,i,j)[1]
            if np.linalg.norm( np.array(temp1,temp2) ) > MW2max:
                MW2max = np.linalg.norm( np.array(temp1,temp2) )

    list_weight2 = []
    for i in range(N):
        for j in range(N):
            temp1 = Te(rho,i,j)[0]
            temp2 = Te(rho,i,j)[1]
            a = n2( temp1,temp2 )
            list_weight2.append( (500*a/MW2max)+10)

    # plot transport optimal par Sinkhorn

    #On a un problème au frontière TKS ne marche pas

    #Grille -5;5  sans le bord
    list2 = []
    for i in range(N):
        for j in range(N):
            list2.append( [s1[ i],s2[j]] )

    abscisses2 = [ t[0] for t in list2 ]
    ordonnées2 = [ t[1] for t in list2 ]

    list_star2 = []
    for i in range(N):
        for j in range(N):
            temp1 = TKS(s1[ i],s2[j],rho, epsilon)[0]
            temp2 = TKS(s1[ i],s2[j],rho, epsilon)[1]
            test = np.array(temp1,temp2)
            list_star2.append([temp1,temp2])

    abscisses_star2 = [ t[0] for t in list_star2 ]
    ordonnées_star2 = [ t[1] for t in list_star2 ]

    Estar = 1

    for i in range(N):
        for j in range(N):
            for k in range(N):
                for l in range(N-1):
                    a0 = abscisses_star2[ i]
                    a1 = abscisses_star2[j]
                    b0 = ordonnées_star2[k]
                    b1 = ordonnées_star2[l]
                    a = np.array(a0,a1)
                    b = np.array(b0,b1)
                    if np.linalg.norm(a-b) > Estar:
                        Estar = np.linalg.norm(a-b)

        # Modification de l'image
    listTOS = []
    for i in range((N)**2):
        temp1 = list_star2[ i][0]
        temp2 = list_star2[ i][1]
        listTOS.append( [(np.sqrt(200)/Estar)*temp1+12,(np.sqrt(200)/Estar)*temp2+12] )

    abscissesTOS = [ t[0] for t in listTOS ]
    ordonnéesTOS = [ t[1] for t in listTOS ]

    list_weight_star2 = []

    MW2max_star = 1
    for i in range(N):
        for j in range(N):
            a0 = abscisses_star2[ i]
            b0 = ordonnées_star2[j]
            if np.linalg.norm( np.array(a0,b0) ) > MW2max_star:
                MW2max_star = np.linalg.norm( np.array(a0,b0) )

    for i in range(N):
        for j in range(N):
            temp1 = TKS(s1[ i],s2[j],rho,epsilon)[0]
            temp2 = TKS(s1[1],s2[j],rho,epsilon)[1]
            a = n2( temp1,temp2 )
            list_weight_star2.append( 500*a/MW2max_star+10 )



    # plot (color : (Red,Green,Bleu)RGBA)
    colors = len(list1)*[(0.5,0.5,1.0)]

    plt.figure()
    plt.subplot(1,2,1)
    plt.scatter(abscisses1, ordonnées1, s=list_weight1, c=colors, marker='o')
    plt.scatter(abscissesT0, ordonnéesT0, s=list_weight2, c='r', marker='o')
    plt.title("Transport de Knothe(rho = {}) entre Gaus1 et Gaus2".format(rho))
    lines = zip(list1, listTO)
    lc = mc.LineCollection(lines, colors= (0.2, 0.9, 0.1))
    ax = plt.gca()
    ax.add_collection(lc)
    plt.axis('off')

    plt.subplot(1,2,2)

    colors_star = len(list2)*[(0.5,0.5,1.0)]

    plt.scatter(abscisses2, ordonnées2, s=list_weight1, c=colors_star, marker='o')
    plt.scatter(abscissesTOS, ordonnéesTOS, s=list_weight_star2, c='r', marker='o')
    plt.title("Transport de Knothe(rho = {}) entre Gaus1 et Gaus2 par Sinkhorn(epsilon = {})".format(rho,epsilon))
    lines_star = zip(list2, listTOS)
    lc_star = mc.LineCollection(lines_star, colors= (0.2, 0.9, 0.1))
    ax = plt.gca()
    ax.add_collection(lc_star)
    plt.axis('off')

    plt.show()
epsilon = 0.0001
rho = epsilon**2
plot(rho,epsilon)

Du coup si vous avez des remarques pour améliorer mon code c'est avec plaisir.

[large][Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD][/large]
[Pour un long programme, le mieux serait de le joindre sous forme de fichier. AD]

Gentil · January 2020

Du coup je réitère ma question de départ, comment relier deux nuages de points sachant que la distance maximale entre les points respectifs des deux nuages peut avoir une échelle très différente. Par exemple x10 ou plus.

Gentil · January 2020

Je rappelle qu'il y a un problème avec Latex pour copier mon code faites "citer" là ça ira je pense.

[Si tu lisais les indications qu'on te donne ! :-X AD]
[Mettre une ' ' entre [ et i sinon, l'afficheur du forum le prend pour une bannière BBcode de mise en italique. AD]
[Pour un long programme, le mieux serait de le joindre sous forme de fichier. AD]

brian · January 2020

J'espère que tu as bien lu les commentaires d'AD ci-dessus. 8-)

Désolé, le temps dont je dispose ne me permet pas de lire sérieusement un gros morceau de code comme ça, d'autant que je ne connais absolument rien à la partie scientifique de ton programme. Ce que je peux dire à vue de nez :

1. Si tu recherches de bonnes performances, tu n'utilises sans doute pas assez la vectorisation, mais la priorité est d'abord d'avoir un code qui marche (même sur un petit nombre de points) avant d'optimiser.

2. La matrice que tu as donnée :
\[ \begin{bmatrix} 3 & 1 \\ 1 & 4 \end{bmatrix} \]
ne ressemble pas à une matrice de dilatation, en effet. Mais n'ayant pas vu comment tu l'utilises, je ne peux guère commenter davantage...

3. J'ai regardé le passage suivant :

    abscisses3 = [ t[0] for t in list3 ]
    ordonnées3 = [ t[1] for t in list3 ]

    # Les pointS des images sont trop espacéS ; on ne voit rien, donc je les rapproche manuellement par dilatation
    E = 1

    for i in range(1,N):
        for j in range(1,N):
            for k in range(1,N):
                for l in range(1,N):
                    a0 = abscisses3[ i]
                    a1 = abscisses3[j]
                    b0 = ordonnées3[k]
                    b1 = ordonnées3[l]
                    a = np.array(a0,a1)
                    b = np.array(b0,b1)
                    if np.linalg.norm(a-b) > E:
                        E = np.linalg.norm(a-b)

        # Modification de l'image
    listTO = []
    for i in range(N**2):
        temp1 = list3[ i][0]
        temp2 = list3[ i][1]
        listTO.append( [(np.sqrt(200)/E)*temp1+12,(np.sqrt(200)/E)*temp2+12] )

Je crois que tu essaies de mettre dans E le maximum des distances entre deux points d'un nuage puis d'utiliser ça pour faire une mise à l'échelle, mais le premier paquet de boucles me laisse perplexe (je suis peut-être un peu bête) :

a) Que signifie 'np.array(a0,a1)' ? Je ne connais pas cette syntaxe ; de plus, quand je teste dans le REPL des trucs comme 'np.array(3, 5)', je me prends des 'TypeError: data type not understood'.

b) np.linalg.norm(a-b), c'est bien une norme appliquée à un truc (je ne sais pas quoi à cause du a) qui, en gros, est formé de différences $\text{abscisse} - \text{ordonnée}$ ? Je ne comprends pas.

c) Les range(1,N) vont de 1 à N-1, est-ce bien voulu ?

d) Si tu es préoccupé par la vitesse, vu le niveau d'imbrication des boucles, il faut probablement utiliser moins de variables intermédiaires peu utiles (a0, a1, b0, b1)... et vectoriser. Par exemple, le + 12 fait à chaque élément séparément, bof bof. Mais pas de précipitation : la priorité numéro 1, c'est de faire marcher le code comme souhaité avant de songer à optimiser.

Gentil · January 2020

Bonjour Brian et AD.

Veillez m'excuser AD, la prochaine j'essayerai de joindre un fichier.

Je comprend que vous ne pouvez pas tout lire Brian, en tout cas vos commentaire m'aide beaucoup et je vous remercie.

Gentil · January 2020

Vous avez tout à fait raison, j'essaie de mettre dans E le maximum des distances entre deux points d'un nuage cible puis d'utiliser ça pour faire une mise à l'échelle.

a) Que signifie 'np.array(a0,a1)' ? np.array permet d'utiliser des tableaux numpy. Ca permet de voir les tableaux comme des vecteurs afin, par suite de calculer leurs normes.

Je pense que j'ai bien corrigé mon code. Le seul problème à présent c'est que sur l'image suivante les deux nuages cibles ne sont pas les mêmes.
Attention peut-être que cela est du à mon algo. Ou alors peut être que c'est du à mon code, ce qui est pas surprenant et dans ce cas j'aimerais bien le corrigé. 95836

Gentil · January 2020

J'aimerais bien vous envoyez mon code mais comment faire ?

Gentil · January 2020

Après il faudrait effectivement que j’accélère l'algo car pour N >= 20 ça rame déjà.

brian · January 2020

@Gentil

Merci, mais je sais à quoi sert numpy.array(). J'ai posé des questions précises sur un bout de code précis que tu as toi-même posté ici. Une partie de ce code est :

a0 = abscisses3[ i]
a1 = abscisses3[j]
b0 = ordonnées3[k]
b1 = ordonnées3[l]
a = np.array(a0,a1)

a0 et a1 ont bien l'air d'être des abscisses, non ? Eh bien, cet usage ne correspond ni à la documentation de la fonction numpy.array(), ni à ce que l'on peut trouver dans le tutoriel NumPy que j'ai déjà cité dans ce fil. Enfin, pas tout à fait ; il y a bien un endroit du tutoriel qui utilise la fonction comme toi, c'est celui-ci :

A frequent error consists in calling array with multiple numeric arguments, rather than providing a single list of numbers as an argument.
>>> a = np.array(1,2,3,4)    # WRONG
>>> a = np.array([1,2,3,4])  # RIGHT

Tu vois la ligne où il y a marqué WRONG ? Tu comprends pourquoi je t'ai posé la question, pourquoi je t'ai dit que j'ai même testé ça dans le REPL Python et que ça génère une erreur dont j'ai précisé le message exact ?

Gentil · January 2020

Ah oui d'accord, ben j'ai dû sûrement mal m'y prendre. Mon but est simplement de calculer la distance maximale entre les points du nuage de points formé par list3. Comment feriez-vous ?

brian · January 2020

On en a déjà parlé récemment, notamment ici. La racine carrée est inutile pour ce que tu veux faire, puisque c'est une fonction (strictement) croissante sur $\R_{+}$ (edité : pas $\R$ !).

Gentil · February 2020

D'accord merci ...

Nuage de points python

Réponses

Bonjour!

Catégories

In this Discussion

Qui est en ligne 15