Seit Version 2.4 unterstützt Python das Konzept der verschachtelten Funktionen. Das ermöglicht elegante Ausdrucksweisen bei der Programmierung, etwa beim Berechnen der konvexen Hülle für eine Punktemenge.

Beginnen Sie mit dem Programmieren, gehören Prozeduren oder Funktionen als in sich abgeschlossener Programmblock zu den festen Bestandteilen jedes Lernmaterials. Zum vermittelten Grundwissen zählt auch, dass Sie jede Funktion über einen eindeutigen, meist frei wählbaren Namen ansprechen. Mit Parametern verfeinern Sie den Aufruf und bringen zusätzliche Flexibilität in die Ausführung des Codeblocks. Eine Funktion verlassen Sie explizit über die Anweisung return, häufig ergänzt um einen Rückgabewert.

Listing 1 enthält in den ersten vier Zeilen eine einfache Funktion summe() mit zwei Parametern. Sie addiert zwei übergebene Werte und liefert das Ergebnis an den Aufrufer zurück. In den letzten beiden Zeilen sehen Sie zwei Aufrufe der Funktion samt Ausgabe der Berechnung.

def summe(wert1, wert2):
  """wert1 und wert2 addieren"""
  ergebnis = wert1 + wert2
  return ergebnis
print("die Summe aus 5 und 10 ist", summe(5, 10))
print("die Summe aus 15 und 10 ist", summe(15, 10))

Nachdem dem Verlassen der Funktion mittels return räumt die interne Müllabfuhr von Python jedes Mal die Umgebung dafür – Stack und Variablen – wieder auf. Alle bisherigen Zwischenstände in der Funktion gehen damit verloren. Rufen Sie sie später ein weiteres Mal auf, erstellt Python für die Funktion eine neue, leere Umgebung.

Verschachtelte Funktionen

Innere oder verschachtelte Funktionen [1] kommen ins Spiel, wenn Sie den Abarbeitungszustand einer Funktion bewahren möchten, um später wieder darauf aufzubauen. Dazu deklarieren Sie die entsprechende(n) Funktion(en) vollständig in einer anderen Funktion. Sie verwenden dazu dieselbe Schreibweise wie bei regulären Funktionen. Auch der Aufruf der Funktion unterscheidet sich nicht vom bisher Gewohnten.

Der Unterschied liegt im Gültigkeitsbereich: Innere Funktionen sind nur innerhalb der umgebenden Funktion gültig und lassen sich nur dort aufrufen. Sie sind also nur lokal zur Funktion sichtbar und nicht global im gesamten Python-Programm. Die innere(n) Funktion(en) erhalten vollen Zugriff auf alle Variablen und Datenstrukturen der umgebenden Funktion.

Verwendung findet dieses Sprachelement beim Abschotten nach außen sowie bei Closures [2] und Decorators [3]. Letzteres ist eine andere (vereinfachte) Schreibweise von Funktionen als Parameter. Die beiden Web-Frameworks Django [4] und Flask [5] nutzen dieses Feature intern sehr intensiv.

Closures

Das Sprachelement Closure gehört zur funktionalen Programmierung [6]. Es beschreibt das Merken des Zustands eines Funktionsaufrufs. Beim Verlassen einer Funktion merkt sich der Interpreter von Python die Werte der genutzten Variablen in der Funktion. Rufen Sie die Funktion erneut auf, erinnert sich Python an diesen Zustand und setzt die Aufrufe mit den vorher gemerkten Werten fort. Das klingt im ersten Moment seltsam, lässt sich aber oft sehr praktisch einsetzen, wie ein Automat, den Sie in einen bestimmten Zustand versetzen.

Listing 2 zeigt eine Closure, wobei zum gemerkten Ausgangswert in der dritten Zeile jeweils 5 addiert wird. Der Aufruf der inneren Funktion erhoehung() samt Übergabe der Parameter erfolgt in Zeile 5.

def summe (w1):        # Start äußere Funktion
  def erhoehung (w2):  # Start innere Funktion
    ergebnis = w1 + w2 # neue Werte
    return ergebnis    # Rückgabe Resultat
  return erhoehung     # Ende der Closure
# Berechnungsbeispiel
a = summe(3)
b = summe(5)
print(a(5))            # gibt 3 + 5 = 8 aus
print(b(5))            # gibt 5 + 5 = 10 aus

Anhand eines ausführlichen Praxisbeispiels verdeutlichen wir nachfolgend, wie Sie innere Funktionen verwenden. Dazu betrachten wir die Bestimmung der konvexen Hülle [7] für eine Menge von Punkten.

Konvexe Hülle

