0031939: Coding - correction of spelling errors in comments [part 3]
[occt.git] / 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