Bases: eoldaslib.eoldas_Operator.Operator
Methods
setup for the differential operator H(x)
A slightly modified J as its efficient to precalculate things for this model
J = 0.5 * x.T D1.T gamma^2 D1 x
A slightly modified J as its efficient to precalculate things for this model
J’ = D.T gamma^2 D x
Calculation of J’‘
We already have the differntial operator self.linear.D1 and self.gamma after we call self.J_prime()
Here, J’’ = D1.T gamma^2 D1
J’ is of shape (nobs,nstates) which is the same as the shape of x
D1 is of shape (nobs,nobs) which needs to be expanded to (nobs,nstates,nobs,nstates)
This is called on initialisation, after data have been read in
Here, we load parameters specifically associated with the model H(x).
In the case of this differential operator, there are:
- model_order : order of the differential operator
- (integer)
- wraparound : edge conditions
- Can be:
- periodic none reflexive
- lag : The (time/space) lag at which
the finite difference is calculated in the differential operator here. If this is 1, then we take the difference between each sample point and its neighbour. This is what we normally use. The main purpose of this mechanism is to allow differences at multiple lags to be calculated (fuller autocorrelation function constraints as in kriging)
Multiple lags can be specified (which you could use to perform kriging), in which case lag weight should also be specified.
- lag_weight : The weight associated with each lag. This will
- generally be decreasing with increasing lag for a ‘usual’ autocorrelation function. There is no point specifying this if only a single lag is specified as the function is normalised.
If the conditions are specified as periodic the period of the function can also be specified, e.g. for time varying data, you could specify 365 for the periodic period.
These are specified in the configuration file as
operator.modelt.rt_model.model_order operator.modelt.rt_model.wraparound operator.modelt.rt_model.lag operator.modelt.rt_model.lag_weight
The default values (set here) are 1, ‘none’, 1 and 1 respectively.
To specify the period for periodic specify e.g.:
[operator.modelt.rt_model] wraparound=periodic,365
The default period is set to 0, which implies that it is periodic on whatever the data extent is.
Or for multiple lags:
[operator.modelt.rt_model] lag=1,2,3,4,5 lag_weight=1,0.7,0.5,0.35,0.2
NB this lag mechanism has not yet been fully tested and should be used with caution. It is intended more as a placeholder for future developments.
Finally, we can also decide to work with inverse gamma (i.e. an uncertainty-based measure)
This is achieved by setting the flag
operator.modelt.rt_model.inverse_gamma=True
This flag should be set if you intend to estimate gamma in the Data Assimilation. Again, the is experimental and should be used with caution.
Here , we use preload_prepare to make sure the x & any y data are gridded for this operator. This greatly simplifies the application of the differential operator.
This method is called before any data are loaded, so ensures they are loaded as a grid.
This method sets up the matrices required for the model.
This operator is written so that it can apply smoothing in different dimensions. This is controlled by the model state vector.
The names of the states are stored in self.x_meta.location and the associated location information in self.x_meta.location. So, we look through these looking for matches, e.g. ‘row’ in location and ‘gamma_row’ in names would mean that we want to apply the model over the row dimension. There should be only one gamma term in the state vectors for this operator. If you give more than one, only the last one will be used.
NOT YET IMPLEMENTED: The model can be applied to multiple dimensions by specifying e.g. gamma_time_row. If you want separate gammas for e.g. time and row, then you should use separate operators. If gamma_roccol is specified, then the model applies to Euclidean distance in row/col space.
Formally, the problem can be stated most simply as a matrix D so that gamma D x is the rate of change of x with respect to the target location variable (time, row, col etc). The job of this method then is to form and store D.
The main complication to this is we have to split up x into those terms that we will apply D to (x2 here) and separately pull out the gamma terms. The resultant matrix D then needs to be re-formed so as to apply to the whole vector x, rather than just x2. We do this with masks.
On input, x is a 1D vector.