Bei einer konvexen Hülle [8] für eine Menge von Punkten handelt es sich um einen Begriff beziehungsweise ein Verfahren aus der angewandten Geometrie. Eine konvexe Hülle [9] bezeichnet eine Menge von Linien, die alle äußeren Punkte der Punktemenge miteinander verbinden. Bildlich gesprochen, spannen Sie ein Gummiband um eine Menge von Punkten, wie Abbildung 1 veranschaulicht.

Abbildung 1: Konvexe Hülle von Punkten mit Markern.

Auf den Alltag angewendet, könnte eine passende Aufgabe lauten: Berechnen Sie die Länge des Zauns, den Sie benötigen, um ein Grundstück vollständig einzuzäunen. Alternativ könnte es darum gehen, die minimal notwendige Verpackung für ein nicht rechteckiges Objekt herauszubekommen. Weitere Beispielanwendungen sind das Ermitteln der Distanz, die ein Hütehund zurücklegt, wenn er einmal eine Schafherde umrundet, oder der Mindestlänge eines Seils zum Umspannen eines bestimmten Gegenstands wie eines Felsbrocken.

Die ermittelten Punkte bilden dann die Basis für weitere Schritte, beispielsweise die Berechnung der Gesamtfläche zwischen den Punkten (Summe der Teilflächen auf der Basis von Dreiecken) oder zur Hindernis- und Kollisionserkennung bei Fahrzeugen (Berechnung der Distanz von einem bestimmten Punkt zu allen anderen Punkten).

Implementierungen

Um eine konvexe Hülle zu ermitteln, gibt es mehrere Verfahren [10] mit neckischen Namen wie Graham Scan, Gift Wrapping, Monotone Chain, Jarvis und Quick Hull. Sie unterscheiden sich hinsichtlich ihrer Arbeitsweise – beispielsweise iterativ von links nach rechts oder nach dem Teile-und-herrsche-Prinzip – und der sich daraus ergebenden Laufzeit.

In den drei populären Python-Bibliotheken SciPy [11], Scikit-geometry  [12] und GeoPandas [13] finden Sie dazu passende Implementierungen. Für diesen Artikel vergleichen wir die Variante aus SciPy mit unserer eigenen Implementierung von Quick Hull. Quick Hull ähnelt stark dem Sortierverfahren Quicksort [14], da es ebenfalls das Prinzip Teile-und-herrsche (“divide and conquer”) verwendet. Dabei zerlegt man ein größeres Problem in kleinere Probleme, die man daraufhin zuerst löst. Die Teilergebnisse werden schließlich Schritt für Schritt zu einem Gesamtergebnis zusammengefasst.

Angewendet auf die vorliegende Punktemenge heißt das, dass eine Gerade sie zunächst in zwei Portionen (“divide”) aufteilt. Die Gerade verläuft dabei durch den kleinsten, sich am weitesten links unten befindenden sowie den größten, sich am weitesten rechts oben liegenden Punkt.

Nun bestimmen Sie die konvexe Hülle für alle Punkte links und rechts der Geraden (“conquer”), wobei Sie zuerst die linke Seite und danach die rechte Seite berechnen. Am Schluss vereinigen Sie die beiden separat berechneten Punktemengen wieder miteinander. Jeder Punkt aus den Teilergebnissen kommt dabei nur einmal im Endergebnis vor.

Dazu kommen in unserer Implementierung fünf Funktionen zum Einsatz: findDistance() (Listing 3) berechnet den senkrechten Abstand zwischen einem Punkt testpoint und der Linie, die durch die Punkte p1 und p2 gebildet wird.

Die Funktion findPosition() (Listing 4) bestimmt, ob sich ein Punkt testpoint(x,y) links oder rechts der Geraden befindet, die durch p1(x,y) und p2(x,y) verläuft.

In Listing 5 ermittelt isPointInTriangle(), ob sich ein Punkt p(x,y) innerhalb eines Dreiecks abc befindet, oder nicht. Zudem berechnet quickHull() (Listing 6) die konvexe Hülle über eine Punktemenge. Bei calculateHull() (Listing 7) handelt es sich um die innere Funktion in quickHull, die die konvexe Hülle über die Punktemenge tatsächlich liefert.

Das Trio findDistance(), findPosition() und isPointInTriangle() fungiert als Hilfsfunktionen für die Funktion quickHull(). In Letzterer befindet sich die innere Funktion calculateHull().

