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