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 | ||
6da30ff1 | 107 | math_Powell::~math_Powell() |
108 | { | |
109 | } | |
7fd59977 | 110 | |
111 | void math_Powell::Perform(math_MultipleVarFunction& F, | |
112 | const math_Vector& StartingPoint, | |
113 | const math_Matrix& StartingDirections) { | |
114 | ||
115 | ||
116 | Done = Standard_False; | |
7fd59977 | 117 | Standard_Integer i, ibig, j; |
118 | Standard_Real t, fptt, del; | |
119 | Standard_Integer n = TheLocation.Length(); | |
120 | math_Vector pt(1,n); | |
121 | math_Vector ptt(1,n); | |
122 | math_Vector xit(1,n); | |
123 | math_Vector Temp1(1, n); | |
124 | math_Vector Temp2(1, n); | |
125 | math_Vector Temp3(1, n); | |
126 | DirFunctionBis F_Dir(Temp1, Temp2, Temp3, F); | |
127 | ||
128 | TheLocation = StartingPoint; | |
129 | TheDirections = StartingDirections; | |
130 | pt = TheLocation; //sauvegarde du point initial | |
131 | ||
132 | ||
133 | for(Iter = 1; Iter<= Itermax; Iter++) { | |
96a95605 | 134 | F.Value(TheLocation, PreviousMinimum); |
7fd59977 | 135 | ibig = 0; |
136 | del = 0.0; | |
137 | for (i = 1; i <= n; i++){ | |
138 | for(j =1; j<= n; j++) xit(j) = TheDirections(j,i); | |
96a95605 | 139 | F.Value(TheLocation, fptt); |
7fd59977 | 140 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, |
141 | TheMinimum, F_Dir); | |
142 | ||
143 | if (!IsGood) { | |
144 | Done = Standard_False; | |
145 | TheStatus = math_DirectionSearchError; | |
146 | return; | |
147 | } | |
148 | ||
149 | if (fabs(fptt - TheMinimum)> del) { | |
150 | del = fabs(fptt- TheMinimum); | |
151 | ibig = i; | |
152 | } | |
153 | } | |
154 | ||
155 | if (IsSolutionReached(F)) { | |
156 | //Termination criterion | |
157 | State = F.GetStateNumber(); | |
158 | Done = Standard_True; | |
159 | TheStatus = math_OK; | |
160 | return; | |
161 | } | |
162 | ||
163 | if (Iter == Itermax) { | |
164 | Done = Standard_False; | |
165 | TheStatus = math_TooManyIterations; | |
166 | return; | |
167 | } | |
168 | ||
169 | ptt = 2.0 * TheLocation - pt; | |
170 | xit = TheLocation - pt; | |
171 | pt = TheLocation; | |
172 | ||
173 | // Valeur de la fonction au point extrapole: | |
174 | ||
96a95605 | 175 | F.Value(ptt, fptt); |
7fd59977 | 176 | |
177 | if (fptt < PreviousMinimum) { | |
178 | t = 2.0 *(PreviousMinimum -2.0*TheMinimum +fptt)* | |
179 | SQR(PreviousMinimum-TheMinimum -del)-del* | |
180 | SQR(PreviousMinimum-fptt); | |
181 | if (t <0.0) { | |
182 | //Minimisation along the direction | |
183 | Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, | |
184 | TheMinimum, F_Dir); | |
185 | if(!IsGood) { | |
186 | Done = Standard_False; | |
187 | TheStatus = math_FunctionError; | |
188 | return; | |
189 | } | |
190 | ||
191 | for(j =1; j <= n; j++) { | |
192 | TheDirections(j, ibig)=xit(j); | |
193 | } | |
194 | } | |
195 | } | |
196 | } | |
197 | } | |
198 | ||
199 | Standard_Boolean math_Powell::IsSolutionReached( | |
200 | // math_MultipleVarFunction& F) { | |
201 | math_MultipleVarFunction& ) { | |
202 | ||
203 | return 2.0*fabs(PreviousMinimum - TheMinimum) <= | |
204 | XTol*(fabs(PreviousMinimum)+fabs(TheMinimum) + EPSZ); | |
205 | } | |
206 | ||
207 | ||
208 | ||
209 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
210 | const math_Vector& StartingPoint, | |
211 | const math_Matrix& StartingDirections, | |
212 | const Standard_Real Tolerance, | |
213 | const Standard_Integer NbIterations, | |
214 | const Standard_Real ZEPS) : | |
215 | TheLocation(1, F.NbVariables()), | |
216 | TheDirections(1, F.NbVariables(), | |
217 | 1, F.NbVariables()) { | |
218 | ||
219 | XTol = Tolerance; | |
220 | EPSZ = ZEPS; | |
221 | Itermax = NbIterations; | |
222 | Perform(F, StartingPoint, StartingDirections); | |
223 | } | |
224 | ||
225 | math_Powell::math_Powell(math_MultipleVarFunction& F, | |
226 | const Standard_Real Tolerance, | |
227 | const Standard_Integer NbIterations, | |
228 | const Standard_Real ZEPS) : | |
229 | TheLocation(1, F.NbVariables()), | |
230 | TheDirections(1, F.NbVariables(), | |
231 | 1, F.NbVariables()) { | |
232 | ||
233 | XTol = Tolerance; | |
234 | EPSZ = ZEPS; | |
235 | Itermax = NbIterations; | |
6da30ff1 | 236 | } |
7fd59977 | 237 | |
238 | void math_Powell::Dump(Standard_OStream& o) const { | |
239 | ||
240 | o << "math_Powell resolution:"; | |
241 | if(Done) { | |
242 | o << " Status = Done \n"; | |
243 | o << " Location Vector = "<< TheLocation << "\n"; | |
244 | o << " Minimum value = " << TheMinimum <<"\n"; | |
245 | o << " Number of iterations = " << Iter <<"\n"; | |
246 | } | |
247 | else { | |
248 | o << " Status = not Done because " << (Standard_Integer)TheStatus << "\n"; | |
249 | } | |
250 | } |