0027300: Boolean operation produces invalid shape in terms of "bopargcheck" command
[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         
132     
133     MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
134     if ( MinEigenValue < CTol) {
135        Convex = Standard_False;
136        if (NoConvexTreatement) {
137            Standard_Real Delta = CTol+0.1*Abs(MinEigenValue) -MinEigenValue ;
138            for (ii=1; ii<=TheGradient.Length(); ii++) {
139                TheHessian(ii, ii) += Delta;
140              }
141          }
142        else {
143          Done = Standard_False;
144          TheStatus = math_FunctionError;
145          return;
146        }
147      }
148
149     // Schemas de Newton
150
151     math_Gauss LU(TheHessian, CTol/100);
152     if ( !LU.IsDone()) {
153       Done = Standard_False;
154       TheStatus = math_DirectionSearchError;
155       return;
156     }
157     LU.Solve(TheGradient, TheStep);
158
159     if (myIsBoundsDefined)
160     {
161       // Project point on bounds or nullify TheStep coords if point lies on boudary.
162
163       *suivant = *precedent - TheStep;
164       Standard_Real aMult = RealLast();
165       for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
166       {
167         if (suivant->Value(anIdx) < myLeft(anIdx))
168         {
169           Standard_Real aValue = Abs(precedent->Value(anIdx) - myLeft(anIdx)) / Abs(TheStep(anIdx));
170           aMult = Min (aValue, aMult);
171         }
172
173         if (suivant->Value(anIdx) > myRight(anIdx))
174         {
175           Standard_Real aValue = Abs(precedent->Value(anIdx) - myRight(anIdx)) / Abs(TheStep(anIdx));
176           aMult = Min (aValue, aMult);
177         }
178       }
179
180       if (aMult != RealLast())
181       {
182         if (aMult > Precision::PConfusion())
183         {
184           // Project point into param space.
185           TheStep *= aMult;
186         }
187         else
188         {
189           // Old point on border and new point out of border:
190           // Nullify corresponding TheStep indexes.
191           for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
192           {
193             if (Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() ||
194                 Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion())
195             {
196               TheStep(anIdx) = 0.0;
197             }
198           }
199         }
200       }
201     }
202
203     Standard_Boolean hasProblem = Standard_False;
204     do
205     {
206       *suivant = *precedent - TheStep;
207
208       //  Gestion de la convergence
209       hasProblem = !(F.Value(*suivant, TheMinimum));
210
211       if (hasProblem)
212       {
213         TheStep /= 2.0;
214       }
215     } while (hasProblem);
216
217     if (IsConverged()) { NbConv++; }
218     else               { NbConv=0; }
219
220     //  Controle et corrections.
221
222     VItere = TheMinimum;
223     TheMinimum = PreviousMinimum;
224     Nreduction =0;
225     while (VItere > VPrecedent && Nreduction < 10) {
226         TheStep *= 0.4;   
227         *suivant = *precedent - TheStep;
228         F.Value(*suivant, VItere);
229         Nreduction++;
230     }
231
232     if (VItere <= VPrecedent) {
233        auxiliaire =  precedent;
234        precedent = suivant;
235        suivant = auxiliaire;
236        PreviousMinimum = VPrecedent;
237        TheMinimum = VItere;
238        Ok = (nbiter < Itermax);
239        if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations;
240     } 
241     else {       
242        Ok = Standard_False;
243        TheStatus = math_DirectionSearchError;
244     }
245   }
246  TheLocation = *precedent;    
247 }
248
249 //=======================================================================
250 //function : Dump
251 //purpose  : 
252 //=======================================================================
253 void math_NewtonMinimum::Dump(Standard_OStream& o) const 
254 {
255   o<< "math_Newton Optimisation: ";
256   o << " Done   ="  << Done << endl; 
257   o << " Status = " << (Standard_Integer)TheStatus << endl;
258   o << " Location Vector = " << Location() << endl;
259   o << " Minimum value = "<< Minimum()<< endl;
260   o << " Previous value = "<< PreviousMinimum << endl;
261   o << " Number of iterations = " <<NbIterations() << endl;
262   o << " Convexity = " << Convex << endl;
263   o << " Eigen Value = " << MinEigenValue << endl;
264 }
265