0022904: Clean up sccsid variables
[occt.git] / src / math / math_Powell.cxx
1 //#ifndef DEB
2 #define No_Standard_RangeError
3 #define No_Standard_OutOfRange
4 #define No_Standard_DimensionError
5 //#endif
6
7 #include <math_Powell.ixx>
8 #include <math_BracketMinimum.hxx>
9 #include <math_BrentMinimum.hxx>
10 #include <math_Function.hxx>
11 #include <math_MultipleVarFunction.hxx>
12
13
14 static Standard_Real sqrarg;
15 #define SQR(a) (sqrarg=(a), sqrarg*sqrarg)
16
17 class DirFunctionBis : public math_Function {
18
19   math_Vector *P0;
20   math_Vector *Dir;
21   math_Vector *P;
22   math_MultipleVarFunction *F;
23
24   public :
25     DirFunctionBis(math_Vector& V1,
26                 math_Vector& V2,
27                 math_Vector& V3,
28                 math_MultipleVarFunction& f);
29
30     void Initialize(const math_Vector& p0, const math_Vector& dir);
31
32     virtual Standard_Boolean Value(const Standard_Real x, Standard_Real& fval);
33 };
34
35  DirFunctionBis::DirFunctionBis(math_Vector& V1,
36                           math_Vector& V2,
37                           math_Vector& V3,
38                           math_MultipleVarFunction& f) {
39    P0 = &V1;
40    Dir = &V2;
41    P = &V3;
42    F = &f;
43  }
44
45
46 void DirFunctionBis::Initialize(const math_Vector& p0,
47                              const math_Vector& dir) {
48
49   *P0 = p0;
50   *Dir = dir;
51 }
52
53 Standard_Boolean DirFunctionBis::Value(const Standard_Real x, Standard_Real& fval) {
54   *P = *Dir;
55   P->Multiply(x);
56   P->Add(*P0);
57   F->Value(*P, fval);
58   return Standard_True;
59 }
60
61
62 static Standard_Boolean MinimizeDirection(math_Vector& P,
63                                  math_Vector& Dir,
64                                  Standard_Real& Result,
65                                  DirFunctionBis& F) {
66
67   Standard_Real ax;
68   Standard_Real xx;
69   Standard_Real bx;
70
71
72   F.Initialize(P, Dir);
73
74   math_BracketMinimum Bracket(F, 0.0, 1.0);
75   if (Bracket.IsDone()) {
76     Bracket.Values(ax, xx, bx);
77     math_BrentMinimum Sol(F, ax, xx, bx, 1.0e-10, 100);
78     if (Sol.IsDone()) {
79       Standard_Real Scale = Sol.Location();
80       Result = Sol.Minimum();
81       Dir.Multiply(Scale);
82       P.Add(Dir);
83       return Standard_True;
84     }
85   }
86   return Standard_False;
87 }
88
89
90
91 void math_Powell::Perform(math_MultipleVarFunction& F,
92                           const math_Vector& StartingPoint,
93                           const math_Matrix& StartingDirections) {
94
95
96   Done = Standard_False;
97   Standard_Boolean Ok;
98   Standard_Integer i, ibig, j;
99   Standard_Real t, fptt, del;
100   Standard_Integer n = TheLocation.Length();
101   math_Vector pt(1,n);
102   math_Vector ptt(1,n);
103   math_Vector xit(1,n);
104   math_Vector Temp1(1, n);
105   math_Vector Temp2(1, n);
106   math_Vector Temp3(1, n);
107   DirFunctionBis F_Dir(Temp1, Temp2, Temp3, F);
108
109   TheLocation = StartingPoint;
110   TheDirections = StartingDirections;
111   pt = TheLocation;  //sauvegarde du point initial
112
113
114   for(Iter = 1; Iter<= Itermax; Iter++) {
115     Ok = F.Value(TheLocation, PreviousMinimum);
116     ibig = 0;
117     del = 0.0;
118     for (i = 1; i <= n; i++){
119       for(j =1; j<= n; j++) xit(j) = TheDirections(j,i);
120       Ok = F.Value(TheLocation, fptt); 
121       Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit, 
122                                          TheMinimum, F_Dir);
123
124       if (!IsGood) {
125         Done = Standard_False;
126         TheStatus = math_DirectionSearchError;
127         return;
128       }
129
130       if (fabs(fptt - TheMinimum)> del) {
131         del = fabs(fptt- TheMinimum);
132         ibig = i;
133       }
134     }
135
136     if (IsSolutionReached(F)) {
137       //Termination criterion
138       State = F.GetStateNumber();
139       Done = Standard_True;
140       TheStatus = math_OK;
141       return;
142     }
143
144     if (Iter == Itermax) {
145       Done = Standard_False;
146       TheStatus = math_TooManyIterations;
147       return;
148     }
149
150       ptt = 2.0 * TheLocation - pt;
151       xit = TheLocation - pt;
152       pt = TheLocation;
153     
154     // Valeur de la fonction au point extrapole:
155
156     Ok = F.Value(ptt, fptt);
157
158     if (fptt < PreviousMinimum) {
159       t = 2.0 *(PreviousMinimum -2.0*TheMinimum +fptt)*
160         SQR(PreviousMinimum-TheMinimum -del)-del*
161           SQR(PreviousMinimum-fptt);
162       if (t <0.0) {
163         //Minimisation along the direction
164         Standard_Boolean IsGood = MinimizeDirection(TheLocation, xit,
165                                            TheMinimum, F_Dir);
166         if(!IsGood) {
167           Done = Standard_False;
168           TheStatus = math_FunctionError;
169           return;
170         }
171         
172         for(j =1; j <= n; j++) {
173           TheDirections(j, ibig)=xit(j);
174         }
175       }
176     }
177   }
178 }
179                                            
180 Standard_Boolean math_Powell::IsSolutionReached(
181 //                       math_MultipleVarFunction& F) {
182                          math_MultipleVarFunction& ) {
183   
184   return 2.0*fabs(PreviousMinimum - TheMinimum) <= 
185     XTol*(fabs(PreviousMinimum)+fabs(TheMinimum) + EPSZ);
186 }
187
188
189
190 math_Powell::math_Powell(math_MultipleVarFunction& F,
191                          const math_Vector& StartingPoint,
192                          const math_Matrix& StartingDirections,
193                          const Standard_Real Tolerance,
194                          const Standard_Integer NbIterations,
195                          const Standard_Real ZEPS) :
196                            TheLocation(1, F.NbVariables()),
197                            TheDirections(1, F.NbVariables(),
198                                          1, F.NbVariables()) {
199
200     XTol = Tolerance;
201     EPSZ = ZEPS;
202     Itermax = NbIterations;
203     Perform(F, StartingPoint, StartingDirections);
204   }
205
206 math_Powell::math_Powell(math_MultipleVarFunction& F,
207                          const Standard_Real Tolerance,
208                          const Standard_Integer NbIterations,
209                          const Standard_Real ZEPS) :
210                            TheLocation(1, F.NbVariables()),
211                            TheDirections(1, F.NbVariables(),
212                                          1, F.NbVariables()) {
213
214     XTol = Tolerance;
215     EPSZ = ZEPS;
216     Itermax = NbIterations;
217   }
218
219 void math_Powell::Delete()
220 {}
221
222 void math_Powell::Dump(Standard_OStream& o) const {
223
224   o << "math_Powell resolution:";
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 }