b311480e |
1 | // Created on: 1996-05-03 |
2 | // Created by: Philippe MANGIN |
3 | // Copyright (c) 1996-1999 Matra Datavision |
973c2be1 |
4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
5 | // |
973c2be1 |
6 | // This file is part of Open CASCADE Technology software library. |
b311480e |
7 | // |
d5f74e42 |
8 | // This library is free software; you can redistribute it and/or modify it under |
9 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
10 | // by the Free Software Foundation, with special exception defined in the file |
11 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
12 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
13 | // |
973c2be1 |
14 | // Alternatively, this file may be used under the terms of Open CASCADE |
15 | // commercial license or contractual agreement. |
7fd59977 |
16 | |
0797d9d3 |
17 | //#ifndef OCCT_DEBUG |
7fd59977 |
18 | #define No_Standard_RangeError |
19 | #define No_Standard_OutOfRange |
20 | #define No_Standard_DimensionError |
7fd59977 |
21 | |
42cf5bc1 |
22 | //#endif |
7fd59977 |
23 | |
24 | #include <math_Gauss.hxx> |
25 | #include <math_Jacobi.hxx> |
42cf5bc1 |
26 | #include <math_MultipleVarFunctionWithHessian.hxx> |
27 | #include <math_NewtonMinimum.hxx> |
28 | #include <Precision.hxx> |
29 | #include <Standard_DimensionError.hxx> |
30 | #include <StdFail_NotDone.hxx> |
7fd59977 |
31 | |
859a47c3 |
32 | //======================================================================= |
33 | //function : math_NewtonMinimum |
34 | //purpose : Constructor |
35 | //======================================================================= |
36 | math_NewtonMinimum::math_NewtonMinimum( |
37 | const math_MultipleVarFunctionWithHessian& theFunction, |
38 | const Standard_Real theTolerance, |
39 | const Standard_Integer theNbIterations, |
40 | const Standard_Real theConvexity, |
41 | const Standard_Boolean theWithSingularity |
42 | ) |
43 | : TheStatus (math_NotBracketed), |
44 | TheLocation(1, theFunction.NbVariables()), |
45 | TheGradient(1, theFunction.NbVariables()), |
46 | TheStep (1, theFunction.NbVariables(), 10.0 * theTolerance), |
47 | TheHessian (1, theFunction.NbVariables(), 1, theFunction.NbVariables()), |
48 | PreviousMinimum (0.0), |
49 | TheMinimum (0.0), |
50 | MinEigenValue (0.0), |
51 | XTol (theTolerance), |
52 | CTol (theConvexity), |
53 | nbiter (0), |
54 | NoConvexTreatement(theWithSingularity), |
55 | Convex (Standard_True), |
91806b90 |
56 | myIsBoundsDefined(Standard_False), |
57 | myLeft(1, theFunction.NbVariables(), 0.0), |
58 | myRight(1, theFunction.NbVariables(), 0.0), |
859a47c3 |
59 | Done (Standard_False), |
60 | Itermax (theNbIterations) |
7fd59977 |
61 | { |
7fd59977 |
62 | } |
63 | |
859a47c3 |
64 | //======================================================================= |
65 | //function : ~math_NewtonMinimum |
66 | //purpose : Destructor |
67 | //======================================================================= |
6da30ff1 |
68 | math_NewtonMinimum::~math_NewtonMinimum() |
69 | { |
70 | } |
71 | |
91806b90 |
72 | //======================================================================= |
73 | //function : SetBoundary |
74 | //purpose : Set boundaries for conditional optimization |
75 | //======================================================================= |
76 | void math_NewtonMinimum::SetBoundary(const math_Vector& theLeftBorder, |
77 | const math_Vector& theRightBorder) |
78 | { |
79 | myLeft = theLeftBorder; |
80 | myRight = theRightBorder; |
81 | myIsBoundsDefined = Standard_True; |
82 | } |
83 | |
859a47c3 |
84 | //======================================================================= |
85 | //function : Perform |
86 | //purpose : |
87 | //======================================================================= |
7fd59977 |
88 | void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, |
859a47c3 |
89 | const math_Vector& StartingPoint) |
7fd59977 |
90 | { |
91 | math_Vector Point1 (1, F.NbVariables()); |
92 | Point1 = StartingPoint; |
93 | math_Vector Point2(1, F.NbVariables()); |
94 | math_Vector* precedent = &Point1; |
95 | math_Vector* suivant = &Point2; |
96 | math_Vector* auxiliaire = precedent; |
97 | |
98 | Standard_Boolean Ok = Standard_True; |
99 | Standard_Integer NbConv = 0, ii, Nreduction; |
100 | Standard_Real VPrecedent, VItere; |
101 | |
102 | Done = Standard_True; |
103 | TheStatus = math_OK; |
104 | nbiter = 0; |
105 | |
106 | while ( Ok && (NbConv < 2) ) { |
107 | nbiter++; |
108 | |
109 | // Positionnement |
110 | |
111 | Ok = F.Values(*precedent, VPrecedent, TheGradient, TheHessian); |
112 | if (!Ok) { |
113 | Done = Standard_False; |
114 | TheStatus = math_FunctionError; |
115 | return; |
116 | } |
117 | if (nbiter==1) { |
118 | PreviousMinimum = VPrecedent; |
119 | TheMinimum = VPrecedent; |
120 | } |
121 | |
122 | // Traitement de la non convexite |
123 | |
124 | math_Jacobi CalculVP(TheHessian); |
125 | if ( !CalculVP.IsDone() ) { |
126 | Done = Standard_False; |
127 | TheStatus = math_FunctionError; |
128 | return; |
129 | } |
130 | |
7fd59977 |
131 | MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min()); |
f79b19a1 |
132 | if ( MinEigenValue < CTol) |
133 | { |
134 | Convex = Standard_False; |
135 | if (NoConvexTreatement && // Treatment is allowed. |
136 | Abs (MinEigenValue) > CTol) // Treatment will have effect. |
137 | { |
138 | Standard_Real Delta = CTol + 0.1 * Abs(MinEigenValue) - MinEigenValue; |
139 | for (ii=1; ii<=TheGradient.Length(); ii++) |
140 | TheHessian(ii, ii) += Delta; |
141 | } |
142 | else |
143 | { |
144 | Done = Standard_False; |
145 | TheStatus = math_FunctionError; |
146 | return; |
147 | } |
148 | } |
7fd59977 |
149 | |
150 | // Schemas de Newton |
151 | |
152 | math_Gauss LU(TheHessian, CTol/100); |
153 | if ( !LU.IsDone()) { |
154 | Done = Standard_False; |
155 | TheStatus = math_DirectionSearchError; |
156 | return; |
157 | } |
7fd59977 |
158 | LU.Solve(TheGradient, TheStep); |
91806b90 |
159 | |
160 | if (myIsBoundsDefined) |
161 | { |
162 | // Project point on bounds or nullify TheStep coords if point lies on boudary. |
163 | |
164 | *suivant = *precedent - TheStep; |
165 | Standard_Real aMult = RealLast(); |
166 | for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++) |
167 | { |
168 | if (suivant->Value(anIdx) < myLeft(anIdx)) |
169 | { |
170 | Standard_Real aValue = Abs(precedent->Value(anIdx) - myLeft(anIdx)) / Abs(TheStep(anIdx)); |
171 | aMult = Min (aValue, aMult); |
172 | } |
173 | |
174 | if (suivant->Value(anIdx) > myRight(anIdx)) |
175 | { |
176 | Standard_Real aValue = Abs(precedent->Value(anIdx) - myRight(anIdx)) / Abs(TheStep(anIdx)); |
177 | aMult = Min (aValue, aMult); |
178 | } |
179 | } |
180 | |
181 | if (aMult != RealLast()) |
182 | { |
183 | if (aMult > Precision::PConfusion()) |
184 | { |
185 | // Project point into param space. |
186 | TheStep *= aMult; |
187 | } |
188 | else |
189 | { |
190 | // Old point on border and new point out of border: |
191 | // Nullify corresponding TheStep indexes. |
192 | for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++) |
193 | { |
194 | if (Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() || |
195 | Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion()) |
196 | { |
197 | TheStep(anIdx) = 0.0; |
198 | } |
199 | } |
200 | } |
201 | } |
202 | } |
203 | |
4bbaf12b |
204 | Standard_Boolean hasProblem = Standard_False; |
205 | do |
206 | { |
207 | *suivant = *precedent - TheStep; |
208 | |
209 | // Gestion de la convergence |
210 | hasProblem = !(F.Value(*suivant, TheMinimum)); |
211 | |
212 | if (hasProblem) |
213 | { |
214 | TheStep /= 2.0; |
215 | } |
216 | } while (hasProblem); |
7fd59977 |
217 | |
218 | if (IsConverged()) { NbConv++; } |
219 | else { NbConv=0; } |
220 | |
221 | // Controle et corrections. |
222 | |
223 | VItere = TheMinimum; |
224 | TheMinimum = PreviousMinimum; |
225 | Nreduction =0; |
226 | while (VItere > VPrecedent && Nreduction < 10) { |
227 | TheStep *= 0.4; |
228 | *suivant = *precedent - TheStep; |
229 | F.Value(*suivant, VItere); |
230 | Nreduction++; |
231 | } |
232 | |
233 | if (VItere <= VPrecedent) { |
234 | auxiliaire = precedent; |
235 | precedent = suivant; |
236 | suivant = auxiliaire; |
237 | PreviousMinimum = VPrecedent; |
238 | TheMinimum = VItere; |
239 | Ok = (nbiter < Itermax); |
240 | if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations; |
241 | } |
242 | else { |
243 | Ok = Standard_False; |
244 | TheStatus = math_DirectionSearchError; |
245 | } |
246 | } |
247 | TheLocation = *precedent; |
248 | } |
249 | |
859a47c3 |
250 | //======================================================================= |
251 | //function : Dump |
252 | //purpose : |
253 | //======================================================================= |
7fd59977 |
254 | void math_NewtonMinimum::Dump(Standard_OStream& o) const |
7fd59977 |
255 | { |
859a47c3 |
256 | o<< "math_Newton Optimisation: "; |
257 | o << " Done =" << Done << endl; |
258 | o << " Status = " << (Standard_Integer)TheStatus << endl; |
259 | o << " Location Vector = " << Location() << endl; |
260 | o << " Minimum value = "<< Minimum()<< endl; |
261 | o << " Previous value = "<< PreviousMinimum << endl; |
262 | o << " Number of iterations = " <<NbIterations() << endl; |
263 | o << " Convexity = " << Convex << endl; |
264 | o << " Eigen Value = " << MinEigenValue << endl; |
7fd59977 |
265 | } |
266 | |