a3776e4ab1def533157409a9302c9b4d974ceabb
[occt.git] / src / math / math_FRPR.cxx
1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 //#ifndef OCCT_DEBUG
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
19 //#endif
20
21 #include <math_FRPR.ixx>
22
23 #include <math_BracketMinimum.hxx>
24 #include <math_BrentMinimum.hxx>
25 #include <math_Function.hxx>
26 #include <math_MultipleVarFunction.hxx>
27 #include <math_MultipleVarFunctionWithGradient.hxx>
28
29 // l'utilisation de math_BrentMinumim pur trouver un minimum dans une direction
30 // donnee n'est pas du tout optimale. voir peut etre interpolation cubique
31 // classique et aussi essayer "recherche unidimensionnelle economique"
32 // PROGRAMMATION MATHEMATIQUE (theorie et algorithmes) tome1 page 82.
33
34 class DirFunctionTer : public math_Function {
35
36      math_Vector *P0;
37      math_Vector *Dir;
38      math_Vector *P;
39      math_MultipleVarFunction *F;
40
41 public :
42
43      DirFunctionTer(math_Vector& V1, 
44                  math_Vector& V2,
45                  math_Vector& V3,
46                  math_MultipleVarFunction& f);
47
48      void Initialize(const math_Vector& p0, const math_Vector& dir);
49
50      virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval);
51 };
52
53      DirFunctionTer::DirFunctionTer(math_Vector& V1, 
54                               math_Vector& V2,
55                               math_Vector& V3,
56                               math_MultipleVarFunction& f) {
57         
58         P0  = &V1;
59         Dir = &V2;
60         P   = &V3;
61         F   = &f;
62      }
63
64      void DirFunctionTer::Initialize(const math_Vector& p0, 
65                                   const math_Vector& dir) {
66
67         *P0 = p0;
68         *Dir = dir;
69      }
70
71      Standard_Boolean DirFunctionTer::Value(const Standard_Real x, Standard_Real& fval) {
72
73         *P = *Dir;
74         P->Multiply(x);
75         P->Add(*P0);        
76         F->Value(*P, fval);
77         return Standard_True;
78      }
79
80 static Standard_Boolean MinimizeDirection(math_Vector& P,
81                                  math_Vector& Dir,
82                                  Standard_Real& Result,
83                                  DirFunctionTer& F) {
84
85      Standard_Real ax, xx, bx;
86
87      F.Initialize(P, Dir);
88      math_BracketMinimum Bracket(F, 0.0, 1.0);
89      if(Bracket.IsDone()) {
90        Bracket.Values(ax, xx, bx);
91        math_BrentMinimum Sol(F, ax, xx, bx, 1.0e-10, 100);
92        if(Sol.IsDone()) {
93          Standard_Real Scale = Sol.Location();
94          Result = Sol.Minimum();
95          Dir.Multiply(Scale);
96          P.Add(Dir);
97          return Standard_True;
98        }
99      }
100      return Standard_False;
101   }
102
103
104 void  math_FRPR::Perform(math_MultipleVarFunctionWithGradient& F,
105                          const math_Vector& StartingPoint) {
106   
107        Standard_Boolean Good;
108        Standard_Integer n = TheLocation.Length();
109        Standard_Integer j, its;
110        Standard_Real gg, gam, dgg;
111
112        math_Vector g(1, n), h(1, n);
113  
114        math_Vector Temp1(1, n);
115        math_Vector Temp2(1, n);
116        math_Vector Temp3(1, n);
117        DirFunctionTer F_Dir(Temp1, Temp2, Temp3, F);
118
119        TheLocation = StartingPoint;
120        Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
121        if(!Good) { 
122          Done = Standard_False;
123          TheStatus = math_FunctionError;
124          return;
125        }
126
127        g = -TheGradient;
128        h = g;
129        TheGradient = g;
130
131        for(its = 1; its <= Itermax; its++) {
132          Iter = its;
133          
134          Standard_Boolean IsGood = MinimizeDirection(TheLocation,
135                                             TheGradient, TheMinimum, F_Dir);
136          if(!IsGood) {
137            Done = Standard_False;
138            TheStatus = math_DirectionSearchError;
139            return;
140          }
141          if(IsSolutionReached(F)) {
142            Done = Standard_True;
143            State = F.GetStateNumber();
144            TheStatus = math_OK;
145            return;
146          }
147          Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
148          if(!Good) { 
149            Done = Standard_False;
150            TheStatus = math_FunctionError;
151            return;
152          }
153
154          dgg =0.0;
155          gg = 0.0;
156          
157          for(j = 1; j<= n; j++) {
158            gg += g(j)*g(j);
159 //         dgg += TheGradient(j)*TheGradient(j);  //for Fletcher-Reeves
160            dgg += (TheGradient(j)+g(j)) * TheGradient(j);  //for Polak-Ribiere
161          }
162
163          if (gg == 0.0) {
164            //Unlikely. If gradient is exactly 0 then we are already done.
165            Done = Standard_False;
166            TheStatus = math_FunctionError;
167            return;
168          }
169          
170          gam = dgg/gg;
171          g = -TheGradient;
172          TheGradient = g + gam*h;
173          h = TheGradient;
174        }
175        Done = Standard_False;
176        TheStatus = math_TooManyIterations;
177        return;
178
179      }
180
181
182
183     Standard_Boolean math_FRPR::IsSolutionReached(
184 //                           math_MultipleVarFunctionWithGradient& F) {
185                            math_MultipleVarFunctionWithGradient& ) {
186
187        return (2.0 * fabs(TheMinimum - PreviousMinimum)) <= 
188               XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
189     }
190
191     math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F,
192                          const math_Vector& StartingPoint, 
193                          const Standard_Real        Tolerance,
194                          const Standard_Integer     NbIterations,
195                          const Standard_Real        ZEPS) 
196                          : TheLocation(1, StartingPoint.Length()),
197                            TheGradient(1, StartingPoint.Length()) {
198
199        XTol = Tolerance;
200        EPSZ = ZEPS;
201        Itermax = NbIterations;
202        Perform(F, StartingPoint);
203     }
204                              
205
206     math_FRPR::math_FRPR(math_MultipleVarFunctionWithGradient& F,
207                          const Standard_Real        Tolerance,
208                          const Standard_Integer     NbIterations,
209                          const Standard_Real        ZEPS) 
210                          : TheLocation(1, F.NbVariables()),
211                            TheGradient(1, F.NbVariables()) {
212
213        XTol = Tolerance;
214        EPSZ = ZEPS;
215        Itermax = NbIterations;
216     }
217
218
219     math_FRPR::~math_FRPR()
220     {
221     }
222
223     void math_FRPR::Dump(Standard_OStream& o) const {
224
225        o << "math_FRPR ";
226        if(Done) {
227          o << " Status = Done \n";
228          o << " Location Vector = "<< TheLocation << "\n";
229          o << " Minimum value = " << TheMinimum <<"\n";
230          o << " Number of iterations = " << Iter <<"\n";
231        }
232        else {
233          o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
234        }
235     }
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250