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-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 //#ifndef DEB
  23 #define No_Standard_RangeError
  24 #define No_Standard_OutOfRange
  25 #define No_Standard_DimensionError
  26 //#endif
  27
  28 #include <math_NewtonMinimum.ixx>
  29
  30 #include <math_Gauss.hxx>
  31 #include <math_Jacobi.hxx>
  32
  33 //============================================================================
  34 math_NewtonMinimum::math_NewtonMinimum(math_MultipleVarFunctionWithHessian& F,
  35                                        const math_Vector& StartingPoint,
  36                                        const Standard_Real Tolerance,
  37                                        const Standard_Integer NbIterations,
  38                                        const Standard_Real Convexity,
  39                                        const Standard_Boolean WithSingularity)
  40 //============================================================================
  41                   : TheLocation(1, F.NbVariables()),
  42                     TheGradient(1, F.NbVariables()),
  43                     TheStep(1, F.NbVariables(), 10*Tolerance),
  44                     TheHessian(1, F.NbVariables(), 1, F.NbVariables() )
  45 {
  46        XTol = Tolerance;
  47        CTol = Convexity;
  48        Itermax = NbIterations;
  49        NoConvexTreatement = WithSingularity;
  50        Convex = Standard_True;
  51 //       Done = Standard_True;
  52 //       TheStatus = math_OK;
  53        Perform ( F, StartingPoint);
  54 }
  55
  56 //============================================================================
  57 math_NewtonMinimum::math_NewtonMinimum(math_MultipleVarFunctionWithHessian& F,
  58                                        const Standard_Real Tolerance,
  59                                        const Standard_Integer NbIterations,
  60                                        const Standard_Real Convexity,
  61                                        const Standard_Boolean WithSingularity)
  62 //============================================================================
  63                   : TheLocation(1, F.NbVariables()),
  64                     TheGradient(1, F.NbVariables()),
  65                     TheStep(1, F.NbVariables(), 10*Tolerance),
  66                     TheHessian(1, F.NbVariables(), 1, F.NbVariables() )
  67 {
  68        XTol = Tolerance;
  69        CTol = Convexity;
  70        Itermax = NbIterations;
  71        NoConvexTreatement = WithSingularity;
  72        Convex = Standard_True;
  73        Done = Standard_False;
  74        TheStatus = math_NotBracketed;
  75 }
  76 //============================================================================
  77 void math_NewtonMinimum::Delete()
  78 {}
  79 //============================================================================
  80 void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
  81                                  const math_Vector& StartingPoint)
  82 //============================================================================
  83 {
  84   math_Vector Point1 (1, F.NbVariables());
  85   Point1 =  StartingPoint;
  86   math_Vector Point2(1, F.NbVariables());
  87   math_Vector* precedent = &Point1;
  88   math_Vector* suivant = &Point2;
  89   math_Vector* auxiliaire = precedent;
  90
  91   Standard_Boolean Ok = Standard_True;
  92   Standard_Integer NbConv = 0, ii, Nreduction;
  93   Standard_Real    VPrecedent, VItere;
  94
  95   Done = Standard_True;
  96   TheStatus = math_OK;
  97   nbiter = 0;
  98
  99   while ( Ok && (NbConv < 2) ) {
 100     nbiter++;
 101
 102     // Positionnement
 103
 104     Ok = F.Values(*precedent, VPrecedent, TheGradient, TheHessian);
 105     if (!Ok) {
 106        Done = Standard_False;
 107        TheStatus = math_FunctionError;
 108        return;
 109     }
 110     if (nbiter==1) {
 111       PreviousMinimum =  VPrecedent;
 112       TheMinimum =  VPrecedent;
 113     }
 114
 115     // Traitement de la non convexite
 116
 117     math_Jacobi CalculVP(TheHessian);
 118     if ( !CalculVP.IsDone() ) {
 119        Done = Standard_False;
 120        TheStatus = math_FunctionError;
 121        return;
 122     }
 123
 124
 125
 126     MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
 127     if ( MinEigenValue < CTol) {
 128        Convex = Standard_False;
 129        if (NoConvexTreatement) {
 130            Standard_Real Delta = CTol+0.1*Abs(MinEigenValue) -MinEigenValue ;
 131            for (ii=1; ii<=TheGradient.Length(); ii++) {
 132                TheHessian(ii, ii) += Delta;
 133              }
 134          }
 135        else {
 136            TheStatus = math_FunctionError;
 137            return;
 138        }
 139      }
 140
 141     // Schemas de Newton
 142
 143     math_Gauss LU(TheHessian, CTol/100);
 144     if ( !LU.IsDone()) {
 145       Done = Standard_False;
 146       TheStatus = math_DirectionSearchError;
 147       return;
 148     }
 149
 150     LU.Solve(TheGradient, TheStep);
 151     *suivant = *precedent - TheStep;
 152
 153
 154     //  Gestion de la convergence
 155
 156     F.Value(*suivant, TheMinimum);
 157
 158     if (IsConverged()) { NbConv++; }
 159     else               { NbConv=0; }
 160
 161     //  Controle et corrections.
 162
 163     VItere = TheMinimum;
 164     TheMinimum = PreviousMinimum;
 165     Nreduction =0;
 166     while (VItere > VPrecedent && Nreduction < 10) {
 167         TheStep *= 0.4;
 168         *suivant = *precedent - TheStep;
 169         F.Value(*suivant, VItere);
 170         Nreduction++;
 171     }
 172
 173     if (VItere <= VPrecedent) {
 174        auxiliaire =  precedent;
 175        precedent = suivant;
 176        suivant = auxiliaire;
 177        PreviousMinimum = VPrecedent;
 178        TheMinimum = VItere;
 179        Ok = (nbiter < Itermax);
 180        if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations;
 181     }
 182     else {
 183        Ok = Standard_False;
 184        TheStatus = math_DirectionSearchError;
 185     }
 186   }
 187  TheLocation = *precedent;
 188 }
 189
 190 //============================================================================
 191 Standard_Boolean math_NewtonMinimum::IsConverged() const
 192 //============================================================================
 193 {
 194   return ( (TheStep.Norm() <= XTol ) ||
 195          ( Abs(TheMinimum-PreviousMinimum) <= XTol*Abs(PreviousMinimum) ));
 196 }
 197
 198 //============================================================================
 199 void math_NewtonMinimum::Dump(Standard_OStream& o) const
 200 //============================================================================
 201 {
 202        o<< "math_Newton Optimisation: ";
 203          o << " Done   ="  << Done << endl;
 204          o << " Status = " << (Standard_Integer)TheStatus << endl;
 205          o <<" Location Vector = " << Location() << endl;
 206          o <<" Minimum value = "<< Minimum()<< endl;
 207          o <<" Previous value = "<< PreviousMinimum << endl;
 208          o <<" Number of iterations = " <<NbIterations() << endl;
 209          o <<" Convexity = " << Convex << endl;
 210          o <<" Eigen Value = " << MinEigenValue << endl;
 211 }
 212