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;
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).
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.
Brents1DMinimization()
Construct with default parameters.
RCP< const ParameterList > getValidParameters() const
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.