0025720: Incorrect code of math classes can lead to unpredicted behavior of algorithms
[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 //#endif
22
23 #include <math_NewtonMinimum.ixx>
24
25 #include <math_Gauss.hxx>
26 #include <math_Jacobi.hxx>
27
28 //=======================================================================
29 //function : math_NewtonMinimum
30 //purpose  : Constructor
31 //=======================================================================
32 math_NewtonMinimum::math_NewtonMinimum(
33   const math_MultipleVarFunctionWithHessian& theFunction,
34   const Standard_Real                        theTolerance,
35   const Standard_Integer                     theNbIterations,
36   const Standard_Real                        theConvexity,
37   const Standard_Boolean                     theWithSingularity
38   )
39 : TheStatus  (math_NotBracketed),
40   TheLocation(1, theFunction.NbVariables()),
41   TheGradient(1, theFunction.NbVariables()),
42   TheStep    (1, theFunction.NbVariables(), 10.0 * theTolerance),
43   TheHessian (1, theFunction.NbVariables(), 1, theFunction.NbVariables()),
44   PreviousMinimum   (0.0),
45   TheMinimum        (0.0),
46   MinEigenValue     (0.0),
47   XTol              (theTolerance),
48   CTol              (theConvexity),
49   nbiter            (0),
50   NoConvexTreatement(theWithSingularity),
51   Convex            (Standard_True),
52   Done              (Standard_False),
53   Itermax           (theNbIterations)
54 {
55 }
56
57 //=======================================================================
58 //function : ~math_NewtonMinimum
59 //purpose  : Destructor
60 //=======================================================================
61 math_NewtonMinimum::~math_NewtonMinimum()
62 {
63 }
64
65 //=======================================================================
66 //function : Perform
67 //purpose  : 
68 //=======================================================================
69 void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
70                                  const math_Vector& StartingPoint)
71 {
72   math_Vector Point1 (1, F.NbVariables());
73   Point1 =  StartingPoint;
74   math_Vector Point2(1, F.NbVariables());
75   math_Vector* precedent = &Point1;
76   math_Vector* suivant = &Point2;
77   math_Vector* auxiliaire = precedent;
78
79   Standard_Boolean Ok = Standard_True;
80   Standard_Integer NbConv = 0, ii, Nreduction;
81   Standard_Real    VPrecedent, VItere; 
82
83   Done = Standard_True;
84   TheStatus = math_OK;
85   nbiter = 0;
86
87   while ( Ok && (NbConv < 2) ) {
88     nbiter++;
89
90     // Positionnement
91
92     Ok = F.Values(*precedent, VPrecedent, TheGradient, TheHessian);
93     if (!Ok) {
94        Done = Standard_False;
95        TheStatus = math_FunctionError;
96        return;
97     }
98     if (nbiter==1) {
99       PreviousMinimum =  VPrecedent;
100       TheMinimum =  VPrecedent;
101     }
102     
103     // Traitement de la non convexite
104
105     math_Jacobi CalculVP(TheHessian);
106     if ( !CalculVP.IsDone() ) {
107        Done = Standard_False;
108        TheStatus = math_FunctionError;
109        return;
110     }
111
112         
113     
114     MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
115     if ( MinEigenValue < CTol) {
116        Convex = Standard_False;
117        if (NoConvexTreatement) {
118            Standard_Real Delta = CTol+0.1*Abs(MinEigenValue) -MinEigenValue ;
119            for (ii=1; ii<=TheGradient.Length(); ii++) {
120                TheHessian(ii, ii) += Delta;
121              }
122          }
123        else {
124            TheStatus = math_FunctionError;
125            return;
126        }
127      }
128
129     // Schemas de Newton
130
131     math_Gauss LU(TheHessian, CTol/100);
132     if ( !LU.IsDone()) {
133       Done = Standard_False;
134       TheStatus = math_DirectionSearchError;
135       return;
136     }
137    
138     LU.Solve(TheGradient, TheStep);
139     Standard_Boolean hasProblem = Standard_False;
140     do
141     {
142       *suivant = *precedent - TheStep;
143
144       //  Gestion de la convergence
145       hasProblem = !(F.Value(*suivant, TheMinimum));
146
147       if (hasProblem)
148       {
149         TheStep /= 2.0;
150       }
151     } while (hasProblem);
152
153     if (IsConverged()) { NbConv++; }
154     else               { NbConv=0; }
155
156     //  Controle et corrections.
157
158     VItere = TheMinimum;
159     TheMinimum = PreviousMinimum;
160     Nreduction =0;
161     while (VItere > VPrecedent && Nreduction < 10) {
162         TheStep *= 0.4;   
163         *suivant = *precedent - TheStep;
164         F.Value(*suivant, VItere);
165         Nreduction++;
166     }
167
168     if (VItere <= VPrecedent) {
169        auxiliaire =  precedent;
170        precedent = suivant;
171        suivant = auxiliaire;
172        PreviousMinimum = VPrecedent;
173        TheMinimum = VItere;
174        Ok = (nbiter < Itermax);
175        if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations;
176     } 
177     else {       
178        Ok = Standard_False;
179        TheStatus = math_DirectionSearchError;
180     }
181   }
182  TheLocation = *precedent;    
183 }
184
185 //=======================================================================
186 //function : Dump
187 //purpose  : 
188 //=======================================================================
189 void math_NewtonMinimum::Dump(Standard_OStream& o) const 
190 {
191   o<< "math_Newton Optimisation: ";
192   o << " Done   ="  << Done << endl; 
193   o << " Status = " << (Standard_Integer)TheStatus << endl;
194   o << " Location Vector = " << Location() << endl;
195   o << " Minimum value = "<< Minimum()<< endl;
196   o << " Previous value = "<< PreviousMinimum << endl;
197   o << " Number of iterations = " <<NbIterations() << endl;
198   o << " Convexity = " << Convex << endl;
199   o << " Eigen Value = " << MinEigenValue << endl;
200 }
201