Commit | Line | Data |
---|---|---|
b311480e | 1 | // Copyright (c) 1997-1999 Matra Datavision |
2 | // Copyright (c) 1999-2012 OPEN CASCADE SAS | |
3 | // | |
4 | // The content of this file is subject to the Open CASCADE Technology Public | |
5 | // License Version 6.5 (the "License"). You may not use the content of this file | |
6 | // except in compliance with the License. Please obtain a copy of the License | |
7 | // at http://www.opencascade.org and read it completely before using this file. | |
8 | // | |
9 | // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its | |
10 | // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. | |
11 | // | |
12 | // The Original Code and all software distributed under the License is | |
13 | // distributed on an "AS IS" basis, without warranty of any kind, and the | |
14 | // Initial Developer hereby disclaims all such warranties, including without | |
15 | // limitation, any warranties of merchantability, fitness for a particular | |
16 | // purpose or non-infringement. Please see the License for the specific terms | |
17 | // and conditions governing the rights and limitations under the License. | |
18 | ||
7fd59977 | 19 | //#ifndef DEB |
20 | #define No_Standard_RangeError | |
21 | #define No_Standard_OutOfRange | |
22 | #define No_Standard_DimensionError | |
23 | //#endif | |
24 | ||
25 | #include <math_Powell.ixx> | |
26 | #include <math_BracketMinimum.hxx> | |
27 | #include <math_BrentMinimum.hxx> | |
28 | #include <math_Function.hxx> | |
29 | #include <math_MultipleVarFunction.hxx> | |
30 | ||
31 | ||
1ef32e96 RL |
32 | namespace { |
33 | static inline Standard_Real SQR (const Standard_Real a) | |
34 | { | |
35 | return a * a; | |
36 | } | |
37 | } | |
7fd59977 | 38 | |
39 | class DirFunctionBis : public math_Function { | |
40 | ||
41 | math_Vector *P0; | |
42 | math_Vector *Dir; | |
43 | math_Vector *P; | |
44 | math_MultipleVarFunction *F; | |
45 | ||
46 | public : | |
47 | DirFunctionBis(math_Vector& V1, | |
48 | math_Vector& V2, | |
49 | math_Vector& V3, | |
50 | math_MultipleVarFunction& f); | |
51 | ||
52 | void Initialize(const math_Vector& p0, const math_Vector& dir); | |
53 | ||
54 | virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval); | |
55 | }; | |
56 | ||
57 | DirFunctionBis::DirFunctionBis(math_Vector& V1, | |
58 | math_Vector& V2, | |
59 | math_Vector& V3, | |
60 | math_MultipleVarFunction& f) { | |
61 | P0 = &V1; | |
62 | Dir = &V2; | |
63 | P = &V3; | |
64 | F = &f; | |
65 | } | |
66 | ||
67 | ||
68 | void DirFunctionBis::Initialize(const math_Vector& p0, | |
69 | const math_Vector& dir) { | |
70 | ||
71 | *P0 = p0; | |
72 | *Dir = dir; | |
73 | } | |
74 | ||
75 | Standard_Boolean DirFunctionBis::Value(const Standard_Real x, Standard_Real& fval) { | |
76 | *P = *Dir; | |
77 | P->Multiply(x); | |
78 | P->Add(*P0); | |
79 | F->Value(*P, fval); | |
80 | return Standard_True; | |
81 | } | |
82 | ||
83 | ||
84 | static Standard_Boolean MinimizeDirection(math_Vector& P, | |
85 | math_Vector& Dir, | |
86 | Standard_Real& Result, | |
87 | DirFunctionBis& F) { | |
88 | ||
89 | Standard_Real ax; | |
90 | Standard_Real xx; | |
91 | Standard_Real bx; | |
92 | ||
93 | ||
94 | F.Initialize(P, Dir); | |
95 | ||
96 | math_BracketMinimum Bracket(F, 0.0, 1.0); | |
97 | if (Bracket.IsDone()) { | |
98 | Bracket.Values(ax, xx, bx); | |
99 | math_BrentMinimum Sol(F, ax, xx, bx, 1.0e-10, 100); | |
100 | if (Sol.IsDone()) { | |
101 | Standard_Real Scale = Sol.Location(); | |
102 | Result = Sol.Minimum(); | |
103 | Dir.Multiply(Scale); | |
104 | P.Add(Dir); | |
105 | return Standard_True; | |
106 | } | |
107 | } | |
108 | return Standard_False; | |
109 | } | |
110 | ||
111 | ||
112 | ||
113 | void math_Powell::Perform(math_MultipleVarFunction& F, | |
114 | const math_Vector& StartingPoint, | |
115 | const math_Matrix& StartingDirections) { | |
116 | ||
117 | ||
118 | Done = Standard_False; | |
119 | Standard_Boolean Ok; | |
120 | Standard_Integer i, ibig, j; | |
121 | Standard_Real t, fptt, del; | |
122 | Standard_Integer n = TheLocation.Length(); | |
123 | math_Vector pt(1,n); | |
124 | math_Vector ptt(1,n); | |
125 | math_Vector xit(1,n); | |
126 | math_Vector Temp1(1, n); | |
127 | math_Vector Temp2(1, n); | |
128 | math_Vector Temp3(1, n); | |
129 | DirFunctionBis F_Dir(Temp1, Temp2, Temp3, F); | |
130 | ||
131 | TheLocation = StartingPoint; | |
132 | TheDirections = StartingDirections; | |
133 | pt = TheLocation; //sauvegarde du point initial | |
134 | ||
135 | ||
136 | for(Iter = 1; Iter<= Itermax; Iter++) { | |
137 | Ok = F.Value(TheLocation, PreviousMinimum); | |
138 | ibig = 0; | |
139 | del = 0.0; | |
140 | for (i = 1; i <= n; i++){ | |
141 | for(j =1; j<= n; j++) xit(j) = TheDirections(j,i); | |
142 | Ok = F.Value(TheLocation, fptt); | |
143 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, | |
144 | TheMinimum, F_Dir); | |
145 | ||
146 | if (!IsGood) { | |
147 | Done = Standard_False; | |
148 | TheStatus = math_DirectionSearchError; | |
149 | return; | |
150 | } | |
151 | ||
152 | if (fabs(fptt - TheMinimum)> del) { | |
153 | del = fabs(fptt- TheMinimum); | |
154 | ibig = i; | |
155 | } | |
156 | } | |
157 | ||
158 | if (IsSolutionReached(F)) { | |
159 | //Termination criterion | |
160 | State = F.GetStateNumber(); | |
161 | Done = Standard_True; | |
162 | TheStatus = math_OK; | |
163 | return; | |
164 | } | |
165 | ||
166 | if (Iter == Itermax) { | |
167 | Done = Standard_False; | |
168 | TheStatus = math_TooManyIterations; | |
169 | return; | |
170 | } | |
171 | ||
172 | ptt = 2.0 * TheLocation - pt; | |
173 | xit = TheLocation - pt; | |
174 | pt = TheLocation; | |
175 | ||
176 | // Valeur de la fonction au point extrapole: | |
177 | ||
178 | Ok = F.Value(ptt, fptt); | |
179 | ||
180 | if (fptt < PreviousMinimum) { | |
181 | t = 2.0 *(PreviousMinimum -2.0*TheMinimum +fptt)* | |
182 | SQR(PreviousMinimum-TheMinimum -del)-del* | |
183 | SQR(PreviousMinimum-fptt); | |
184 | if (t <0.0) { | |
185 | //Minimisation along the direction | |
186 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, | |
187 | TheMinimum, F_Dir); | |
188 | if(!IsGood) { | |
189 | Done = Standard_False; | |
190 | TheStatus = math_FunctionError; | |
191 | return; | |
192 | } | |
193 | ||
194 | for(j =1; j <= n; j++) { | |
195 | TheDirections(j, ibig)=xit(j); | |
196 | } | |
197 | } | |
198 | } | |
199 | } | |
200 | } | |
201 | ||
202 | Standard_Boolean math_Powell::IsSolutionReached( | |
203 | // math_MultipleVarFunction& F) { | |
204 | math_MultipleVarFunction& ) { | |
205 | ||
206 | return 2.0*fabs(PreviousMinimum - TheMinimum) <= | |
207 | XTol*(fabs(PreviousMinimum)+fabs(TheMinimum) + EPSZ); | |
208 | } | |
209 | ||
210 | ||
211 | ||
212 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
213 | const math_Vector& StartingPoint, | |
214 | const math_Matrix& StartingDirections, | |
215 | const Standard_Real Tolerance, | |
216 | const Standard_Integer NbIterations, | |
217 | const Standard_Real ZEPS) : | |
218 | TheLocation(1, F.NbVariables()), | |
219 | TheDirections(1, F.NbVariables(), | |
220 | 1, F.NbVariables()) { | |
221 | ||
222 | XTol = Tolerance; | |
223 | EPSZ = ZEPS; | |
224 | Itermax = NbIterations; | |
225 | Perform(F, StartingPoint, StartingDirections); | |
226 | } | |
227 | ||
228 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
229 | const Standard_Real Tolerance, | |
230 | const Standard_Integer NbIterations, | |
231 | const Standard_Real ZEPS) : | |
232 | TheLocation(1, F.NbVariables()), | |
233 | TheDirections(1, F.NbVariables(), | |
234 | 1, F.NbVariables()) { | |
235 | ||
236 | XTol = Tolerance; | |
237 | EPSZ = ZEPS; | |
238 | Itermax = NbIterations; | |
239 | } | |
240 | ||
241 | void math_Powell::Delete() | |
242 | {} | |
243 | ||
244 | void math_Powell::Dump(Standard_OStream& o) const { | |
245 | ||
246 | o << "math_Powell resolution:"; | |
247 | if(Done) { | |
248 | o << " Status = Done \n"; | |
249 | o << " Location Vector = "<< TheLocation << "\n"; | |
250 | o << " Minimum value = " << TheMinimum <<"\n"; | |
251 | o << " Number of iterations = " << Iter <<"\n"; | |
252 | } | |
253 | else { | |
254 | o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n"; | |
255 | } | |
256 | } |