src/daComposant/daAlgorithms/Atoms/mmqr.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2022 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     MMQR
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import sys, math, numpy
  29 from daCore.PlatformInfo import PlatformInfo
  30 mpr = PlatformInfo().MachinePrecision()
  31 mfp = PlatformInfo().MaximumPrecision()
  32
  33 # ==============================================================================
  34 def mmqr(
  35         func     = None,
  36         x0       = None,
  37         fprime   = None,
  38         bounds   = None,
  39         quantile = 0.5,
  40         maxfun   = 15000,
  41         toler    = 1.e-06,
  42         y        = None,
  43         ):
  44     """
  45     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
  46     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
  47     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
  48     """
  49     #
  50     # Récupération des données et informations initiales
  51     # --------------------------------------------------
  52     variables = numpy.ravel( x0 )
  53     mesures   = numpy.ravel( y )
  54     increment = sys.float_info[0]
  55     p         = variables.size
  56     n         = mesures.size
  57     quantile  = float(quantile)
  58     #
  59     # Calcul des paramètres du MM
  60     # ---------------------------
  61     tn      = float(toler) / n
  62     e0      = -tn / math.log(tn)
  63     epsilon = (e0-tn)/(1+math.log(e0))
  64     #
  65     # Calculs d'initialisation
  66     # ------------------------
  67     residus  = mesures - numpy.ravel( func( variables ) )
  68     poids    = 1./(epsilon+numpy.abs(residus))
  69     veps     = 1. - 2. * quantile - residus * poids
  70     lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
  71     iteration = 0
  72     #
  73     # Recherche itérative
  74     # -------------------
  75     while (increment > toler) and (iteration < maxfun) :
  76         iteration += 1
  77         #
  78         Derivees  = numpy.array(fprime(variables))
  79         Derivees  = Derivees.reshape(n,p) # ADAO & check shape
  80         DeriveesT = Derivees.transpose()
  81         M         =   numpy.dot( DeriveesT , (numpy.array(p*[poids,]).T * Derivees) )
  82         SM        =   numpy.transpose(numpy.dot( DeriveesT , veps ))
  83         step      = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
  84         #
  85         variables = variables + step
  86         if bounds is not None:
  87             # Attention : boucle infinie à éviter si un intervalle est trop petit
  88             while( (variables < numpy.ravel(numpy.asarray(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asarray(bounds)[:,1])).any() ):
  89                 step      = step/2.
  90                 variables = variables - step
  91         residus   = mesures - numpy.ravel( func(variables) )
  92         surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
  93         #
  94         while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
  95             step      = step/2.
  96             variables = variables - step
  97             residus   = mesures - numpy.ravel( func(variables) )
  98             surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
  99         #
 100         increment     = lastsurrogate-surrogate
 101         poids         = 1./(epsilon+numpy.abs(residus))
 102         veps          = 1. - 2. * quantile - residus * poids
 103         lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
 104     #
 105     # Mesure d'écart
 106     # --------------
 107     Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
 108     #
 109     return variables, Ecart, [n,p,iteration,increment,0]
 110
 111 # ==============================================================================
 112 if __name__ == "__main__":
 113     print('\n AUTODIAGNOSTIC\n')