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