src/math/math_NewtonMinimum.cxx

   1 // Created on: 1996-05-03
   2 // Created by: Philippe MANGIN
   3 // Copyright (c) 1996-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //#ifndef OCCT_DEBUG
  18 #define No_Standard_RangeError
  19 #define No_Standard_OutOfRange
  20 #define No_Standard_DimensionError
  21
  22 //#endif
  23
  24 #include <math_Gauss.hxx>
  25 #include <math_Jacobi.hxx>
  26 #include <math_MultipleVarFunctionWithHessian.hxx>
  27 #include <math_NewtonMinimum.hxx>
  28 #include <Precision.hxx>
  29 #include <Standard_DimensionError.hxx>
  30 #include <StdFail_NotDone.hxx>
  31
  32 //=======================================================================
  33 //function : math_NewtonMinimum
  34 //purpose  : Constructor
  35 //=======================================================================
  36 math_NewtonMinimum::math_NewtonMinimum(
  37   const math_MultipleVarFunctionWithHessian& theFunction,
  38   const Standard_Real                        theTolerance,
  39   const Standard_Integer                     theNbIterations,
  40   const Standard_Real                        theConvexity,
  41   const Standard_Boolean                     theWithSingularity
  42   )
  43 : TheStatus  (math_NotBracketed),
  44   TheLocation(1, theFunction.NbVariables()),
  45   TheGradient(1, theFunction.NbVariables()),
  46   TheStep    (1, theFunction.NbVariables(), 10.0 * theTolerance),
  47   TheHessian (1, theFunction.NbVariables(), 1, theFunction.NbVariables()),
  48   PreviousMinimum   (0.0),
  49   TheMinimum        (0.0),
  50   MinEigenValue     (0.0),
  51   XTol              (theTolerance),
  52   CTol              (theConvexity),
  53   nbiter            (0),
  54   NoConvexTreatement(theWithSingularity),
  55   Convex            (Standard_True),
  56   myIsBoundsDefined(Standard_False),
  57   myLeft(1, theFunction.NbVariables(), 0.0),
  58   myRight(1, theFunction.NbVariables(), 0.0),
  59   Done              (Standard_False),
  60   Itermax           (theNbIterations)
  61 {
  62 }
  63
  64 //=======================================================================
  65 //function : ~math_NewtonMinimum
  66 //purpose  : Destructor
  67 //=======================================================================
  68 math_NewtonMinimum::~math_NewtonMinimum()
  69 {
  70 }
  71
  72 //=======================================================================
  73 //function : SetBoundary
  74 //purpose  : Set boundaries for conditional optimization
  75 //=======================================================================
  76 void math_NewtonMinimum::SetBoundary(const math_Vector& theLeftBorder,
  77                                      const math_Vector& theRightBorder)
  78 {
  79   myLeft = theLeftBorder;
  80   myRight = theRightBorder;
  81   myIsBoundsDefined = Standard_True;
  82 }
  83
  84 //=======================================================================
  85 //function : Perform
  86 //purpose  :
  87 //=======================================================================
  88 void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
  89                                  const math_Vector& StartingPoint)
  90 {
  91   math_Vector Point1 (1, F.NbVariables());
  92   Point1 =  StartingPoint;
  93   math_Vector Point2(1, F.NbVariables());
  94   math_Vector* precedent = &Point1;
  95   math_Vector* suivant = &Point2;
  96   math_Vector* auxiliaire = precedent;
  97
  98   Standard_Boolean Ok = Standard_True;
  99   Standard_Integer NbConv = 0, ii, Nreduction;
 100   Standard_Real    VPrecedent, VItere;
 101
 102   Done = Standard_True;
 103   TheStatus = math_OK;
 104   nbiter = 0;
 105
 106   while ( Ok && (NbConv < 2) ) {
 107     nbiter++;
 108
 109     // Positionnement
 110
 111     Ok = F.Values(*precedent, VPrecedent, TheGradient, TheHessian);
 112     if (!Ok) {
 113        Done = Standard_False;
 114        TheStatus = math_FunctionError;
 115        return;
 116     }
 117     if (nbiter==1) {
 118       PreviousMinimum =  VPrecedent;
 119       TheMinimum =  VPrecedent;
 120     }
 121
 122     // Traitement de la non convexite
 123
 124     math_Jacobi CalculVP(TheHessian);
 125     if ( !CalculVP.IsDone() ) {
 126        Done = Standard_False;
 127        TheStatus = math_FunctionError;
 128        return;
 129     }
 130
 131     MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
 132     if ( MinEigenValue < CTol)
 133     {
 134       Convex = Standard_False;
 135       if (NoConvexTreatement && // Treatment is allowed.
 136           Abs (MinEigenValue) > CTol) // Treatment will have effect.
 137       {
 138         Standard_Real Delta = CTol + 0.1 * Abs(MinEigenValue) - MinEigenValue;
 139         for (ii=1; ii<=TheGradient.Length(); ii++)
 140           TheHessian(ii, ii) += Delta;
 141       }
 142       else
 143       {
 144         Done = Standard_False;
 145         TheStatus = math_FunctionError;
 146         return;
 147       }
 148     }
 149
 150     // Schemas de Newton
 151
 152     math_Gauss LU(TheHessian, CTol/100);
 153     if ( !LU.IsDone()) {
 154       Done = Standard_False;
 155       TheStatus = math_DirectionSearchError;
 156       return;
 157     }
 158     LU.Solve(TheGradient, TheStep);
 159
 160     if (myIsBoundsDefined)
 161     {
 162       // Project point on bounds or nullify TheStep coords if point lies on boundary.
 163
 164       *suivant = *precedent - TheStep;
 165       Standard_Real aMult = RealLast();
 166       for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
 167       {
 168         const Standard_Real anAbsStep = Abs(TheStep(anIdx));
 169         if (anAbsStep < gp::Resolution())
 170           continue;
 171
 172         if (suivant->Value(anIdx) < myLeft(anIdx))
 173         {
 174           Standard_Real aValue = Abs(precedent->Value(anIdx) - myLeft(anIdx)) / anAbsStep;
 175           aMult = Min (aValue, aMult);
 176         }
 177
 178         if (suivant->Value(anIdx) > myRight(anIdx))
 179         {
 180           Standard_Real aValue = Abs(precedent->Value(anIdx) - myRight(anIdx)) / anAbsStep;
 181           aMult = Min (aValue, aMult);
 182         }
 183       }
 184
 185       if (aMult != RealLast())
 186       {
 187         if (aMult > Precision::PConfusion())
 188         {
 189           // Project point into param space.
 190           TheStep *= aMult;
 191         }
 192         else
 193         {
 194           // Old point on border and new point out of border:
 195           // Nullify corresponding TheStep indexes.
 196           for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
 197           {
 198             if ((Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() && TheStep(anIdx) < 0.0) ||
 199                 (Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion() && TheStep(anIdx) > 0.0) )
 200             {
 201               TheStep(anIdx) = 0.0;
 202             }
 203           }
 204         }
 205       }
 206     }
 207
 208     Standard_Boolean hasProblem = Standard_False;
 209     do
 210     {
 211       *suivant = *precedent - TheStep;
 212
 213       //  Gestion de la convergence
 214       hasProblem = !(F.Value(*suivant, TheMinimum));
 215
 216       if (hasProblem)
 217       {
 218         TheStep /= 2.0;
 219       }
 220     } while (hasProblem);
 221
 222     if (IsConverged()) { NbConv++; }
 223     else               { NbConv=0; }
 224
 225     //  Controle et corrections.
 226
 227     VItere = TheMinimum;
 228     TheMinimum = PreviousMinimum;
 229     Nreduction =0;
 230     while (VItere > VPrecedent && Nreduction < 10) {
 231         TheStep *= 0.4;
 232         *suivant = *precedent - TheStep;
 233         F.Value(*suivant, VItere);
 234         Nreduction++;
 235     }
 236
 237     if (VItere <= VPrecedent) {
 238        auxiliaire =  precedent;
 239        precedent = suivant;
 240        suivant = auxiliaire;
 241        PreviousMinimum = VPrecedent;
 242        TheMinimum = VItere;
 243        Ok = (nbiter < Itermax);
 244        if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations;
 245     }
 246     else {
 247        Ok = Standard_False;
 248        TheStatus = math_DirectionSearchError;
 249     }
 250   }
 251  TheLocation = *precedent;
 252 }
 253
 254 //=======================================================================
 255 //function : Dump
 256 //purpose  :
 257 //=======================================================================
 258 void math_NewtonMinimum::Dump(Standard_OStream& o) const
 259 {
 260   o<< "math_Newton Optimisation: ";
 261   o << " Done   ="  << Done << std::endl;
 262   o << " Status = " << (Standard_Integer)TheStatus << std::endl;
 263   o << " Location Vector = " << Location() << std::endl;
 264   o << " Minimum value = "<< Minimum()<< std::endl;
 265   o << " Previous value = "<< PreviousMinimum << std::endl;
 266   o << " Number of iterations = " <<NbIterations() << std::endl;
 267   o << " Convexity = " << Convex << std::endl;
 268   o << " Eigen Value = " << MinEigenValue << std::endl;
 269 }
 270