A Graph a Day

Vor einiger Zeit habe ich @randomGraphs geschrieben: Ein Twitterbot, der einen Zufallsgraphen pro Tag tweetet.

Die meisten Graphtypen, die er darstellen kann stammen aus der NetworkX Bibliothek oder sind reale Netzwerke. Ein paar Proximity Graphs habe ich selbst geschrieben. Die Darstellung und gegebenenfalls das Layout übernimmt Cytoscape oder graph-tool (dessen Autor diesem Bot folgt).

Bei diesem Projekt habe ich exzessiv Gebrauch von Pythons Decorator und Introspection gemacht, sodass man, um einen neuen Graphtyp einzuführen nur eine Methode schreiben muss, die eine Graph-Datenstruktur zurück gibt. Einstellungen, welche Darstellungen erlaubt sind, werden per decorator getätigt und alle Methoden werden per Introspection automatisch zum Pool hinzugefügt, aus dem der Zufallsgenerator zieht.

Eine typische Methode sieht etwa so aus.

@synonym("Barabasi Albert")
@synonym("preferential attachment")
@style(styles_all)
@layout(["kamada-kawai", "force-directed", "sfdp", "fruchterman_reingold", "arf", "radial_tree"])
def generateBarabasiAlbert(self, N=None, m=None, **kwargs):
    if N is None: N = random.randint(4, 400)
    if m is None: m = random.randint(1, 5)

    G = gen.barabasi_albert_graph(N, m)  # gen is networkx Generator
    details = dict(name="Barabási-Albert Graph", N=N, m=m, seed=self.seed,
                   template="{name}, N = {N}, m = {m}")

    return G, details

Und liefert für \(N=226, m=1\) und das radial_tree Layout beispielsweise diesen Graph. Die Größe der Knoten wird hier von der Betweenness Centrality bestimmt.

Graph

Die @synonym Decorators ermöglichen die zweite Funktion des Bots, denn er tweetet nicht nur einmal am Tag einen zufälligen Graphen, sondern reagiert auch auf Mentions. Falls in der Mention der Name der Methode oder eines der per @synonym registrierten Worte auftaucht, antwortet er mit einem Bild des entsprechenden Graphen. Dank fuzzywuzzy ist es sogar resistent gegen Tippfehler.

Twitter unterstützt leider keine Vektorgrafiken und wandelt Bilder gerne in stark komprimierte .jpg, was gerade bei diesen Graphen zu störenden Artefakten führt. Dagegen hilft es, wenn ich einen Rand aus transparenten Pixeln dem Bild hinzufüge. Das führt dazu, dass Twitter .jpg nicht als geeignetes Format ansieht und die Bilder im verlustfreien .png ausliefert.

convert -alpha on -channel RGBA -bordercolor "rgba(0,0,0,0)" -border "1x1" input.png output.png

Graph

Der komplette Quellcode ist auf Github.

SimulatedSort

Simulated Annealing ist eine Optimierungsmethode, die von natürlichen Kristallisationsprozessen inspiriert ist. Man startet in der Schmelze bei hohen Temperaturen und lässt es dann abkühlen, sodass die Atome sich in einem Zustand minimaler Energie anordnen, dem Kristallgitter. Wenn man also für ein Optimierungsproblem die zu optimierende Größe als Energie ansieht, und man eine Lösung durch eine kleine Änderung in eine andere Lösung verwandeln kann, kann man mit dieser Methode eine Näherung für das Optimum finden.

Wenn wir also eine Sequenz \(S\) von \(N\) Zahlen sortieren wollen, können wir die Summe der Differenzen zwischen benachbarten Zahlen als Energie betrachten, denn die ist minimal in einer sortierten Liste.

\begin{equation} \mathcal{H} = \sum_{i=1}^{N-1} \left| S_i - S_{i+1} \right| \end{equation}

Um eine Lösung in eine andere zu verwandeln, reicht es zwei Elemente der Sequenz zu tauschen.

Der Kern von Simulated Annealing ist der Metropolis Algorithmus.

  1. Starte bei einer hohen Temperatur \(T\).
  2. Berechne die Energie \(\mathcal{H}(S)\) der aktuellen Konfiguration \(S\).
  3. Erzeuge eine neue Konfiguration \(R\) durch eine kleine Änderunge von \(S\).
  4. Akzeptiere \(R\) mit der Wahscheinlichkeit
    $$p_\mathrm{acc} = \min\left[1 ,\exp(-(\mathcal{H}(R) - \mathcal{H}(S))/T) \right],$$
    sodass eine „sortiertere“ Sequenz immer akzeptiert wird und eine „unsortiertere“ vor allem bei hohen Temperaturen. Wenn \(R\) akzeptiert wird, gilt \(S:=R\), ansonsten wir die alte Konfiguration \(S\) weiter benutzt.
  5. Reduziere die Temperatur (beispielsweise durch Multiplikation mit einer Zahl etwas kleiner als 1) und breche ab, wenn die Zieltemperatur erreicht ist. Ansonsten beginne wieder bei Punkt 2.

Genug der Theorie: In einem Gist auf GitHub präsentiere ich ein schnell terminierendes Sortierprogramm, das zwar nicht immer eine sortierte Liste findet, aber zumindest eine Näherung! Es ist also Bogosort in mehr als nur einer Hinsicht überlegen!

