Merge remote-tracking branch 'origin/main'

GHA_triaxial/gha2_ES.py

@@ -2,11 +2,12 @@ import numpy as np
 from Hansen_ES_CMA import escma
 from ellipsoide import EllipsoidTriaxial
 from numpy.typing import NDArray
+from typing import Tuple
 import plotly.graph_objects as go
 from GHA_triaxial.gha2_num import gha2_num
-from GHA_triaxial.utils import sigma2alpha
+from GHA_triaxial.utils import sigma2alpha, pq_ell
+
 
-ell_ES: EllipsoidTriaxial = None
 P_left: NDArray = None
 P_right: NDArray = None
 
@@ -20,6 +21,7 @@ def Sehne(P1: NDArray, P2: NDArray) -> float:
     """
     R12 = P2-P1
     s = np.linalg.norm(R12)
+
     return s
 
 
@@ -31,45 +33,42 @@ def midpoint_fitness(x: tuple) -> float:
     :param x: enthält die Startwerte von u und v
     :return: Fitnesswert (f)
     """
-    global ell_ES, P_left, P_right
+    global P_left, P_right
 
     u, v = x
-    P_middle = ell_ES.para2cart(u, v)
+    P_middle = ell.para2cart(u, v)
     d1 = Sehne(P_left, P_middle)
     d2 = Sehne(P_middle, P_right)
     base = d1 + d2
 
     # midpoint penalty (dimensionslos)
-    # relative Differenz, skaliert stabil über verschiedene Segmentlängen
+    # relative Differenz, skaliert über verschiedene Segmentlängen
     denom = max(base, 1e-9)
     pen_equal = ((d1 - d2) / denom) ** 2
     w_equal = 10.0
 
     f = base + denom * w_equal * pen_equal
+
     return f
 
 
-def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, maxIter: int = 10000, all_points: bool = False):
+def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, all_points: bool = True) -> Tuple[float, float, float, NDArray]:
     """
     Berechnen der 2. GHA mithilfe der CMA-ES.
     Die CMA-ES optimiert sukzessive den Mittelpunkt zwischen Start- und Zielpunkt. Der Abbruch der Berechnung erfolgt, wenn alle Segmentlängen <= maxSegLen sind.
     Die Distanzen zwischen den einzelnen Punkten werden als direkte 3D-Distanzen berechnet und aufaddiert.
-    :param ell: Parameter des triaxialen Ellipsoids
+    :param ell: Ellipsoid
     :param P0: Startpunkt
     :param Pk: Zielpunkt
     :param maxSegLen: maximale Segmentlänge
-    :param maxIter: maximale Durchläufe der Mittelpunktsgenerierung
+    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird
-    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben werden kann
+    :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge
-    :return: Richtungswinkel des Start- und Zielpunktes und Gesamtlänge
     """
-    global ell_ES
-    ell_ES = ell
     R0 = (ell.ax + ell.ay + ell.b) / 3
     if maxSegLen is None:
         maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius
 
     sigma_uv_nom = 1e-3 * (maxSegLen / R0)  # ~1e-5
-
     points: list[NDArray] = [P0, Pk]
     startIter = 0
     level = 0
@@ -82,7 +81,6 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
 
         level += 1
         new_points: list[NDArray] = [points[0]]
-
         for i in range(len(points) - 1):
             A = points[i]
             B = points[i+1]
@@ -92,22 +90,22 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
             if dAB > maxSegLen:
                 global P_left, P_right
                 P_left, P_right = A, B
-                Au, Av = ell_ES.cart2para(A)
+                Au, Av = ell.cart2para(A)
-                Bu, Bv = ell_ES.cart2para(B)
+                Bu, Bv = ell.cart2para(B)
                 u0 = (Au + Bu) / 2
                 v0 = Av + 0.5 * np.arctan2(np.sin(Bv - Av), np.cos(Bv - Av))
                 xmean = [u0, v0]
 
                 sigmaStep = sigma_uv_nom * (Sehne(A, B) / maxSegLen)
 
-                u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep)
+                u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) # Aufruf CMA-ES
 
                 P_next = ell.para2cart(u, v)
                 new_points.append(P_next)
                 startIter += 1
+                maxIter = 10000
                 if startIter > maxIter:
                     raise RuntimeError("Abbruch: maximale Iterationen überschritten.")
-
             new_points.append(B)
 
         points = new_points
@@ -117,13 +115,13 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
     totalLen = float(np.sum(np.linalg.norm(P_all[1:] - P_all[:-1], axis=1)))
 
     if len(points) >= 3:
-        p0i = ell_ES.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0))
+        p0i = ell.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0))
         sigma0 = (p0i - P0) / np.linalg.norm(p0i - P0)
-        alpha0 = sigma2alpha(ell_ES, sigma0, P0)
+        alpha0 = sigma2alpha(ell, sigma0, P0)
 
-        p1i = ell_ES.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2]))
+        p1i = ell.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2]))
         sigma1 = (Pk - p1i) / np.linalg.norm(Pk - p1i)
-        alpha1 = sigma2alpha(ell_ES, sigma1, Pk)
+        alpha1 = sigma2alpha(ell, sigma1, Pk)
     else:
         alpha0 = None
         alpha1 = None
@@ -166,7 +164,7 @@ if __name__ == '__main__':
 
     beta0, lamb0 = (0.2, 0.1)
     P0 = ell.ell2cart(beta0, lamb0)
-    beta1, lamb1 = (0.7, 0.3)
+    beta1, lamb1 = (0.3, 0.2)
     P1 = ell.ell2cart(beta1, lamb1)
 
     alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000, all_points=True)
@@ -175,8 +173,10 @@ if __name__ == '__main__':
         points_num.append(ell.ell2cart(beta, lamb))
     points_num = np.array(points_num)
 
-    alpha0, alpha1, s, points = gha2_ES(ell, P0, P1, all_points=True)
+    alpha0, alpha1, s, points = gha2_ES(ell, P0, P1)
     print(s_num)
     print(s)
+    print(alpha0)
+    print(alpha1)
     print(s - s_num)
     show_points(points, points_num, P0, P1)