From a2a380b643be9798cf206dcee83433eb56b5ed0e Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Tue, 16 Jul 2019 14:47:41 +0200 Subject: [PATCH] Various file reorganization and call fixes --- .../daAlgorithms/EnsembleKalmanFilter.py | 4 +- .../daAlgorithms/QuantileRegression.py | 5 +- src/daComposant/daAlgorithms/mmqr.py | 108 ---- src/daComposant/daCore/BasicObjects.py | 111 ++-- src/daComposant/daCore/NumericObjects.py | 508 ++++++++++++++++++ 5 files changed, 573 insertions(+), 163 deletions(-) delete mode 100644 src/daComposant/daAlgorithms/mmqr.py create mode 100644 src/daComposant/daCore/NumericObjects.py diff --git a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py index 13b5f4f..1f89857 100644 --- a/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py +++ b/src/daComposant/daAlgorithms/EnsembleKalmanFilter.py @@ -125,8 +125,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): or self._toStore("APosterioriCovariance"): BI = B.getI() RI = R.getI() - BIdemi = B.choleskyI() - RIdemi = R.choleskyI() + # BIdemi = B.choleskyI() + # RIdemi = R.choleskyI() # # Initialisation # -------------- diff --git a/src/daComposant/daAlgorithms/QuantileRegression.py b/src/daComposant/daAlgorithms/QuantileRegression.py index 19d637b..4b0dd22 100644 --- a/src/daComposant/daAlgorithms/QuantileRegression.py +++ b/src/daComposant/daAlgorithms/QuantileRegression.py @@ -21,7 +21,7 @@ # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D import logging -from daCore import BasicObjects +from daCore import BasicObjects, NumericObjects import numpy # ============================================================================== @@ -144,8 +144,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # Minimisation de la fonctionnelle # -------------------------------- if self._parameters["Minimizer"] == "MMQR": - import mmqr - Minimum, J_optimal, Informations = mmqr.mmqr( + Minimum, J_optimal, Informations = NumericObjects.mmqr( func = CostFunction, x0 = Xini, fprime = GradientOfCostFunction, diff --git a/src/daComposant/daAlgorithms/mmqr.py b/src/daComposant/daAlgorithms/mmqr.py deleted file mode 100644 index ce8c7c5..0000000 --- a/src/daComposant/daAlgorithms/mmqr.py +++ /dev/null @@ -1,108 +0,0 @@ -# -*- coding: utf-8 -*- -# -# Copyright (C) 2008-2019 EDF R&D -# -# This library is free software; you can redistribute it and/or -# modify it under the terms of the GNU Lesser General Public -# License as published by the Free Software Foundation; either -# version 2.1 of the License. -# -# This library is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -# Lesser General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this library; if not, write to the Free Software -# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -# -# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com -# -# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D - -__doc__ = """ - Implémentation informatique de l'algorithme MMQR, basée sur la publication : - David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm", - Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000. -""" -__author__ = "Jean-Philippe ARGAUD" - -import sys, math -from numpy import array, matrix, asarray, asmatrix -from numpy import sum, dot, linalg, ravel, max, min, hstack, argmin, argmax - -# ============================================================================== -def mmqr( - func = None, - x0 = None, - fprime = None, - bounds = None, - quantile = 0.5, - maxfun = 15000, - toler = 1.e-06, - y = None, - ): - # - # Recuperation des donnees et informations initiales - # -------------------------------------------------- - variables = asmatrix(ravel( x0 )) - mesures = asmatrix(ravel( y )).T - increment = sys.float_info[0] - p = len(variables.flat) - n = len(mesures.flat) - quantile = float(quantile) - # - # Calcul des parametres du MM - # --------------------------- - tn = float(toler) / n - e0 = -tn / math.log(tn) - epsilon = (e0-tn)/(1+math.log(e0)) - # - # Calculs d'initialisation - # ------------------------ - residus = ravel( mesures - func( variables ) ) - poids = asarray( 1./(epsilon+abs(residus)) ) - veps = 1. - 2. * quantile - residus * poids - lastsurrogate = -sum(residus*veps) - (1.-2.*quantile)*sum(residus) - iteration = 0 - # - # Recherche iterative - # ------------------- - while (increment > toler) and (iteration < maxfun) : - iteration += 1 - # - 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, rcond=-1 )[0] - # - variables = variables + step - if bounds is not None: - while( (variables < ravel(asmatrix(bounds)[:,0])).any() or (variables > ravel(asmatrix(bounds)[:,1])).any() ): - step = step/2. - variables = variables - step - residus = ravel( mesures - func(variables) ) - surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus) - # - while ( (surrogate > lastsurrogate) and ( max(list(abs(step))) > 1.e-16 ) ) : - step = step/2. - variables = variables - step - residus = ravel( mesures-func(variables) ) - surrogate = sum(residus**2 * poids) + (4.*quantile-2.) * sum(residus) - # - increment = lastsurrogate-surrogate - poids = 1./(epsilon+abs(residus)) - veps = 1. - 2. * quantile - residus * poids - lastsurrogate = -sum(residus * veps) - (1.-2.*quantile)*sum(residus) - # - # Mesure d'écart : q*Sum(residus)-sum(residus negatifs) - # ---------------- - Ecart = quantile * sum(residus) - sum( residus[residus<0] ) - # - return variables, Ecart, [n,p,iteration,increment,0] - -# ============================================================================== -if __name__ == "__main__": - print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index 592c56f..5ac2e62 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -32,11 +32,8 @@ import logging import copy import numpy from functools import partial -from daCore import Persistence -from daCore import PlatformInfo -from daCore import Interfaces +from daCore import Persistence, PlatformInfo, Interfaces from daCore import Templates -from daCore.Interfaces import ImportFromScript, ImportFromFile # ============================================================================== class CacheManager(object): @@ -422,16 +419,16 @@ class FullOperator(object): if asScript is not None: __Matrix, __Function = None, None if asMatrix: - __Matrix = ImportFromScript(asScript).getvalue( self.__name ) + __Matrix = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) elif asOneFunction: - __Function = { "Direct":ImportFromScript(asScript).getvalue( "DirectOperator" ) } + __Function = { "Direct":Interfaces.ImportFromScript(asScript).getvalue( "DirectOperator" ) } __Function.update({"useApproximatedDerivatives":True}) __Function.update(__Parameters) elif asThreeFunctions: __Function = { - "Direct" :ImportFromScript(asScript).getvalue( "DirectOperator" ), - "Tangent":ImportFromScript(asScript).getvalue( "TangentOperator" ), - "Adjoint":ImportFromScript(asScript).getvalue( "AdjointOperator" ), + "Direct" :Interfaces.ImportFromScript(asScript).getvalue( "DirectOperator" ), + "Tangent":Interfaces.ImportFromScript(asScript).getvalue( "TangentOperator" ), + "Adjoint":Interfaces.ImportFromScript(asScript).getvalue( "AdjointOperator" ), } __Function.update(__Parameters) else: @@ -491,8 +488,8 @@ class FullOperator(object): if "withLenghtOfRedundancy" not in __Function: __Function["withLenghtOfRedundancy"] = -1 if "NumberOfProcesses" not in __Function: __Function["NumberOfProcesses"] = None if "withmfEnabled" not in __Function: __Function["withmfEnabled"] = inputAsMF - from daNumerics.ApproximatedDerivatives import FDApproximation - FDA = FDApproximation( + from daCore import NumericObjects + FDA = NumericObjects.FDApproximation( Function = __Function["Direct"], centeredDF = __Function["CenteredFiniteDifference"], increment = __Function["DifferentialIncrement"], @@ -520,7 +517,7 @@ class FullOperator(object): self.__FO["Adjoint"] = Operator( fromMatrix = __matrice.T, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF ) del __matrice else: - raise ValueError("Improperly defined operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.") + raise ValueError("The %s object is improperly defined or undefined, it requires at minima either a matrix, a Direct operator for approximate derivatives or a Tangent/Adjoint operators pair. Please check your operator input."%self.__name) # if __appliedInX is not None: self.__FO["AppliedInX"] = {} @@ -541,11 +538,11 @@ class FullOperator(object): def __repr__(self): "x.__repr__() <==> repr(x)" - return repr(self.__V) + return repr(self.__FO) def __str__(self): "x.__str__() <==> str(x)" - return str(self.__V) + return str(self.__FO) # ============================================================================== class Algorithm(object): @@ -613,6 +610,8 @@ class Algorithm(object): self.__required_parameters = {} self.__required_inputs = {"RequiredInputValues":{"mandatory":(), "optional":()}} self.__variable_names_not_public = {"nextStep":False} # Duplication dans AlgorithmAndParameters + self.__canonical_parameter_name = {} # Correspondance "lower"->"correct" + self.__canonical_stored_name = {} # Correspondance "lower"->"correct" # self.StoredVariables = {} self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations") @@ -653,6 +652,9 @@ class Algorithm(object): self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState") self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum") self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles") + # + for k in self.StoredVariables: + self.__canonical_stored_name[k.lower()] = k def _pre_run(self, Parameters, Xb=None, Y=None, R=None, B=None, Q=None ): "Pré-calcul" @@ -750,13 +752,16 @@ class Algorithm(object): des classes de persistance. """ if key is not None: - return self.StoredVariables[key] + return self.StoredVariables[self.__canonical_stored_name[key.lower()]] else: return self.StoredVariables def __contains__(self, key=None): "D.__contains__(k) -> True if D has a key k, else False" - return key in self.StoredVariables + if key is None or key.lower() not in self.__canonical_stored_name: + return False + else: + return self.__canonical_stored_name[key.lower()] in self.StoredVariables def keys(self): "D.keys() -> list of D's keys" @@ -767,8 +772,8 @@ class Algorithm(object): def pop(self, k, d): "D.pop(k[,d]) -> v, remove specified key and return the corresponding value" - if hasattr(self, "StoredVariables"): - return self.StoredVariables.pop(k, d) + if hasattr(self, "StoredVariables") and k.lower() in self.__canonical_stored_name: + return self.StoredVariables.pop(self.__canonical_stored_name[k.lower()], d) else: try: msg = "'%s'"%k @@ -803,6 +808,7 @@ class Algorithm(object): "listval" : listval, "message" : message, } + self.__canonical_parameter_name[name.lower()] = name logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name)) def getRequiredParameters(self, noDetails=True): @@ -819,11 +825,12 @@ class Algorithm(object): """ Renvoie la valeur d'un paramètre requis de manière contrôlée """ - default = self.__required_parameters[name]["default"] - typecast = self.__required_parameters[name]["typecast"] - minval = self.__required_parameters[name]["minval"] - maxval = self.__required_parameters[name]["maxval"] - listval = self.__required_parameters[name]["listval"] + __k = self.__canonical_parameter_name[name.lower()] + default = self.__required_parameters[__k]["default"] + typecast = self.__required_parameters[__k]["typecast"] + minval = self.__required_parameters[__k]["minval"] + maxval = self.__required_parameters[__k]["maxval"] + listval = self.__required_parameters[__k]["listval"] # if value is None and default is None: __val = None @@ -832,19 +839,23 @@ class Algorithm(object): else: __val = typecast( default ) else: if typecast is None: __val = value - else: __val = typecast( value ) + else: + try: + __val = typecast( value ) + except: + raise ValueError("The value '%s' for the parameter named '%s' can not be correctly evaluated with type '%s'."%(value, __k, typecast)) # if minval is not None and (numpy.array(__val, float) < minval).any(): - raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval)) + raise ValueError("The parameter named '%s' of value '%s' can not be less than %s."%(__k, __val, minval)) if maxval is not None and (numpy.array(__val, float) > maxval).any(): - raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval)) + raise ValueError("The parameter named '%s' of value '%s' can not be greater than %s."%(__k, __val, maxval)) if listval is not None: if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple): for v in __val: if v not in listval: - raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval)) + raise ValueError("The value '%s' is not allowed for the parameter named '%s', it has to be in the list %s."%(v, __k, listval)) elif __val not in listval: - raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval)) + raise ValueError("The value '%s' is not allowed for the parameter named '%s', it has to be in the list %s."%( __val, __k,listval)) return __val def requireInputArguments(self, mandatory=(), optional=()): @@ -890,7 +901,7 @@ class AlgorithmAndParameters(object): self.updateParameters( asDict, asScript ) # if asAlgorithm is None and asScript is not None: - __Algo = ImportFromScript(asScript).getvalue( "Algorithm" ) + __Algo = Interfaces.ImportFromScript(asScript).getvalue( "Algorithm" ) else: __Algo = asAlgorithm # @@ -908,7 +919,7 @@ class AlgorithmAndParameters(object): ): "Mise a jour des parametres" if asDict is None and asScript is not None: - __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" ) + __Dict = Interfaces.ImportFromScript(asScript).getvalue( self.__name, "Parameters" ) else: __Dict = asDict # @@ -1254,12 +1265,12 @@ class RegulationAndParameters(object): self.__P = {} # if asAlgorithm is None and asScript is not None: - __Algo = ImportFromScript(asScript).getvalue( "Algorithm" ) + __Algo = Interfaces.ImportFromScript(asScript).getvalue( "Algorithm" ) else: __Algo = asAlgorithm # if asDict is None and asScript is not None: - __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" ) + __Dict = Interfaces.ImportFromScript(asScript).getvalue( self.__name, "Parameters" ) else: __Dict = asDict # @@ -1321,7 +1332,7 @@ class DataObserver(object): elif (asTemplate is not None) and (asTemplate in Templates.ObserverTemplates): __FunctionText = Templates.ObserverTemplates[asTemplate] elif asScript is not None: - __FunctionText = ImportFromScript(asScript).getstring() + __FunctionText = Interfaces.ImportFromScript(asScript).getstring() else: __FunctionText = "" __Function = ObserverF(__FunctionText) @@ -1378,10 +1389,10 @@ class State(object): contenant des valeurs en colonnes, elles-mêmes nommées "colNames" (s'il n'y a pas de nom de colonne indiquée, on cherche une colonne nommée "name"), on récupère les colonnes et on les range ligne après - ligne (colMajor=False) ou colonne après colonne (colMajor=True). La - variable résultante est de type "asVector" (par défaut) ou - "asPersistentVector" selon que l'une de ces variables est placée à - "True". + ligne (colMajor=False, par défaut) ou colonne après colonne + (colMajor=True). La variable résultante est de type "asVector" (par + défaut) ou "asPersistentVector" selon que l'une de ces variables est + placée à "True". """ self.__name = str(name) self.__check = bool(toBeChecked) @@ -1394,25 +1405,25 @@ class State(object): if asScript is not None: __Vector, __Series = None, None if asPersistentVector: - __Series = ImportFromScript(asScript).getvalue( self.__name ) + __Series = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) else: - __Vector = ImportFromScript(asScript).getvalue( self.__name ) + __Vector = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) elif asDataFile is not None: __Vector, __Series = None, None if asPersistentVector: if colNames is not None: - __Series = ImportFromFile(asDataFile).getvalue( colNames )[1] + __Series = Interfaces.ImportFromFile(asDataFile).getvalue( colNames )[1] else: - __Series = ImportFromFile(asDataFile).getvalue( [self.__name,] )[1] - if bool(colMajor) and not ImportFromFile(asDataFile).getformat() == "application/numpy.npz": + __Series = Interfaces.ImportFromFile(asDataFile).getvalue( [self.__name,] )[1] + if bool(colMajor) and not Interfaces.ImportFromFile(asDataFile).getformat() == "application/numpy.npz": __Series = numpy.transpose(__Series) - elif not bool(colMajor) and ImportFromFile(asDataFile).getformat() == "application/numpy.npz": + elif not bool(colMajor) and Interfaces.ImportFromFile(asDataFile).getformat() == "application/numpy.npz": __Series = numpy.transpose(__Series) else: if colNames is not None: - __Vector = ImportFromFile(asDataFile).getvalue( colNames )[1] + __Vector = Interfaces.ImportFromFile(asDataFile).getvalue( colNames )[1] else: - __Vector = ImportFromFile(asDataFile).getvalue( [self.__name,] )[1] + __Vector = Interfaces.ImportFromFile(asDataFile).getvalue( [self.__name,] )[1] if bool(colMajor): __Vector = numpy.ravel(__Vector, order = "F") else: @@ -1442,7 +1453,7 @@ class State(object): self.shape = (self.shape[0],1) self.size = self.shape[0] * self.shape[1] else: - raise ValueError("The %s object is improperly defined, it requires at minima either a vector, a list/tuple of vectors or a persistent object. Please check your vector input."%self.__name) + raise ValueError("The %s object is improperly defined or undefined, it requires at minima either a vector, a list/tuple of vectors or a persistent object. Please check your vector input."%self.__name) # if scheduledBy is not None: self.__T = scheduledBy @@ -1514,13 +1525,13 @@ class Covariance(object): if asScript is not None: __Matrix, __Scalar, __Vector, __Object = None, None, None, None if asEyeByScalar: - __Scalar = ImportFromScript(asScript).getvalue( self.__name ) + __Scalar = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) elif asEyeByVector: - __Vector = ImportFromScript(asScript).getvalue( self.__name ) + __Vector = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) elif asCovObject: - __Object = ImportFromScript(asScript).getvalue( self.__name ) + __Object = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) else: - __Matrix = ImportFromScript(asScript).getvalue( self.__name ) + __Matrix = Interfaces.ImportFromScript(asScript).getvalue( self.__name ) else: __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject # diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py new file mode 100644 index 0000000..0fc96ef --- /dev/null +++ b/src/daComposant/daCore/NumericObjects.py @@ -0,0 +1,508 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2008-2019 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D + +__doc__ = """ + Définit les versions approximées des opérateurs tangents et adjoints. +""" +__author__ = "Jean-Philippe ARGAUD" + +import os, numpy, time, copy, types, sys, math, logging +from daCore.BasicObjects import Operator +from daCore.PlatformInfo import PlatformInfo +mpr = PlatformInfo().MachinePrecision() +mfp = PlatformInfo().MaximumPrecision() +# logging.getLogger().setLevel(logging.DEBUG) + +# ============================================================================== +def ExecuteFunction( paire ): + assert len(paire) == 2, "Incorrect number of arguments" + X, funcrepr = paire + __X = numpy.asmatrix(numpy.ravel( X )).T + __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"]) + __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), []) + __fonction = getattr(__module,funcrepr["__userFunction__name"]) + sys.path = __sys_path_tmp ; del __sys_path_tmp + __HX = __fonction( __X ) + return numpy.ravel( __HX ) + +# ============================================================================== +class FDApproximation(object): + """ + Cette classe sert d'interface pour définir les opérateurs approximés. A la + création d'un objet, en fournissant une fonction "Function", on obtient un + objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et + "AdjointOperator". On contrôle l'approximation DF avec l'incrément + multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe + "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF + centrées si le booléen "centeredDF" est vrai. + """ + def __init__(self, + Function = None, + centeredDF = False, + increment = 0.01, + dX = None, + avoidingRedundancy = True, + toleranceInRedundancy = 1.e-18, + lenghtOfRedundancy = -1, + mpEnabled = False, + mpWorkers = None, + mfEnabled = False, + ): + if mpEnabled: + try: + import multiprocessing + self.__mpEnabled = True + except ImportError: + self.__mpEnabled = False + else: + self.__mpEnabled = False + self.__mpWorkers = mpWorkers + if self.__mpWorkers is not None and self.__mpWorkers < 1: + self.__mpWorkers = None + logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers)) + # + if mfEnabled: + self.__mfEnabled = True + else: + self.__mfEnabled = False + logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,)) + # + if avoidingRedundancy: + self.__avoidRC = True + self.__tolerBP = float(toleranceInRedundancy) + self.__lenghtRJ = int(lenghtOfRedundancy) + self.__listJPCP = [] # Jacobian Previous Calculated Points + self.__listJPCI = [] # Jacobian Previous Calculated Increment + self.__listJPCR = [] # Jacobian Previous Calculated Results + self.__listJPPN = [] # Jacobian Previous Calculated Point Norms + self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms + else: + self.__avoidRC = False + # + if self.__mpEnabled: + if isinstance(Function,types.FunctionType): + logging.debug("FDA Calculs en multiprocessing : FunctionType") + self.__userFunction__name = Function.__name__ + try: + mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename']) + except: + mod = os.path.abspath(Function.__globals__['__file__']) + if not os.path.isfile(mod): + raise ImportError("No user defined function or method found with the name %s"%(mod,)) + self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','') + self.__userFunction__path = os.path.dirname(mod) + del mod + self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled ) + self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct + elif isinstance(Function,types.MethodType): + logging.debug("FDA Calculs en multiprocessing : MethodType") + self.__userFunction__name = Function.__name__ + try: + mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename']) + except: + mod = os.path.abspath(Function.__func__.__globals__['__file__']) + if not os.path.isfile(mod): + raise ImportError("No user defined function or method found with the name %s"%(mod,)) + self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','') + self.__userFunction__path = os.path.dirname(mod) + del mod + self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled ) + self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct + else: + raise TypeError("User defined function or method has to be provided for finite differences approximation.") + else: + self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled ) + self.__userFunction = self.__userOperator.appliedTo + # + self.__centeredDF = bool(centeredDF) + if abs(float(increment)) > 1.e-15: + self.__increment = float(increment) + else: + self.__increment = 0.01 + if dX is None: + self.__dX = None + else: + self.__dX = numpy.asmatrix(numpy.ravel( dX )).T + logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC) + if self.__avoidRC: + logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP) + + # --------------------------------------------------------- + def __doublon__(self, e, l, n, v=None): + __ac, __iac = False, -1 + for i in range(len(l)-1,-1,-1): + if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]: + __ac, __iac = True, i + if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac)) + break + return __ac, __iac + + # --------------------------------------------------------- + def DirectOperator(self, X ): + """ + Calcul du direct à l'aide de la fonction fournie. + """ + logging.debug("FDA Calcul DirectOperator (explicite)") + if self.__mfEnabled: + _HX = self.__userFunction( X, argsAsSerie = True ) + else: + _X = numpy.asmatrix(numpy.ravel( X )).T + _HX = numpy.ravel(self.__userFunction( _X )) + # + return _HX + + # --------------------------------------------------------- + def TangentMatrix(self, X ): + """ + Calcul de l'opérateur tangent comme la Jacobienne par différences finies, + c'est-à-dire le gradient de H en X. On utilise des différences finies + directionnelles autour du point X. X est un numpy.matrix. + + Différences finies centrées (approximation d'ordre 2): + 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation + dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et + on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi = + H( X_moins_dXi ) + 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par + le pas 2*dXi + 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne + + Différences finies non centrées (approximation d'ordre 1): + 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la + composante X[i] pour composer X_plus_dXi, et on calcule la réponse + HX_plus_dXi = H( X_plus_dXi ) + 2/ On calcule la valeur centrale HX = H(X) + 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par + le pas dXi + 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne + + """ + logging.debug("FDA Début du calcul de la Jacobienne") + logging.debug("FDA Incrément de............: %s*X"%float(self.__increment)) + logging.debug("FDA Approximation centrée...: %s"%(self.__centeredDF)) + # + if X is None or len(X)==0: + raise ValueError("Nominal point X for approximate derivatives can not be None or void (X=%s)."%(str(X),)) + # + _X = numpy.asmatrix(numpy.ravel( X )).T + # + if self.__dX is None: + _dX = self.__increment * _X + else: + _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T + # + if (_dX == 0.).any(): + moyenne = _dX.mean() + if moyenne == 0.: + _dX = numpy.where( _dX == 0., float(self.__increment), _dX ) + else: + _dX = numpy.where( _dX == 0., moyenne, _dX ) + # + __alreadyCalculated = False + if self.__avoidRC: + __bidon, __alreadyCalculatedP = self.__doublon__(_X, self.__listJPCP, self.__listJPPN, None) + __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None) + if __alreadyCalculatedP == __alreadyCalculatedI > -1: + __alreadyCalculated, __i = True, __alreadyCalculatedP + logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i) + # + if __alreadyCalculated: + logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i) + _Jacobienne = self.__listJPCR[__i] + else: + logging.debug("FDA Calcul Jacobienne (explicite)") + if self.__centeredDF: + # + if self.__mpEnabled and not self.__mfEnabled: + funcrepr = { + "__userFunction__path" : self.__userFunction__path, + "__userFunction__modl" : self.__userFunction__modl, + "__userFunction__name" : self.__userFunction__name, + } + _jobs = [] + for i in range( len(_dX) ): + _dXi = _dX[i] + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dXi + _X_moins_dXi = numpy.array( _X.A1, dtype=float ) + _X_moins_dXi[i] = _X[i] - _dXi + # + _jobs.append( (_X_plus_dXi, funcrepr) ) + _jobs.append( (_X_moins_dXi, funcrepr) ) + # + import multiprocessing + self.__pool = multiprocessing.Pool(self.__mpWorkers) + _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs ) + self.__pool.close() + self.__pool.join() + # + _Jacobienne = [] + for i in range( len(_dX) ): + _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) ) + # + elif self.__mfEnabled: + _xserie = [] + for i in range( len(_dX) ): + _dXi = _dX[i] + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dXi + _X_moins_dXi = numpy.array( _X.A1, dtype=float ) + _X_moins_dXi[i] = _X[i] - _dXi + # + _xserie.append( _X_plus_dXi ) + _xserie.append( _X_moins_dXi ) + # + _HX_plusmoins_dX = self.DirectOperator( _xserie ) + # + _Jacobienne = [] + for i in range( len(_dX) ): + _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) ) + # + else: + _Jacobienne = [] + for i in range( _dX.size ): + _dXi = _dX[i] + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dXi + _X_moins_dXi = numpy.array( _X.A1, dtype=float ) + _X_moins_dXi[i] = _X[i] - _dXi + # + _HX_plus_dXi = self.DirectOperator( _X_plus_dXi ) + _HX_moins_dXi = self.DirectOperator( _X_moins_dXi ) + # + _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) ) + # + else: + # + if self.__mpEnabled and not self.__mfEnabled: + funcrepr = { + "__userFunction__path" : self.__userFunction__path, + "__userFunction__modl" : self.__userFunction__modl, + "__userFunction__name" : self.__userFunction__name, + } + _jobs = [] + _jobs.append( (_X.A1, funcrepr) ) + for i in range( len(_dX) ): + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dX[i] + # + _jobs.append( (_X_plus_dXi, funcrepr) ) + # + import multiprocessing + self.__pool = multiprocessing.Pool(self.__mpWorkers) + _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs ) + self.__pool.close() + self.__pool.join() + # + _HX = _HX_plus_dX.pop(0) + # + _Jacobienne = [] + for i in range( len(_dX) ): + _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) ) + # + elif self.__mfEnabled: + _xserie = [] + _xserie.append( _X.A1 ) + for i in range( len(_dX) ): + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dX[i] + # + _xserie.append( _X_plus_dXi ) + # + _HX_plus_dX = self.DirectOperator( _xserie ) + # + _HX = _HX_plus_dX.pop(0) + # + _Jacobienne = [] + for i in range( len(_dX) ): + _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) ) + # + else: + _Jacobienne = [] + _HX = self.DirectOperator( _X ) + for i in range( _dX.size ): + _dXi = _dX[i] + _X_plus_dXi = numpy.array( _X.A1, dtype=float ) + _X_plus_dXi[i] = _X[i] + _dXi + # + _HX_plus_dXi = self.DirectOperator( _X_plus_dXi ) + # + _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) ) + # + # + _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T + if self.__avoidRC: + if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size + while len(self.__listJPCP) > self.__lenghtRJ: + self.__listJPCP.pop(0) + self.__listJPCI.pop(0) + self.__listJPCR.pop(0) + self.__listJPPN.pop(0) + self.__listJPIN.pop(0) + self.__listJPCP.append( copy.copy(_X) ) + self.__listJPCI.append( copy.copy(_dX) ) + self.__listJPCR.append( copy.copy(_Jacobienne) ) + self.__listJPPN.append( numpy.linalg.norm(_X) ) + self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) ) + # + logging.debug("FDA Fin du calcul de la Jacobienne") + # + return _Jacobienne + + # --------------------------------------------------------- + def TangentOperator(self, paire ): + """ + Calcul du tangent à l'aide de la Jacobienne. + """ + if self.__mfEnabled: + assert len(paire) == 1, "Incorrect lenght of arguments" + _paire = paire[0] + assert len(_paire) == 2, "Incorrect number of arguments" + else: + assert len(paire) == 2, "Incorrect number of arguments" + _paire = paire + X, dX = _paire + _Jacobienne = self.TangentMatrix( X ) + if dX is None or len(dX) == 0: + # + # Calcul de la forme matricielle si le second argument est None + # ------------------------------------------------------------- + if self.__mfEnabled: return [_Jacobienne,] + else: return _Jacobienne + else: + # + # Calcul de la valeur linéarisée de H en X appliqué à dX + # ------------------------------------------------------ + _dX = numpy.asmatrix(numpy.ravel( dX )).T + _HtX = numpy.dot(_Jacobienne, _dX) + if self.__mfEnabled: return [_HtX.A1,] + else: return _HtX.A1 + + # --------------------------------------------------------- + def AdjointOperator(self, paire ): + """ + Calcul de l'adjoint à l'aide de la Jacobienne. + """ + if self.__mfEnabled: + assert len(paire) == 1, "Incorrect lenght of arguments" + _paire = paire[0] + assert len(_paire) == 2, "Incorrect number of arguments" + else: + assert len(paire) == 2, "Incorrect number of arguments" + _paire = paire + X, Y = _paire + _JacobienneT = self.TangentMatrix( X ).T + if Y is None or len(Y) == 0: + # + # Calcul de la forme matricielle si le second argument est None + # ------------------------------------------------------------- + if self.__mfEnabled: return [_JacobienneT,] + else: return _JacobienneT + else: + # + # Calcul de la valeur de l'adjoint en X appliqué à Y + # -------------------------------------------------- + _Y = numpy.asmatrix(numpy.ravel( Y )).T + _HaY = numpy.dot(_JacobienneT, _Y) + if self.__mfEnabled: return [_HaY.A1,] + else: return _HaY.A1 + +# ============================================================================== +def mmqr( + func = None, + x0 = None, + fprime = None, + bounds = None, + quantile = 0.5, + maxfun = 15000, + toler = 1.e-06, + y = None, + ): + """ + Implémentation informatique de l'algorithme MMQR, basée sur la publication : + David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm", + Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000. + """ + # + # Recuperation des donnees et informations initiales + # -------------------------------------------------- + variables = numpy.ravel( x0 ) + mesures = numpy.ravel( y ) + increment = sys.float_info[0] + p = variables.size + n = mesures.size + quantile = float(quantile) + # + # Calcul des parametres du MM + # --------------------------- + tn = float(toler) / n + e0 = -tn / math.log(tn) + epsilon = (e0-tn)/(1+math.log(e0)) + # + # Calculs d'initialisation + # ------------------------ + residus = mesures - numpy.ravel( func( variables ) ) + poids = 1./(epsilon+numpy.abs(residus)) + veps = 1. - 2. * quantile - residus * poids + lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus) + iteration = 0 + # + # Recherche iterative + # ------------------- + while (increment > toler) and (iteration < maxfun) : + iteration += 1 + # + Derivees = numpy.array(fprime(variables)) + Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS + DeriveesT = Derivees.transpose() + M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) ) + SM = numpy.transpose(numpy.dot( DeriveesT , veps )) + step = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0] + # + variables = variables + step + if bounds is not None: + while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ): + step = step/2. + variables = variables - step + residus = mesures - numpy.ravel( func(variables) ) + surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus) + # + while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) : + step = step/2. + variables = variables - step + residus = mesures - numpy.ravel( func(variables) ) + surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus) + # + increment = lastsurrogate-surrogate + poids = 1./(epsilon+numpy.abs(residus)) + veps = 1. - 2. * quantile - residus * poids + lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus) + # + # Mesure d'écart + # -------------- + Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] ) + # + return variables, Ecart, [n,p,iteration,increment,0] + +# ============================================================================== +if __name__ == "__main__": + print('\n AUTODIAGNOSTIC\n') -- 2.39.2