[Homologie] Comptez le nombre de trous dans les données avec Python

Jupyter Notebook est un outil pratique qui vous permet de vérifier les résultats des calculs en cours de route lors de l'écriture d'un programme. Un autre avantage de l'utilisation de Google Colaboratory est que vous pouvez calculer quel que soit votre environnement.

L'une des méthodes d'analyse des données qui est devenue un sujet brûlant ces jours-ci est l'analyse des données par étapes. Il s'agit d'une méthode d'analyse qui capture les données sous forme de formulaire et compte le nombre d'entités telles que les trous. Tout d'abord, comme son nom l'indique, nous discuterons de la topologie.

Qu'est-ce que la topologie?

Dans wikipedia, «La topologie est une propriété (propriété de phase ou phase» qui peut être conservée même si elle est continuellement déformée (étirée ou pliée mais pas coupée ou collée) sous une forme (forme ou «espace»)). Il se concentre sur (invariant). C'est dit.

Beaucoup d'entre vous ont peut-être entendu dire que les «tasses à café» et les «beignets», qui sont souvent utilisés pour expliquer la topologie, sont les mêmes du point de vue de la topologie.

Cité de la géométrie de phase de wikipedia

A partir de cette figure, on peut voir que la "tasse à café" et le "beignet" peuvent être déformés par pliage et étirement sans couper ni coller. Cela signifie que le nombre de trous, qui est un invariant de phase, a une chose en commun. Dans cet article, nous avons calculé l'homologie, qui compte le nombre de trous dans une figure à deux dimensions.

Procédure de calcul d'homologie

Le calcul spécifique de l'homologie est Entrée: coordonnées des données, rayon des données Sortie: nombre de trous dans les données Je vais le faire comme.

Entrer des données

Saisissez les données que vous souhaitez calculer. Cette fois, le calcul a été effectué à l'aide de données créées aléatoirement.

import os
import random as rd
import numpy as np
import copy
#Création de données de nuage de points
#Dimension des données
dim = 2
#Le nombre de données
data_size = 500
#Rayon complexe
sim_r = 0.08
def make_point_cloud(dim,size):
    arr = np.random.rand(size, dim).astype(np.float64)
    return arr
data = make_point_cloud(dim,data_size)

Ensuite, afin de calculer le nombre de trous, le corps connecté dans lequel les données sont connectées est calculé.

Calcul de la concaténation

Puisque les données sont juste pour s'asseoir, la figure ci-dessous montre les données avec un rayon.

def distance(data1,data2):
    r = data1 - data2;
    r = np.sqrt(np.sum(np.power(r,2.0)))
    return r
def make_path_matrix(data,r):
    size ,data_dim = data.shape
    path_matrix = [[0 for j in range(size)]for i in range(size)]
    for i in range(size):
        for j in range(size):
            if distance(data[i],data[j]) <= 2*r:
                path_matrix[i][j] = 1;
    return path_matrix
#Créer une matrice adjacente
pathM = make_path_matrix(data,sim_r)
#Indicateur d'arrivée Deta non atteint: 0 atteint: 1
flaglist = [0 for i in range(data_size)]
#Une fonction qui recherche récursivement une matrice adjacente
def visit(i,path):
    flaglist[i] = 1
    for j in range(data_size):
        if(path[i][j] == 1 and flaglist[j] == 0):
            visit(j,path)
#Corps connecté
K = []
for i in range(data_size):
    #Enregistrez l'un des conjugués
    comp = [];
    #Enregistrer la flaglist avant l'exploration
    flaglistsave = copy.deepcopy(flaglist);
    #Vérifiez s'il a été fouillé
    if flaglist[i] == 0:
        #Rechercher sans recherche
        visit(i,pathM);
        #Consultez la liste des indicateurs de recherche
        for flagindex in range(data_size):
            #Confirmer les points recherchés et nouvellement recherchés
            if flaglist[flagindex] == 1 and flaglistsave[flagindex] == 0:
                #Ajouter à la concaténation
                comp.append(flagindex);
        #Ajouté à l'ensemble des concaténations
        K.append(comp);

En donnant aux données un rayon, on peut voir que les cercles peuvent être divisés en connexions connectées. Ensuite, nous étudierons la plus petite unité contenue dans le conjugué.

Créer une seule unité

Une seule unité est la plus petite unité qui constitue un corps connecté, et une seule unité est un point, une ligne ou une surface.

Lorsqu'il n'y a que des points, on l'appelle 0 unité, une ligne s'appelle 1 unité et une surface triangulaire s'appelle 2 unité. Cette fois, le nombre de trous en 2 dimensions est limité à 2 unités, mais en 3 dimensions, le tétraèdre devient 3 unités, ce qui est défini en dimensions supérieures.

Dans la configuration à 2 unités, comme le montre la figure ci-dessous, il y a des cas où chaque morceau du triangle constitue 1 unité mais ne constitue pas 2 unités. ~~ Il semble que vous puissiez utiliser le composé alpha, qui a une vitesse de calcul rapide, mais je ne sais pas, donc je l'ai calculé avec le composé de contrôle ~~

