Merge remote-tracking branch 'origin/main'

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)\
-        -> Tuple[NDArray, float, NDArray]:
+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]:
     """
     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])
 
+    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,7 +24,25 @@ def Sehne(P1: NDArray, P2: NDArray) -> float:
     return s
 
 
-def midpoint_fitness(x: tuple) -> float:
+
+
+
+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.
+    Die Distanzen zwischen den einzelnen Punkten werden als direkte 3D-Distanzen berechnet und aufaddiert.
+    :param ell: Ellipsoid
+    :param P0: Startpunkt
+    :param Pk: Zielpunkt
+    :param maxSegLen: maximale Segmentlänge
+    :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)
@@ -33,7 +50,7 @@ def midpoint_fitness(x: tuple) -> float:
         :param x: enthält die Startwerte von u und v
         :return: Fitnesswert (f)
         """
-    global P_left, P_right
+        nonlocal P_left, P_right, ell
 
         u, v = x
         P_middle = ell.para2cart(u, v)
@@ -51,19 +68,6 @@ def midpoint_fitness(x: tuple) -> float:
 
         return f
 
-
-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: Ellipsoid
-    :param P0: Startpunkt
-    :param Pk: Zielpunkt
-    :param maxSegLen: maximale Segmentlänge
-    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird
-    :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge
-    """
     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