ReClaM: BFGS.cpp Source File

00001 //===========================================================================
00073 //===========================================================================
00074 
00075 
00076 #include <cmath>
00077 #include <iostream>
00078 #include <float.h>
00079 #include <SharkDefs.h>
00080 #include <ReClaM/BFGS.h>
00081 
00082 
00083 #define DLINMIN
00084 
00085 
00086 //===========================================================================
00116 double BFGS::optimize(
00117     Model& model,
00118     ErrorFunction& errorfunction,
00119     const Array<double>& input,
00120     const Array<double>& target)
00121 {
00122     bfgs(model, errorfunction, input, target, FRET);
00123     return FRET;
00124 }
00125 
00126 double BFGS::optimize(
00127     Model& model,
00128     ErrorFunction& errorfunction,
00129     const Array<double>& input,
00130     const Array<double>& target,
00131     double gtol,
00132     int iterMax)
00133 {
00134     unsigned int iter;
00135     bfgs(model, errorfunction, input, target, gtol, iter, FRET, iterMax);
00136     return FRET;
00137 }
00138 
00139 
00140 void BFGS::bfgs(
00141     Model& model,
00142     ErrorFunction& errorfunction,
00143     const Array<double> &in,
00144     const Array<double> &target,
00145     double gtol,
00146     unsigned& iter,
00147     double&   fret,
00148     unsigned iterMax)
00149 {
00150     unsigned i, j, n = BFGS_n;
00151     DEDW.resize(n, false);
00152 //    part of Init
00153 
00154     if (firstIteration)
00155     {
00156         fret = errorfunction.errorDerivative(model, in, target, DEDW);
00157         BFGS_g0 = DEDW;
00158         firstIteration = false;
00159     }
00160     Array<double> p(n);
00161     for (i = 0;i < n;i++)p(i) = model.getParameter(i);
00162 
00163     double   d = 0., err, scale;
00164     double   test, temp;
00165 
00166     Array< double > p1(n);
00167     Array< double > g0(n);
00168     Array< double > g1(n);
00169     Array< double > s(n);
00170 
00171     Array< double > H(n, n);
00172 
00173     Array< double > gamma(n);
00174     Array< double > delta(n);
00175 
00176     iter = 0;
00177 
00178     double r = errorfunction.errorDerivative(model, in, target, DEDW);
00179     g0 = DEDW;
00180 
00181     for (err = 0., i = 0; i < n; ++i)
00182         err += g0(i) * g0(i);
00183     if (err < gtol*gtol)
00184     {
00185         fret = r;
00186         for (i = 0;i < n;i++)model.setParameter(i, p(i));
00187         return;
00188     }
00189 
00190     for (i = 0; i < n; ++i)
00191     {
00192         for (j = 0; j < n; ++j)
00193             H(i, j) = 0.;
00194         H(i, i) = 1.;
00195     }
00196 
00197     for (; ;)
00198     {
00199         for (i = 0; i < n; ++i)
00200             for (s(i) = 0., j = 0; j < n; ++j)
00201                 s(i) -= H(i, j) * g0(j);
00202 
00203         //
00204         // line search
00205         //
00206 
00207 
00208         switch (lineSearchType)
00209         {
00210         case Dlinmin:
00211 //  dlinmin( p1, s, fret, in, target);
00212             dlinmin(model, errorfunction, DEDW, p1, s, fret, in, target);
00213             break;
00214         case Linmin:
00215 //    linmin( p1, s, fret, in, target);
00216             linmin(model, errorfunction, DEDW, p1, s, fret, in, target);
00217             break;
00218         case Cblnsrch:
00219 //  cblnsrch(model, errorfunction, CG_fret, dedw, CG_xi, CG_fret, input, target);
00220             cblnsrch(model, errorfunction, fret, DEDW, BFGS_p1, fret, in, target);
00221 //  cblnsrch( BFGS_p, errorfuncion.errorDerrivative(input,target,DEDW), BFGS_g0, BFGS_s, BFGS_p1, fret, in, target);
00222             break;
00223         }
00224 
00225         for (i = 0;i < n;i++)model.setParameter(i, p1(i));
00226         errorfunction.errorDerivative(model, in, target, DEDW);
00227         g1 = DEDW;
00228         for (err = 0., i = 0; i < n; ++i)
00229             err += g1(i) * g1(i);
00230         if (++iter >= iterMax || err <= gtol*gtol)
00231             break;
00232 
00233         //
00234         // test for convergence
00235         //
00236         test = 0.;
00237         for (i = 0; i < n; i++)
00238         {
00239             temp = fabs(p1(i) - p(i)) / Shark::max((p(i)), 1.);
00240             if (temp > test)
00241                 test = temp;
00242         }
00243         if (test < gtol)
00244         {
00245             return;
00246         }
00247 
00248         for (d = 0., i = 0; i < n; ++i)
00249             d += (gamma(i) = g1(i) - g0(i)) * (delta(i) = p1(i) - p(i));
00250 
00251         for (i = 0; i < n; ++i)
00252             for (s(i) = 0., j = 0; j < n; ++j)
00253                 s(i) += H(i, j) * gamma(j);
00254 
00255         if (d < 1e-8)
00256         {
00257             for (i = 0; i < n; ++i)
00258             {
00259                 for (j = 0; j < n; ++j)
00260                     H(i, j) = 0.;
00261                 H(i, i) = 1.;
00262             }
00263         }
00264         else
00265         {
00266             for (scale = 0., i = 0; i < n; ++i)
00267                 scale += gamma(i) * s(i);
00268             scale = (scale / d + 1) / d;
00269 
00270             for (i = 0; i < n; ++i)
00271             {
00272                 g0(i) = g1(i);
00273                 p(i) = p1(i);
00274                 for (j = 0; j < n; ++j)
00275                     H(i, j) += scale * delta(i) * delta(j)
00276                                - (s(i) * delta(j) + s(j) * delta(i)) / d;
00277             }
00278         }
00279     }
00280 
00281     for (i = 0; i < n; ++i)model.setParameter(i, p1(i));
00282 
00283 }
00284 
00285 
00286 //===========================================================================
00321 unsigned i;
00322 void BFGS::init(Model& model)
00323 {
00324     initBfgs(model);
00325 }
00326 
00327 
00328 void BFGS::initBfgs(Model& model, LineSearch ls,
00329                     double ax, double bx, double lambda)
00330 {
00331     firstIteration = true;
00332     BFGS_n = model.getParameterDimension();
00333     unsigned i, j;
00334     BFGS_p.resize(BFGS_n, false);
00335     for (i = 0;i < BFGS_n;i++)BFGS_p(i) = model.getParameter(i);
00336     lineSearchType = ls;
00337     LS_lambda      = lambda;
00338     LS_ax          = ax;
00339     LS_bx          = bx;
00340 
00341 
00342 
00343 
00344     H.resize(BFGS_n, BFGS_n, false);
00345     BFGS_p1.resize(BFGS_n, false);
00346     BFGS_g0.resize(BFGS_n, false);
00347     BFGS_g1.resize(BFGS_n, false);
00348     BFGS_s.resize(BFGS_n, false);
00349 
00350     BFGS_gamma.resize(BFGS_n, false);
00351     BFGS_delta.resize(BFGS_n, false);
00352 
00353     for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p(i));
00354     for (i = 0; i < BFGS_n; ++i)
00355     {
00356         for (j = 0; j < BFGS_n; ++j) H(i, j) = 0.;
00357         H(i, i) = 1.;
00358     }
00359 }
00360 
00361 void BFGS::reset()
00362 {
00363     unsigned i, j;
00364     BFGS_g0 = DEDW;
00365     BFGS_p1 = 0.; BFGS_g1 = 0.; BFGS_s = 0.; BFGS_gamma = 0.; BFGS_delta = 0.;
00366     for (i = 0; i < BFGS_n; ++i)
00367     {
00368         for (j = 0; j < BFGS_n; ++j)
00369             H(i, j) = 0.;
00370         H(i, i) = 1.;
00371     }
00372 }
00373 
00374 
00375 //===========================================================================
00404 void BFGS::bfgs(Model& model,
00405                 ErrorFunction& errorfunction,
00406                 const Array<double> &in,
00407                 const Array<double> &target,
00408                 double& fret)
00409 {
00410     unsigned i, j;
00411     DEDW.resize(BFGS_n, false);
00412     //    part of Init
00413 
00414     if (firstIteration)
00415     {
00416         fret = errorfunction.errorDerivative(model, in, target, DEDW);
00417         BFGS_g0 = DEDW;
00418         firstIteration = false;
00419     }
00420     for (i = 0; i < BFGS_n; ++i)
00421     {
00422         for (BFGS_s(i) = 0., j = 0; j < BFGS_n; ++j)
00423             BFGS_s(i) -= H(i, j) * BFGS_g0(j);
00424     }
00425 
00426     //
00427     // line search
00428     //
00429 
00430     BFGS_p1 = BFGS_p;
00431 
00432     switch (lineSearchType)
00433     {
00434     case Dlinmin:
00435 //  dlinmin( p1, s, fret, in, target);
00436         dlinmin(model, errorfunction, DEDW, BFGS_p1, BFGS_s, fret, in, target);
00437         break;
00438     case Linmin:
00439 //    linmin( p1, s, fret, in, target);
00440         linmin(model, errorfunction, DEDW, BFGS_p1, BFGS_s, fret, in, target);
00441         break;
00442     case Cblnsrch:
00443 //  cblnsrch(model, errorfunction, CG_fret, dedw, CG_xi, CG_fret, input, target);
00444         cblnsrch(model, errorfunction, fret, DEDW, BFGS_p1, fret, in, target);
00445 //  cblnsrch( BFGS_p, errorfuncion.errorDerrivative(input,target,DEDW), BFGS_g0, BFGS_s, BFGS_p1, fret, in, target);
00446         break;
00447     }
00448 
00449     for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p1(i));
00450 
00451     errorfunction.errorDerivative(model, in, target, DEDW);
00452     BFGS_g1 = DEDW;
00453 
00454     for (BFGS_d = 0., i = 0; i < BFGS_n; ++i)
00455         BFGS_d += (BFGS_gamma(i) = BFGS_g1(i) - BFGS_g0(i)) * (BFGS_delta(i) = BFGS_p1(i) - BFGS_p(i));
00456 
00457     for (i = 0; i < BFGS_n; ++i)
00458         for (BFGS_s(i) = 0., j = 0; j < BFGS_n; ++j)
00459             BFGS_s(i) += H(i, j) * BFGS_gamma(j);
00460 
00461     if (BFGS_d < 1e-8)
00462     {
00463         for (i = 0; i < BFGS_n; ++i)
00464         {
00465             for (j = 0; j < BFGS_n; ++j)
00466                 H(i, j) = 0.;
00467             H(i, i) = 1.;
00468         }
00469     }
00470     else
00471     {
00472         for (BFGS_scale = 0., i = 0; i < BFGS_n; ++i)
00473             BFGS_scale += BFGS_gamma(i) * BFGS_s(i);
00474         BFGS_scale = (BFGS_scale / BFGS_d + 1) / BFGS_d;
00475 
00476         for (i = 0; i < BFGS_n; ++i)
00477         {
00478             BFGS_g0(i) = BFGS_g1(i);
00479             BFGS_p(i) = BFGS_p1(i);
00480             for (j = 0; j < BFGS_n; ++j)
00481                 H(i, j) += BFGS_scale * BFGS_delta(i) * BFGS_delta(j)
00482                            - (BFGS_s(i) * BFGS_delta(j) + BFGS_s(j) * BFGS_delta(i)) / BFGS_d;
00483         }
00484     }
00485 
00486     for (i = 0; i < BFGS_n; ++i)
00487         BFGS_p(i) = BFGS_p1(i);
00488 
00489     for (i = 0;i < BFGS_n;i++)model.setParameter(i, BFGS_p(i));
00490 }
00491 
00492 
00494 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
00495 
00496 
00497 
00498 void BFGS::dlinmin(
00499     Model& model,
00500     ErrorFunction& errorfunction,
00501     Array<double>& dedw,
00502     Array< double >& p,
00503     const Array< double >& xi,
00504     double& fret,
00505     const Array< double >&in,
00506     const Array< double >&target)
00507 {
00508 
00509     SIZE_CHECK(p.ndim() == 1 && p.samedim(xi))
00510     const double   GOLD   = 1.618034;
00511     const double   GLIMIT = 100.;
00512     const double   TINY   = 1.0e-20;
00513     const unsigned ITMAX  = 100;
00514     const double   ZEPS   = 1.0e-10;
00515     const double   TOL    = 2.0e-4;
00516     bool     ok1, ok2;
00517     unsigned i, iter;
00518     double   fa, fb, fc, fp, ax, bx, cx;
00519     double   ulim, dum;
00520     double   a, b, d = 0, e, fu, fv, fw, fx, q, r, tol1, tol2, u, v, _w, x, xm;
00521     double   du, dv, dw, dx, d1, d2, u1, u2, olde;
00522 
00523 
00524 
00525     unsigned        n = BFGS_n;
00526     Array< double > xt(n);
00527     Array< double > generalDerivative(n);
00528 
00529     ax = LS_ax;
00530     bx = LS_bx;
00531 
00532     //===================================================================
00533 
00534     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*ax);
00535     fp = fa = errorfunction.error(model, in, target);
00536 
00537     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*bx);
00538     fb = errorfunction.error(model, in, target);
00539 
00540     if (fb > fa)
00541     {
00542 
00543         dum = ax;
00544         ax = bx;
00545         bx = dum;
00546         dum = fb;
00547         fb = fa;
00548         fa = dum;
00549     }
00550 
00551     cx = bx + GOLD * (bx - ax);
00552 
00553     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*cx);
00554     fc = errorfunction.error(model, in, target);
00555 
00556     while (fb > fc)
00557     {
00558         r = (bx - ax) * (fb - fc);
00559         q = (bx - cx) * (fb - fa);
00560         u = bx - ((bx - cx) * q - (bx - ax) * r) /
00561             (2. * SIGN(Shark::max(fabs(q - r), TINY), q - r));
00562         ulim = bx + GLIMIT * (cx - bx);
00563         if ((bx - u) *(u - cx) > 0.)
00564         {
00565             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00566 
00567             fu = errorfunction.error(model, in, target);
00568 
00569             if (fu < fc)
00570             {
00571                 ax = bx;
00572                 bx = u;
00573                 fa = fb;
00574                 fb = fu;
00575                 break;
00576             }
00577             else if (fu > fb)
00578             {
00579                 cx = u;
00580                 fc = fu;
00581                 break;
00582             }
00583             u = cx + GOLD * (cx - bx);
00584 
00585             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00586             fu = errorfunction.error(model, in, target);
00587 
00588         }
00589         else if ((cx - u) *(u - ulim) > 0.)
00590         {
00591             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00592             fu = errorfunction.error(model, in, target);
00593 
00594             if (fu < fc)
00595             {
00596                 bx = cx;
00597                 cx = u;
00598                 u  = cx + GOLD * (cx - bx);
00599                 fb = fc;
00600                 fc = fu;
00601                 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00602 
00603                 fu = errorfunction.error(model, in, target);
00604             }
00605         }
00606         else if ((u - ulim) *(ulim - cx) >= 0.)
00607         {
00608             u = ulim;
00609             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00610             fu = errorfunction.error(model, in, target);
00611         }
00612         else
00613         {
00614             u = cx + GOLD * (cx - bx);
00615             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00616             fu = errorfunction.error(model, in, target);
00617         }
00618         ax = bx;
00619         bx = cx;
00620         cx = u;
00621         fa = fb;
00622         fb = fc;
00623         fc = fu;
00624     }
00625     //=======================================================================
00626 
00627     e = 0.;
00628     if (ax < cx)
00629     {
00630         a = ax;
00631         b = cx;
00632     }
00633     else
00634     {
00635         a = cx;
00636         b = ax;
00637     }
00638     x = _w = v = bx;
00639     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*x);
00640     fw = fv = fx = errorfunction.errorDerivative(model, in, target, dedw);
00641     generalDerivative = dedw;
00642 
00643     for (dx = 0., i = 0; i < n; ++i)
00644         dx += xi(i) * generalDerivative(i);
00645     dw = dv = dx;
00646 
00647     for (iter = 0; iter < ITMAX; iter++)
00648     {
00649         xm = 0.5 * (a + b);
00650         tol2 = 2. * (tol1 = TOL * fabs(x) + ZEPS);
00651         if (fabs(x - xm) <= (tol2 - 0.5 *(b - a)))
00652         {
00653             break;
00654         }
00655         if (fabs(e) > tol1)
00656         {
00657             d1 = 2. * (b - a);
00658             d2 = d1;
00659             if (dw != dx)
00660                 d1 = (_w - x) * dx / (dx - dw);
00661             if (dv != dx)
00662                 d2 = (v - x) * dx / (dx - dv);
00663             u1   = x + d1;
00664             u2   = x + d2;
00665             ok1  = (a - u1) * (u1 - b) > 0. && dx * d1 <= 0.;
00666             ok2  = (a - u2) * (u2 - b) > 0. && dx * d2 <= 0.;
00667             olde = e;
00668             e    = d;
00669             if (ok1 || ok2)
00670             {
00671                 if (ok1 && ok2)
00672                     d = (fabs(d1) < fabs(d2) ? d1 : d2);
00673                 else if (ok1)
00674                     d = d1;
00675                 else
00676                     d = d2;
00677                 if (fabs(d) <= fabs(0.5 * olde))
00678                 {
00679                     u = x + d;
00680                     if (u - a < tol2 || b - u < tol2)
00681                         d = SIGN(tol1, xm - x);
00682                 }
00683                 else
00684                 {
00685                     d = 0.5 * (e = (dx >= 0. ? a - x : b - x));
00686                 }
00687             }
00688             else
00689             {
00690                 d = 0.5 * (e = (dx >= 0. ? a - x : b - x));
00691             }
00692         }
00693         else
00694         {
00695             d = 0.5 * (e = (dx >= 0. ? a - x : b - x));
00696         }
00697         if (fabs(d) >= tol1)
00698         {
00699             u = x + d;
00700             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00701             fu = errorfunction.errorDerivative(model, in, target, dedw);
00702             generalDerivative = dedw;
00703         }
00704         else
00705         {
00706             u = x + SIGN(tol1, d);
00707             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00708             fu = errorfunction.errorDerivative(model, in, target, dedw);
00709             generalDerivative = dedw;
00710             if (fu > fx)
00711                 break;
00712         }
00713         for (du = 0., i = 0; i < n; ++i)
00714             du += xi(i) * generalDerivative(i);
00715         if (fu <= fx)
00716         {
00717             if (u >= x)
00718                 a = x;
00719             else
00720                 b = x;
00721             v  = _w;
00722             _w  = x;
00723             x  = u;
00724             fv = fw;
00725             fw = fx;
00726             fx = fu;
00727             dv = dw;
00728             dw = dx;
00729             dx = du;
00730         }
00731         else
00732         {
00733             if (u < x)
00734                 a = u;
00735             else
00736                 b = u;
00737             if (fu <= fw || _w == x)
00738             {
00739                 v  = _w;
00740                 _w  = u;
00741                 fv = fw;
00742                 fw = fu;
00743                 dv = dw;
00744                 dw = du;
00745             }
00746             else if (fu < fv || v == x || v == _w)
00747             {
00748                 v  = u;
00749                 fv = fu;
00750                 dv = du;
00751             }
00752         }
00753     }
00754     if (iter >= ITMAX)
00755         throw SHARKEXCEPTION("BFGS: too many iterations in routine dbrent");
00756     if (fx < fp)
00757     {
00758         fret = fx;
00759         for (i = 0; i < n; ++i)
00760             p(i) += xi(i) * x;
00761     }
00762     else
00763         fret = fp;
00764 }
00765 
00766 // linmin.cpp
00767 
00768 void BFGS::linmin(
00769     Model& model,
00770     ErrorFunction& errorfunction,
00771     Array<double>& dedw,
00772     Array< double >& p,
00773     const Array< double >& xi,
00774     double& fret,
00775     const Array< double >&in,
00776     const Array< double >&target)
00777 {
00778 
00779     SIZE_CHECK(p.ndim() == 1 && p.samedim(xi))
00780     const double   GOLD   = 1.618034;
00781     const double   GLIMIT = 100.;
00782     const double   TINY   = 1.0e-20;
00783     const unsigned ITMAX  = 100;
00784     const double   CGOLD  = 0.3819660;
00785     const double   ZEPS   = 1.0e-10;
00786     const double   TOL    = 2.0e-4;
00787 
00788     unsigned i, iter;
00789     double   fa, fb, fc, ax, bx, cx;
00790     double   ulim, dum;
00791     double   a, b, d = 0., e, etemp, fu, fv, fw, fx, o, q, r, tol1, tol2, u, v, _w, x, xm;
00792     unsigned        n = BFGS_n;
00793     Array< double > xt(n);
00794 
00795     ax = LS_ax;
00796     bx = LS_bx;
00797 
00798     //===================================================================
00799 
00800     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*ax);
00801     double initial = fa = errorfunction.error(model, in, target);
00802 
00803     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*bx);
00804     fb = errorfunction.error(model, in, target);
00805 
00806     if (fb > fa)
00807     {
00808         dum = ax;
00809         ax  = bx;
00810         bx  = dum;
00811         dum = fb;
00812         fb  = fa;
00813         fa  = dum;
00814     }
00815     cx = bx + GOLD * (bx - ax);
00816     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*cx);
00817     fc = errorfunction.error(model, in, target);
00818     while (fb > fc)
00819     {
00820         r = (bx - ax) * (fb - fc);
00821         q = (bx - cx) * (fb - fa);
00822         u = bx - ((bx - cx) * q - (bx - ax) * r) /
00823             (2. * SIGN(Shark::max(fabs(q - r), TINY), q - r));
00824         ulim = bx + GLIMIT * (cx - bx);
00825         if ((bx - u) *(u - cx) > 0.)
00826         {
00827             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00828             fu = errorfunction.error(model, in, target);
00829             if (fu < fc)
00830             {
00831                 ax = bx;
00832                 bx = u;
00833                 fa = fb;
00834                 fb = fu;
00835                 break;
00836             }
00837             else if (fu > fb)
00838             {
00839                 cx = u;
00840                 fc = fu;
00841                 break;
00842             }
00843             u = cx + GOLD * (cx - bx);
00844             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00845             fu = errorfunction.error(model, in, target);
00846         }
00847         else if ((cx - u) *(u - ulim) > 0.)
00848         {
00849             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00850             fu = errorfunction.error(model, in, target);
00851             if (fu < fc)
00852             {
00853                 bx = cx;
00854                 cx = u;
00855                 u  = cx + GOLD * (cx - bx);
00856                 fb = fc;
00857                 fc = fu;
00858                 for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00859                 fu = errorfunction.error(model, in, target);
00860             }
00861         }
00862         else if ((u - ulim) *(ulim - cx) >= 0.)
00863         {
00864             u = ulim;
00865             for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00866             fu = errorfunction.error(model, in, target);
00867         }
00868         else
00869         {
00870             u = cx + GOLD * (cx - bx);
00871             for (i = 0; i < n; ++i)  model.setParameter(i, p(i) + xi(i)*u);
00872             fu = errorfunction.error(model, in, target);
00873         }
00874         ax = bx;
00875         bx = cx;
00876         cx = u;
00877         fa = fb;
00878         fb = fc;
00879         fc = fu;
00880     }
00881     //=======================================================================
00882 
00883     e = 0.;
00884     if (ax < cx)
00885     {
00886         a = ax;
00887         b = cx;
00888     }
00889     else
00890     {
00891         a = cx;
00892         b = ax;
00893     }
00894     x = _w = v = bx;
00895     for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*x);
00896     fw = fv = fx = errorfunction.error(model, in, target);
00897 
00898     for (iter = 0; iter < ITMAX; iter++)
00899     {
00900         xm = 0.5 * (a + b);
00901         tol2 = 2. * (tol1 = TOL * fabs(x) + ZEPS);
00902         if (fabs(x - xm) <= (tol2 - 0.5 *(b - a)))
00903         {
00904             break;
00905         }
00906         if (fabs(e) > tol1)
00907         {
00908             r = (x - _w) * (fx - fv);
00909             q = (x - v) * (fx - fw);
00910             o = (x - v) * q - (x - _w) * r;
00911             q = 2. * (q - r);
00912             if (q > 0.)
00913                 o = -o;
00914             q = fabs(q);
00915             etemp = e;
00916             e = d;
00917             if (fabs(o) >= fabs(0.5 * q * etemp) ||
00918                     o <= q *(a - x) ||
00919                     o >= q *(b - x))
00920                 d = CGOLD * (e = (x >= xm ? a - x : b - x));
00921             else
00922             {
00923                 d = o / q;
00924                 u = x + d;
00925                 if (u - a < tol2 || b - u < tol2)
00926                     d = SIGN(tol1, xm - x);
00927             }
00928         }
00929         else
00930             d = CGOLD * (e = (x >= xm ? a - x : b - x));
00931         u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
00932         for (i = 0; i < n; ++i) model.setParameter(i, p(i) + xi(i)*u);
00933         fu = errorfunction.error(model, in, target);
00934         if (fu <= fx)
00935         {
00936             if (u >= x)
00937                 a = x;
00938             else
00939                 b = x;
00940             v = _w;
00941             _w = x;
00942             x = u;
00943             fv = fw;
00944             fw = fx;
00945             fx = fu;
00946         }
00947         else
00948         {
00949             if (u < x)
00950                 a = u;
00951             else
00952                 b = u;
00953             if (fu <= fw || _w == x)
00954             {
00955                 v = _w;
00956                 _w = u;
00957                 fv = fw;
00958                 fw = fu;
00959             }
00960             else if (fu <= fv || v == x || v == _w)
00961             {
00962                 v = u;
00963                 fv = fu;
00964             }
00965         }
00966     }
00967 
00968     if (iter >= ITMAX)
00969         throw SHARKEXCEPTION("BFGS: too many iterations in brent");
00970     fret = fx;
00971     //=======================================================================
00972 
00973     if (fx < initial)
00974     {
00975         fret = fx;
00976         for (i = 0; i < n; ++i)
00977             p(i) += xi(i) * x;
00978     }
00979     else
00980     {
00981         fret = initial;
00982     }
00983 }
00984 //======================================================================
00985 // Definition of the cubic line search class
00986 // Armijo and Goldstein's line search algorithm
00987 
00988 // author:  Doug Hart, Adapted to the COOOL library by Wenceslau Gouveia
00989 
00990 // Modified to fit into new classes.  H. Lydia Deng, 02/21/94, 03/15/94
00991 
00992 //========================================================================
00993 
00994 
00995 /*
00996 
00997  * cubic line search
00998 
00999  *
01000 
01001  * The parameter lambda controls the accuraccy of the line search.
01002 
01003  * lambda = .25 is a good choice.
01004 
01005  */
01006 
01007 void BFGS::cblnsrch(
01008     Model& model,
01009     ErrorFunction& errorfunction,
01010     double fold,
01011     Array< double >& g,
01012     Array< double >& p,
01013     double& f,
01014     const Array< double >&in,
01015     const Array< double >&target)
01016 {
01017     //
01018     // the maximum of performed iterations should be choosen small enough
01019     // to ensure that neither (3^iterMax) nor (0.1^(2*iterMax)) causes a
01020     // numerical overflow or underflow.
01021     //
01022     const unsigned iterMax = 20;
01023     unsigned iterNum;
01024     unsigned i, n;
01025     bool     tst = false;
01026     double   slope;
01027     double   alpha, alpha2;
01028     double   alpha_prev, alpha_prev2;
01029     double   alpha_tmp = 0;
01030     double   f1, f2, fprev = 0;
01031     double   a, b, c;
01032     double   cm11, cm12, cm21, cm22;
01033     double   disc;
01034     double   lambda = LS_lambda;
01035 
01036     n = BFGS_n;
01037     /*
01038 
01039      * dot product of search direction and gradient
01040 
01041      */
01042 
01043     for (slope = 0., i = 0; i < n; i++)
01044         slope += g(i) * p(i);
01045 
01046     iterNum = 0;            /* iteration counter */
01047     alpha   = 1.;           /* updating step */
01048 
01049     /*
01050      * updating
01051      */
01052 
01053     for (i = 0; i < n; i++)
01054         model.setParameter(i, model.getParameter(i) + alpha * p(i));
01055     f = errorfunction.error(model, in, target);
01056     iterNum++;
01057     /*
01058      * Implementing Goldstein's test for alpha too small
01059      */
01060     while (f < fold + (1. - lambda)*alpha*slope && iterNum < iterMax)
01061     {
01062         alpha *= 3;
01063         for (i = 0; i < n; i++)
01064             model.setParameter(i, model.getParameter(i) + alpha * p(i));
01065         f = errorfunction.error(model, in, target);
01066         iterNum++;
01067     }
01068     if (iterNum >= iterMax)
01069         throw SHARKEXCEPTION("BFGS: alpha overflowed!");
01070 
01071     /*
01072      * Armijo's test for alpha too large
01073      */
01074 
01075     alpha_prev = alpha; /* H.L. Deng, 6/13/95 */
01076     while (f > fold + lambda*alpha*slope && iterNum < iterMax)
01077     {
01078         alpha2 = alpha * alpha;
01079         f1 = f - fold - slope * alpha;
01080         if (tst == false)
01081         {
01082             alpha_tmp = -slope * alpha2 / (f1 * 2.); /* tentative alpha */
01083             tst = true;
01084         }
01085         else
01086         {
01087             alpha_prev2 = alpha_prev * alpha_prev;
01088             f2   = fprev - fold - alpha_prev * slope;
01089             c    = 1. / (alpha - alpha_prev);
01090             cm11 = 1. / alpha2;
01091             cm12 = -1. / alpha_prev2;
01092             cm21 = -alpha_prev / alpha2;
01093             cm22 = alpha / alpha_prev2;
01094             a    = c * (cm11 * f1 + cm12 * f2);
01095             b    = c * (cm21 * f1 + cm22 * f2);
01096             disc = b * b - 3. * a * slope;
01097 
01098             if ((fabs(a) > FLT_MIN) && (disc > FLT_MIN))
01099                 alpha_tmp = (-b + sqrt(disc)) / (3. * a);
01100             else
01101                 alpha_tmp = slope * alpha2 / (2. * f1);
01102             if (alpha_tmp >= .5 * alpha)
01103                 alpha_tmp = .5 * alpha;
01104         }
01105         alpha_prev = alpha;
01106         fprev = f;
01107         if (alpha_tmp < .1 * alpha)
01108             alpha *= .1;
01109         else
01110             alpha = alpha_tmp;
01111         for (i = 0; i < n; i++)
01112             model.setParameter(i, model.getParameter(i) + alpha * p(i));
01113         f = errorfunction.error(model, in, target);
01114         iterNum++;
01115     }
01116     if (iterNum >= iterMax)
01117         throw SHARKEXCEPTION("BFGS: alpha underflowed!");
01118 }
01119