#2 Vérifiez la configuration d'une seule unité et améliorez
def checksim2(point,r):
    a = point[0];
    b = point[1];
    c = point[2];
    da = distance(b,c)
    db = distance(a,c)
    dc = distance(b,a)
    #Vérifiez si chacun des trois côtés constitue une unité
    sim1flag = da <(2*r) and db<(2*r) and dc<(2*r)
    #Officiel de Heron
    ss = (da + db + dc)/2
    #print(ss * (ss - da)*(ss - db)*(ss - dc))
    S = np.sqrt(ss * (ss - da)*(ss - db)*(ss - dc))
    MaxS = np.power(r,2) * np.sin(np.pi/6)*np.cos(np.pi/6)/2*3
    #Confirmation de la zone
    Sflag = S<MaxS
    #Jugement du cercle externe
    cosC = (np.power(da,2) + np.power(db,2) - np.power(dc,2))/(2*da*db)
    sinC =np.sqrt(1 - np.power(cosC,2))
    outerflag = 2*sinC*sim_r > dc
    if (Sflag and sim1flag) or (outerflag and sim1flag):
        return 1
    else:
        return 0
#Nombre de connecteurs
betti0 = len(K)
#1 Liste unique
sim1list = [[] for i in range(betti0)]
#2 Liste unique
sim2list = [[] for i in range(betti0)]

#Extraire une concaténation avec une seule unité 1 et 2 unités
for i in range(betti0):
    if len(K[i]) <4 and len(K[i]) != 1:
        if len(K[i]) == 2:
            sim1list[i] = copy.deepcopy([K[i]])
        if len(K[i]) == 3 and checksim2([data[K[i][0]],data[K[i][1]],data[K[i][2]]],sim_r) == 1:
            sim2list[i] = copy.deepcopy([K[i]])
#Extraire 1 unité du conjugué
for b in range(betti0):
    if len(K[b]) > 2:
        for i in range(1,len(K[b]),1):
            for j in range(0,i,1):
                if pathM[K[b][i]][K[b][j]] == 1:
                    sim1list[b].append([K[b][i],K[b][j]])
#Extraire 2 unités du conjugué
for b in range(betti0):
    if len(K[b]) > 3:
        for i in range(2,len(K[b]),1):
            for j in range(1,i,1):
                for k in range(0,j,1):
                    l = [data[K[b][i]],data[K[b][j]],data[K[b][k]]]
                    if checksim2(l,sim_r) == 1:
                        sim2list[b].append([K[b][i],K[b][j],K[b][k]])

En comptant les unités individuelles dans la concaténation, un ensemble contenant un nombre fini d'unités uniques est créé. Ensuite, une cartographie pour convertir une seule unité en une seule unité d'une autre dimension et une cartographie des limites seront décrites.

Cartographie des limites

La cartographie des limites est une cartographie qui transforme une seule unité en une petite unité unique d'une dimension.

Comme le montre la figure, il s'agit d'un mappage qui convertit de 2 unités en somme alternée de 1 unité. La cartographie concrète est la suivante

\langle 1,2,3 \rangle = \langle 2,3 \rangle - \langle 1,3 \rangle + \langle 1,2 \rangle

Le nombre supprimé de l'unité d'origine est un mappage qui prend la somme en changeant le signe du nombre de l'avant avec des nombres pairs ou impairs. La cartographie des limites qui passe de 2 unités à 1 unité contenue dans le conjugué s'appelle la cartographie quadratique des limites, et la cartographie des limites qui transfère de 1 unité à 0 unités est appelée la cartographie des limites principale. Le mappage des limites du second ordre et le mappage des limites du premier ordre obtenus sont représentés comme une matrice d'expression. Ensuite, le calcul pour compter le nombre de trous de cet agoniste de frontière sera décrit.

#Trier la concaténation, 1 unité, 2 unités
cmp = [sorted(l, reverse=True) for l in K]
K = copy.deepcopy(cmp)
del cmp
cmp = [[ sorted(s, reverse=True) for s in l]for l in sim1list]
sim1list = copy.deepcopy(cmp)
del cmp
cmp = [[ sorted(s, reverse=True) for s in l]for l in sim2list]
sim2list = copy.deepcopy(cmp)
del cmp
#Agoniste de frontière primaire
boundary1 = [[] for _ in range(betti0)]
for b in range(betti0):
    if len(sim1list[b]) > 0:
        M = np.ones((len(K[b]),len(sim1list[b])),dtype=np.int8)
        for i in range(len(K[b])):
            for j in range(len(sim1list[b])):
                if len(K[b]) > 1:
                    if sim1list[b][j][0] == K[b][i]:
                        M[i][j] = 1
                    else:
                        if sim1list[b][j][1] == K[b][i]:
                            M[i][j] = -1
                        else:
                            M[i][j] = 0
                else:
                    if sim1list[b][j][0] == K[b]:
                        M[i][j] = 1
                    else:
                        if sim1list[b][j][1] == K[b]:
                            M[i][j] = -1
                        else:
                            M[i][j] = 0
        boundary1[b] = copy.deepcopy(M)