Wer braucht da noch \(\mathcal{O}(N \log(N))\) Sortier-Algorithmen?!

DGLshow

Nachdem ich so vielen Differenzialgleichungssystemen [1, 2, 3, 4] begegnet bin, die sich nicht analytisch lösen lassen, habe ich mir ein Programm zur numerischen Lösung und Visualisierung derselben geschrieben.

Die grundlegende Idee zur numerischen Lösung von Differentialgleichungen ist es, die Zeit in diskreten Schritten \(\tau\) vergehen zu lassen. Nach jedem Schritt wird der Zustand so geändert, als ob sich während des Zeitschrittes nichts geändert hätte und die „Kräfte“ werden entsprechend der Bewegungsgleichungen neu berechnet. Für infinitesimal kleine \(\tau \to \mathrm{d}t\) ist diese Methode schließlich exakt.

Im einfachsten Fall, dem Euler Verfahren, sähe das für ein einfaches Fadenpendel nach \(k\) Zeitschritten so aus

\begin{align*} \ddot\vartheta_{k+1} &= - mgl \sin(\vartheta_k)\\ \dot\vartheta_{k+1} &= \tau \ddot\vartheta_{k} + \dot\vartheta_{k}\\ \vartheta_{k+1} &= \tau \dot\vartheta_{k} + \vartheta_{k} \end{align*}

Unglücklicherweise hat dieses Verfahren ernsthafte Probleme mit der Energieerhaltung und braucht sehr kleine \(\tau\) für brauchbare Ergebnisse. Es gibt deutlich ausgefeiltere Methoden, wie den klassischen Runge-Kutta Algorithmus. Es gibt Methoden, den Zeitschritt \(\tau\) während der Simulation adaptiv anzupassen, um nur wenig Rechenaufwand in den wenig fehleranfälligen Phasen zu verbringen. Es gibt spezialisierte Methoden, die sehr gut für bestimmte Bewegungsgleichungen funktionieren, wie Velocity-Verlet, der oft für Molekulardynamiksimulationen eingesetzt wird.

Chaotische Systeme haben in der Regel etwas kompliziertere Bewegungsgleichungen. Das oben abgebildete Doppelpendel etwa wird, wie ich in einem anderen Post beschrieben habe durch folgendes Ungetüm beschrieben.

\begin{align*} \ddot\theta_1 &= \frac{m_2 \cos(\theta_1 - \theta_2) (l_1 \sin(\theta_1 - \theta_2) \dot\theta_1^2 - g \sin(\theta_2)) + m_2 l_2 \sin(\theta_1 - \theta_2) \dot\theta_2^2 + (m_1 + m_2) g \sin(\theta_1)}{m_2 l_1 \cos^2(\theta_1 - \theta_2) - (m_1+m_2) l_1} \\ \ddot\theta_2 &= \frac{m_2 l_2 \cos(\theta_1 - \theta_2) \sin(\theta_1 - \theta_2) \dot\theta_2^2 + (m_1+m_2) l_1 \sin(\theta_1 - \theta_2) \dot\theta_1^2 + (m_1+m_2) g \cos(\theta_1 - \theta_2) \sin(\theta_1) - (m_1+m_2) g \sin(\theta_2)}{(m_1+m_2) l_2 - m_2 l_2 \cos^2(\theta_1 - \theta_2)} \end{align*}

Anfangs empfiehlt es sich also etwas einfacheres und vertrauteres zu lösen, wie den Lorenz-Attraktor

\begin{align*} \dot{X} &= a(Y - X) \\ \dot{Y} &= X(b - Z) - Y \\ \dot{Z} &= XY - cZ \\ \end{align*}

Oder das Dreikörperproblem

\begin{align*} \ddot{\vec{x}_1} &= -\frac{Gm_2}{\left(x_1 - x_2\right)^3} (\vec{x}_1 - \vec{x}_2) - \frac{Gm_3}{\left(x_1 - x_3\right)^3} (\vec{x}_1 - \vec{x}_3)\\ \ddot{\vec{x}_2} &= -\frac{Gm_1}{\left(x_2 - x_1\right)^3} (\vec{x}_2 - \vec{x}_1) - \frac{Gm_3}{\left(x_2 - x_3\right)^3} (\vec{x}_2 - \vec{x}_3)\\ \ddot{\vec{x}_3} &= -\frac{Gm_1}{\left(x_3 - x_1\right)^3} (\vec{x}_3 - \vec{x}_1) - \frac{Gm_2}{\left(x_3 - x_2\right)^3} (\vec{x}_3 - \vec{x}_2)\\ \end{align*}

Da man das 3-Körperproblem trivial auf ein \(N\)-Körperproblem erweitern kann, habe ich hier ein „Sonnensystem“ bzw. Bohrsches „Atom“-modell simuliert.

Sonnensystem

Um die obigen (bewegten) Bilder zu erzeugen und um ein bewegtes Doppelpendel für meinen Schreibtisch zu haben, — wennauch nur auf einem Bildschirm — habe ich in C++ einen adaptiven Runge-Kutta-4 Löser geschrieben, der mit den Qt Zeichenprimitiven animiert wird.

Auch wenn der Code nicht sehr aufgeräumt ist und Startwerte im Quellcode angepasst werden müssen, sind die Quellen auf GitHub: github.com/surt91/DGLshow.