import numpy as np
import os
import itertools
from scipy.spatial.distance import cdist
[docs]
class conesRegion():
"""
A class to define a Data Assimilation region
"""
def __init__(self):
self.id = int()
self.ranks = []
self.globalCells = []
self.localCells = []
self.cell_coordinates = []
self.state = np.empty((0, 0), dtype="float")
self.stateNvar = int()
self.params = np.empty([0, 0], dtype="float")
self.sampling = np.empty([0, 0], dtype="float")
self.obsList = list()
self.current_time = float()
self.inflationType = str()
self.stateInflation = float()
self.parametersInflation = float()
self.obsType = "instantaneous"
self.obsAverageWIndow = int()
self.detlaTSim = float()
[docs]
def setRanks(self):
"""
Set the number of ranks that spans the region
:returns: None
"""
for i, a in enumerate(self.globalCells):
if len(a) > 0:
self.ranks.append(i)
return
[docs]
def addObs(self, obs):
"""
Method to include an observation to the region
:param obs: The observation
:type obs: conesStaticObservation
:returns: None
"""
self.obsList.append(obs)
return
[docs]
def setState(self, state):
"""
Set the state vector for the region
:param state: State Matrix
:type state: numpy ndarray
:returns: None
"""
self.state = state
return
[docs]
def splitState(self):
"""
Split the region's stacked state matrix back into one array per state variable
:returns: List of the state arrays, one per state variable
:rtype: list
"""
splitState = np.split(self.state, self.stateNvar)
return splitState
[docs]
def filterState(self, state, stateCells):
"""
From the complete state matrix, this function filters the values belonging to the region
:param state: The raw state matrix gathered
:type state: numpy ndarray
:param stateCells: the list of cells id of the state
:type stateCells: list
:returns: None
"""
regionCells = [x for xs in self.globalCells for x in xs]
stateCells = [x for xs in stateCells for x in xs]
idx = np.where(np.isin(stateCells, regionCells))[0]
self.state = np.zeros((state.shape[0], len(idx), state[0].shape[1]))
stateNvar = 0
for ii, var in enumerate(state):
self.state[ii] = var[idx]
stateNvar = ii
self.stateNvar = stateNvar+1
self.state = np.vstack(self.state)
return
[docs]
def filterSampling(self, sampling, sampling_order):
"""
From the complete sampling matrix, this function sorts and filters the sampled values for the region
:param sampling: The raw sampling matrix gathered
:type sampling: numpy ndarray
:param sampling_order: A list of observation id ordered by gather
:type sampling_order: int
:returns: None
"""
n_obs = len(self.obsList)
n_var = sum(1 for s in self.obsList[0].var_list if not s.endswith("_std"))
filtered_sampling = np.zeros((n_obs*n_var, sampling.shape[1]))
global_order_map = {obs_id: i for i, obs_id in enumerate(sampling_order)}
for i_local, obs in enumerate(self.obsList):
i_global = global_order_map[obs.id]
start_local = i_local * n_var
end_local = (i_local + 1) * n_var
start_global = i_global * n_var
end_global = (i_global + 1) * n_var
filtered_sampling[start_local:end_local, :] = sampling[start_global:end_global, :]
self.sampling = filtered_sampling
return
[docs]
def setParams(self, params):
"""
Set the parameters associated to the region
:param params: array of parameters
:type params: numpy ndarray
:returns: None
"""
self.params = params
return
[docs]
def set_current_time(self, time):
"""
Set the current simulation time associated with the region's data
:param time: The current simulation time
:type time: float
:returns: None
"""
self.current_time = time
return
[docs]
def setInflation(self, conesSettings):
"""
Set the hyperparameters related to Inflation
:param conesSettings: the cones object containing the settings of the calculation
:type conesSettings: conesSettings class
:return: None
"""
self.inflationType = conesSettings.inflationType
self.stateInflation = conesSettings.stateInflation
self.parametersInflation = conesSettings.parametersInflation
[docs]
def setObsSettings(self, conesSettings):
"""
Copy the observation-operator settings from the cones settings object
"""
self.obsType = getattr(conesSettings, "obsType", "instantaneous")
self.obsAverageWindow = getattr(conesSettings, "obsAverageWindow", conesSettings.obsWindow)
self.deltaTSim = conesSettings.deltaTSim
[docs]
def set_cells_coordinates(self, mesh_dir):
"""
Read the region's cell center coordinates from an OpenFOAM 'C' field,
used for state covariance localisation
:param mesh_dir: Path to the OpenFOAM case directory
:type mesh_dir: str
:raises FileNotFoundError: if the cell-centres file was not generated
:returns: None
"""
file = os.path.join(mesh_dir, '0/C')
if not os.path.exists(file):
raise FileNotFoundError(
f"Cells coordinate file {file} not found. Please, run postProcess -func writeCellCentres before"
)
with open(file, 'r') as f:
lines = f.readlines()
ncells = int(lines[19].strip())
coord = np.zeros((ncells, 3))
for ir, rank in enumerate(self.globalCells):
self.cell_coordinates.append([])
for ic, cell in enumerate(rank):
for j in range(0, 3):
coord[ic][j] = float(lines[21+cell].strip().replace('(', '').replace(')', '').split(" ")[j])
self.cell_coordinates[ir].append(coord[ic])
return
[docs]
def get_distance_matrix(self):
"""
Compute the euclidean distance matrix between the region's state cells
and its observations, used for state covaraince localisation
:returns: Distance matrix, shaped (n_state_cells, n_observations)
:rtype: numpy ndarray
"""
s = list(itertools.chain.from_iterable(self.cell_coordinates))
o = [obs.coordinates for obs in self.obsList]
d = cdist(s, o, metric='euclidean')
return d
[docs]
def report(self):
"""
Prints a report of the region to stdout
"""
print("========REGION REPORT========")
print("Region", self.id)
print("Cells per proc:", [len(x) for x in self.localCells])
print("Observations", len(self.obsList))
print("Ranks", self.ranks)
print("State", self.state)
print("Params", self.params)
# for obs in self.obsList:
# obs.report()
print("============================")
[docs]
def full_report(self):
"""
Prints an extended report of the region to stdout
"""
print("========REGION REPORT========")
print("Region", self.id)
print("Global Cell list per proc:",
len([cell for cells in self.globalCells for cell in cells]), self.globalCells)
print("Local Cell list per proc", self.localCells)
print("Observations", len(self.obsList))
print("Ranks", self.ranks)
print("State", self.state)
print("Params", self.params)
for obs in self.obsList:
obs.report()
print("============================")