src/daComposant/daCore/NumericObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2019 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 __doc__ = """
  24     Définit les versions approximées des opérateurs tangents et adjoints.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, numpy, time, copy, types, sys, math, logging
  29 from daCore.BasicObjects import Operator
  30 from daCore.PlatformInfo import PlatformInfo
  31 mpr = PlatformInfo().MachinePrecision()
  32 mfp = PlatformInfo().MaximumPrecision()
  33 # logging.getLogger().setLevel(logging.DEBUG)
  34
  35 # ==============================================================================
  36 def ExecuteFunction( paire ):
  37     assert len(paire) == 2, "Incorrect number of arguments"
  38     X, funcrepr = paire
  39     __X = numpy.asmatrix(numpy.ravel( X )).T
  40     __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
  41     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  42     __fonction = getattr(__module,funcrepr["__userFunction__name"])
  43     sys.path = __sys_path_tmp ; del __sys_path_tmp
  44     __HX  = __fonction( __X )
  45     return numpy.ravel( __HX )
  46
  47 # ==============================================================================
  48 class FDApproximation(object):
  49     """
  50     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  51     création d'un objet, en fournissant une fonction "Function", on obtient un
  52     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  53     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  54     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  55     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  56     centrées si le booléen "centeredDF" est vrai.
  57     """
  58     def __init__(self,
  59             Function              = None,
  60             centeredDF            = False,
  61             increment             = 0.01,
  62             dX                    = None,
  63             avoidingRedundancy    = True,
  64             toleranceInRedundancy = 1.e-18,
  65             lenghtOfRedundancy    = -1,
  66             mpEnabled             = False,
  67             mpWorkers             = None,
  68             mfEnabled             = False,
  69             ):
  70         if mpEnabled:
  71             try:
  72                 import multiprocessing
  73                 self.__mpEnabled = True
  74             except ImportError:
  75                 self.__mpEnabled = False
  76         else:
  77             self.__mpEnabled = False
  78         self.__mpWorkers = mpWorkers
  79         if self.__mpWorkers is not None and self.__mpWorkers < 1:
  80             self.__mpWorkers = None
  81         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
  82         #
  83         if mfEnabled:
  84             self.__mfEnabled = True
  85         else:
  86             self.__mfEnabled = False
  87         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
  88         #
  89         if avoidingRedundancy:
  90             self.__avoidRC = True
  91             self.__tolerBP = float(toleranceInRedundancy)
  92             self.__lenghtRJ = int(lenghtOfRedundancy)
  93             self.__listJPCP = [] # Jacobian Previous Calculated Points
  94             self.__listJPCI = [] # Jacobian Previous Calculated Increment
  95             self.__listJPCR = [] # Jacobian Previous Calculated Results
  96             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
  97             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
  98         else:
  99             self.__avoidRC = False
 100         #
 101         if self.__mpEnabled:
 102             if isinstance(Function,types.FunctionType):
 103                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 104                 self.__userFunction__name = Function.__name__
 105                 try:
 106                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 107                 except:
 108                     mod = os.path.abspath(Function.__globals__['__file__'])
 109                 if not os.path.isfile(mod):
 110                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 111                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 112                 self.__userFunction__path = os.path.dirname(mod)
 113                 del mod
 114                 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 115                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 116             elif isinstance(Function,types.MethodType):
 117                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 118                 self.__userFunction__name = Function.__name__
 119                 try:
 120                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 121                 except:
 122                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 123                 if not os.path.isfile(mod):
 124                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 125                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 126                 self.__userFunction__path = os.path.dirname(mod)
 127                 del mod
 128                 self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 129                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 130             else:
 131                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 132         else:
 133             self.__userOperator = Operator( fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 134             self.__userFunction = self.__userOperator.appliedTo
 135         #
 136         self.__centeredDF = bool(centeredDF)
 137         if abs(float(increment)) > 1.e-15:
 138             self.__increment  = float(increment)
 139         else:
 140             self.__increment  = 0.01
 141         if dX is None:
 142             self.__dX     = None
 143         else:
 144             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
 145         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
 146         if self.__avoidRC:
 147             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
 148
 149     # ---------------------------------------------------------
 150     def __doublon__(self, e, l, n, v=None):
 151         __ac, __iac = False, -1
 152         for i in range(len(l)-1,-1,-1):
 153             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
 154                 __ac, __iac = True, i
 155                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
 156                 break
 157         return __ac, __iac
 158
 159     # ---------------------------------------------------------
 160     def DirectOperator(self, X ):
 161         """
 162         Calcul du direct à l'aide de la fonction fournie.
 163         """
 164         logging.debug("FDA Calcul DirectOperator (explicite)")
 165         if self.__mfEnabled:
 166             _HX = self.__userFunction( X, argsAsSerie = True )
 167         else:
 168             _X = numpy.asmatrix(numpy.ravel( X )).T
 169             _HX = numpy.ravel(self.__userFunction( _X ))
 170         #
 171         return _HX
 172
 173     # ---------------------------------------------------------
 174     def TangentMatrix(self, X ):
 175         """
 176         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 177         c'est-à-dire le gradient de H en X. On utilise des différences finies
 178         directionnelles autour du point X. X est un numpy.matrix.
 179
 180         Différences finies centrées (approximation d'ordre 2):
 181         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 182            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 183            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 184            H( X_moins_dXi )
 185         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 186            le pas 2*dXi
 187         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 188
 189         Différences finies non centrées (approximation d'ordre 1):
 190         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 191            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 192            HX_plus_dXi = H( X_plus_dXi )
 193         2/ On calcule la valeur centrale HX = H(X)
 194         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 195            le pas dXi
 196         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 197
 198         """
 199         logging.debug("FDA Début du calcul de la Jacobienne")
 200         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 201         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 202         #
 203         if X is None or len(X)==0:
 204             raise ValueError("Nominal point X for approximate derivatives can not be None or void (X=%s)."%(str(X),))
 205         #
 206         _X = numpy.asmatrix(numpy.ravel( X )).T
 207         #
 208         if self.__dX is None:
 209             _dX  = self.__increment * _X
 210         else:
 211             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 212         #
 213         if (_dX == 0.).any():
 214             moyenne = _dX.mean()
 215             if moyenne == 0.:
 216                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 217             else:
 218                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 219         #
 220         __alreadyCalculated  = False
 221         if self.__avoidRC:
 222             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 223             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 224             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 225                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 226                 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
 227         #
 228         if __alreadyCalculated:
 229             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 230             _Jacobienne = self.__listJPCR[__i]
 231         else:
 232             logging.debug("FDA   Calcul Jacobienne (explicite)")
 233             if self.__centeredDF:
 234                 #
 235                 if self.__mpEnabled and not self.__mfEnabled:
 236                     funcrepr = {
 237                         "__userFunction__path" : self.__userFunction__path,
 238                         "__userFunction__modl" : self.__userFunction__modl,
 239                         "__userFunction__name" : self.__userFunction__name,
 240                     }
 241                     _jobs = []
 242                     for i in range( len(_dX) ):
 243                         _dXi            = _dX[i]
 244                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 245                         _X_plus_dXi[i]  = _X[i] + _dXi
 246                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 247                         _X_moins_dXi[i] = _X[i] - _dXi
 248                         #
 249                         _jobs.append( (_X_plus_dXi,  funcrepr) )
 250                         _jobs.append( (_X_moins_dXi, funcrepr) )
 251                     #
 252                     import multiprocessing
 253                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 254                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 255                     self.__pool.close()
 256                     self.__pool.join()
 257                     #
 258                     _Jacobienne  = []
 259                     for i in range( len(_dX) ):
 260                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 261                     #
 262                 elif self.__mfEnabled:
 263                     _xserie = []
 264                     for i in range( len(_dX) ):
 265                         _dXi            = _dX[i]
 266                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 267                         _X_plus_dXi[i]  = _X[i] + _dXi
 268                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 269                         _X_moins_dXi[i] = _X[i] - _dXi
 270                         #
 271                         _xserie.append( _X_plus_dXi )
 272                         _xserie.append( _X_moins_dXi )
 273                     #
 274                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 275                      #
 276                     _Jacobienne  = []
 277                     for i in range( len(_dX) ):
 278                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 279                     #
 280                 else:
 281                     _Jacobienne  = []
 282                     for i in range( _dX.size ):
 283                         _dXi            = _dX[i]
 284                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 285                         _X_plus_dXi[i]  = _X[i] + _dXi
 286                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 287                         _X_moins_dXi[i] = _X[i] - _dXi
 288                         #
 289                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 290                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 291                         #
 292                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 293                 #
 294             else:
 295                 #
 296                 if self.__mpEnabled and not self.__mfEnabled:
 297                     funcrepr = {
 298                         "__userFunction__path" : self.__userFunction__path,
 299                         "__userFunction__modl" : self.__userFunction__modl,
 300                         "__userFunction__name" : self.__userFunction__name,
 301                     }
 302                     _jobs = []
 303                     _jobs.append( (_X.A1, funcrepr) )
 304                     for i in range( len(_dX) ):
 305                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 306                         _X_plus_dXi[i] = _X[i] + _dX[i]
 307                         #
 308                         _jobs.append( (_X_plus_dXi, funcrepr) )
 309                     #
 310                     import multiprocessing
 311                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 312                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 313                     self.__pool.close()
 314                     self.__pool.join()
 315                     #
 316                     _HX = _HX_plus_dX.pop(0)
 317                     #
 318                     _Jacobienne = []
 319                     for i in range( len(_dX) ):
 320                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 321                     #
 322                 elif self.__mfEnabled:
 323                     _xserie = []
 324                     _xserie.append( _X.A1 )
 325                     for i in range( len(_dX) ):
 326                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 327                         _X_plus_dXi[i] = _X[i] + _dX[i]
 328                         #
 329                         _xserie.append( _X_plus_dXi )
 330                     #
 331                     _HX_plus_dX = self.DirectOperator( _xserie )
 332                     #
 333                     _HX = _HX_plus_dX.pop(0)
 334                     #
 335                     _Jacobienne = []
 336                     for i in range( len(_dX) ):
 337                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 338                    #
 339                 else:
 340                     _Jacobienne  = []
 341                     _HX = self.DirectOperator( _X )
 342                     for i in range( _dX.size ):
 343                         _dXi            = _dX[i]
 344                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 345                         _X_plus_dXi[i]  = _X[i] + _dXi
 346                         #
 347                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 348                         #
 349                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 350                 #
 351             #
 352             _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
 353             if self.__avoidRC:
 354                 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
 355                 while len(self.__listJPCP) > self.__lenghtRJ:
 356                     self.__listJPCP.pop(0)
 357                     self.__listJPCI.pop(0)
 358                     self.__listJPCR.pop(0)
 359                     self.__listJPPN.pop(0)
 360                     self.__listJPIN.pop(0)
 361                 self.__listJPCP.append( copy.copy(_X) )
 362                 self.__listJPCI.append( copy.copy(_dX) )
 363                 self.__listJPCR.append( copy.copy(_Jacobienne) )
 364                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 365                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 366         #
 367         logging.debug("FDA Fin du calcul de la Jacobienne")
 368         #
 369         return _Jacobienne
 370
 371     # ---------------------------------------------------------
 372     def TangentOperator(self, paire ):
 373         """
 374         Calcul du tangent à l'aide de la Jacobienne.
 375         """
 376         if self.__mfEnabled:
 377             assert len(paire) == 1, "Incorrect lenght of arguments"
 378             _paire = paire[0]
 379             assert len(_paire) == 2, "Incorrect number of arguments"
 380         else:
 381             assert len(paire) == 2, "Incorrect number of arguments"
 382             _paire = paire
 383         X, dX = _paire
 384         _Jacobienne = self.TangentMatrix( X )
 385         if dX is None or len(dX) == 0:
 386             #
 387             # Calcul de la forme matricielle si le second argument est None
 388             # -------------------------------------------------------------
 389             if self.__mfEnabled: return [_Jacobienne,]
 390             else:                return _Jacobienne
 391         else:
 392             #
 393             # Calcul de la valeur linéarisée de H en X appliqué à dX
 394             # ------------------------------------------------------
 395             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 396             _HtX = numpy.dot(_Jacobienne, _dX)
 397             if self.__mfEnabled: return [_HtX.A1,]
 398             else:                return _HtX.A1
 399
 400     # ---------------------------------------------------------
 401     def AdjointOperator(self, paire ):
 402         """
 403         Calcul de l'adjoint à l'aide de la Jacobienne.
 404         """
 405         if self.__mfEnabled:
 406             assert len(paire) == 1, "Incorrect lenght of arguments"
 407             _paire = paire[0]
 408             assert len(_paire) == 2, "Incorrect number of arguments"
 409         else:
 410             assert len(paire) == 2, "Incorrect number of arguments"
 411             _paire = paire
 412         X, Y = _paire
 413         _JacobienneT = self.TangentMatrix( X ).T
 414         if Y is None or len(Y) == 0:
 415             #
 416             # Calcul de la forme matricielle si le second argument est None
 417             # -------------------------------------------------------------
 418             if self.__mfEnabled: return [_JacobienneT,]
 419             else:                return _JacobienneT
 420         else:
 421             #
 422             # Calcul de la valeur de l'adjoint en X appliqué à Y
 423             # --------------------------------------------------
 424             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 425             _HaY = numpy.dot(_JacobienneT, _Y)
 426             if self.__mfEnabled: return [_HaY.A1,]
 427             else:                return _HaY.A1
 428
 429 # ==============================================================================
 430 def mmqr(
 431         func     = None,
 432         x0       = None,
 433         fprime   = None,
 434         bounds   = None,
 435         quantile = 0.5,
 436         maxfun   = 15000,
 437         toler    = 1.e-06,
 438         y        = None,
 439         ):
 440     """
 441     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
 442     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
 443     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
 444     """
 445     #
 446     # Recuperation des donnees et informations initiales
 447     # --------------------------------------------------
 448     variables = numpy.ravel( x0 )
 449     mesures   = numpy.ravel( y )
 450     increment = sys.float_info[0]
 451     p         = variables.size
 452     n         = mesures.size
 453     quantile  = float(quantile)
 454     #
 455     # Calcul des parametres du MM
 456     # ---------------------------
 457     tn      = float(toler) / n
 458     e0      = -tn / math.log(tn)
 459     epsilon = (e0-tn)/(1+math.log(e0))
 460     #
 461     # Calculs d'initialisation
 462     # ------------------------
 463     residus  = mesures - numpy.ravel( func( variables ) )
 464     poids    = 1./(epsilon+numpy.abs(residus))
 465     veps     = 1. - 2. * quantile - residus * poids
 466     lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
 467     iteration = 0
 468     #
 469     # Recherche iterative
 470     # -------------------
 471     while (increment > toler) and (iteration < maxfun) :
 472         iteration += 1
 473         #
 474         Derivees  = numpy.array(fprime(variables))
 475         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
 476         DeriveesT = Derivees.transpose()
 477         M         =   numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
 478         SM        =   numpy.transpose(numpy.dot( DeriveesT , veps ))
 479         step      = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
 480         #
 481         variables = variables + step
 482         if bounds is not None:
 483             while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
 484                 step      = step/2.
 485                 variables = variables - step
 486         residus   = mesures - numpy.ravel( func(variables) )
 487         surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
 488         #
 489         while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
 490             step      = step/2.
 491             variables = variables - step
 492             residus   = mesures - numpy.ravel( func(variables) )
 493             surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
 494         #
 495         increment     = lastsurrogate-surrogate
 496         poids         = 1./(epsilon+numpy.abs(residus))
 497         veps          = 1. - 2. * quantile - residus * poids
 498         lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
 499     #
 500     # Mesure d'écart
 501     # --------------
 502     Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
 503     #
 504     return variables, Ecart, [n,p,iteration,increment,0]
 505
 506 # ==============================================================================
 507 if __name__ == "__main__":
 508     print('\n AUTODIAGNOSTIC\n')