Nice

GHA_triaxial/gha1_ES.py

@@ -209,8 +209,8 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
     return beta_best, omega_best, P_best, alpha_end
 
 
-def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000)\
+def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000, all_points: boolean = False)\
-        -> Tuple[NDArray, float, NDArray]:
+        -> Tuple[NDArray, float, NDArray] | Tuple[NDArray, float]:
     """
     Aufruf der 1. GHA mittels CMA-ES
     :param ell: Ellipsoid
@@ -219,6 +219,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
     :param alpha0: Azimut Startkoordinate
     :param s_total: Gesamtstrecke
     :param maxSegLen: maximale Segmentlänge
+    :param all_points: Alle Punkte ausgeben?
     :return: Zielpunkt Pk, Azimut am Zielpunkt und Punktliste
     """
     beta = float(beta0)
@@ -236,7 +237,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
     while s_acc < s_total - 1e-9:
         step += 1
         ds = min(maxSegLen, s_total - s_acc)
-        print(f"[GHA1-ES] Step {step}/{nsteps_est}  ds={ds:.3f} m  s_acc={s_acc:.3f} m  beta={beta:.6f}  omega={omega:.6f}  alpha={alpha:.6f}")
+        # print(f"[GHA1-ES] Step {step}/{nsteps_est}  ds={ds:.3f} m  s_acc={s_acc:.3f} m  beta={beta:.6f}  omega={omega:.6f}  alpha={alpha:.6f}")
 
         beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_i=alpha, ds=ds, gamma0=gamma0,
                                                     ell=ell, maxSegLen=maxSegLen)
@@ -248,7 +249,10 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
     Pk = P_all[-1]
     alpha1 = float(alpha_end[-1])
 
-    return Pk, alpha1, np.array(P_all)
+    if all_points:
+        return Pk, alpha1, np.array(P_all)
+    else:
+        return Pk, alpha1
 
 
 if __name__ == "__main__":

GHA_triaxial/gha1_ana.py

@@ -1,10 +1,11 @@
+import math
 from math import comb
 from typing import Tuple
 
 import numpy as np
 from numpy import sin, cos, arctan2
 from numpy._typing import NDArray
-from scipy.special import factorial as fact
+import winkelumrechnungen as wu
 
 from ellipsoide import EllipsoidTriaxial
 from GHA_triaxial.utils import pq_para
@@ -64,7 +65,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
             x_m.append(x_(m))
             y_m.append(y_(m))
             z_m.append(z_(m))
-        fact_m = fact(m)
+        fact_m = math.factorial(m)
 
         # 22-24
         a_m.append(x_m[m] / fact_m)
@@ -112,7 +113,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
     return p1, alpha1
 
 
-def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 2) -> Tuple[NDArray, float]:
+def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]:
     """
     :param ell: Ellipsoid
     :param point: Punkt in kartesischen Koordinaten
@@ -134,3 +135,11 @@ def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, ma
         raise Exception("Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
 
     return point_end, alpha_end
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    p0 = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(20))
+    p1, alpha1 = gha1_ana(ell, p0, wu.deg2rad(36), 200000, 70)
+    print(p1, wu.rad2gms(alpha1))
+

GHA_triaxial/gha2_ES.py

@@ -8,8 +8,7 @@ from GHA_triaxial.gha2_num import gha2_num
 from GHA_triaxial.utils import sigma2alpha, pq_ell
 
 
-P_left: NDArray = None
+
-P_right: NDArray = None
 
 
 def Sehne(P1: NDArray, P2: NDArray) -> float:
@@ -25,34 +24,10 @@ def Sehne(P1: NDArray, P2: NDArray) -> float:
     return s
 
 
-def midpoint_fitness(x: tuple) -> float:
-    """
-    Fitness für einen Mittelpunkt P_middle zwischen P_left und P_right auf dem triaxialen Ellipsoid:
-    - Minimiert d(P_left, P_middle) + d(P_middle, P_right)
-    - Erzwingt d(P_left,P_middle) ≈ d(P_middle,P_right)  (echter Mittelpunkt im Sinne der Polygonkette)
-    :param x: enthält die Startwerte von u und v
-    :return: Fitnesswert (f)
-    """
-    global P_left, P_right
-
-    u, v = x
-    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 ü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, all_points: bool = True) -> Tuple[float, float, float, NDArray]:
+
+def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, all_points: bool = False) -> Tuple[float, float, float, NDArray] | Tuple[float, float, float]:
     """
     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.
