44 #ifndef GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP 45 #define GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP 48 #include "GlobiPack_Brents1DMinimization_decl.hpp" 49 #include "Teuchos_TabularOutputter.hpp" 58 template<
typename Scalar>
60 :rel_tol_(Brents1DMinimizationUtils::rel_tol_default),
61 bracket_tol_(Brents1DMinimizationUtils::bracket_tol_default),
62 max_iters_(Brents1DMinimizationUtils::max_iters_default)
69 template<
typename Scalar>
72 typedef ScalarTraits<Scalar> ST;
73 namespace BMU = Brents1DMinimizationUtils;
74 using Teuchos::getParameter;
75 paramList->validateParametersAndSetDefaults(*this->getValidParameters());
76 rel_tol_ = getParameter<double>(*paramList, BMU::rel_tol_name);
77 bracket_tol_ = getParameter<double>(*paramList, BMU::bracket_tol_name);
78 max_iters_ = getParameter<int>(*paramList, BMU::max_iters_name);
79 TEUCHOS_ASSERT_INEQUALITY( rel_tol_, >, ST::zero() );
80 TEUCHOS_ASSERT_INEQUALITY( bracket_tol_, >, ST::zero() );
81 TEUCHOS_ASSERT_INEQUALITY( max_iters_, >=, 0 );
82 setMyParamList(paramList);
86 template<
typename Scalar>
89 namespace BMU = Brents1DMinimizationUtils;
90 static RCP<const ParameterList> validPL;
91 if (is_null(validPL)) {
92 RCP<Teuchos::ParameterList>
93 pl = Teuchos::rcp(
new Teuchos::ParameterList());
94 pl->set( BMU::rel_tol_name, BMU::rel_tol_default );
95 pl->set( BMU::bracket_tol_name, BMU::bracket_tol_default );
96 pl->set( BMU::max_iters_name, BMU::max_iters_default );
106 template<
typename Scalar>
112 const Ptr<int> &numIters
117 using Teuchos::TabularOutputter;
118 typedef Teuchos::TabularOutputter TO;
119 typedef ScalarTraits<Scalar> ST;
120 using Teuchos::OSTab;
123 #endif // TEUCHOS_DEBUG 128 TEUCHOS_TEST_FOR_EXCEPT(is_null(pointMiddle));
129 TEUCHOS_ASSERT_INEQUALITY(pointLower.
alpha, <, pointMiddle->alpha);
130 TEUCHOS_ASSERT_INEQUALITY(pointMiddle->alpha, <, pointUpper.
alpha);
131 TEUCHOS_ASSERT_INEQUALITY(pointLower.
phi, !=, PE1D::valNotGiven());
132 TEUCHOS_ASSERT_INEQUALITY(pointMiddle->phi, !=, PE1D::valNotGiven());
133 TEUCHOS_ASSERT_INEQUALITY(pointUpper.
phi, !=, PE1D::valNotGiven());
134 #endif // TEUCHOS_DEBUG 136 const RCP<Teuchos::FancyOStream> out = this->getOStream();
138 *out <<
"\nStarting Brent's 1D minimization algorithm ...\n\n";
140 TabularOutputter tblout(out);
142 tblout.pushFieldSpec(
"itr", TO::INT);
143 tblout.pushFieldSpec(
"alpha_a", TO::DOUBLE);
144 tblout.pushFieldSpec(
"alpha_min", TO::DOUBLE);
145 tblout.pushFieldSpec(
"alpha_b", TO::DOUBLE);
146 tblout.pushFieldSpec(
"phi(alpha_min)", TO::DOUBLE);
147 tblout.pushFieldSpec(
"alpha_b - alpha_a", TO::DOUBLE);
148 tblout.pushFieldSpec(
"alpha_min - alpha_avg", TO::DOUBLE);
149 tblout.pushFieldSpec(
"tol", TO::DOUBLE);
151 tblout.outputHeader();
153 const Scalar INV_GOLD2=0.3819660112501051518;
154 const Scalar TINY = ST::squareroot(ST::eps());
156 const Scalar alpha_l = pointLower.
alpha, phi_l = pointLower.
phi;
157 Scalar &alpha_m = pointMiddle->alpha, &phi_m = pointMiddle->phi;
158 const Scalar alpha_u = pointUpper.
alpha, phi_u = pointUpper.
phi;
160 Scalar d = ST::nan();
161 Scalar e = ST::nan();
162 Scalar u = ST::nan();
164 Scalar phi_w = min(phi_l, phi_u);
166 Scalar alpha_v = ST::nan();
167 Scalar alpha_w = ST::nan();
168 Scalar phi_v = ST::nan();
181 Scalar alpha_min = alpha_m;
182 Scalar phi_min = phi_m;
183 Scalar alpha_a = alpha_l;
184 Scalar alpha_b = alpha_u;
186 bool foundMin =
false;
190 for ( ; iteration <= max_iters_; ++iteration) {
193 e = 2.0 * (alpha_b - alpha_a);
195 const Scalar alpha_avg = 0.5 *(alpha_a + alpha_b);
196 const Scalar tol1 = rel_tol_ * ST::magnitude(alpha_min) + TINY;
197 const Scalar tol2 = 2.0 * tol1;
199 const Scalar step_diff = alpha_min - alpha_avg;
200 const Scalar step_diff_tol = (tol2 + bracket_tol_ * (alpha_b - alpha_a));
208 tblout.outputField(iteration);
209 tblout.outputField(alpha_a);
210 tblout.outputField(alpha_min);
211 tblout.outputField(alpha_b);
212 tblout.outputField(phi_min);
213 tblout.outputField(alpha_b - alpha_a);
214 tblout.outputField(step_diff);
215 tblout.outputField(step_diff_tol);
225 ST::magnitude(step_diff) <= step_diff_tol
237 if (ST::magnitude(e) > tol1 || iteration < 2) {
239 const Scalar r = (alpha_min - alpha_w) * (phi_min - phi_v);
240 Scalar q = (alpha_min - alpha_v) * (phi_min - phi_w);
241 Scalar p = (alpha_min - alpha_v) * q - (alpha_min - alpha_w) * r;
245 q = ST::magnitude(q);
246 const Scalar etemp = e;
249 if ( ST::magnitude(p) >= ST::magnitude(0.5 * q * etemp)
250 || p <= q * (alpha_a - alpha_min)
251 || p >= q * (alpha_b - alpha_min)
254 e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
260 if (u - alpha_a < tol2 || alpha_b - u < tol2)
262 d = ( alpha_avg - alpha_min > ST::zero()
263 ? ST::magnitude(tol1)
264 : -ST::magnitude(tol1) );
270 e = (alpha_min >= alpha_avg ? alpha_a - alpha_min : alpha_b - alpha_min);
275 u = ( ST::magnitude(d) >= tol1
277 : alpha_min + (d >= 0 ? ST::magnitude(tol1) : -ST::magnitude(tol1))
280 const Scalar phi_eval_u = computeValue<Scalar>(phi, u);
282 if (phi_eval_u <= phi_min) {
294 phi_min = phi_eval_u;
304 if (phi_eval_u <= phi_w || alpha_w == alpha_min) {
310 else if (phi_eval_u <= phi_v || alpha_v == alpha_min || alpha_v == alpha_w) {
320 if (!is_null(numIters))
321 *numIters = iteration;
324 *out <<
"\nFound the minimum alpha="<<alpha_m<<
", phi(alpha)="<<phi_m<<
"\n";
327 *out <<
"\nExceeded maximum number of iterations!\n";
340 #endif // GLOBIPACK_BRENTS_1D_MINIMIZATION_DEF_HPP Scalar phi
The value of the merit function phi(alpha).
Brents1DMinimization()
Construct with default parameters.
Scalar alpha
The value of the unknown alpha.
Represents the evaluation point of the merit function phi(alpha) and/or is derivative Dphi(alpha)...
void setParameterList(RCP< ParameterList > const ¶mList)
Base class for 1D merit fucntions used in globalization methods.
RCP< const ParameterList > getValidParameters() const
bool approxMinimize(const MeritFunc1DBase< Scalar > &phi, const PointEval1D< Scalar > &pointLower, const Ptr< PointEval1D< Scalar > > &pointMiddle, const PointEval1D< Scalar > &pointUpper, const Ptr< int > &numIters=Teuchos::null) const
Approximatly mimimize a 1D function.