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_BFGS.ixx>
8
9#include <math_BracketMinimum.hxx>
10#include <math_BrentMinimum.hxx>
11#include <math_FunctionWithDerivative.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
27class DirFunction : public math_FunctionWithDerivative {
28
29 math_Vector *P0;
30 math_Vector *Dir;
31 math_Vector *P;
32 math_Vector *G;
34
35public :
36
37 DirFunction(math_Vector& V1,
38 math_Vector& V2,
39 math_Vector& V3,
40 math_Vector& V4,
42
43 void Initialize(const math_Vector& p0, const math_Vector& dir) const;
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,
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
74 }
75
76
77 Standard_Boolean DirFunction::Value(const Standard_Real x,
78 Standard_Real& fval)
79 {
80 *P = *Dir;
81 P->Multiply(x);
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);
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);
104 Standard_Real fval;
105 F->Values(*P, fval, *G);
106 D = (*G).Multiplied(*Dir);
107 return Standard_True;
108 }
109
110
111static 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;
7fd59977 124 dy1 = Gr*Dir;
125 if (dy1 != 0) {
126 Hnr1 = Dir.Norm2();
127 alfa = 0.7*(-F0)/dy1;
128 lambda = 0.015/Sqrt(Hnr1);
129 }
130 else lambda = 1.0;
131 if (lambda > alfa) lambda = alfa;
132 F.Value(lambda, F1);
133 math_BracketMinimum Bracket(F, 0.0, lambda, F0, F1);
134 if(Bracket.IsDone()) {
135 Bracket.Values(ax, xx, bx);
136 Bracket.FunctionValues(Fax, Fxx, Fbx);
137
138 Standard_Integer niter = 100;
139 Standard_Real tol = 1.e-03;
140 math_BrentMinimum Sol(tol, Fxx, niter, 1.e-08);
141 Sol.Perform(F, ax, xx, bx);
142 if(Sol.IsDone()) {
143 Standard_Real Scale = Sol.Location();
144 Result = Sol.Minimum();
145 Dir.Multiply(Scale);
147 return Standard_True;
148 }
149 }
150 return Standard_False;
151 }
152
153
155 const math_Vector& StartingPoint) {
156
157 Standard_Boolean Good;
158 Standard_Integer n = TheLocation.Length();
159 Standard_Integer j, i;
161
162 math_Vector xi(1, n), dg(1, n), hdg(1, n);
163 math_Matrix hessin(1, n, 1, n);
164 hessin.Init(0.0);
165
166 math_Vector Temp1(1, n);
167 math_Vector Temp2(1, n);
168 math_Vector Temp3(1, n);
169 math_Vector Temp4(1, n);
170 DirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
171
172 TheLocation = StartingPoint;
173 Good = F.Values(TheLocation, PreviousMinimum, TheGradient);
174 if(!Good) {
175 Done = Standard_False;
176 TheStatus = math_FunctionError;
177 return;
178 }
179 for(i = 1; i <= n; i++) {
180 hessin(i, i) = 1.0;
182 }
183
184
185 for(nbiter = 1; nbiter <= Itermax; nbiter++) {
186 TheMinimum = PreviousMinimum;
187 Standard_Boolean IsGood = MinimizeDirection(TheLocation, TheMinimum,
189 xi, TheMinimum, F_Dir);
190 if(!IsGood) {
191 Done = Standard_False;
192 TheStatus = math_DirectionSearchError;
193 return;
194 }
195 if(IsSolutionReached(F)) {
196 Done = Standard_True;
197 TheStatus = math_OK;
198 return;
199 }
200 if (nbiter == Itermax) {
201 Done = Standard_False;
202 TheStatus = math_TooManyIterations;
203 return;
204 }
205 PreviousMinimum = TheMinimum;
206
208
209 Good = F.Values(TheLocation, TheMinimum, TheGradient);
210 if(!Good) {
211 Done = Standard_False;
212 TheStatus = math_FunctionError;
213 return;
214 }
215
216 for(i = 1; i <= n; i++) {
217 dg(i) = TheGradient(i) - dg(i);
218 }
219 for(i = 1; i <= n; i++) {
220 hdg(i) = 0.0;
221 for (j = 1; j <= n; j++)
222 hdg(i) += hessin(i, j) * dg(j);
223 }
224 fac = fae = 0.0;
225 for(i = 1; i <= n; i++) {
226 fac += dg(i) * xi(i);
227 fae += dg(i) * hdg(i);
228 }
229 fac = 1.0 / fac;
230 fad = 1.0 / fae;
231
232 for(i = 1; i <= n; i++)
233 dg(i) = fac * xi(i) - fad * hdg(i);
234
235 for(i = 1; i <= n; i++)
236 for(j = 1; j <= n; j++)
237 hessin(i, j) += fac * xi(i) * xi(j)
238 - fad * hdg(i) * hdg(j) + fae * dg(i) * dg(j);
239
240 for(i = 1; i <= n; i++) {
241 xi(i) = 0.0;
242 for (j = 1; j <= n; j++) xi(i) -= hessin(i, j) * TheGradient(j);
243 }
244 }
245 Done = Standard_False;
246 TheStatus = math_TooManyIterations;
247 return;
248 }
249
250 Standard_Boolean math_BFGS::IsSolutionReached(
252
253 return 2.0 * fabs(TheMinimum - PreviousMinimum) <=
254 XTol * (fabs(TheMinimum) + fabs(PreviousMinimum) + EPSZ);
255 }
256
258 const math_Vector& StartingPoint,
259 const Standard_Real Tolerance,
260 const Standard_Integer NbIterations,
261 const Standard_Real ZEPS)
262 : TheLocation(1, StartingPoint.Length()),
264
265 XTol = Tolerance;
266 EPSZ = ZEPS;
267 Itermax = NbIterations;
268 Perform(F, StartingPoint);
269 }
270
272 const Standard_Real Tolerance,
273 const Standard_Integer NbIterations,
274 const Standard_Real ZEPS)
275 : TheLocation(1, F.NbVariables()),
277
278 XTol = Tolerance;
279 EPSZ = ZEPS;
280 Itermax = NbIterations;
281 }
282
283
284 void math_BFGS::Delete()
285 {}
286
287 void math_BFGS::Dump(Standard_OStream& o) const {
288
289 o<< "math_BFGS resolution: ";
290 if(Done) {
291 o << " Status = Done \n";
292 o <<" Location Vector = " << Location() << "\n";
293 o <<" Minimum value = "<< Minimum()<<"\n";
294 o <<" Number of iterations = "<<NbIterations() <<"\n";;
295 }
296 else {
297 o<< " Status = not Done because " << (Standard_Integer)TheStatus << "\n";
298 }
299 }
300
301
302