GRF package¶

Submodules¶

GRF.GRF module¶

Gaussian Random Field module handles the data assimilation and the cost valley computation.

Objectives:

Construct the Gaussian Random Field (GRF) kernel.
Update the prior mean and covariance matrix.
Compute the cost valley fields.
Assimilate in-situ data.

Methodology:

Construct the GRF kernel.

1.1. Construct the distance matrix using
\[d_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}\]

1.2. Construct the covariance matrix.
\[\Sigma_{ij} = \sigma^2 (1 + \eta d_{ij}) \exp(-\eta d_{ij})\]
Update the prior mean and covariance matrix.
2.1. Update the prior mean. 2.2. Update the prior covariance matrix.
Compute the cost valley fields.
3.1. Compute the expected improvement of best value (EIBV). 3.2. Compute the inverse of the variance reduction (IVR). 3.3. Compute the cost valley fields by summing weighted EIBV and weighted IVR fields.

class GRF.GRF.GRF(sigma: float = 1.0, nugget: float = 0.4, approximate_eibv: bool = True)¶

Bases: object

__Sigma¶: Cost valley

__ar1_corr¶: Empirical parameters

__construct_grf_field() → None¶

Construct distance matrix and thus Covariance matrix for the kernel.

Methodology:

Construct the distance matrix using
\[d_{ij} = \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}\]
Construct the covariance matrix.
\[\Sigma_{ij} = \sigma^2 (1 + \eta d_{ij}) \exp(-\eta d_{ij})\]

__construct_prior_mean() → None¶

Construct prior mean for the kernel.

Methodology:

Construct the prior mean using the SINMOD dataset.
Interpolate the prior mean onto the grid.

Returns: None

__get_eibv(mu: numpy.ndarray, sigma_diag: numpy.ndarray, vr_diag: numpy.ndarray) → float¶

Calculate the eibv using the analytical formula with a bivariate cumulative dentisty function.

Parameters

mu – n x 1 dimension
sigma_diag – n x 1 dimension
vr_diag – n x 1 dimension

Methodology:

Calculate the probability of exceedance of the threshold using a bivariate cumulative dentisty function.
\[p = \Phi(\frac{\theta - \mu}{\sigma}) - \Phi(\frac{\theta - \mu}{\sigma}) \Phi(\frac{\theta - \mu}{\sigma})\]
Calculate eibv by summing up the product of p*(1-p).

Returns: information based on variance reduction, dimension: n x 1
Return type: eibv

Examples

>>> grf = GRF()
>>> eibv = grf.__get_eibv(grf.__mu, grf.__sigma_diag, grf.__vr_diag)

__get_ibv(mu: numpy.ndarray, sigma_diag: numpy.ndarray) → numpy.ndarray¶

Calculate the ibv using the approximate formula with a univariate cumulative dentisty function.

Parameters

mu – n x 1 dimension
sigma_diag – n x 1 dimension

Methodology:

Calculate the probability of exceedance of the threshold.
Calculate ibv by summing up the product of probability of exceedance and (1 - probability of exceedance).

Returns: information based on variance reduction, dimension: n x 1
Return type: ibv

Examples

>>> grf = GRF()
>>> ibv = grf.__get_ibv(grf.__mu, grf.__sigma_diag)

__tau¶: Conditional field

__update(ind_measured: numpy.ndarray, salinity_measured: numpy.ndarray) → None¶

Update GRF kernel based on sampled data.

Parameters

ind_measured – indices where the data is assimilated.
salinity_measured – measurements at sampeld locations, dimension: m x 1

Methodology:

Loop through each measurement and construct the measurement matrix F.
Construct the measurement noise matrix R.
Update the kernel mean and covariance matrix using
\[ \begin{align}\begin{aligned}\mu = \mu + \Sigma F^T (F \Sigma F^T + R)^{-1} (y - F \mu)\\\Sigma = \Sigma - \Sigma F^T (F \Sigma F^T + R)^{-1} F \Sigma\end{aligned}\end{align} \]

Returns: None

Examples

>>> ind_measured = np.array([0, 1])
>>> salinity_measured = np.array([0, 1])
>>> grf = GRF()
>>> grf.__update(ind_measured, salinity_measured)
>>> grf.get_mu()

__update_temporal(ind_measured: numpy.ndarray, salinity_measured: numpy.ndarray, timestep=0) → None¶

Update GRF kernel based on sampled data.

Parameters

ind_measured – indices where the data is assimilated.
salinity_measured – measurements at sampeld locations, dimension: m x 1
timestep – timestep of the data assimilation, default is 0, which means that the data assimilation is performed at the same time step as the model.

Methodology:

Loop through each measurement and construct the measurement matrix F.
Construct the measurement noise matrix R.
Update the kernel mean and covariance matrix using
\[ \begin{align}\begin{aligned}\mu_{t|t-1} = \mu + \rho * (\mu_{t-1|t-1} - \mu)\\\Sigma_{t|t-1} = \rho^2 * \Sigma_{t-1|t-1} + (1 - \rho^2) * \Sigma\\G_t = \Sigma_{t|t-1} F^T (F \Sigma_{t|t-1} F^T + R)^{-1}\\\mu_{t|t} = \mu_{t|t-1} + G_t (y - F \mu_{t|t-1})\\\Sigma_{t|t} = \Sigma_{t|t-1} - G_t F \Sigma_{t|t-1}\end{aligned}\end{align} \]

