Merge remote-tracking branch 'origin/main'

Campusnetz.ipynb

@@ -23,7 +23,8 @@
     "import Stochastisches_Modell\n",
     "from Stochastisches_Modell import StochastischesModell\n",
     "import Export\n",
-    "import Netzqualität_Genauigkeit"
+    "import Netzqualität_Genauigkeit\n",
+    "import Datumsfestlegung"
    ],
    "outputs": [],
    "execution_count": 1
@@ -934,6 +935,24 @@
    "id": "122dca077d1d267c",
    "outputs": [],
    "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "outputs": [],
+   "execution_count": null,
+   "source": [
+    "# Datumsfestlegung: Bitte geben Sie nachfolgend die Koordinatenkomponenten an, die das Datum definieren sollen\n",
+    "\n",
+    "auswahl = [\n",
+    "    (\"101\",\"X\"), (\"101\",\"Y\"),          # Punkt 101 nur Lage\n",
+    "    (\"205\",\"X\"), (\"205\",\"Y\"), (\"205\",\"Z\"),  # Punkt 205 voll\n",
+    "    (\"330\",\"Z\")                        # Punkt 330 nur Höhe\n",
+    "]\n",
+    "\n",
+    "aktive_unbekannte_indices = Datumsfestlegung.datumskomponenten(auswahl, liste_punktnummern)"
+   ],
+   "id": "c37670a07848d977"
   }
  ],
  "metadata": {

Datumsfestlegung.py

@@ -1,81 +1,154 @@
 import sympy as sp
-from typing import List, Iterable, Tuple
+from typing import Iterable, List, Sequence, Tuple, Optional
 
+
+class Datumsfestlegung:
+
+    @staticmethod
+    def datumskomponenten(
+        auswahl: Iterable[Tuple[str, str]],
+        liste_punktnummern: Sequence[str],
+        *,
+        layout: str = "XYZ"
+    ) -> List[int]:
+        punkt2pos = {str(p): i for i, p in enumerate(liste_punktnummern)}
+
+        layout = layout.upper()
+        if layout != "XYZ":
+            raise ValueError("Nur layout='XYZ' unterstützt (wie bei euch).")
+        comp2off = {"X": 0, "Y": 1, "Z": 2}
+
+        aktive: List[int] = []
+        for pt, comp in auswahl:
+            spt = str(pt)
+            c = comp.upper()
+            if spt not in punkt2pos:
+                raise KeyError(f"Punkt '{pt}' nicht in liste_punktnummern.")
+            if c not in comp2off:
+                raise ValueError(f"Komponente '{comp}' ungültig. Nur X,Y,Z.")
+            p = punkt2pos[spt]
+            aktive.append(3 * p + comp2off[c])
+
+        # Duplikate entfernen
+        out, seen = [], set()
+        for i in aktive:
+            if i not in seen:
+                seen.add(i)
+                out.append(i)
+        return out
+
+    @staticmethod
+    def auswahlmatrix_E(u: int, aktive_unbekannte_indices: Iterable[int]) -> sp.Matrix:
+        E = sp.zeros(u, u)
+        for idx in aktive_unbekannte_indices:
+            i = int(idx)
+            if not (0 <= i < u):
+                raise IndexError(f"Aktiver Index {i} außerhalb [0,{u-1}]")
+            E[i, i] = 1
+        return E
+
+    @staticmethod
     def raenderungsmatrix_G(
         x0: sp.Matrix,
-    idx_X: List[int],
+        liste_punktnummern: Sequence[str],
-    idx_Y: List[int],
+        *,
-    idx_Z: List[int],
         mit_massstab: bool = True,
+        layout: str = "XYZ",
     ) -> sp.Matrix:
+        if x0.cols != 1:
+            raise ValueError("x0 muss Spaltenvektor sein.")
+        layout = layout.upper()
+        if layout != "XYZ":
+            raise ValueError("Nur layout='XYZ' unterstützt (wie bei euch).")
 
+        nP = len(liste_punktnummern)
         u = x0.rows
         d = 7 if mit_massstab else 6
         G = sp.zeros(u, d)
 
-    # --- Translationen ---
+        for p in range(nP):
-    for i in idx_X:
+            ix, iy, iz = 3*p, 3*p+1, 3*p+2
-        G[i, 0] = 1
+            xi, yi, zi = x0[ix, 0], x0[iy, 0], x0[iz, 0]
-    for i in idx_Y:
-        G[i, 1] = 1
-    for i in idx_Z:
-        G[i, 2] = 1
 
-    # --- Rotationen ---
+            # Translationen
-    # Rotation um X-Achse
+            G[ix, 0] = 1
-    for iy, iz in zip(idx_Y, idx_Z):
+            G[iy, 1] = 1
-        zi = x0[iz, 0]
+            G[iz, 2] = 1
-        yi = x0[iy, 0]
-        G[iy, 3] = -zi
-        G[iz, 3] = yi
 
-    # Rotation um Y-Achse
+            # Rotationen
-    for ix, iz in zip(idx_X, idx_Z):
+            G[iy, 3] = -zi; G[iz, 3] =  yi     # Rx
-        zi = x0[iz, 0]
+            G[ix, 4] =  zi; G[iz, 4] = -xi     # Ry
-        xi = x0[ix, 0]
+            G[ix, 5] = -yi; G[iy, 5] =  xi     # Rz
-        G[ix, 4] = zi
-        G[iz, 4] = -xi
 
-    # Rotation um Z-Achse
+            # Maßstab
-    for ix, iy in zip(idx_X, idx_Y):
-        yi = x0[iy, 0]
-        xi = x0[ix, 0]
-        G[ix, 5] = -yi
-        G[iy, 5] = xi
-
-    # --- Maßstab ---
             if mit_massstab:
-        for ix, iy, iz in zip(idx_X, idx_Y, idx_Z):
-            xi = x0[ix, 0]
-            yi = x0[iy, 0]
-            zi = x0[iz, 0]
                 G[ix, 6] = xi
                 G[iy, 6] = yi
                 G[iz, 6] = zi
+
         return G
 
-
+    @staticmethod
-def auswahlmatrix_E(u: int, aktive_unbekannte_indices: Iterable[int]) -> sp.Matrix:
-    E = sp.zeros(u, u)
-    for idx in aktive_unbekannte_indices:
-        E[int(idx), int(idx)] = 1
-    return E
-
-
-def teilspurminimierung_Gi(G: sp.Matrix, E: sp.Matrix) -> sp.Matrix:
-    Gi = E * G
-    return Gi
-
-
     def berechne_dx_geraendert(N: sp.Matrix, n: sp.Matrix, Gi: sp.Matrix) -> sp.Matrix:
-    u = N.rows
+        if N.rows != N.cols:
-    d = Gi.shape[1]
+            raise ValueError("N muss quadratisch sein.")
+        if n.cols != 1:
+            raise ValueError("n muss Spaltenvektor sein.")
+        if Gi.rows != N.rows:
+            raise ValueError("Gi hat falsche Zeilenzahl.")
 
+        u = N.rows
+        d = Gi.cols
         K = N.row_join(Gi)
         K = K.col_join(Gi.T.row_join(sp.zeros(d, d)))
-
         rhs = n.col_join(sp.zeros(d, 1))
-
         sol = K.LUsolve(rhs)
-    dx = sol[:u, :]
+        return sol[:u, :]
-    return dx
+
+    @staticmethod
+    def weiches_datum(
+        A: sp.Matrix,
+        dl: sp.Matrix,
+        Q_ll: sp.Matrix,
+        x0: sp.Matrix,
+        anschluss_indices: Sequence[int],
+        anschluss_werte: sp.Matrix,
+        Sigma_AA: Optional[sp.Matrix] = None,
+    ) -> Tuple[sp.Matrix, sp.Matrix, sp.Matrix]:
+        if dl.cols != 1 or x0.cols != 1:
+            raise ValueError("dl und x0 müssen Spaltenvektoren sein.")
+        if A.rows != dl.rows:
+            raise ValueError("A.rows muss dl.rows entsprechen.")
+        if A.cols != x0.rows:
+            raise ValueError("A.cols muss x0.rows entsprechen.")
+        if Q_ll.rows != Q_ll.cols or Q_ll.rows != A.rows:
+            raise ValueError("Q_ll muss (n×n) sein und zu A.rows passen.")
+
+        u = A.cols
+        idx = [int(i) for i in anschluss_indices]
+        m = len(idx)
+
+        if anschluss_werte.cols != 1 or anschluss_werte.rows != m:
+            raise ValueError("anschluss_werte muss (m×1) sein.")
+        if Sigma_AA is None:
+            Sigma_AA = sp.eye(m)
+        if Sigma_AA.rows != m or Sigma_AA.cols != m:
+            raise ValueError("Sigma_AA muss (m×m) sein.")
+
+        A_A = sp.zeros(m, u)
+        for r, j in enumerate(idx):
+            if not (0 <= j < u):
+                raise IndexError(f"Anschluss-Index {j} außerhalb [0,{u-1}]")
+            A_A[r, j] = 1
+
+        x0_A = sp.Matrix([[x0[j, 0]] for j in idx])
+        dl_A = anschluss_werte - x0_A
+
+        A_ext = A.col_join(A_A)
+        dl_ext = dl.col_join(dl_A)
+
+        Q_ext = sp.zeros(Q_ll.rows + m, Q_ll.cols + m)
+        Q_ext[:Q_ll.rows, :Q_ll.cols] = Q_ll
+        Q_ext[Q_ll.rows:, Q_ll.cols:] = Sigma_AA
+
+        return A_ext, dl_ext, Q_ext

Parameterschaetzung.py

@@ -1,29 +1,62 @@
-from Datumsfestlegung import *
 from Stochastisches_Modell import StochastischesModell
+from Netzqualität_Genauigkeit import Genauigkeitsmaße
+from Datumsfestlegung import Datumsfestlegung
 import sympy as sp
 import Export
-import Netzqualität_Genauigkeit
 
 
-def ausgleichung_global(A, dl, stoch_modell: StochastischesModell):
+def ausgleichung_global(
+    A: sp.Matrix,
+    dl: sp.Matrix,
+    Q_ll: sp.Matrix,
+    x0: sp.Matrix,
+    idx_X, idx_Y, idx_Z,
+    anschluss_indices,
+    anschluss_werte,
+    Sigma_AA,
+):
+    # 1) Datumsfestlegung (weiches Datum) System erweitern
+    A_ext, dl_ext, Q_ext = Datumsfestlegung.weiches_datum(
+        A=A,
+        dl=dl,
+        Q_ll=Q_ll,
+        x0=x0,
+        anschluss_indices=anschluss_indices,
+        anschluss_werte=anschluss_werte,
+        Sigma_AA=Sigma_AA,
+    )
 
-    #Q_ll, P = stoch_modell.berechne_Qll_P()        #Kofaktormatrix und P-Matrix
+    # 2) Gewichtsmatrix P
-    P = sp.eye(A.shape[0])
+    P = StochastischesModell.berechne_P(Q_ext)
-    N = A.T * P * A                                 #Normalgleichungsmatrix N
-    Q_xx = N.inv()                                  #Kofaktormatrix der Unbekannten Qxx
-    n = A.T * P * dl                                #Absolutgliedvektor n
 
-    dx = N.LUsolve(n)                               #Zuschlagsvektor dx
+    # 3) Normalgleichungsmatrix N und Absolutgliedvektor n
+    N = A_ext.T * P * A_ext
+    n = A_ext.T * P * dl_ext
 
-    v = dl - A * dx                                 #Residuenvektor v
+    # 4) Zuschlagsvektor dx
+    dx = N.LUsolve(n)
 
+    # 5) Residuenvektor v
+    v = dl - A * dx
+
+    # 6) Kofaktormatrix der Unbekannten Q_xx
+    Q_xx = StochastischesModell.berechne_Q_xx(N)
+
+    # 7) Kofaktormatrix der Beobachtungen Q_ll_dach
     Q_ll_dach = A * Q_xx * A.T
-    Q_vv = stoch_modell.berechne_Qvv(A, P, Q_xx)          #Kofaktormatrix der Verbesserungen Qvv
-    R = stoch_modell.berechne_R(Q_vv, P)                  #Redundanzmatrix R
-    r = stoch_modell.berechne_r(R)                        #Redundanzanteile als Vektor r
-    redundanzanteile = A.shape[0] - A.shape[1] #n-u+d
-    soaposteriori = Netzqualität_Genauigkeit.Genauigkeitsmaße.s0apost(v, P, redundanzanteile)
 
+    # 8) Kofaktormatrix der Verbesserungen Q_vv
+    Q_vv = StochastischesModell.berechne_Qvv(A, P, Q_xx)
+
+    # 9) Redundanzmatrix R und Redundanzanteile r
+    R = StochastischesModell.berechne_R(Q_vv, P)                  #Redundanzmatrix R
+    r = StochastischesModell.berechne_r(R)                        #Redundanzanteile als Vektor r
+    redundanzanteile = A.shape[0] - A.shape[1] #n-u+d
+
+    # 10) s0 a posteriori
+    soaposteriori = Genauigkeitsmaße.s0apost(v, P, redundanzanteile)
+
+    # 11) Ausgabe
     dict_ausgleichung = {
         "dx": dx,
         "v": v,
@@ -38,68 +71,65 @@ def ausgleichung_global(A, dl, stoch_modell: StochastischesModell):
     }
 
     Export.Export.ausgleichung_to_datei(r"Zwischenergebnisse\Ausgleichung_Iteration0.csv", dict_ausgleichung)
-
     return dict_ausgleichung, dx
 
 
 
+
 def ausgleichung_lokal(
     A: sp.Matrix,
     dl: sp.Matrix,
-    stoch_modell: StochastischesModell,
+    Q_ll: sp.Matrix,
     x0: sp.Matrix,
     idx_X, idx_Y, idx_Z,
     aktive_unbekannte_indices,
     mit_massstab: bool = True,
 ):
-    # 1) Gewichte
+    # 1) Gewichtsmatrix P
-    Q_ll, P = stoch_modell.berechne_Qll_P()
+    P = StochastischesModell.berechne_P(Q_ll)
-    # Debug-Option:
-    # P = sp.eye(A.rows)
 
-    # 2) Normalgleichungen
+    # 2) Normalgleichungsmatrix N und Absolutgliedvektor n
     N = A.T * P * A
     n = A.T * P * dl
 
-    # 3) Datum (G, E, Gi)
+    # 3) Datumsfestlegung (Teilspurminimierung)
-    G = raenderungsmatrix_G(x0, idx_X, idx_Y, idx_Z, mit_massstab=mit_massstab)
+    G = Datumsfestlegung.raenderungsmatrix_G(x0, liste_punktnummern, mit_massstab=mit_massstab)
-    E = auswahlmatrix_E(u=A.cols, aktive_unbekannte_indices=aktive_unbekannte_indices)
+    aktive = Datumsfestlegung.datumskomponenten(auswahl, liste_punktnummern)
+    E = Datumsfestlegung.auswahlmatrix_E(u=A.cols, aktive_unbekannte_indices=aktive)
     Gi = E * G
 
-    # 4) Geränderte Lösung (dx)
+    # 4) Zuschlagsvektor dx
-    dx = berechne_dx_geraendert(N, n, Gi)
+    dx = Datumsfestlegung.berechne_dx_geraendert(N, n, Gi)
 
-    # 5) Residuen
+    # 5) Residuenvektor v
     v = dl - A * dx
 
-    # 6) KORREKTE Q_xx für gerändertes Problem:
+    # 6) Kofaktormatrix der Unbekannten Q_xx
-    # Q_xx = N^{-1} - N^{-1}Gi (Gi^T N^{-1} Gi)^{-1} Gi^T N^{-1}
+    N_inv = N.inv()
-    # numerisch besser via LUsolve statt inv:
-    N_inv = N.inv()  # wenn N groß ist, kann man das unten auch ohne inv machen (siehe Hinweis)
     N_inv_G = N_inv * Gi
     S = Gi.T * N_inv_G
     S_inv = S.inv()
     Q_xx = N_inv - N_inv_G * S_inv * N_inv_G.T
 
-    # 7) Q_lhat_lhat und Q_vv
+    # 7) Kofaktormatrix der Beobachtungen Q_ll_dach
     Q_lhat_lhat = A * Q_xx * A.T
+
+    # 8) Kofaktormatrix der Verbesserungen Q_vv
     Q_vv = P.inv() - Q_lhat_lhat
 
-    # 8) Redundanzmatrix und -anteile
+    # 9) Redundanzmatrix R, Redundanzanteile r, Redundanz
     R = Q_vv * P
     r_vec = sp.Matrix(R.diagonal())
-
-    # 9) Freiheitsgrade (Redundanz gesamt)
     n_beob = A.rows
     u = A.cols
     d = Gi.shape[1]
     r_gesamt = n_beob - u + d
 
-    # 10) sigma0 a posteriori
+    # 10) s0 a posteriori
-    omega = float((v.T * P * v)[0, 0])
+    sigma0_apost = Genauigkeitsmaße.s0apost(v, P, r_gesamt)
-    sigma0_hat = (omega / float(r_gesamt)) ** 0.5
 
-    return {
+    # 11) Ausgabe
+    dict_ausgleichung_lokal = {
         "dx": dx,
         "v": v,
         "Q_ll": Q_ll,
@@ -111,7 +141,10 @@ def ausgleichung_lokal(
         "R": R,
         "r": r_vec,
         "r_gesamt": r_gesamt,
-        "sigma0_hat": sigma0_hat,
+        "sigma0_apost": sigma0_apost,
         "G": G,
         "Gi": Gi,
     }
+
+    Export.Export.ausgleichung_to_datei(r"Zwischenergebnisse\Ausgleichung_Iteration0_lokal.csv", dict_ausgleichung_lokal)
+    return dict_ausgleichung_lokal, dx

Stochastisches_Modell.py

@@ -167,11 +167,18 @@ class StochastischesModell:
 
         return Qll_numerisch
 
-    def berechne_P(Q_ll):
+    def berechne_P(Q_ll: sp.Matrix) -> sp.Matrix:
         P = Q_ll.inv()
         return P
 
 
+    def berechne_Q_xx(N: sp.Matrix) -> sp.Matrix:
+        if N.rows != N.cols:
+            raise ValueError("N muss eine quadratische Matrix sein")
+        Q_xx = N.inv()
+        return Q_xx
+
+
     def berechne_Qvv(self, A: sp.Matrix, P: sp.Matrix, Q_xx: sp.Matrix) -> sp.Matrix:
         Q_vv = P.inv() - A * Q_xx * A.T
         return Q_vv                                 #Kofaktormatrix der Beobachtungsresiduen