@@ -64,6 +39,35 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
     :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird
     :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge
     """
+    P_left: NDArray = None
+    P_right: NDArray = None
+
+    def midpoint_fitness(x: tuple) -> float:
+        """
+        Fitness für einen Mittelpunkt P_middle zwischen P_left und P_right auf dem triaxialen Ellipsoid:
+        - Minimiert d(P_left, P_middle) + d(P_middle, P_right)
+        - Erzwingt d(P_left,P_middle) ≈ d(P_middle,P_right)  (echter Mittelpunkt im Sinne der Polygonkette)
+        :param x: enthält die Startwerte von u und v
+        :return: Fitnesswert (f)
+        """
+        nonlocal P_left, P_right, ell
+
+        u, v = x
+        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 ü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
+
     R0 = (ell.ax + ell.ay + ell.b) / 3
     if maxSegLen is None:
         maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius
@@ -85,10 +89,10 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
             A = points[i]
             B = points[i+1]
             dAB = Sehne(A, B)
-            print(dAB)
+            # print(dAB)
 
             if dAB > maxSegLen:
-                global P_left, P_right
+                # global P_left, P_right
                 P_left, P_right = A, B
                 Au, Av = ell.cart2para(A)
                 Bu, Bv = ell.cart2para(B)
@@ -109,7 +113,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
             new_points.append(B)
 
         points = new_points
-        print(f"[Level {level}] Punkte: {len(points)} | max Segment: {max_len:.3f} m")
+        # print(f"[Level {level}] Punkte: {len(points)} | max Segment: {max_len:.3f} m")
 
     P_all = np.vstack(points)
     totalLen = float(np.sum(np.linalg.norm(P_all[1:] - P_all[:-1], axis=1)))

Hansen_ES_CMA.py

@@ -64,7 +64,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
 
     gen = 0
 
-    print(f'    [CMA-ES] Start: lambda = {lambda_}, sigma ={round(sigma, 6)}, stopeval = {stopeval}')
+    # print(f'    [CMA-ES] Start: lambda = {lambda_}, sigma ={round(sigma, 6)}, stopeval = {stopeval}')
 
     while counteval < stopeval:
         gen += 1
@@ -95,14 +95,15 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
 
 
         if gen == 1 or gen%50==0:
-            print(f'    [CMA-ES] Gen {gen}, best = {round(fbest, 6)}, sigma = {sigma:.3g}')
+            # print(f'    [CMA-ES] Gen {gen}, best = {round(fbest, 6)}, sigma = {sigma:.3g}')
+            pass
 
         if noImproveGen >= maxNoImproveGen:
-            print(f'    [CMA-ES] Abbruch: keine Verbesserung > {round(absTolImprove, 3)} in {maxNoImproveGen} Generationen.')
+            # print(f'    [CMA-ES] Abbruch: keine Verbesserung > {round(absTolImprove, 3)} in {maxNoImproveGen} Generationen.')
             break
 
         if sigma < sigmaImprove:
-            print(f'    [CMA-ES] Abbruch: sigma zu klein {sigma:.3g}')
+            # print(f'    [CMA-ES] Abbruch: sigma zu klein {sigma:.3g}')
             break
 
 
@@ -140,7 +141,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
         # Escape flat fitness, or better terminate?
         if arfitness[0] == arfitness[int(np.ceil(0.7 * lambda_)) - 1]:
             sigma = sigma * np.exp(0.2 + cs / damps)
-            print('    [CMA-ES] stopfitness erreicht.')
+            # print('    [CMA-ES] stopfitness erreicht.')
             #print("warning: flat fitness, consider reformulating the objective")
             break
 
@@ -150,7 +151,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
     #print(f"{counteval}: {arfitness[0]}")
     xmin = arx[:, arindex[0]]
     bestValue = arfitness[0]
-    print(f'    [CMA-ES] Ende: Gen = {gen}, best = {round(bestValue, 6)}')
+    # print(f'    [CMA-ES] Ende: Gen = {gen}, best = {round(bestValue, 6)}')
     return xmin