import numpy as np
from ..conesFunctions import printCones
from .conesEnKF import conesEnKF
from scipy import stats
from sklearn.preprocessing import StandardScaler
rng = np.random.default_rng(seed=1)
[docs]
class conesMGEnKF(conesEnKF):
"""
Multi-Grid Ensemble Kalman Filter.
Extends conesEnKF (coarse ensemble, base class) with the state, parameters
and sampling of a single fine-mesh member, updated with the same Kalman
gain computed from the coarse ensemvle.
"""
def __init__(self, region):
"""
:param region: conesMGRegion where the analysis phase will be performed
:type region: conesMGRegion
"""
super().__init__(region)
self.set_x_anomaly()
self.set_R()
self.set_y_matrix()
self.stateFine = region.stateFine if hasattr(region, "stateFine") else np.empty((0, 1))
self.paramsFine = region.paramFine
self.xFine = np.vstack((self.stateFine.T.reshape(-1, 1), self.paramsFine.reshape(-1, 1)))
self.upxFine = np.empty(1, dtype="float")
self.upStateFine = np.empty(1, dtype="float")
self.upParamsFine = np.empty(1, dtype="float")
self.samplingFine = region.samplingFine
[docs]
def updateFine(self):
"""
Updates the state and parameters matrix for the fine member
:returns: None
"""
printCones("xFine =", self.xFine.shape)
printCones("K shape =", self.K.shape)
printCones("y shape =", np.mean(self.y, axis=1).shape)
printCones("sampFine shape =", self.samplingFine.shape)
self.upxFine = self.xFine + self.K @ (np.mean(self.y, axis=1).reshape(-1, 1) - self.samplingFine)
printCones("upxFine shape =", self.upxFine.shape)
return
[docs]
def updateStateFine(self):
"""
Sets the updated state array for the fine Member
:returns: None
"""
self.upStateFine = self.upxFine[:self.state.shape[0]]
return
[docs]
def updateParamsFine(self):
"""
Sets the updated parameters array for the fine Member
"""
self.upParamsFine = self.upxFine[self.state.shape[0]:]
return
[docs]
def mergeUpdate(self):
"""
Concatenate the updated fine member with the updated coarse ensemble
into a single upx array
:returns: None
"""
self.upx = np.hstack((self.upxFine, self.upx))
return
[docs]
def small_mgenkf_report(self):
"""
Prints to stdout a smaller report of the Multi-Grid EnKF
:returns: None
"""
print("\n================ MGENKF REPORT ==================\n")
print("Region: ", self.region.id)
print("Parameters: ", self.paramsFine, self.params)
print("sampling: ", self.samplingFine, self.sampling)
print("Kalman Gain: ", self.K)
print("y", self.y)
print("y-Hx", np.mean(self.y, axis=1).reshape(-1, 1) - self.samplingFine, self.y - self.H)
print("Updated Parameters", self.upParams)
print("avg = ", np.mean(self.upParams))
print("std = ", np.std(self.upParams))
print("RMSE", np.linalg.norm(np.abs(self.y - self.sampling))/np.linalg.norm(self.y))
print("\n===============================================\n")