From a813117f2a16a45ee04bf5f99ce4452798c73875 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Fri, 3 Aug 2012 14:04:13 +0200 Subject: [PATCH] Improving array/matrix treatment and correcting commentaries --- src/daComposant/daAlgorithms/Blue.py | 4 ++-- .../daAlgorithms/LinearLeastSquares.py | 2 ++ .../daAlgorithms/QuantileRegression.py | 9 ++++++++- src/daComposant/daCore/BasicObjects.py | 2 +- src/daComposant/daNumerics/mmqr.py | 18 ++++++++++-------- 5 files changed, 23 insertions(+), 12 deletions(-) diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py index 1f22b5f..42a2e05 100644 --- a/src/daComposant/daAlgorithms/Blue.py +++ b/src/daComposant/daAlgorithms/Blue.py @@ -88,8 +88,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): d = Y - HXb logging.debug("%s Innovation d = %s"%(self._name, d)) # - # Calcul de la matrice de gain dans l'espace le plus petit et de l'analyse - # ------------------------------------------------------------------------ + # Calcul de la matrice de gain et de l'analyse + # -------------------------------------------- if Y.size <= Xb.size: if self._parameters["R_scalar"] is not None: R = self._parameters["R_scalar"] * numpy.eye(len(Y), dtype=numpy.float) diff --git a/src/daComposant/daAlgorithms/LinearLeastSquares.py b/src/daComposant/daAlgorithms/LinearLeastSquares.py index bf8b289..5b8863d 100644 --- a/src/daComposant/daAlgorithms/LinearLeastSquares.py +++ b/src/daComposant/daAlgorithms/LinearLeastSquares.py @@ -52,6 +52,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): else: raise ValueError("Observation error covariance matrix has to be properly defined!") # + # Calcul de la matrice de gain et de l'analyse + # -------------------------------------------- K = (Ha * RI * Hm ).I * Ha * RI Xa = K * Y logging.debug("%s Analyse Xa = %s"%(self._name, Xa)) diff --git a/src/daComposant/daAlgorithms/QuantileRegression.py b/src/daComposant/daAlgorithms/QuantileRegression.py index 78a602f..8e93044 100644 --- a/src/daComposant/daAlgorithms/QuantileRegression.py +++ b/src/daComposant/daAlgorithms/QuantileRegression.py @@ -57,6 +57,12 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): typecast = float, message = "Maximum de variation de la fonction d'estimation lors de l'arrêt", ) + self.defineRequiredParameter( + name = "StoreInternalVariables", + default = False, + typecast = bool, + message = "Stockage des variables internes ou intermédiaires du calcul", + ) def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None): """ @@ -105,7 +111,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): logging.debug("%s CostFunction Jb = %s"%(self._name, Jb)) logging.debug("%s CostFunction Jo = %s"%(self._name, Jo)) logging.debug("%s CostFunction J = %s"%(self._name, J)) - self.StoredVariables["CurrentState"].store( _X.A1 ) + if self._parameters["StoreInternalVariables"]: + self.StoredVariables["CurrentState"].store( _X.A1 ) self.StoredVariables["CostFunctionJb"].store( Jb ) self.StoredVariables["CostFunctionJo"].store( Jo ) self.StoredVariables["CostFunctionJ" ].store( J ) diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index 650ad22..16d44b3 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -99,7 +99,7 @@ class Operator: if self.__Matrix is not None: return self.__Matrix elif ValueForMethodForm is not None: - return self.__Method( ValueForMethodForm ) + return self.__Method( (ValueForMethodForm, None) ) else: raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.") diff --git a/src/daComposant/daNumerics/mmqr.py b/src/daComposant/daNumerics/mmqr.py index dea513b..fa8cb80 100644 --- a/src/daComposant/daNumerics/mmqr.py +++ b/src/daComposant/daNumerics/mmqr.py @@ -40,10 +40,11 @@ def mmqr( # Recuperation des donnees et informations initiales # -------------------------------------------------- variables = asmatrix(x0).A1 - mesures = asmatrix(asmatrix(y).A1).T + mesures = asmatrix(y).flatten().T increment = sys.float_info[0] p = len(variables.flat) n = len(mesures.flat) + quantile = float(quantile) # # Calcul des parametres du MM # --------------------------- @@ -54,7 +55,7 @@ def mmqr( # Calculs d'initialisation # ------------------------ residus = asmatrix( mesures - func( variables ) ).A1 - poids = 1./(epsilon+abs(residus)) + poids = asarray( 1./(epsilon+abs(residus)) ) veps = 1. - 2. * quantile - residus * poids lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus) iteration = 0 @@ -64,14 +65,15 @@ def mmqr( while (increment > toler) and (iteration < maxfun) : iteration += 1 # - Derivees = fprime(variables) - DeriveesT = matrix(Derivees).T - M = - dot( DeriveesT , (array(matrix(p*[poids]).T)*array(Derivees)) ) + Derivees = array(fprime(variables)) + Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS + DeriveesT = array(matrix(Derivees).T) + M = - dot( DeriveesT , (array(matrix(p*[poids,]).T)*Derivees) ) SM = dot( DeriveesT , veps ).T - step = linalg.lstsq( M, SM )[0].A1 + step = linalg.lstsq( M, SM )[0] # - variables = asarray(variables) + asarray(step) - residus = ( mesures - func(variables) ).A1 + variables = variables + step + residus = asmatrix( mesures - func(variables) ).A1 surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus) # while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) : -- 2.39.2