0022904: Clean up sccsid variables
[occt.git] / src / math / math_BFGS.cxx
1 //#ifndef DEB
2 #define No_Standard_RangeError
3 #define No_Standard_OutOfRange
4 #define No_Standard_DimensionError
5 //#endif
6
7 #include <math_BFGS.ixx>
8
9 #include <math_BracketMinimum.hxx>
10 #include <math_BrentMinimum.hxx>
11 #include <math_FunctionWithDerivative.hxx>
12 #include <math_MultipleVarFunctionWithGradient.hxx>
13 #include <math_Matrix.hxx>
14
15 #define R 0.61803399
16 #define C (1.0-R)
17 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
18 #define SIGN(a, b) ((b) > 0.0 ? fabs(a) : -fabs(a))
19 #define MOV3(a,b,c, d, e, f) (a)=(d); (b)= (e); (c)=(f);
20
21
22 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
23 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
24 // classique et aussi essayer "recherche unidimensionnelle economique"
25 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
26
27 class DirFunction : public math_FunctionWithDerivative {
28
29      math_Vector *P0;
30      math_Vector *Dir;
31      math_Vector *P;
32      math_Vector *G;
33      math_MultipleVarFunctionWithGradient *F;
34
35 public :
36
37      DirFunction(math_Vector& V1, 
38                  math_Vector& V2,
39                  math_Vector& V3,
40                  math_Vector& V4,
41                  math_MultipleVarFunctionWithGradient& f) ;
42
43      void Initialize(const math_Vector& p0, const math_Vector& dir) const;
44      void TheGradient(math_Vector& Grad);
45      virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
46      virtual Standard_Boolean Values(const Standard_Real x, Standard_Real& fval, Standard_Real& D) ;
47      virtual Standard_Boolean Derivative(const Standard_Real x, Standard_Real& D) ;
48
49      
50 };
51
52      DirFunction::DirFunction(math_Vector& V1, 
53                               math_Vector& V2,
54                               math_Vector& V3,
55                               math_Vector& V4,
56                               math_MultipleVarFunctionWithGradient& f) {
57         
58         P0  = &V1;
59         Dir = &V2;
60         P   = &V3;
61         F   = &f;
62         G =   &V4;
63      }
64
65      void DirFunction::Initialize(const math_Vector& p0, 
66                                   const math_Vector& dir)  const{
67
68         *P0 = p0;
69         *Dir = dir;
70      }
71
72      void DirFunction::TheGradient(math_Vector& Grad) {
73        Grad = *G;
74      }
75
76
77      Standard_Boolean DirFunction::Value(const Standard_Real x, 
78                                          Standard_Real& fval) 
79      {
80         *P = *Dir;
81         P->Multiply(x);
82         P->Add(*P0);        
83         F->Value(*P, fval);
84         return Standard_True;
85      }
86
87      Standard_Boolean DirFunction::Values(const Standard_Real x, 
88                                           Standard_Real& fval, 
89                                           Standard_Real& D) 
90      {
91         *P = *Dir;
92         P->Multiply(x);
93         P->Add(*P0);  
94         F->Values(*P, fval, *G);
95         D = (*G).Multiplied(*Dir);
96         return Standard_True;
97      }
98      Standard_Boolean DirFunction::Derivative(const Standard_Real x, 
99                                               Standard_Real& D) 
100      {
101         *P = *Dir;
102         P->Multiply(x);
103         P->Add(*P0);        
104         Standard_Real fval;
105         F->Values(*P, fval, *G);
106         D = (*G).Multiplied(*Dir);
107         return Standard_True;
108      }
109
110
111 static Standard_Boolean MinimizeDirection(math_Vector&   P,
112                                           Standard_Real  F0,
113                                           math_Vector&   Gr,
114                                           math_Vector&   Dir,
115                                           Standard_Real& Result,
116                                           DirFunction&   F) {
117
118
119
120      Standard_Real ax, xx, bx, Fax, Fxx, Fbx, F1;
121      F.Initialize(P, Dir);
122
123      Standard_Real dy1, Hnr1, lambda, alfa=0;
124 #ifdef DEB
125      Standard_Integer n = 
126 #endif
127        Dir.Length();
128      dy1 = Gr*Dir;
129      if (dy1 != 0) {
130        Hnr1 = Dir.Norm2();
131        alfa = 0.7*(-F0)/dy1;
132        lambda = 0.015/Sqrt(Hnr1);
133      }
134      else lambda = 1.0;
135      if (lambda > alfa) lambda = alfa;
136      F.Value(lambda, F1);
137      math_BracketMinimum Bracket(F, 0.0, lambda, F0, F1);
138      if(Bracket.IsDone()) {
139        Bracket.Values(ax, xx, bx);
140        Bracket.FunctionValues(Fax, Fxx, Fbx);
141
142        Standard_Integer niter = 100;
143        Standard_Real tol = 1.e-03;
144        math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
145        Sol.Perform(F, ax, xx, bx);
146        if(Sol.IsDone()) {
147          Standard_Real Scale = Sol.Location();
148          Result = Sol.Minimum();
149          Dir.Multiply(Scale);
150          P.Add(Dir);
151          return Standard_True;
152        }
153      }
154      return Standard_False;
155    }
156
157
158 void  math_BFGS::Perform(math_MultipleVarFunctionWithGradient& F,
159                          const math_Vector& StartingPoint) {
160   
161        Standard_Boolean Good;
162        Standard_Integer n = TheLocation.Length();
163        Standard_Integer j, i;
164        Standard_Real fae, fad, fac;
165
166        math_Vector xi(1, n), dg(1, n), hdg(1, n);
167        math_Matrix hessin(1, n, 1, n);
168        hessin.Init(0.0);
169
170        math_Vector Temp1(1, n);
171        math_Vector Temp2(1, n);
172        math_Vector Temp3(1, n);
173        math_Vector Temp4(1, n);
174        DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
175
176        TheLocation = StartingPoint;
177        Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
178        if(!Good) { 
179          Done = Standard_False;
180          TheStatus = math_FunctionError;
181          return;
182        }
183        for(i = 1; i <= n; i++) {
184          hessin(i, i) = 1.0;
185          xi(i) = -TheGradient(i);
186        }
187
188
189        for(nbiter = 1; nbiter <= Itermax; nbiter++) {
190          TheMinimum = PreviousMinimum;
191          Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum,
192                                                      TheGradient,
193                                                      xi, TheMinimum, F_Dir);
194          if(!IsGood) {
195            Done = Standard_False;
196            TheStatus = math_DirectionSearchError;
197            return;
198          }
199          if(IsSolutionReached(F)) {
200            Done = Standard_True;
201            TheStatus = math_OK;
202            return;
203          }
204          if (nbiter == Itermax) {
205            Done = Standard_False;
206            TheStatus = math_TooManyIterations;
207            return;
208          }
209          PreviousMinimum = TheMinimum;
210
211          dg = TheGradient;
212
213          Good = F.Values(TheLocation, TheMinimum, TheGradient);
214          if(!Good) { 
215            Done = Standard_False;
216            TheStatus = math_FunctionError;
217            return;
218          }
219
220          for(i = 1; i <= n; i++) {
221            dg(i) = TheGradient(i) - dg(i);
222          }
223          for(i = 1; i <= n; i++) {
224            hdg(i) = 0.0;
225            for (j = 1; j <= n; j++) 
226              hdg(i) += hessin(i, j) * dg(j);
227          }
228          fac = fae = 0.0;
229          for(i = 1; i <= n; i++) {
230            fac += dg(i) * xi(i);
231            fae += dg(i) * hdg(i);
232          }
233          fac = 1.0 / fac;
234          fad = 1.0 / fae;
235
236          for(i = 1; i <= n; i++) 
237            dg(i) = fac * xi(i) - fad * hdg(i);
238          
239          for(i = 1; i <= n; i++) 
240            for(j = 1; j <= n; j++)
241              hessin(i, j) += fac * xi(i) * xi(j)
242                - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
243
244          for(i = 1; i <= n; i++) {
245            xi(i) = 0.0;
246            for (j = 1; j <= n; j++) xi(i) -= hessin(i, j) * TheGradient(j);
247          }  
248        }
249        Done = Standard_False;
250        TheStatus = math_TooManyIterations;
251        return;
252    }
253
254     Standard_Boolean math_BFGS::IsSolutionReached(
255                            math_MultipleVarFunctionWithGradient&) const {
256
257        return 2.0 * fabs(TheMinimum - PreviousMinimum) <= 
258               XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
259     }
260
261     math_BFGS::math_BFGS(math_MultipleVarFunctionWithGradient& F,
262                          const math_Vector& StartingPoint, 
263                          const Standard_Real        Tolerance,
264                          const Standard_Integer     NbIterations,
265                          const Standard_Real        ZEPS) 
266                          : TheLocation(1, StartingPoint.Length()),
267                            TheGradient(1, StartingPoint.Length()) {
268
269        XTol = Tolerance;
270        EPSZ = ZEPS;
271        Itermax = NbIterations;
272        Perform(F, StartingPoint);
273     }
274                              
275     math_BFGS::math_BFGS(math_MultipleVarFunctionWithGradient& F,
276                          const Standard_Real        Tolerance,
277                          const Standard_Integer     NbIterations,
278                          const Standard_Real        ZEPS) 
279                          : TheLocation(1, F.NbVariables()),
280                            TheGradient(1, F.NbVariables()) {
281
282        XTol = Tolerance;
283        EPSZ = ZEPS;
284        Itermax = NbIterations;
285     }
286
287
288     void math_BFGS::Delete()
289     {}
290
291     void math_BFGS::Dump(Standard_OStream& o) const {
292
293        o<< "math_BFGS resolution: ";
294        if(Done) {
295          o << " Status = Done \n";
296          o <<" Location Vector = " << Location() << "\n";
297          o <<" Minimum value = "<< Minimum()<<"\n";
298          o <<" Number of iterations = "<<NbIterations() <<"\n";;
299        }
300        else {
301          o<< " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
302        }
303     }
304
305
306