Die Funktion findDistance() ermittelt den geometrischen Abstand eines Punkts testpoint zu einer Geraden, auf der die beiden Punkte p1 und p2 liegen. Der Abstand ist 0, sofern sich der Punkt direkt auf der Geraden befindet, ansonsten hat er einen positiven Wert. Wo der Punkt testpoint liegt, klärt findPosition(). Ein Rückgabewert von 1 steht für links der Geraden, 0 für auf der Geraden und -1 für rechts der Geraden.

def findDistance(p1, p2, testpoint):
  """
  Berechne den senkrechten Abstand zwischen einem
  Punkt testpoint und der Geraden, die durch die
  Punkte p1 und p2 gebildet wird.
  """
  numerator = abs((p2[1]-p1[1])*testpoint[0] - (p2[0]-p1[0])*testpoint[1] + p2[0]*p1[1] - (p2[1]*p1[0])
  denominator = math.sqrt((p2[1]-p1[1]) * (p2[1]-p1[1]) + (p2[0]-p1[0])*(p2[0]-p1[0]))
  try:
    result = numerator / denominator
  except:
    result = 0
  return result

def findPosition(p1, p2, testpoint):
  """
  Berechne, ob ein Punkt testpoint auf der rechten
  oder linken Seite einer Geraden liegt, die durch
  die Punkte p1 und p2 gebildet wird.
  """
  result = (p2[0] - p1[0]) * (testpoint[1] - p1[1]) - (p2[1] - p1[1]) * (testpoint[0] - p1[0])
  # Auswertung
  if result == 0:
    return 0
  if result < 0:
    return -1
  return 1

Die Funktion isPointInTriangle() berechnet die baryzentrischen Koordinaten [15] des Punkts testpoint mit Bezug auf das Dreieck. Daraus ergibt sich, ob sich der Punkt testpoint innerhalb des Dreiecks befindet, oder nicht. Innerhalb des Dreiecks liegende Punkte werden später aus der Berechnung der konvexen Hülle herausgenommen. Damit verringert sich die insgesamt erforderliche Rechenzeit.

def isPointInTriangle(p1, p2, p3, testpoint):
  """
  Bestimme, ob der Punkt testpoint innerhalb des
  Dreiecks (p1,p2,p3) liegt
  """
  # definiere den Standardrückgabewert: False
  result = False
   # Berechne die baryzentrischen Koordinaten des
   # Punkts testpoint mit Bezug auf das Dreieck
   denominator = ((p2[1] - p3[1]) * (p1[0] - p3[0]) + (p3[0] - p2[0]) * (p1[1] - p3[1]))
   if denominator > 0:
     a = ((p2[1] - p3[1]) * (testpoint[0] - p3[0]) + (p3[0] - p2[0]) * (testpoint[1] - p3[1])) / denominator
     b = ((p3[1] - p1[1]) * (testpoint[0] - p3[0]) + (p1[0] - p3[0]) * (testpoint[1] - p3[1])) / denominator
     c = 1 - a - b
     # Prüfe, ob alle baryzentrischen Koordinaten
     # nicht-negativ sind
     if a >= 0 and b >= 0 and c >= 0:
       result = True
  return result

Die Funktion quickHull() (Listing 6) berechnet die konvexe Hülle für die Punktemenge auf der Basis des Quick-Hull-Verfahrens. Die Grundlage dazu bilden die der Funktion als Parameter übergebenen Punkte. Sind es weniger als drei Punkte, braucht es keine weitere, aufwendige Berechnung – die konvexe Hülle liegt bereits vor. Aus den maximal zwei Punkten wird die sortierte Schnittmenge errechnet und als Ergebnis zurückgegeben.

Bei mehr als zwei Punkten ermitteln Sie daraus im ersten Schritt den kleinsten und größten Punkt. Durch die beiden Punkte legen Sie eine Gerade und teilen die verbleibenden Punkte anhand ihrer Position zur Geraden ein – entweder links oder rechts davon. Mithilfe der Funktion calculateHull() (Listing 6, Zeile 10) berechnen Sie zuerst die konvexe Hülle für alle links liegenden Punkte und danach für alle rechts liegenden Punkte. Beide Ergebnisse werden nachfolgend kombiniert, indem die kleinste Schnittmenge beider konvexen Hüllen ermittelt wird.

def quickHull(points = []):
  """
  berechne die konvexe Hülle auf der Basis des
  Quick Hull-Verfahrens
  """
  # definiere Basiswert: leer
  hull = []
  # prüfe, ob weniger als drei Punkte vorliegen
  if len(points) < 3:
    return sorted(set(points)) # nichts zu tun
  def calculateHull(p1, p2, pointsList):
    # siehe Listing 7
    return result
   # berechne den kleinsten und größten Punkt
   minimalPoint = min(points)
   maximalPoint = max(points)
   # bestimme die Punkte, die sich links oder
   # rechts der Geraden befinden
   leftPoints = []  # Annahme: leere Menge
   rightPoints = [] # Annahme: leere Menge
   for singlePoint in points:
     position = findPosition(minimalPoint, maximalPoint, singlePoint)
     if position > 0:
       leftPoints.append(singlePoint)
     if position < 0:
       rightPoints.append(singlePoint)
   # berechne konvexe Hülle für linke Punkte
   part1 = calculateHull(minimalPoint, maximalPoint, leftPoints)
   # berechne konvexe Hülle für rechte Punkte
   part2 = calculateHull(minimalPoint, maximalPoint, rightPoints)
   # Menge aus dem kleinsten und größten Punkt
   part3 = [minimalPoint, maximalPoint]
   # kleinste Schnittmenge aller Punkte
   hull = sorted(set(part1 + part2 + part3))
   return hull

Die innere Funktion calculateHull() ist nach einem ähnlichen Prinzip aufgebaut wie quickHull(). Zunächst wird geprüft, ob die Liste der Punkte leer ist. Falls nicht, ermittelt die Funktion findDistance() aus der Punktemenge den am weitesten entfernten Punkt bezogen auf die Gerade aus den beiden Punkten p1 und p2. Nun folgt der Test, welche der verbleibenden Punkte sich innerhalb des Dreiecks aus den drei Punkten p1, p2 und dem am weitesten entfernt liegenden Punkt befinden.

Zum Einsatz kommt dabei die bereits oben beschriebene Hilfsfunktion isPointInTriangle() (Listing 7, Zeile 14). Die errechnete Punktemenge wird wieder in zwei Teile zerlegt: bezüglich p1 und dem am weitesten entfernten Punkt sowie dem am weitesten entfernten Punkt und p2. In einem nachfolgenden Rekursionsschritt rufen Sie die Funktion calculateHull() über die beiden Teilmengen auf. Die sich daraus ergebenden Ergebnisse werden miteinander kombiniert, woraus sich am Ende die konvexe Hülle über die gesamte Punktemenge ergibt.

def calculateHull(p1, p2, pointsList):
  """
  berechne die konvexe Hülle für eine Punktemenge
  bzgl. einer Geraden
  """
  # prüfe auf eine leere Liste von Punkten
  if pointsList == []:
    return []
  # berechne am weitesten entfernten Punkt
  farthestPoint = max(pointsList, key=lambda p: findDistance(p1, p2, p))
  # berechne die übrigen Punkte
  nextPoints = []
  for testpoint in pointsList:
    if testpoint != farthestPoint:
      result = isPointInTriangle(p1, p2, farthestPoint, testpoint)
      if not result:
        nextPoints.append(testpoint)
  if not len(nextPoints):
    result = sorted([p1, p2, farthestPoint])
  else:
    # berechne die konvexe Hülle
    # - zwischen Punkt 1 und dem am
    #   weitesten entfernten Punkt
    part1 = calculateHull(p1, farthestPoint, nextPoints)
    # - zwischen am weitesten entfernten
    #   Punkt und Punkt 2
    part2 = calculateHull(farthestPoint, p2, nextPoints)
    # kombiniere beide Ergebnisse
    result = sorted(set(part1 + part2))
  return result

Abschließend fehlt noch das Hauptprogramm, in dem Sie zunächst die Punktemenge definieren (Listing 8). Hier sind das 500 zufällige Punkte mit ganzzahligen Werten (Zeile 5). Über sie berechnen Sie danach die konvexe Hülle mithilfe der zuvor besprochenen Funktion quickHull() aus Listing 6. Dazu dient der Funktionsaufruf in Zeile 8 von Listing 8. Anschließend geben Sie das Ergebnis aus (Zeile 9).

import math
import random
if __name__ == '__main__':
  # define 500 random data points
  points = [(int(random.random()*50), int(random.random()*50)) for value in range(500)]
  print(f"Points: {points}")
  hull = quickHull(points)
  print(f"Output: {hull}")

Rufen Sie Listing 8 auf, erhalten Sie eine Ausgabe wie in Listing 9. Das Ergebnis ist eine aufsteigend sortierte Punkteliste, wobei jeder Punkt aus einem Tupel aus zwei Ganzzahlen besteht.

$ python3 calculate-convex-hull.py
Output: [(0, 5), (0, 8), (0, 14), (0, 25), (0, 26), (0, 29), (0, 30), (0, 32), (0, 37), (0, 41), (1, 0), (1, 49), (3, 3), (3, 5), (3, 49), (4, 8), (5, 5), (5, 6), (6, 2), (6, 4), (6, 9), (6, 49), (7, 8), (8, 8), (8, 10), (8, 49), (10, 1), (11, 2), (12, 0), (13, 15), (14, 1), (15, 1), (16, 0), (17, 0), (17, 49), (19, 2), (20, 22), (21, 21), (23, 1), (24, 0), (24, 24), (26, 2), (31, 30), (32, 0), (34, 1), (35, 1), (35, 49), (39, 38), (39, 48), (40, 49), (41, 0), (41, 48), (41, 49), (43, 49), (45, 43), (46, 0), (47, 2), (47, 4), (47, 6), (47, 10), (47, 11), (47, 14), (47, 15), (48, 2), (48, 24), (48, 25), (48, 48), (48, 49), (49, 7), (49, 13), (49, 37), (49, 46), (49, 47)]

Bibliotheksfunktionen

Wie oben erwähnt benötigen Sie noch einen Vergleichswert. Dazu nutzen Sie die Klasse ConvexHull aus Pythons SciPy-Bibliothek. In den ersten Zeilen von Listing 10 importieren Sie zunächst die Klasse und dann die Bibliothek NumPy. Die verwenden Sie, um den Zufallszahlengenerator zu starten (Zeile 6) und damit 500 Paare von Zufallswerten zu erzeugen (Zeile 8).

In Zeile 10 berechnen Sie mithilfe der Klasse ConvexHull die konvexe Hülle über die 500 Punkte. Das Ergebnis – also die Punkte, die die konvexe Hülle bilden – geben Sie danach über die Methode simplices aus (Zeile 12). Die errechneten Punkte sehen Sie untereinander aufgeführt in Listing 11.

# ConvexHull aus SciPy
from scipy.spatial import ConvexHull
# NumPy
import numpy as np
# Zufallsgenerator starten
rng = np.random.default_rng()
# 500 Zufallswerte erzeugen
points = rng.random((500, 2))
# konvexe Hülle berechnen
hull = ConvexHull(points)
# berechnete Punkte ausgeben
print(hull.simplices)

$ python3 calculate-convex-hull-scipy.py
[[196  25]
 [196  69]
 [255  32]
 [255  23]
 [ 54  32]
 [ 54  80]
 [178  69]
 [178  80]
 [127 187]
 [228 187]
 [495  25]
 [495 127]
 [117 228]
 [399  23]
 [399 117]]

Messen und Vergleichen

Messen Sie die Ausführungszeit beider Implementierungen mit dem Unix-Werkzeug time, errechnet es stets zwischen 0.1 und 0.2 Sekunden. In der Tabelle “Laufzeit (Sekunden)” finden Sie weitere Vergleichsmessungen für 100, 250 und 1000 zufällige Werte. Daraus können Sie ablesen, dass die Laufzeit beider Verfahren bei einer geringen Werteanzahl nicht signifikant voneinander abweicht, später aber deutlich zunimmt.

Darüber hinaus haben wir bei unseren Tests beobachtet, dass Sie bereits ab 250 Datenpunkten deutlich mehr Zeit einplanen müssen – vor allem, wenn Sie zusätzlich Optimierungen wie Threads und Parallelisierung nutzen. Derzeit liegen der Grund für das Verhalten und die Stelle des Flaschenhalses noch im Dunkeln. Etwas Licht bringt das Werkzeug Scalene [16], das dabei hilft, die internen Aufrufe zu analysieren. Unter Debian GNU/Linux steht es über das Paket Python3-scalene [17] bereit. Abbildung 2 zeigt die Analyse von Listing 8, die die am häufigsten aufgerufenen Programmzeilen identifiziert.

Abbildung 2: Nach der Analyse von Scalene gehören in Listing 8 die Zeilen 29 und 30 zu den am häufigsten aufgerufenen.

Ausblick: Threads

Die bisherige Implementierung funktioniert und liefert in überschaubarer Zeit korrekte Ergebnisse. Je mehr Datenpunkte vorliegen, umso mehr Rechenzeit müssen Sie ansetzen.

Der Python-Interpreter spannt ohne explizite weitere Angaben nur einen einzigen Prozessor beziehungsweise Kern ein. Möchten Sie alle Kerne einer modernen CPU ausnutzen, müssen Sie den Programmcode anpassen [18]. Das gelingt über Parallelisierung einzelner Abschnitte und das Aufteilen in Prozesse oder Threads. Ohne eine Korrektur des Programmcodes klappt das allerdings selten. Die entsprechende Anpassung an die parallele Verarbeitung behandeln wir in einem Folgebeitrag. (csi)