0023952: Improving thread-safety of intersections, approximations and other modeling...
[occt.git] / src / math / math_Powell.cxx
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
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
32 namespace {
33 static inline Standard_Real SQR (const Standard_Real a)
34 {
35     return a * a;
36 }
37 }
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 }