pipt.update_schemes.es

  1"""
  2ES type schemes
  3"""
  4from pipt.update_schemes.enkf import enkf_approx
  5from pipt.update_schemes.enkf import enkf_full
  6from pipt.update_schemes.enkf import enkf_subspace
  7
  8import numpy as np
  9from copy import deepcopy
 10from pipt.misc_tools import analysis_tools as at
 11
 12
 13class esMixIn():
 14    """
 15    This is the straightforward ES analysis scheme. We treat this as a all-data-at-once EnKF step, hence the
 16    calc_analysis method here is identical to that in the `enkf` class. Since, for the moment, ASSIMINDEX is parsed in a
 17    specific manner (or more precise, single rows and columns in the PIPT init. file is parsed to a 1D list), a
 18    `Simultaneous` 'loop' had to be implemented, and `es` will use this to do the inversion. Maybe in the future, we can
 19    make the `enkf` class do simultaneous updating also. The consequence of all this is that we inherit BOTH `enkf` and
 20    `Simultaneous` classes, which is convenient. The `Simultaneous` class is inherited to set up the correct inversion
 21    structure and `enkf` is inherited to get `calc_analysis`, so we do not have to implement it again.
 22    """
 23
 24    def __init__(self, keys_da, keys_fwd, sim):
 25        """
 26        The class is initialized by passing the PIPT init. file upwards in the hierarchy to be read and parsed in
 27        `pipt.input_output.pipt_init.ReadInitFile`.
 28
 29        Parameters
 30        ----------
 31        init_file: str
 32            PIPT init. file containing info. to run the inversion algorithm
 33        """
 34        # Pass init. file to Simultaneous parent class (Python searches parent classes from left to right).
 35        super().__init__(keys_da, keys_fwd, sim)
 36
 37        if self.restart is False:
 38            # At the moment, the iterative loop is threated as an iterative smoother an thus we check if assim. indices
 39            # are given as in the Simultaneous loop.
 40            self.check_assimindex_simultaneous()
 41
 42            # Extract no. assimilation steps from MDA keyword in DATAASSIM part of init. file and set this equal to
 43            # the number of iterations pluss one. Need one additional because the iter=0 is the prior run.
 44            self.max_iter = 2
 45
 46    def check_convergence(self):
 47        """
 48        Calculate the "convergence" of the method. Important to
 49        """
 50        self.prev_data_misfit = self.prior_data_misfit
 51        # only calulate for the final (posterior) estimate
 52        if self.iteration == len(self.keys_da['assimindex']):
 53            assim_index = [self.keys_da['obsname'], list(
 54                np.concatenate(self.keys_da['assimindex']))]
 55            list_datatypes = self.list_datatypes
 56            obs_data_vector, pred_data = at.aug_obs_pred_data(self.obs_data, self.pred_data, assim_index,
 57                                                              list_datatypes)
 58
 59            data_misfit = at.calc_objectivefun(
 60                self.full_real_obs_data, pred_data, self.full_cov_data)
 61            self.data_misfit = np.mean(data_misfit)
 62            self.data_misfit_std = np.std(data_misfit)
 63
 64        else:  # sequential updates not finished. Misfit is not relevant
 65            self.data_misfit = self.prior_data_misfit
 66
 67        # Logical variables for conv. criteria
 68        why_stop = {'rel_data_misfit': 1 - (self.data_misfit / self.prev_data_misfit),
 69                    'data_misfit': self.data_misfit,
 70                    'prev_data_misfit': self.prev_data_misfit}
 71
 72        if self.data_misfit == self.prev_data_misfit:
 73            self.logger.info(
 74                f'ES update {self.iteration} complete!')
 75            self.current_state = deepcopy(self.state)
 76        else:
 77            if self.data_misfit < self.prior_data_misfit:
 78                self.logger.info(
 79                    f'ES update complete! Objective function decreased from {self.prior_data_misfit:0.1f} to {self.data_misfit:0.1f}.')
 80            else:
 81                self.logger.info(
 82                    f'ES update complete! Objective function increased from {self.prior_data_misfit:0.1f} to {self.data_misfit:0.1f}.')
 83        # Return conv = False, why_stop var.
 84        return False, True, why_stop
 85
 86
 87class es_approx(esMixIn, enkf_approx):
 88    """
 89    Mixin of ES class and approximate update
 90    """
 91    pass
 92
 93
 94class es_full(esMixIn, enkf_full):
 95    """
 96    mixin of ES class and full update.
 97    Note that since we do not iterate there is no difference between is full and approx.
 98    """
 99    pass
100
101
102class es_subspace(esMixIn, enkf_subspace):
103    """
104    mixin of ES class and subspace update.
105    """
106    pass