From 018c6093cf0f1f41fb2a99f3c32cfa4c7cdb9d12 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Thu, 30 Sep 2021 20:40:33 +0200 Subject: [PATCH] Minor internal clarification and performance improvements --- src/daComposant/daAlgorithms/FunctionTest.py | 4 +- src/daComposant/daAlgorithms/ObserverTest.py | 2 +- src/daComposant/daCore/Aidsm.py | 4 +- src/daComposant/daCore/BasicObjects.py | 2 +- src/daComposant/daCore/NumericObjects.py | 76 ++++++++++---------- 5 files changed, 44 insertions(+), 44 deletions(-) diff --git a/src/daComposant/daAlgorithms/FunctionTest.py b/src/daComposant/daAlgorithms/FunctionTest.py index c7d622b..8da7096 100644 --- a/src/daComposant/daAlgorithms/FunctionTest.py +++ b/src/daComposant/daAlgorithms/FunctionTest.py @@ -97,7 +97,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): msgs += (" -----------------------------\n") msgs += (" Characteristics of input vector X, internally converted:\n") msgs += (" Type...............: %s\n")%type( Xn ) - msgs += (" Lenght of vector...: %i\n")%max(numpy.matrix( Xn ).shape) + msgs += (" Lenght of vector...: %i\n")%max(numpy.ravel( Xn ).shape) msgs += (" Minimum value......: %."+str(_p)+"e\n")%numpy.min( Xn ) msgs += (" Maximum value......: %."+str(_p)+"e\n")%numpy.max( Xn ) msgs += (" Mean of vector.....: %."+str(_p)+"e\n")%numpy.mean( Xn, dtype=mfp ) @@ -131,7 +131,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): msgs = ("===> Information after evaluation:\n") msgs += ("\n Characteristics of simulated output vector Y=H(X), to compare to others:\n") msgs += (" Type...............: %s\n")%type( Yn ) - msgs += (" Lenght of vector...: %i\n")%max(numpy.matrix( Yn ).shape) + msgs += (" Lenght of vector...: %i\n")%max(numpy.ravel( Yn ).shape) msgs += (" Minimum value......: %."+str(_p)+"e\n")%numpy.min( Yn ) msgs += (" Maximum value......: %."+str(_p)+"e\n")%numpy.max( Yn ) msgs += (" Mean of vector.....: %."+str(_p)+"e\n")%numpy.mean( Yn, dtype=mfp ) diff --git a/src/daComposant/daAlgorithms/ObserverTest.py b/src/daComposant/daAlgorithms/ObserverTest.py index 39446d2..323370a 100644 --- a/src/daComposant/daAlgorithms/ObserverTest.py +++ b/src/daComposant/daAlgorithms/ObserverTest.py @@ -68,7 +68,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): self.StoredVariables["SigmaObs2"].store( 1. ) self.StoredVariables["SigmaBck2"].store( 1. ) self.StoredVariables["MahalanobisConsistency"].store( 1. ) - self.StoredVariables["SimulationQuantiles"].store( numpy.matrix((__YY,__YY,__YY)) ) + self.StoredVariables["SimulationQuantiles"].store( numpy.array((__YY,__YY,__YY)) ) self.StoredVariables["SimulatedObservationAtBackground"].store( __YY ) self.StoredVariables["SimulatedObservationAtCurrentState"].store( __YY ) self.StoredVariables["SimulatedObservationAtOptimum"].store( __YY ) diff --git a/src/daComposant/daCore/Aidsm.py b/src/daComposant/daCore/Aidsm.py index 1ca48cd..79ddef2 100644 --- a/src/daComposant/daCore/Aidsm.py +++ b/src/daComposant/daCore/Aidsm.py @@ -763,9 +763,9 @@ class Aidsm(object): self.__adaoObject["AlgorithmParameters"].executePythonScheme( self.__adaoObject ) if "UserPostAnalysis" in self.__adaoObject and len(self.__adaoObject["UserPostAnalysis"])>0: self.__objname = self.__retrieve_objname() - __Upa = [str(val).replace("ADD.","self.").replace(self.__objname+".","self.") for val in self.__adaoObject["UserPostAnalysis"]] + __Upa = map(str, self.__adaoObject["UserPostAnalysis"]) __Upa = eval("\n".join(__Upa)) - exec(__Upa, {}, {'self':self}) + exec(__Upa, {}, {'self':self, 'ADD':self, 'case':self, 'adaopy':self, self.__objname:self}) return 0 def __executeYACSScheme(self, FileName=None): diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index bb0bec8..aba702c 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -226,7 +226,7 @@ class Operator(object): else: if self.__Matrix is not None: self.__addOneMatrixCall() - _xv = numpy.matrix(numpy.ravel(xv)).T + _xv = numpy.ravel(xv).reshape((-1,1)) _hv = self.__Matrix * _xv else: self.__addOneMethodCall() diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index c5b6306..85ecff2 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -890,7 +890,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): XEtnnp = [] for point in range(nbSpts): if selfA._parameters["EstimationOf"] == "State": - XEtnnpi = numpy.asmatrix(numpy.ravel( Mm( (Xnp[:,point], Un) ) )).T + XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1)) if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon ! Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape XEtnnpi = XEtnnpi + Cm * Un @@ -899,10 +899,10 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): elif selfA._parameters["EstimationOf"] == "Parameters": # --- > Par principe, M = Id, Q = 0 XEtnnpi = Xnp[:,point] - XEtnnp.append( XEtnnpi ) - XEtnnp = numpy.hstack( XEtnnp ) + XEtnnp.append( numpy.ravel(XEtnnpi).reshape((-1,1)) ) + XEtnnp = numpy.concatenate( XEtnnp, axis=1 ) # - Xncm = numpy.matrix( XEtnnp.getA()*numpy.array(Wm) ).sum(axis=1) + Xncm = ( XEtnnp * Wm ).sum(axis=1) # if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection": Xncm = ApplyBounds( Xncm, selfA._parameters["Bounds"] ) @@ -910,14 +910,14 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): if selfA._parameters["EstimationOf"] == "State": Pnm = Q elif selfA._parameters["EstimationOf"] == "Parameters": Pnm = 0. for point in range(nbSpts): - Pnm += Wc[i] * (XEtnnp[:,point]-Xncm) * (XEtnnp[:,point]-Xncm).T + Pnm += Wc[i] * ((XEtnnp[:,point]-Xncm).reshape((-1,1)) * (XEtnnp[:,point]-Xncm)) # if selfA._parameters["EstimationOf"] == "Parameters" and selfA._parameters["Bounds"] is not None: Pnmdemi = selfA._parameters["Reconditioner"] * numpy.real(scipy.linalg.sqrtm(Pnm)) else: Pnmdemi = numpy.real(scipy.linalg.sqrtm(Pnm)) # - Xnnp = numpy.hstack([Xncm, Xncm+Gamma*Pnmdemi, Xncm-Gamma*Pnmdemi]) + Xnnp = numpy.hstack([Xncm.reshape((-1,1)), Xncm.reshape((-1,1))+Gamma*Pnmdemi, Xncm.reshape((-1,1))-Gamma*Pnmdemi]) # if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection": for point in range(nbSpts): @@ -926,27 +926,27 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): Ynnp = [] for point in range(nbSpts): if selfA._parameters["EstimationOf"] == "State": - Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], None) ) )).T + Ynnpi = Hm( (Xnnp[:,point], None) ) elif selfA._parameters["EstimationOf"] == "Parameters": - Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], Un) ) )).T - Ynnp.append( Ynnpi ) - Ynnp = numpy.hstack( Ynnp ) + Ynnpi = Hm( (Xnnp[:,point], Un) ) + Ynnp.append( numpy.ravel(Ynnpi).reshape((-1,1)) ) + Ynnp = numpy.concatenate( Ynnp, axis=1 ) # - Yncm = numpy.matrix( Ynnp.getA()*numpy.array(Wm) ).sum(axis=1) + Yncm = ( Ynnp * Wm ).sum(axis=1) # Pyyn = R Pxyn = 0. for point in range(nbSpts): - Pyyn += Wc[i] * (Ynnp[:,point]-Yncm) * (Ynnp[:,point]-Yncm).T - Pxyn += Wc[i] * (Xnnp[:,point]-Xncm) * (Ynnp[:,point]-Yncm).T + Pyyn += Wc[i] * ((Ynnp[:,point]-Yncm).reshape((-1,1)) * (Ynnp[:,point]-Yncm)) + Pxyn += Wc[i] * ((Xnnp[:,point]-Xncm).reshape((-1,1)) * (Ynnp[:,point]-Yncm)) # - _Innovation = Ynpu - Yncm + _Innovation = Ynpu - Yncm.reshape((-1,1)) if selfA._parameters["EstimationOf"] == "Parameters": if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! _Innovation = _Innovation - Cm * Un # Kn = Pxyn * Pyyn.I - Xn = Xncm + Kn * _Innovation + Xn = Xncm.reshape((-1,1)) + Kn * _Innovation Pn = Pnm - Kn * Pyyn * Kn.T # if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection": @@ -988,7 +988,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): or selfA._toStore("CurrentOptimum") \ or selfA._toStore("APosterioriCovariance"): Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) ) - Jo = float( 0.5 * _Innovation.T * RI * _Innovation ) + Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) ) J = Jb + Jo selfA.StoredVariables["CostFunctionJb"].store( Jb ) selfA.StoredVariables["CostFunctionJo"].store( Jo ) @@ -1300,11 +1300,11 @@ def enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterForm # if U is not None: if hasattr(U,"store") and len(U)>1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # @@ -1442,11 +1442,11 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # if U is not None: if hasattr(U,"store") and len(U)>1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # @@ -4001,51 +4001,51 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): XEtnnp = [] for point in range(nbSpts): if selfA._parameters["EstimationOf"] == "State": - XEtnnpi = numpy.asmatrix(numpy.ravel( Mm( (Xnp[:,point], Un) ) )).T + XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1)) if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon ! Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape XEtnnpi = XEtnnpi + Cm * Un elif selfA._parameters["EstimationOf"] == "Parameters": # --- > Par principe, M = Id, Q = 0 XEtnnpi = Xnp[:,point] - XEtnnp.append( XEtnnpi ) - XEtnnp = numpy.hstack( XEtnnp ) + XEtnnp.append( numpy.ravel(XEtnnpi).reshape((-1,1)) ) + XEtnnp = numpy.concatenate( XEtnnp, axis=1 ) # - Xncm = numpy.matrix( XEtnnp.getA()*numpy.array(Wm) ).sum(axis=1) + Xncm = ( XEtnnp * Wm ).sum(axis=1) # if selfA._parameters["EstimationOf"] == "State": Pnm = Q elif selfA._parameters["EstimationOf"] == "Parameters": Pnm = 0. for point in range(nbSpts): - Pnm += Wc[i] * (XEtnnp[:,point]-Xncm) * (XEtnnp[:,point]-Xncm).T + Pnm += Wc[i] * ((XEtnnp[:,point]-Xncm).reshape((-1,1)) * (XEtnnp[:,point]-Xncm)) # Pnmdemi = numpy.real(scipy.linalg.sqrtm(Pnm)) # - Xnnp = numpy.hstack([Xncm, Xncm+Gamma*Pnmdemi, Xncm-Gamma*Pnmdemi]) + Xnnp = numpy.hstack([Xncm.reshape((-1,1)), Xncm.reshape((-1,1))+Gamma*Pnmdemi, Xncm.reshape((-1,1))-Gamma*Pnmdemi]) # Ynnp = [] for point in range(nbSpts): if selfA._parameters["EstimationOf"] == "State": - Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], None) ) )).T + Ynnpi = Hm( (Xnnp[:,point], None) ) elif selfA._parameters["EstimationOf"] == "Parameters": - Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], Un) ) )).T - Ynnp.append( Ynnpi ) - Ynnp = numpy.hstack( Ynnp ) + Ynnpi = Hm( (Xnnp[:,point], Un) ) + Ynnp.append( numpy.ravel(Ynnpi).reshape((-1,1)) ) + Ynnp = numpy.concatenate( Ynnp, axis=1 ) # - Yncm = numpy.matrix( Ynnp.getA()*numpy.array(Wm) ).sum(axis=1) + Yncm = ( Ynnp * Wm ).sum(axis=1) # Pyyn = R Pxyn = 0. for point in range(nbSpts): - Pyyn += Wc[i] * (Ynnp[:,point]-Yncm) * (Ynnp[:,point]-Yncm).T - Pxyn += Wc[i] * (Xnnp[:,point]-Xncm) * (Ynnp[:,point]-Yncm).T + Pyyn += Wc[i] * ((Ynnp[:,point]-Yncm).reshape((-1,1)) * (Ynnp[:,point]-Yncm)) + Pxyn += Wc[i] * ((Xnnp[:,point]-Xncm).reshape((-1,1)) * (Ynnp[:,point]-Yncm)) # - _Innovation = Ynpu - Yncm + _Innovation = Ynpu - Yncm.reshape((-1,1)) if selfA._parameters["EstimationOf"] == "Parameters": if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! _Innovation = _Innovation - Cm * Un # Kn = Pxyn * Pyyn.I - Xn = Xncm + Kn * _Innovation + Xn = Xncm.reshape((-1,1)) + Kn * _Innovation Pn = Pnm - Kn * Pyyn * Kn.T # Xa = Xn # Pointeurs @@ -4084,7 +4084,7 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): or selfA._toStore("CurrentOptimum") \ or selfA._toStore("APosterioriCovariance"): Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) ) - Jo = float( 0.5 * _Innovation.T * RI * _Innovation ) + Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) ) J = Jb + Jo selfA.StoredVariables["CostFunctionJb"].store( Jb ) selfA.StoredVariables["CostFunctionJo"].store( Jo ) -- 2.39.2