src/daComposant/daAlgorithms/EnsembleKalmanFilter.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2018 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 import logging
  24 from daCore import BasicObjects, PlatformInfo
  25 import numpy, math
  26 mfp = PlatformInfo.PlatformInfo().MaximumPrecision()
  27
  28 # Using "Ensemble Kalman filtering", L. HOUTEKAMER and HERSCHEL L. MITCHELL, QJRMS (2005), 131, pp. 3269–3289
  29
  30 # ==============================================================================
  31 class ElementaryAlgorithm(BasicObjects.Algorithm):
  32     def __init__(self):
  33         BasicObjects.Algorithm.__init__(self, "ENSEMBLEKALMANFILTER")
  34         self.defineRequiredParameter(
  35             name     = "NumberOfMembers",
  36             default  = 100,
  37             typecast = int,
  38             message  = "Nombre de membres dans l'ensemble",
  39             minval   = -1,
  40             )
  41         self.defineRequiredParameter(
  42             name     = "Minimizer",
  43             default  = "EnKF",
  44             typecast = str,
  45             message  = "Schéma de mise a jour des informations d'ensemble",
  46             listval  = ["EnKF", "ETKF", "DEnKF"],
  47             )
  48         self.defineRequiredParameter(
  49             name     = "ConstrainedBy",
  50             default  = "EstimateProjection",
  51             typecast = str,
  52             message  = "Prise en compte des contraintes",
  53             listval  = ["EstimateProjection"],
  54             )
  55         self.defineRequiredParameter(
  56             name     = "EstimationOf",
  57             default  = "State",
  58             typecast = str,
  59             message  = "Estimation d'etat ou de parametres",
  60             listval  = ["State", "Parameters"],
  61             )
  62         self.defineRequiredParameter(
  63             name     = "StoreInternalVariables",
  64             default  = False,
  65             typecast = bool,
  66             message  = "Stockage des variables internes ou intermédiaires du calcul",
  67             )
  68         self.defineRequiredParameter(
  69             name     = "StoreSupplementaryCalculations",
  70             default  = [],
  71             typecast = tuple,
  72             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
  73             listval  = ["APosterioriCorrelations", "APosterioriCovariance", "APosterioriStandardDeviations", "APosterioriVariances", "BMA", "CurrentState", "CostFunctionJ", "CostFunctionJb", "CostFunctionJo", "Innovation"]
  74             )
  75         self.defineRequiredParameter( # Pas de type
  76             name     = "Bounds",
  77             message  = "Liste des valeurs de bornes",
  78             )
  79         self.defineRequiredParameter(
  80             name     = "SetSeed",
  81             typecast = numpy.random.seed,
  82             message  = "Graine fixée pour le générateur aléatoire",
  83             )
  84         self.requireInputArguments(
  85             mandatory= ("Xb", "Y", "HO", "R", "B"),
  86             optional = ("U", "EM", "CM", "Q"),
  87             )
  88
  89     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
  90         self._pre_run(Parameters, Xb, Y, R, B, Q)
  91         #
  92         if self._parameters["EstimationOf"] == "Parameters":
  93             self._parameters["StoreInternalVariables"] = True
  94         #
  95         # Opérateurs
  96         # ----------
  97         H = HO["Direct"].appliedControledFormTo
  98         #
  99         if self._parameters["EstimationOf"] == "State":
 100             M = EM["Direct"].appliedControledFormTo
 101         #
 102         if CM is not None and "Tangent" in CM and U is not None:
 103             Cm = CM["Tangent"].asMatrix(Xb)
 104         else:
 105             Cm = None
 106         #
 107         # Nombre de pas identique au nombre de pas d'observations
 108         # -------------------------------------------------------
 109         if hasattr(Y,"stepnumber"):
 110             duration = Y.stepnumber()
 111             __p = numpy.cumprod(Y.shape())[-1]
 112         else:
 113             duration = 2
 114             __p = numpy.array(Y).size
 115         #
 116         # Précalcul des inversions de B et R
 117         # ----------------------------------
 118         if self._parameters["StoreInternalVariables"]:
 119             BI = B.getI()
 120             RI = R.getI()
 121         BIdemi = B.choleskyI()
 122         RIdemi = R.choleskyI()
 123         #
 124         # Initialisation
 125         # --------------
 126         __n = Xb.size
 127         __m = self._parameters["NumberOfMembers"]
 128         Xn = numpy.asmatrix(numpy.dot( Xb.reshape(__n,1), numpy.ones((1,__m)) ))
 129         if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
 130         else:                         Pn = B
 131         if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
 132         else:                         Rn = R
 133         if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
 134         else:                         Qn = Q
 135         #
 136         self.StoredVariables["Analysis"].store( Xb.A1 )
 137         if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 138             self.StoredVariables["APosterioriCovariance"].store( Pn )
 139             covarianceXa = Pn
 140         Xa               = Xb
 141         previousJMinimum = numpy.finfo(float).max
 142         #
 143         # Predimensionnement
 144         Xn_predicted = numpy.asmatrix(numpy.zeros((__n,__m)))
 145         HX_predicted = numpy.asmatrix(numpy.zeros((__p,__m)))
 146         #
 147         for step in range(duration-1):
 148             if hasattr(Y,"store"):
 149                 Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
 150             else:
 151                 Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
 152             #
 153             if U is not None:
 154                 if hasattr(U,"store") and len(U)>1:
 155                     Un = numpy.asmatrix(numpy.ravel( U[step] )).T
 156                 elif hasattr(U,"store") and len(U)==1:
 157                     Un = numpy.asmatrix(numpy.ravel( U[0] )).T
 158                 else:
 159                     Un = numpy.asmatrix(numpy.ravel( U )).T
 160             else:
 161                 Un = None
 162             #
 163             if self._parameters["EstimationOf"] == "State":
 164                 for i in range(__m):
 165                     qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn)).T
 166                     Xn_predicted[:,i] = numpy.asmatrix(numpy.ravel( M((Xn[:,i], Un)) )).T + qi
 167                     HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( H((Xn_predicted[:,i], Un)) )).T
 168                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
 169                     Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
 170                     Xn_predicted = Xn_predicted + Cm * Un
 171             elif self._parameters["EstimationOf"] == "Parameters":
 172                 # --- > Par principe, M = Id, Q = 0
 173                 Xn_predicted = Xn
 174             #
 175             Xfm = numpy.asmatrix(numpy.ravel(Xn_predicted.mean(axis=1, dtype=mfp))).T
 176             Hfm = numpy.asmatrix(numpy.ravel(HX_predicted.mean(axis=1, dtype=mfp))).T
 177             Af  = Xn_predicted - Xfm
 178             Hf  = HX_predicted - Hfm
 179             #
 180             PfHT, HPfHT = 0., 0.
 181             for i in range(__m):
 182                 PfHT  += Af[:,i] * Hf[:,i].T
 183                 HPfHT += Hf[:,i] * Hf[:,i].T
 184             PfHT  = (1./(__m-1)) * PfHT
 185             HPfHT = (1./(__m-1)) * HPfHT
 186             #
 187             K = PfHT * ( R + HPfHT ).I
 188             #
 189             Yo = numpy.asmatrix(numpy.zeros((__p,__m)))
 190             for i in range(__m):
 191                 ri = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__p), Rn)).T
 192                 Yo[:,i] = Ynpu + ri
 193             #
 194             for i in range(__m):
 195                 Xn[:,i] = Xn_predicted[:,i] + K * (Yo[:,i] - HX_predicted[:,i])
 196             #
 197             Xa = Xn.mean(axis=1, dtype=mfp)
 198             self.StoredVariables["Analysis"].store( Xa )
 199             #
 200             del Yo, PfHT, HPfHT
 201             if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 202                 Ht = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
 203                 Ht = Ht.reshape(__p,__n) # ADAO & check shape
 204                 Pf = 0.
 205                 for i in range(__m):
 206                     Pf += Af[:,i] * Af[:,i].T
 207                 Pf = (1./(__m-1)) * Pf
 208                 Pn = (1. - K * Ht) * Pf
 209                 self.StoredVariables["APosterioriCovariance"].store( Pn )
 210         #
 211         self._post_run(HO)
 212         return 0
 213
 214 # ==============================================================================
 215 if __name__ == "__main__":
 216     print('\n AUTODIAGNOSTIC \n')