fa2338/Masterprojekt-Campusnetz
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Pythonfiles

Browse Source
This commit is contained in:
mi2847
2025-12-07 16:07:51 +01:00
parent 350941e690
commit 2a4fa83c47
1 changed files with 36 additions and 46 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
Stochastisches_Modell.py
View File

@@ -4,68 +4,60 @@ from typing import Dict, Tuple
@dataclass @dataclass
class StochastischesModellApriori: 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): def __post_init__(self):
# In sympy-Objekte konvertieren self.sigma_beob = sp.Matrix(list(self.sigma_beob)) #Spaltenvektor
self.sigma_obs = sp.Matrix(list(self.sigma_obs)) # Spaltenvektor self.group_beob = sp.Matrix(list(self.group_beob)) #Spaltenvektor
self.group_ids = sp.Matrix(list(self.group_ids)) # 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.") 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_beob}) #jede Beobachtungsgruppe wird genau einmal berücksichtigt
unique_groups = sorted({int(g) for g in self.group_ids})
for g in unique_groups: for g in unique_groups:
if g not in self.sigma0_sq_groups: if g not in self.sigma0_groups: #Fehlende Gruppen mit σ_0j^2 = 1.0
self.sigma0_sq_groups[g] = 1.0 self.sigma0_groups[g] = 1.0
@property @property
def n_beob(self) -> int:
def n_obs(self) -> int: return int(self.sigma_beob.rows)
return int(self.sigma_obs.rows)
def build_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]: def aufstellen_Qll_P(self) -> Tuple[sp.Matrix, sp.Matrix]:
n = self.n_beob
n = self.n_obs
Q_ll = sp.zeros(n, n) Q_ll = sp.zeros(n, n)
P = sp.zeros(n, n) P = sp.zeros(n, n)
for i in range(n): for i in range(n):
sigma_i = self.sigma_obs[i, 0] sigma_i = self.sigma_beob[i, 0] #σ-Wert der i-ten Beobachtung holen
g = int(self.group_ids[i, 0]) g = int(self.group_beob[i, 0]) #Gruppenzugehörigkeit der Beobachtung bestimmen
sigma0_sq = self.sigma0_sq_groups[g] sigma0_sq = self.sigma0_groups[g] #Den Varianzfaktor der Gruppe holen
q_ii = sigma_i**2 #σ² berechnen
q_ii = sigma_i**2 Q_ll[i, i] = q_ii #Diagonale
Q_ll[i, i] = q_ii P[i, i] = 1 / (sigma0_sq * q_ii) #durch VKS nicht mehr P=Qll^-1
P[i, i] = 1 / (sigma0_sq * q_ii)
return Q_ll, P return Q_ll, P
@staticmethod @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 sqrtP = sp.zeros(n_beob, n_beob) #Wurzel von P (der Diagonale)
n_param = A.cols for i in range(n_beob):
# 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]) 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() 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 = A_tilde.row(i) # 1 × n_param
a_i_row = sp.Matrix([a_i]) # explizit 1×n-Matrix 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_i = 1 - (a_i_row * M * a_i_row.T)[0, 0]
@@ -74,29 +66,29 @@ class StochastischesModellApriori:
return r_vec return r_vec
def varianzkomponenten_schaetzung( def varianzkomponentenschaetzung(
self, self,
v: sp.Matrix, # Residuenvektor (n × 1) v: sp.Matrix, # Residuenvektor (n × 1)
A: sp.Matrix, # Designmatrix A: sp.Matrix, # Designmatrix
) -> Dict[int, float]: ) -> 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.") raise ValueError("Länge von v passt nicht zur Anzahl Beobachtungen im Modell.")
# Aktuelle Gewichte # Aktuelle Gewichte
Q_ll, P = self.build_Qll_P() Q_ll, P = self.aufstellen_Qll_P()
# Redundanzzahlen pro Beobachtung # 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] = {} new_sigma0_sq: Dict[int, float] = {}
# Für jede Gruppe j: # 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: for g in unique_groups:
# Indizes der Beobachtungen in dieser Gruppe # 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: if not idx:
continue continue
@@ -113,10 +105,8 @@ class StochastischesModellApriori:
# als float rausgeben, kann man aber auch symbolisch lassen # als float rausgeben, kann man aber auch symbolisch lassen
new_sigma0_sq[g] = float(sigma_hat_j_sq) new_sigma0_sq[g] = float(sigma_hat_j_sq)
return new_sigma0_sq return new_sigma0_sq
def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None: def update_sigma0(self, new_sigma0_sq: Dict[int, float]) -> None:
for g, val in new_sigma0_sq.items(): for g, val in new_sigma0_sq.items():
self.sigma0_sq_groups[int(g)] = float(val) self.sigma0_groups[int(g)] = float(val)