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