Using low-level forecast/analysis steps

Work in progress

This section should be expanded. Give a complete snippet ?

AugKF

We present here how to use the forecast and analysis steps of the AugKF algorithm.

Forecast/Analyse a CoreState

Principle

An algorithm is implemented through a AugkfAlgo object. It contains all the informations needed for the algorithm (configuration, priors, correlation matrices, drift_matrices). It uses a method parallel_forecast_step to forecast CoreStates that is, under the hood, delegated to a AugkfForecaster that is an object created at the initialisation of AugkfAlgo.

Similarly, the AugkfAlgo also implements analysis_step that is delegated to an AugkfAnalyser object, that also stores the observations needed for the analysis_step.

The easiest way to set up these objects is to use the an create_augkf_algo function from a Config that will take care of initialising all objects and data needed:

from pygeodyn.inout.config import ComputationConfig
from pygeodyn.augkf.algo import create_augkf_algo

config = ComputationConfig(0,'augkf.conf') # Or whatever config path
nb_realisations = 10
nb_times = 3
seed = 1
attributed_models = np.array(list(range(nb_realisations)))
algo = create_augkf_algo(config, nb_realisations, seed, attributed_models)

Then, we can create CoreState objects to perform the forecast/analysis on them.

AR1 Example

Both parallel_forecast_step and analysis_step take a 2D CoreState as input (dim: nb_reals x nb_coefs).

The following code shows an example of a forecast step followed by an analysis step using an AR1 process:

# Build a complete CoreState with 10 realisations, 3 times and the number of coefficients given in the config.

# You can use your own corestate initialization or use the algo.init_corestate() methode provided

###################################################
#PROVIDED#

cs, _, _, _, Z_AR3 = algo.init_corestates()

###################################################
#OWN#

cs = CoreState()
corestate_measures={}
corestate_measures.update({'MF': self.config.Nb,
                            'U': self.config.Nu2,
                            'SV': self.config.Nsv,
                            'ER': self.config.Nsv,
                            'Z': self.config.Nz})

# Do whatever operation to initialise the CoreState at t=0 (ex: set all realisations of B according to a magnetic model)
cs.B[:, 0] = my_awesome_magnetic_model
# Z_AR3 is None if AR1
Z_AR3 = None

###################################################
#check if observation is found

algo.analyser.check_if_analysis_data(0)

###################################################
#FORECAST

# forecast to t=1
cs[:, 1] = algo.forecaster.parallel_forecast_step(cs[:,0], seed, 1)

##################################################
#ANALYSIS

# Now do an analysis step on the forecast CoreState at t=1
# if AR1: requires 1 observation date at t=1 for analysis
analysis_result = algo.analyser.analysis_step(cs[:, 1])

# The result is a CoreState containing all analysed realisations that we can again use to update our CoreState (stored here at t=2)
cs[:, 2] = analysis_result

# Then we can access the data as previously:
cs.B[:, 0, 0] # Initial value of g10 for all realisations
cs.B[:, 1, 0] # Forecast value of g10 for all realisations
cs.U[0, 2] # Analysed values of all core flow coefs for the realisation 0
# and so on...

AR3 Example

parallel_forecast_step takes a 2D CoreState as input (dim: nb_reals x nb_coefs).

analysis_step takes a 3D CoreState as input (dim: nb_reals x 2*dt_a_f ratio + 1 x nb_coefs).

The following code shows an example of a forecast step followed by an analysis step using an AR3 process:

# Build a complete CoreState with 10 realisations, 3 times and the number of coefficients given in the config.

# You can use your own corestate initialization or use the algo.init_corestate() methode provided

###################################################
#PROVIDED#

cs, _, _, _, Z_AR3 = algo.init_corestates()

###################################################
#OWN#

cs = CoreState()
corestate_measures={}
corestate_measures.update({'MF': self.config.Nb,
                            'U': self.config.Nu2,
                            'SV': self.config.Nsv,
                            'ER': self.config.Nsv,
                            'Z': self.config.Nz,
                            'dUdt': self.config.Nu2,
                            'd2Udt2': self.config.Nu2,
                            'dERdt': self.config.Nsv,
                            'd2ERdt2': self.config.Nsv})

# Do whatever operation to initialise the CoreState at t=0 (ex: set all realisations of B according to a magnetic model)
cs.B[:, 0] = my_awesome_magnetic_model
# init Z_AR3 at t=0 (Z_AR3 stores states at 3 last times. dimension: nb_reals x 3 x nb_coefs)
Z_AR3 = my_init_Z_AR3

###################################################
#check if observation is found

algo.analyser.check_if_analysis_data(0)

###################################################
#FORECAST

# forecast to t=2
cs[:, 1], Z_AR3 = algo.forecaster.parallel_forecast_step(cs[:, 0], Z_AR3, seed, 1)
cs[:, 2], Z_AR3 = algo.forecaster.parallel_forecast_step(cs[:, 1], Z_AR3, seed, 1)

##################################################
#ANALYSIS

# Now do an analysis step on the forecast CoreState at t=1
# if AR3: requires 3 observation dates at t=0,1,2 for analysis
analysis_result, Z_AR3 = algo.analyser.analysis_step(cs, Z_AR3)

# The result is a CoreState containing all analysed realisations that we can again use to update our CoreState (stored here at t=2)
cs[:, 2] = analysis_result

# Then we can access the data as previously:
cs.B[:, 0, 0] # Initial value of g10 for all realisations
cs.B[:, 1, 0] # Forecast value of g10 for all realisations
cs.U[0, 2] # Analysed values of all core flow coefs for the realisation 0
# and so on...

More information

See the API reference of the Augkf objects:

Common operations are implemented in the motherclass Computer from which AugkfForecaster and AugkfAnalyser inherit. The version of the algorithm with a PCA is implemented in PCA-specific objects (AugkfAlgoPCA, AugkfForecasterPCA, AugkfAnalyserPCA).