Returns: None

References

[1] https://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf

Examples

>>> ind_measured = np.array([0, 1])
>>> salinity_measured = np.array([0, 1])
>>> grf = GRF()
>>> grf.__update_temporal(ind_measured, salinity_measured)
>>> grf.get_mu()

assimilate_data(dataset: numpy.ndarray) → None¶

Assimilate dataset to GRF kernel.

Parameters: dataset – np.array([x, y, sal])

Methodology:

Construct the distance matrix between gmrf grid and dataset grid.
Average the values to each cell.
Update the kernel mean and covariance matrix.

Returns: None

Examples

>>> dataset = np.array([[0, 0, 0], [1, 1, 1]])
>>> grf = GRF()
>>> grf.assimilate_data(dataset)
>>> grf.get_mu()

assimilate_temporal_data(dataset: numpy.ndarray) → tuple¶

Assimilate temporal dataset to GRF kernel.

Parameters: dataset – np.array([timestamp, x, y, sal])

Methodology:

Construct the distance matrix between grf grid and dataset grid.
Average the values to each cell.
Update the kernel mean and covariance matrix.

Returns: (ind, salinity) for visualising eda plots.

Examples

>>> dataset = np.array([[0, 0, 0, 0], [1, 1, 1, 1]])
>>> grf = GRF()
>>> grf.assimilate_temporal_data(dataset)
>>> grf.get_mu()

get_covariance_matrix() → numpy.ndarray¶

Return Covariance.

Returns: Covariance matrix
Return type: Sigma

Examples

>>> grf = GRF()
>>> grf.get_covariance_matrix()
array([[1.00000000e+00, 9.99999998e-01, 9.99999994e-01],
       [9.99999998e-01, 1.00000000e+00, 9.99999998e-01],
       [9.99999994e-01, 9.99999998e-01, 1.00000000e+00]])

get_ei_field() → tuple¶

Get expected information field.

Methodology:

Loop through each grid point and calculate the EIBV and IVR fields.
Calculate the EI field using
\[EI = weight * EIBV + (1 - weight) * IVR\]

Returns: expected information based on variance reduction, dimension: Ngrid x 1
Return type: eibv_field

Examples

>>> grf = GRF()
>>> eibv_field = grf.get_ei_field()

get_eibv_field() → numpy.ndarray¶

Return the computed eibv field, given which method to be called.

Returns: eibv field
Return type: eibv_field

Examples

>>> grf = GRF()
>>> grf.get_eibv_field()
array([0.1, 0.2, 0.3])

get_ivr_field() → numpy.ndarray¶

Return the computed ivr field, given which method to be called.

Returns: ivr field
Return type: ivr_field

Examples

>>> grf = GRF()
>>> grf.get_ivr_field()
array([0.1, 0.2, 0.3])

get_lateral_range() → float¶

Return lateral range.

Returns: lateral range
Return type: lateral_range

Examples

>>> grf = GRF()
>>> grf.get_lateral_range()
600.0

get_mu() → numpy.ndarray¶

Return mean vector.

Returns: mean vector
Return type: mu

Examples

>>> grf = GRF()
>>> grf.get_mu()
array([0.1, 0.2, 0.3])

get_nugget() → float¶

Return nugget of the field.

Returns: nugget
Return type: nugget

Examples

>>> grf = GRF()
>>> grf.get_nugget()
0.0

get_sigma() → float¶

Return variability of the field.

Returns: space variability
Return type: sigma

Examples

>>> grf = GRF()
>>> grf.get_sigma()
1.0

get_threshold() → float¶

Return threshold.

Returns: threshold
Return type: threshold

Examples

>>> grf = GRF()
>>> grf.get_threshold()
27.0

set_lateral_range(value: float) → None¶

Set lateral range.

Parameters: value – lateral range

Examples

>>> grf = GRF()
>>> grf.set_lateral_range(0.1)

set_mu(value: numpy.ndarray) → None¶

Set mean of the field.

Parameters: value – mean of the field

Examples

>>> grf = GRF()
>>> grf.set_mu(np.array([0.1, 0.2, 0.3]))

set_nugget(value: float) → None¶

Set nugget.

Parameters: value – nugget

Examples

>>> grf = GRF()
>>> grf.set_nugget(0.1)

set_sigma(value: float) → None¶

Set space variability.

Parameters: value – space variability

Examples

>>> grf = GRF()
>>> grf.set_sigma(0.1)

set_threshold(value: float) → None¶

Set threshold.

Parameters: value – threshold

Examples

>>> grf = GRF()
>>> grf.set_threshold(0.1)

GRF package¶

Submodules¶

GRF.GRF module¶

Module contents¶

Table of Contents

Related Topics

This Page