boundary2 = [[] for _ in range(betti0)]
for b in range(betti0):
    if len(sim2list[b]) > 0:
        M = np.ones((len(sim1list[b]),len(sim2list[b])),dtype=np.int8)
        for i in range(len(sim1list[b])):
            for j in range(len(sim2list[b])):
                if [sim2list[b][j][1],sim2list[b][j][2]] == sim1list[b][i]:
                    M[i][j] = 1
                else:
                    if [sim2list[b][j][0],sim2list[b][j][2]] == sim1list[b][i]:
                        M[i][j] = -1
                    else:
                        if [sim2list[b][j][0],sim2list[b][j][1]] == sim1list[b][i]:
                            M[i][j] = 1
                        else:
                            M[i][j] = 0
        boundary2[b] = copy.deepcopy(M)

Comptez le nombre de trous

Le nombre de trous dans les données est dans la zone indiquée dans la figure suivante.

L'image de la cartographie de frontière secondaire $ \ partial_2 $ et le noyau de la cartographie de frontière primaire $ \ partial_1 $ sont chacun dans un groupe convertible, et il semble que le nombre de trous puisse être obtenu à partir de la dimension du groupe de quotient. Je ne connais pas les détails exacts, donc la formule que j'ai calculée cette fois est indiquée ci-dessous.

B_1 = dim(Ker(\partial_1)) - dim(Im(\partial_2))\\
B_1 = (dim(C_1) - rank(\partial_1)) - rank(\partial_2)

$ dim (C_1) $ est le nombre de 1 unité contenue dans la concaténation

#Variable qui compte le nombre de vétérinaires primaires
betti1 = 0
for b in range(betti0):
    rank1 = 0
    rank2 = 0
    #Confirmation de la présence ou de l'absence d'agonistes de frontière primaires
    if len(boundary1[b]) != 0:
        rank1 = GPU_rank(boundary1[b])
    if len(boundary2[b]) != 0:
        rank2 = GPU_rank(boundary2[b])
    betti1 += len(sim1list[b]) - rank1 - rank2
    print("Numéro concaténé",b,"  ",len(sim1list[b]),"-",rank1,"-",rank2,"=",len(sim1list[b]) - rank1 - rank2)
print("Le nombre de vétérinaires primaires est",betti1,"Devient")

Résultat du calcul

Lorsque vous visualisez les données On peut voir qu'il y a cinq trous blancs dans la figure. De plus, avec ces données, le résultat de l'exécution du programme de calcul d'homologie

Numéro concaténé 0414- 99 - 310 = 5
Le nombre de vétérinaires primaires sera de 5

Vous pouvez voir que le nombre de trous a été compté.

Visualisation de données

Le traitement a été utilisé pour la visualisation des données. Le programme est ici.

String lin;
String linR;
int ln;
static final int img_size = 800;
float r;
String lines[];
String linesR[];
void setup() {
  ln = 0;
  //r = 0.1;
  lines = loadStrings("pointcloud.csv");
  linesR = loadStrings("r.csv");
  linR = linesR[0];//Stocker la ligne ln
  String [] coR = split(linR,',');
  if (coR.length == 1){
    r = float(coR[0]);
  }
//Lire le fichier txt
  background(#FFFFFF);
  size(1000,1000);
  println(r);
  for(ln = 0;ln < lines.length;ln++){
    lin = lines[ln];//Stocker la ligne ln
    String [] co = split(lin,',');

    if (co.length == 2){
      float x = float(co[0]);
      float y = float(co[1]);
      fill(#000000);
      noStroke();
      ellipse(x*img_size + 100,
      1000 -(y*img_size + 100),
      2*r*img_size,
      2*r*img_size);
    }
  }
}

void draw() {
  
}

Enregistrez les données de coordonnées sous pointcloud.csv et la valeur de rayon sous r.csv.

La page que j'ai utilisée comme référence

https://cookie-box.hatenablog.com/entry/2016/12/04/132021 http://peng225.hatenablog.com/entry/2017/02/05/235951 http://peng225.hatenablog.com/entry/2017/02/18/133838

finalement

Le temps de calcul étant devenu long, comme les pièces qui composent le complexe, j'aimerais savoir comment le raccourcir. L'homologie persistante est connue comme une analyse de données de phase. À l'origine, j'avais l'intention de le faire, mais je ne pouvais pas écrire de programme concret, alors j'ai écrit un programme d'homologie. Si vous lisez cet article, je vous serais reconnaissant si vous pouviez m'apprendre comment accélérer le calcul de rang des complexes alpha et des matrices, et calculer l'homologie persistante.