[Homologie] Zählen Sie mit Python die Anzahl der Löcher in den Daten

Jupyter Notebook ist ein praktisches Tool, mit dem Sie die Berechnungsergebnisse unterwegs beim Schreiben eines Programms überprüfen können. Ein weiterer Vorteil der Verwendung von Google Colaboratory besteht darin, dass Sie unabhängig von Ihrer Umgebung rechnen können.

Eine der Methoden der Datenanalyse, die heutzutage ein heißes Thema geworden ist, ist die schrittweise Datenanalyse. Dies ist eine Analysemethode, die Daten als Formular erfasst und zählt, wie viele Features wie Löcher darin enthalten sind. Wie der Name schon sagt, werden wir zunächst die Topologie diskutieren.

Was ist Topologie?

In Wikipedia ist "Topologie eine Eigenschaft (Phaseneigenschaft oder Phase"), die beibehalten wird, selbst wenn sie in irgendeiner Form (Form oder "Raum") kontinuierlich verformt (gedehnt oder gebogen, aber nicht ausgeschnitten oder eingefügt) wird. Es konzentriert sich auf (invariant). Es wird gesagt.

Viele von Ihnen haben vielleicht gehört, dass "Kaffeetassen" und "Donuts", die häufig zur Erklärung der Topologie verwendet werden, aus topologischer Sicht identisch sind.

Zitiert aus der Wikipedia-Phasengeometrie

Aus dieser Abbildung können Sie erkennen, dass die "Kaffeetasse" und der "Donut" durch Biegen und Strecken ohne Schneiden und Einfügen verformt werden können. Dies bedeutet, dass die Anzahl der Löcher, die eine Phaseninvariante ist, eines gemeinsam hat. In diesem Artikel haben wir die Homologie berechnet, die die Anzahl der Löcher in einer zweidimensionalen Figur zählt.

Verfahren zur Berechnung der Homologie

Die spezifische Berechnung der Homologie ist Eingabe: Datenkoordinaten, Datenradius Ausgabe: Anzahl der Löcher in den Daten Ich werde es als tun.

Daten eingeben

Geben Sie die Daten ein, die Sie berechnen möchten. Diesmal wurde die Berechnung mit zufällig erstellten Daten durchgeführt.

import os
import random as rd
import numpy as np
import copy
#Punktwolkendatenerstellung
#Datendimension
dim = 2
#Die Anzahl der Daten
data_size = 500
#Komplexer Radius
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)

Um die Anzahl der Löcher zu berechnen, wird als nächstes der verbundene Körper berechnet, in dem die Daten verbunden sind.

Berechnung der Verkettung

Da die Daten nur zum Sitzen bestimmt sind, zeigt die folgende Abbildung die Daten mit einem Radius.

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
#Erstellen Sie eine benachbarte Matrix
pathM = make_path_matrix(data,sim_r)
#Deta Ankunftsflagge nicht erreicht: 0 erreicht: 1
flaglist = [0 for i in range(data_size)]
#Eine Funktion, die rekursiv eine benachbarte Matrix durchsucht
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)
#Verbundener Körper
K = []
for i in range(data_size):
    #Speichern Sie eines der Konjugate
    comp = [];
    #Speichern Sie die Flaggenliste vor der Erkundung
    flaglistsave = copy.deepcopy(flaglist);
    #Überprüfen Sie, ob es gesucht wurde
    if flaglist[i] == 0:
        #Suchen, wenn nicht gesucht
        visit(i,pathM);
        #Überprüfen Sie die Liste der Suchflaggen
        for flagindex in range(data_size):
            #Bestätigen Sie die gesuchten und neu gesuchten Punkte
            if flaglist[flagindex] == 1 and flaglistsave[flagindex] == 0:
                #Zur Verkettung hinzufügen
                comp.append(flagindex);
        #Zum Satz von Verkettungen hinzugefügt
        K.append(comp);

Indem Sie den Daten einen Radius geben, können Sie sehen, dass die Kreise in verbundene Verbindungen unterteilt werden können. Als nächstes untersuchen wir die kleinste Einheit, die im Konjugat enthalten ist.

Eine einzelne Einheit erstellen

Eine einzelne Einheit ist die kleinste Einheit, aus der ein verbundener Körper besteht, und eine einzelne Einheit ist ein Punkt, eine Linie oder eine Oberfläche.

Wenn es nur Punkte gibt, heißt es 0 Einheit, eine Linie heißt 1 Einheit und eine dreieckige Fläche heißt 2 Einheit. Dieses Mal ist es auf 2 Einheiten begrenzt, um die Löcher in 2 Dimensionen zu zählen, aber in 3 Dimensionen wird das Tetraeder zu 3 Einheiten und ist in höheren Dimensionen definiert.

Bei der Konfiguration von 2 Einheiten, wie in der folgenden Abbildung gezeigt, gibt es Fälle, in denen jedes Stück des Dreiecks 1 Einheit darstellt, jedoch nicht 2 Einheiten. ~~ Es scheint, dass Sie die Alpha-Verbindung verwenden können, die eine schnelle Berechnungsgeschwindigkeit hat, aber ich weiß es nicht, also habe ich sie mit der Prüfverbindung ~~ berechnet

