From cd17508c03fcdff25c65c46f31aea582f43ff70f Mon Sep 17 00:00:00 2001 From: Michelle Burfeind Date: Sat, 6 Dec 2025 16:19:34 +0100 Subject: [PATCH 1/3] Pythonfiles --- Stochastisches_Modell.py | 122 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) diff --git a/Stochastisches_Modell.py b/Stochastisches_Modell.py index e69de29..5eabecd 100644 --- a/Stochastisches_Modell.py +++ b/Stochastisches_Modell.py @@ -0,0 +1,122 @@ +import sympy as sp +from dataclasses import dataclass, field +from typing import Dict, Tuple + +@dataclass +class StochastischesModellApriori: + + sigma_obs: Iterable[float] # σ_i + group_ids: Iterable[int] # Gruppenzugehörigkeit der i-ten Beobachtung + sigma0_sq_groups: Dict[int, float] = field(default_factory=dict) + + def __post_init__(self): + # In sympy-Objekte konvertieren + self.sigma_obs = sp.Matrix(list(self.sigma_obs)) # Spaltenvektor + self.group_ids = sp.Matrix(list(self.group_ids)) # Spaltenvektor + + if self.sigma_obs.rows != self.group_ids.rows: + raise ValueError("sigma_obs und group_ids müssen gleich viele Einträge haben.") + + # Fehlende Gruppen mit σ_0j^2 = 1.0 initialisieren + unique_groups = sorted({int(g) for g in self.group_ids}) + for g in unique_groups: + if g not in self.sigma0_sq_groups: + self.sigma0_sq_groups[g] = 1.0 + + @property + + def n_obs(self) -> int: + return int(self.sigma_obs.rows) + + + def build_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]: + + n = self.n_obs + Q_ll = sp.zeros(n, n) + P = sp.zeros(n, n) + + for i in range(n): + sigma_i = self.sigma_obs[i, 0] + g = int(self.group_ids[i, 0]) + sigma0_sq = self.sigma0_sq_groups[g] + + q_ii = sigma_i**2 + Q_ll[i, i] = q_ii + + P[i, i] = 1 / (sigma0_sq * q_ii) + + return Q_ll, P + + @staticmethod + def _redundanz_pro_beobachtung(A: sp.Matrix, P: sp.Matrix) -> sp.Matrix: + + n_obs = P.rows + n_param = A.cols + + # P^(1/2) aufbauen (diagonal, sqrt der Diagonale) + sqrtP = sp.zeros(n_obs, n_obs) + for i in range(n_obs): + sqrtP[i, i] = sp.sqrt(P[i, i]) + + A_tilde = sqrtP * A # Ã + + # M = (Ãᵀ Ã)^(-1) + M = (A_tilde.T * A_tilde).inv() + + r_vec = sp.zeros(n_obs, 1) + + for i in range(n_obs): + a_i = A_tilde.row(i) # 1 × n_param + a_i_row = sp.Matrix([a_i]) # explizit 1×n-Matrix + r_i = 1 - (a_i_row * M * a_i_row.T)[0, 0] + r_vec[i, 0] = r_i + + return r_vec + + + def varianzkomponenten_schaetzung( + self, + v: sp.Matrix, # Residuenvektor (n × 1) + A: sp.Matrix, # Designmatrix + ) -> Dict[int, float]: + + if v.rows != self.n_obs: + raise ValueError("Länge von v passt nicht zur Anzahl Beobachtungen im Modell.") + + # Aktuelle Gewichte + Q_ll, P = self.build_Qll_P() + + # Redundanzzahlen pro Beobachtung + r_vec = self._redundanz_pro_beobachtung(A, P) + + new_sigma0_sq: Dict[int, float] = {} + + # Für jede Gruppe j: + unique_groups = sorted({int(g) for g in self.group_ids}) + + for g in unique_groups: + # Indizes der Beobachtungen in dieser Gruppe + idx = [i for i in range(self.n_obs) if int(self.group_ids[i, 0]) == g] + if not idx: + continue + + # v_j, P_j, r_j extrahieren + v_j = sp.Matrix([v[i, 0] for i in idx]) # (m_j × 1) + P_j = sp.zeros(len(idx), len(idx)) + r_j = 0 + for ii, i in enumerate(idx): + P_j[ii, ii] = P[i, i] + r_j += r_vec[i, 0] + + # σ̂_j^2 = (v_jᵀ P_j v_j) / r_j + sigma_hat_j_sq = (v_j.T * P_j * v_j)[0, 0] / r_j + + # als float rausgeben, kann man aber auch symbolisch lassen + new_sigma0_sq[g] = float(sigma_hat_j_sq) + + return new_sigma0_sq + + def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None: + + for g, val in new_sigma0_sq.items(): + self.sigma0_sq_groups[int(g)] = float(val) \ No newline at end of file From 350941e690f8c6969d076edfbf5536932d09f43b Mon Sep 17 00:00:00 2001 From: Michelle Burfeind Date: Sat, 6 Dec 2025 16:20:13 +0100 Subject: [PATCH 2/3] Pythonfiles --- .idea/.gitignore | 5 +++++ .idea/Masterprojekt-Campusnetz.iml | 8 ++++++++ .idea/inspectionProfiles/profiles_settings.xml | 6 ++++++ .idea/misc.xml | 7 +++++++ .idea/modules.xml | 8 ++++++++ .idea/vcs.xml | 6 ++++++ 6 files changed, 40 insertions(+) create mode 100644 .idea/.gitignore create mode 100644 .idea/Masterprojekt-Campusnetz.iml create mode 100644 .idea/inspectionProfiles/profiles_settings.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/vcs.xml diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..b58b603 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,5 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ diff --git a/.idea/Masterprojekt-Campusnetz.iml b/.idea/Masterprojekt-Campusnetz.iml new file mode 100644 index 0000000..c03f621 --- /dev/null +++ b/.idea/Masterprojekt-Campusnetz.iml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..590a59e --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..7276197 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..35eb1dd --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file From 2a4fa83c473d2fbba6ae0a2d652e700589fc1d99 Mon Sep 17 00:00:00 2001 From: Michelle Burfeind Date: Sun, 7 Dec 2025 16:07:51 +0100 Subject: [PATCH 3/3] Pythonfiles --- Stochastisches_Modell.py | 82 ++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 46 deletions(-) diff --git a/Stochastisches_Modell.py b/Stochastisches_Modell.py index 5eabecd..bd882bd 100644 --- a/Stochastisches_Modell.py +++ b/Stochastisches_Modell.py @@ -4,68 +4,60 @@ from typing import Dict, Tuple @dataclass class StochastischesModellApriori: + sigma_beob: Iterable[float] #σ der einzelnen Beobachtung + group_beob: Iterable[int] #Gruppenzugehörigkeit jeder Beobachtung (Distanz, Richtung, GNSS, Nivellement,...,) + sigma0_groups: Dict[int, float] = field(default_factory=dict) #σ0² für jede Gruppe - sigma_obs: Iterable[float] # σ_i - group_ids: Iterable[int] # Gruppenzugehörigkeit der i-ten Beobachtung - sigma0_sq_groups: Dict[int, float] = field(default_factory=dict) def __post_init__(self): - # In sympy-Objekte konvertieren - self.sigma_obs = sp.Matrix(list(self.sigma_obs)) # Spaltenvektor - self.group_ids = sp.Matrix(list(self.group_ids)) # Spaltenvektor + self.sigma_beob = sp.Matrix(list(self.sigma_beob)) #Spaltenvektor + self.group_beob = sp.Matrix(list(self.group_beob)) #Spaltenvektor - if self.sigma_obs.rows != self.group_ids.rows: + if self.sigma_beob.rows != self.group_beob.rows: raise ValueError("sigma_obs und group_ids müssen gleich viele Einträge haben.") - # Fehlende Gruppen mit σ_0j^2 = 1.0 initialisieren - unique_groups = sorted({int(g) for g in self.group_ids}) + unique_groups = sorted({int(g) for g in self.group_beob}) #jede Beobachtungsgruppe wird genau einmal berücksichtigt for g in unique_groups: - if g not in self.sigma0_sq_groups: - self.sigma0_sq_groups[g] = 1.0 + if g not in self.sigma0_groups: #Fehlende Gruppen mit σ_0j^2 = 1.0 + self.sigma0_groups[g] = 1.0 + @property - - def n_obs(self) -> int: - return int(self.sigma_obs.rows) + def n_beob(self) -> int: + return int(self.sigma_beob.rows) - def build_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]: - - n = self.n_obs + def aufstellen_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]: + n = self.n_beob Q_ll = sp.zeros(n, n) P = sp.zeros(n, n) for i in range(n): - sigma_i = self.sigma_obs[i, 0] - g = int(self.group_ids[i, 0]) - sigma0_sq = self.sigma0_sq_groups[g] - - q_ii = sigma_i**2 - Q_ll[i, i] = q_ii - - P[i, i] = 1 / (sigma0_sq * q_ii) - + sigma_i = self.sigma_beob[i, 0] #σ-Wert der i-ten Beobachtung holen + g = int(self.group_beob[i, 0]) #Gruppenzugehörigkeit der Beobachtung bestimmen + sigma0_sq = self.sigma0_groups[g] #Den Varianzfaktor der Gruppe holen + q_ii = sigma_i**2 #σ² berechnen + Q_ll[i, i] = q_ii #Diagonale + P[i, i] = 1 / (sigma0_sq * q_ii) #durch VKS nicht mehr P=Qll^-1 return Q_ll, P + @staticmethod - def _redundanz_pro_beobachtung(A: sp.Matrix, P: sp.Matrix) -> sp.Matrix: + def redundanz_pro_beobachtung(A: sp.Matrix, P: sp.Matrix) -> sp.Matrix: + n_beob = P.rows #Anzahl der Beobachtungen (Zeilen in P) + n_param = A.cols #Anzahl der Unbekannten (Spalten in A) - n_obs = P.rows - n_param = A.cols - - # P^(1/2) aufbauen (diagonal, sqrt der Diagonale) - sqrtP = sp.zeros(n_obs, n_obs) - for i in range(n_obs): + sqrtP = sp.zeros(n_beob, n_beob) #Wurzel von P (der Diagonale) + for i in range(n_beob): sqrtP[i, i] = sp.sqrt(P[i, i]) - A_tilde = sqrtP * A # Ã + A_tilde = sqrtP * A - # M = (Ãᵀ Ã)^(-1) M = (A_tilde.T * A_tilde).inv() - r_vec = sp.zeros(n_obs, 1) + r_vec = sp.zeros(n_beob, 1) - for i in range(n_obs): + for i in range(n_beob): a_i = A_tilde.row(i) # 1 × n_param a_i_row = sp.Matrix([a_i]) # explizit 1×n-Matrix r_i = 1 - (a_i_row * M * a_i_row.T)[0, 0] @@ -74,29 +66,29 @@ class StochastischesModellApriori: return r_vec - def varianzkomponenten_schaetzung( + def varianzkomponentenschaetzung( self, v: sp.Matrix, # Residuenvektor (n × 1) A: sp.Matrix, # Designmatrix ) -> Dict[int, float]: - if v.rows != self.n_obs: + if v.rows != self.n_beob: raise ValueError("Länge von v passt nicht zur Anzahl Beobachtungen im Modell.") # Aktuelle Gewichte - Q_ll, P = self.build_Qll_P() + Q_ll, P = self.aufstellen_Qll_P() # Redundanzzahlen pro Beobachtung - r_vec = self._redundanz_pro_beobachtung(A, P) + r_vec = self.redundanz_pro_beobachtung(A, P) new_sigma0_sq: Dict[int, float] = {} # Für jede Gruppe j: - unique_groups = sorted({int(g) for g in self.group_ids}) + unique_groups = sorted({int(g) for g in self.group_beob}) for g in unique_groups: # Indizes der Beobachtungen in dieser Gruppe - idx = [i for i in range(self.n_obs) if int(self.group_ids[i, 0]) == g] + idx = [i for i in range(self.n_beob) if int(self.group_beob[i, 0]) == g] if not idx: continue @@ -113,10 +105,8 @@ class StochastischesModellApriori: # als float rausgeben, kann man aber auch symbolisch lassen new_sigma0_sq[g] = float(sigma_hat_j_sq) - return new_sigma0_sq def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None: - for g, val in new_sigma0_sq.items(): - self.sigma0_sq_groups[int(g)] = float(val) \ No newline at end of file + self.sigma0_groups[int(g)] = float(val) \ No newline at end of file