1 // Created on: 1996-05-03
2 // Created by: Philippe MANGIN
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
18 #define No_Standard_RangeError
19 #define No_Standard_OutOfRange
20 #define No_Standard_DimensionError
24 #include <math_Gauss.hxx>
25 #include <math_Jacobi.hxx>
26 #include <math_MultipleVarFunctionWithHessian.hxx>
27 #include <math_NewtonMinimum.hxx>
28 #include <Precision.hxx>
29 #include <Standard_DimensionError.hxx>
30 #include <StdFail_NotDone.hxx>
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
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),
54 NoConvexTreatement(theWithSingularity),
55 Convex (Standard_True),
56 myIsBoundsDefined(Standard_False),
57 myLeft(1, theFunction.NbVariables(), 0.0),
58 myRight(1, theFunction.NbVariables(), 0.0),
59 Done (Standard_False),
60 Itermax (theNbIterations)
64 //=======================================================================
65 //function : ~math_NewtonMinimum
66 //purpose : Destructor
67 //=======================================================================
68 math_NewtonMinimum::~math_NewtonMinimum()
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)
79 myLeft = theLeftBorder;
80 myRight = theRightBorder;
81 myIsBoundsDefined = Standard_True;
84 //=======================================================================
87 //=======================================================================
88 void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
89 const math_Vector& StartingPoint)
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;
98 Standard_Boolean Ok = Standard_True;
99 Standard_Integer NbConv = 0, ii, Nreduction;
100 Standard_Real VPrecedent, VItere;
102 Done = Standard_True;
106 while ( Ok && (NbConv < 2) ) {
111 Ok = F.Values(*precedent, VPrecedent, TheGradient, TheHessian);
113 Done = Standard_False;
114 TheStatus = math_FunctionError;
118 PreviousMinimum = VPrecedent;
119 TheMinimum = VPrecedent;
122 // Traitement de la non convexite
124 math_Jacobi CalculVP(TheHessian);
125 if ( !CalculVP.IsDone() ) {
126 Done = Standard_False;
127 TheStatus = math_FunctionError;
131 MinEigenValue = CalculVP.Values() ( CalculVP.Values().Min());
132 if ( MinEigenValue < CTol)
134 Convex = Standard_False;
135 if (NoConvexTreatement && // Treatment is allowed.
136 Abs (MinEigenValue) > CTol) // Treatment will have effect.
138 Standard_Real Delta = CTol + 0.1 * Abs(MinEigenValue) - MinEigenValue;
139 for (ii=1; ii<=TheGradient.Length(); ii++)
140 TheHessian(ii, ii) += Delta;
144 Done = Standard_False;
145 TheStatus = math_FunctionError;
152 math_Gauss LU(TheHessian, CTol/100);
154 Done = Standard_False;
155 TheStatus = math_DirectionSearchError;
158 LU.Solve(TheGradient, TheStep);
160 if (myIsBoundsDefined)
162 // Project point on bounds or nullify TheStep coords if point lies on boudary.
164 *suivant = *precedent - TheStep;
165 Standard_Real aMult = RealLast();
166 for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
168 const Standard_Real anAbsStep = Abs(TheStep(anIdx));
169 if (anAbsStep < gp::Resolution())
172 if (suivant->Value(anIdx) < myLeft(anIdx))
174 Standard_Real aValue = Abs(precedent->Value(anIdx) - myLeft(anIdx)) / anAbsStep;
175 aMult = Min (aValue, aMult);
178 if (suivant->Value(anIdx) > myRight(anIdx))
180 Standard_Real aValue = Abs(precedent->Value(anIdx) - myRight(anIdx)) / anAbsStep;
181 aMult = Min (aValue, aMult);
185 if (aMult != RealLast())
187 if (aMult > Precision::PConfusion())
189 // Project point into param space.
194 // Old point on border and new point out of border:
195 // Nullify corresponding TheStep indexes.
196 for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
198 if (Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() ||
199 Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion())
201 TheStep(anIdx) = 0.0;
208 Standard_Boolean hasProblem = Standard_False;
211 *suivant = *precedent - TheStep;
213 // Gestion de la convergence
214 hasProblem = !(F.Value(*suivant, TheMinimum));
220 } while (hasProblem);
222 if (IsConverged()) { NbConv++; }
225 // Controle et corrections.
228 TheMinimum = PreviousMinimum;
230 while (VItere > VPrecedent && Nreduction < 10) {
232 *suivant = *precedent - TheStep;
233 F.Value(*suivant, VItere);
237 if (VItere <= VPrecedent) {
238 auxiliaire = precedent;
240 suivant = auxiliaire;
241 PreviousMinimum = VPrecedent;
243 Ok = (nbiter < Itermax);
244 if (!Ok && NbConv < 2) TheStatus = math_TooManyIterations;
248 TheStatus = math_DirectionSearchError;
251 TheLocation = *precedent;
254 //=======================================================================
257 //=======================================================================
258 void math_NewtonMinimum::Dump(Standard_OStream& o) const
260 o<< "math_Newton Optimisation: ";
261 o << " Done =" << Done << std::endl;
262 o << " Status = " << (Standard_Integer)TheStatus << std::endl;
263 o << " Location Vector = " << Location() << std::endl;
264 o << " Minimum value = "<< Minimum()<< std::endl;
265 o << " Previous value = "<< PreviousMinimum << std::endl;
266 o << " Number of iterations = " <<NbIterations() << std::endl;
267 o << " Convexity = " << Convex << std::endl;
268 o << " Eigen Value = " << MinEigenValue << std::endl;