#2 Überprüfen Sie die Konfiguration einer einzelnen Einheit und verbessern Sie sie
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)
    #Überprüfen Sie, ob jede der drei Seiten eine Einheit darstellt
    sim1flag = da <(2*r) and db<(2*r) and dc<(2*r)
    #Reihers Beamter
    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
    #Bereichsbestätigung
    Sflag = S<MaxS
    #Externes Kreisurteil
    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
#Anzahl der Anschlüsse
betti0 = len(K)
#1 Einzelliste
sim1list = [[] for i in range(betti0)]
#2 Einzelliste
sim2list = [[] for i in range(betti0)]

#Extrahieren Sie eine Verkettung mit nur einer Einheit und einer Einheit
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]])
#1 Einheit aus dem Konjugat extrahieren
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]])
#Extrahiere 2 Einheiten aus dem Konjugat
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]])

Durch Zählen der einzelnen Einheiten in der Verkettung wird eine Menge erstellt, die eine endliche Anzahl einzelner Einheiten enthält. Als nächstes wird eine Abbildung zum Konvertieren einer einzelnen Einheit in eine einzelne Einheit einer anderen Dimension und eine Grenzzuordnung beschrieben.

Grenzzuordnung

Eine Grenzzuordnung ist eine Zuordnung, die eine einzelne Einheit in eine kleine einzelne Einheit einer Dimension umwandelt.

Wie in der Abbildung gezeigt, handelt es sich um eine Abbildung, die von 2 Einheiten in die alternierende Summe von 1 Einheit umgerechnet wird. Die konkrete Zuordnung ist wie folgt

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

Die von der ursprünglichen Einheit entfernte Zahl ist eine Zuordnung, bei der die Summe durch Ändern des Vorzeichens der Zahl von vorne durch gerade oder ungerade Zahlen ermittelt wird. Die Grenzabbildung, die sich von 2 Einheiten zu 1 Einheit bewegt, die im Konjugat enthalten ist, wird als quadratische Grenzabbildung bezeichnet, und die Grenzabbildung, die von 1 Einheit zu 0 Einheiten übertragen wird, wird als primäre Grenzabbildung bezeichnet. Die erhaltene Grenzabbildung zweiter Ordnung und die Grenzabbildung erster Ordnung werden als Ausdrucksmatrix dargestellt. Als nächstes wird die Berechnung zum Zählen der Anzahl von Löchern von diesem Grenzagonisten beschrieben.

#Sortierverkettung, 1 Einheit, 2 Einheit
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
#Primärer Grenzagonist
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)

Zählen Sie die Anzahl der Löcher

Die Anzahl der Löcher in den Daten liegt in dem in der folgenden Abbildung gezeigten Bereich.

Das Bild der sekundären Grenzabbildung $ \ Partial_2 $ und der Kern der primären Grenzabbildung $ \ Partial_1 $ befinden sich jeweils in einer konvertierbaren Gruppe, und es scheint, dass die Anzahl der Löcher aus der Dimension der Quotientengruppe erhalten werden kann. Ich kenne die genauen Details nicht, daher wird die Formel, die ich dieses Mal berechnet habe, unten gezeigt.

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

$ dim (C_1) $ ist die Nummer 1 in der Verkettung

#Variable, die die Anzahl der primären Wicken zählt
betti1 = 0
for b in range(betti0):
    rank1 = 0
    rank2 = 0
    #Bestätigung der Anwesenheit oder Abwesenheit von primären Grenzagonisten
    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("Verkettete Nummer",b,"  ",len(sim1list[b]),"-",rank1,"-",rank2,"=",len(sim1list[b]) - rank1 - rank2)
print("Die Anzahl der primären Wicken beträgt",betti1,"Wird")

Berechnungsergebnis

Wenn Sie die Daten visualisieren Es ist zu sehen, dass die Figur fünf weiße Löcher enthält. Mit diesen Daten auch das Ergebnis der Ausführung des Homologie-Berechnungsprogramms

Verkettete Nummer 0 414- 99 - 310 = 5
Die Anzahl der primären Wicken beträgt 5

Sie können sehen, dass die Anzahl der Löcher gezählt wurde.

Datenvisualisierung

Die Verarbeitung wurde zur Datenvisualisierung verwendet. Das Programm ist da.

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];//Speichern Sie die In-Zeile
  String [] coR = split(linR,',');
  if (coR.length == 1){
    r = float(coR[0]);
  }
//TXT-Datei lesen
  background(#FFFFFF);
  size(1000,1000);
  println(r);
  for(ln = 0;ln < lines.length;ln++){
    lin = lines[ln];//Speichern Sie die In-Zeile
    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() {
  
}

Speichern Sie die Koordinatendaten als pointcloud.csv und den Radiuswert als r.csv.

Die Seite, die ich als Referenz verwendet habe

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

Schließlich

Da die Berechnungszeit lang geworden ist, wie z. B. die Teile, aus denen der Komplex besteht, möchte ich herausfinden, wie er verkürzt werden kann. Die persistente Homologie ist berühmt für die Phasendatenanalyse. Ursprünglich wollte ich das tun, aber ich konnte kein konkretes Programm schreiben, also schrieb ich ein Homologieprogramm. Wenn Sie diesen Artikel lesen, würde ich es begrüßen, wenn Sie mir beibringen könnten, wie man die Rangberechnung von Alpha-Komplexen und Matrizen beschleunigt und die persistente Homologie berechnet.