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