Commit | Line | Data |
---|---|---|
b311480e | 1 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 | 2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e | 3 | // |
973c2be1 | 4 | // This file is part of Open CASCADE Technology software library. |
b311480e | 5 | // |
d5f74e42 | 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 | |
973c2be1 | 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. | |
b311480e | 11 | // |
973c2be1 | 12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. | |
b311480e | 14 | |
0797d9d3 | 15 | //#ifndef OCCT_DEBUG |
7fd59977 | 16 | #define No_Standard_RangeError |
17 | #define No_Standard_OutOfRange | |
18 | #define No_Standard_DimensionError | |
19 | //#endif | |
20 | ||
21 | #include <math_Powell.ixx> | |
22 | #include <math_BracketMinimum.hxx> | |
23 | #include <math_BrentMinimum.hxx> | |
24 | #include <math_Function.hxx> | |
25 | #include <math_MultipleVarFunction.hxx> | |
26 | ||
27 | ||
1ef32e96 RL |
28 | namespace { |
29 | static inline Standard_Real SQR (const Standard_Real a) | |
30 | { | |
31 | return a * a; | |
32 | } | |
33 | } | |
7fd59977 | 34 | |
35 | class DirFunctionBis : public math_Function { | |
36 | ||
37 | math_Vector *P0; | |
38 | math_Vector *Dir; | |
39 | math_Vector *P; | |
40 | math_MultipleVarFunction *F; | |
41 | ||
42 | public : | |
43 | DirFunctionBis(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 | DirFunctionBis::DirFunctionBis(math_Vector& V1, | |
54 | math_Vector& V2, | |
55 | math_Vector& V3, | |
56 | math_MultipleVarFunction& f) { | |
57 | P0 = &V1; | |
58 | Dir = &V2; | |
59 | P = &V3; | |
60 | F = &f; | |
61 | } | |
62 | ||
63 | ||
64 | void DirFunctionBis::Initialize(const math_Vector& p0, | |
65 | const math_Vector& dir) { | |
66 | ||
67 | *P0 = p0; | |
68 | *Dir = dir; | |
69 | } | |
70 | ||
71 | Standard_Boolean DirFunctionBis::Value(const Standard_Real x, Standard_Real& fval) { | |
72 | *P = *Dir; | |
73 | P->Multiply(x); | |
74 | P->Add(*P0); | |
75 | F->Value(*P, fval); | |
76 | return Standard_True; | |
77 | } | |
78 | ||
79 | ||
80 | static Standard_Boolean MinimizeDirection(math_Vector& P, | |
81 | math_Vector& Dir, | |
82 | Standard_Real& Result, | |
83 | DirFunctionBis& F) { | |
84 | ||
85 | Standard_Real ax; | |
86 | Standard_Real xx; | |
87 | Standard_Real bx; | |
88 | ||
89 | ||
90 | F.Initialize(P, Dir); | |
91 | ||
92 | math_BracketMinimum Bracket(F, 0.0, 1.0); | |
93 | if (Bracket.IsDone()) { | |
94 | Bracket.Values(ax, xx, bx); | |
95 | math_BrentMinimum Sol(F, ax, xx, bx, 1.0e-10, 100); | |
96 | if (Sol.IsDone()) { | |
97 | Standard_Real Scale = Sol.Location(); | |
98 | Result = Sol.Minimum(); | |
99 | Dir.Multiply(Scale); | |
100 | P.Add(Dir); | |
101 | return Standard_True; | |
102 | } | |
103 | } | |
104 | return Standard_False; | |
105 | } | |
106 | ||
107 | ||
108 | ||
109 | void math_Powell::Perform(math_MultipleVarFunction& F, | |
110 | const math_Vector& StartingPoint, | |
111 | const math_Matrix& StartingDirections) { | |
112 | ||
113 | ||
114 | Done = Standard_False; | |
7fd59977 | 115 | Standard_Integer i, ibig, j; |
116 | Standard_Real t, fptt, del; | |
117 | Standard_Integer n = TheLocation.Length(); | |
118 | math_Vector pt(1,n); | |
119 | math_Vector ptt(1,n); | |
120 | math_Vector xit(1,n); | |
121 | math_Vector Temp1(1, n); | |
122 | math_Vector Temp2(1, n); | |
123 | math_Vector Temp3(1, n); | |
124 | DirFunctionBis F_Dir(Temp1, Temp2, Temp3, F); | |
125 | ||
126 | TheLocation = StartingPoint; | |
127 | TheDirections = StartingDirections; | |
128 | pt = TheLocation; //sauvegarde du point initial | |
129 | ||
130 | ||
131 | for(Iter = 1; Iter<= Itermax; Iter++) { | |
96a95605 | 132 | F.Value(TheLocation, PreviousMinimum); |
7fd59977 | 133 | ibig = 0; |
134 | del = 0.0; | |
135 | for (i = 1; i <= n; i++){ | |
136 | for(j =1; j<= n; j++) xit(j) = TheDirections(j,i); | |
96a95605 | 137 | F.Value(TheLocation, fptt); |
7fd59977 | 138 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, |
139 | TheMinimum, F_Dir); | |
140 | ||
141 | if (!IsGood) { | |
142 | Done = Standard_False; | |
143 | TheStatus = math_DirectionSearchError; | |
144 | return; | |
145 | } | |
146 | ||
147 | if (fabs(fptt - TheMinimum)> del) { | |
148 | del = fabs(fptt- TheMinimum); | |
149 | ibig = i; | |
150 | } | |
151 | } | |
152 | ||
153 | if (IsSolutionReached(F)) { | |
154 | //Termination criterion | |
155 | State = F.GetStateNumber(); | |
156 | Done = Standard_True; | |
157 | TheStatus = math_OK; | |
158 | return; | |
159 | } | |
160 | ||
161 | if (Iter == Itermax) { | |
162 | Done = Standard_False; | |
163 | TheStatus = math_TooManyIterations; | |
164 | return; | |
165 | } | |
166 | ||
167 | ptt = 2.0 * TheLocation - pt; | |
168 | xit = TheLocation - pt; | |
169 | pt = TheLocation; | |
170 | ||
171 | // Valeur de la fonction au point extrapole: | |
172 | ||
96a95605 | 173 | F.Value(ptt, fptt); |
7fd59977 | 174 | |
175 | if (fptt < PreviousMinimum) { | |
176 | t = 2.0 *(PreviousMinimum -2.0*TheMinimum +fptt)* | |
177 | SQR(PreviousMinimum-TheMinimum -del)-del* | |
178 | SQR(PreviousMinimum-fptt); | |
179 | if (t <0.0) { | |
180 | //Minimisation along the direction | |
181 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, | |
182 | TheMinimum, F_Dir); | |
183 | if(!IsGood) { | |
184 | Done = Standard_False; | |
185 | TheStatus = math_FunctionError; | |
186 | return; | |
187 | } | |
188 | ||
189 | for(j =1; j <= n; j++) { | |
190 | TheDirections(j, ibig)=xit(j); | |
191 | } | |
192 | } | |
193 | } | |
194 | } | |
195 | } | |
196 | ||
197 | Standard_Boolean math_Powell::IsSolutionReached( | |
198 | // math_MultipleVarFunction& F) { | |
199 | math_MultipleVarFunction& ) { | |
200 | ||
201 | return 2.0*fabs(PreviousMinimum - TheMinimum) <= | |
202 | XTol*(fabs(PreviousMinimum)+fabs(TheMinimum) + EPSZ); | |
203 | } | |
204 | ||
205 | ||
206 | ||
207 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
208 | const math_Vector& StartingPoint, | |
209 | const math_Matrix& StartingDirections, | |
210 | const Standard_Real Tolerance, | |
211 | const Standard_Integer NbIterations, | |
212 | const Standard_Real ZEPS) : | |
213 | TheLocation(1, F.NbVariables()), | |
214 | TheDirections(1, F.NbVariables(), | |
215 | 1, F.NbVariables()) { | |
216 | ||
217 | XTol = Tolerance; | |
218 | EPSZ = ZEPS; | |
219 | Itermax = NbIterations; | |
220 | Perform(F, StartingPoint, StartingDirections); | |
221 | } | |
222 | ||
223 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
224 | const Standard_Real Tolerance, | |
225 | const Standard_Integer NbIterations, | |
226 | const Standard_Real ZEPS) : | |
227 | TheLocation(1, F.NbVariables()), | |
228 | TheDirections(1, F.NbVariables(), | |
229 | 1, F.NbVariables()) { | |
230 | ||
231 | XTol = Tolerance; | |
232 | EPSZ = ZEPS; | |
233 | Itermax = NbIterations; | |
234 | } | |
235 | ||
236 | void math_Powell::Delete() | |
237 | {} | |
238 | ||
239 | void math_Powell::Dump(Standard_OStream& o) const { | |
240 | ||
241 | o << "math_Powell resolution:"; | |
242 | if(Done) { | |
243 | o << " Status = Done \n"; | |
244 | o << " Location Vector = "<< TheLocation << "\n"; | |
245 | o << " Minimum value = " << TheMinimum <<"\n"; | |
246 | o << " Number of iterations = " << Iter <<"\n"; | |
247 | } | |
248 | else { | |
249 | o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n"; | |
250 | } | |
251 | } |