ADMB Documentation  11.1.2274
 All Classes Files Functions Variables Typedefs Friends Defines
fmmtr1.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fmmtr1.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: Unknown
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  *
00007  * This file was originally written in FORTRAN II by and unknown author.
00008  * In the 1980s, it was ported to C and C++ and extensively modified by
00009  * David Fournier.
00010  *
00011  */
00016 #ifdef __ZTC__
00017   #include <conio.h>
00018 #endif
00019 #include <fvar.hpp>
00020 extern int ctlc_flag;
00021 #ifdef __TURBOC__
00022   #pragma hdrstop
00023   #include <iostream.h>
00024   #include <conio.h>
00025 #endif
00026 #if defined (__WAT32__) || defined(_MSC_VER)
00027   #include <conio.h>
00028 #endif
00029 #ifdef __ZTC__
00030   #include <iostream.hpp>
00031   #include <disp.h>
00032   #define endl "\n"
00033 #endif
00034 #ifdef __SUN__
00035   #include <iostream.h>
00036   #include <signal.h>
00037   #define getch getchar
00038   #ifdef __HP__
00039   extern "C" void onintr(int k);
00040   #endif
00041 #endif
00042 #ifdef __NDPX__
00043   #include <iostream.hxx>
00044 #endif
00045 #if defined (_MSC_VER)
00046   void __cdecl clrscr();
00047 #else
00048   #include <iostream>
00049   #include <signal.h>
00050   #define getch getchar
00051   extern "C" void onintr(int k);
00052   extern "C" void clrscr();
00053 #endif
00054 #include <math.h>
00055 #include <stdlib.h>
00056 #include <stdio.h>
00057 #include <ctype.h>
00058 dvector update(int nvar, int iter, int m, const dvector& g,
00059   const dmatrix& xalpha, dmatrix& y, const dvector& x, const dvector& xold,
00060   const dvector& gold, const dvector& xrho);
00061 double dafsqrt( double x );
00062 
00067 void fmmt1::fmin(const double& _f, const dvector & _x, const dvector& _g)
00068 {
00069   double& f=(double&) _f;
00070   independent_variables & x=(independent_variables &) _x;
00071   dvector& g=(dvector&) _g;
00072     #ifdef DIAG
00073       cout << "On entry to fmin: " << *this << endl;
00074     #endif
00075   if (use_control_c)
00076   {
00077 #if !defined (_MSC_VER)
00078     #if defined( __SUN__) && !(defined __GNU__)
00079       #if defined( __HP__)
00080         if (ireturn <= 0 )
00081         {
00082           signal(SIGINT, &onintr);
00083         }
00084       #else
00085         if (ireturn <= 0 )
00086         {
00087           signal(SIGINT, (SIG_PF)&onintr);
00088         }
00089       #endif
00090     #endif
00091  #endif
00092   }
00093   if (use_control_c)
00094   {
00095     #if defined( __GNU__)
00096       if (ireturn <= 0 )
00097       {
00098         signal(SIGINT, &onintr);
00099       }
00100     #endif
00101   }
00102     #ifdef __ZTC__
00103       if (ireturn <= 0 )
00104       {
00105         if (disp_inited == 0)
00106         {
00107           disp_open();
00108           disp_usebios();
00109         }
00110       }
00111     #endif
00112       if (ireturn >= 3)
00113       {
00114          derch(f, x, g, n, ireturn);
00115          return;
00116       }
00117       if (ireturn == 1) goto call1;
00118       if (ireturn == 2) goto call2;
00119       ihflag=0;
00120      if (n==0)
00121      {
00122        cerr << "Error -- the number of active parameters"
00123          " fmin must be > 0\n";
00124        ad_exit(1);
00125      }
00126      if (x.indexmin() !=1)
00127      {
00128        cerr << "Error -- minimum valid index"
00129          " for independent_variables in fmin must be 1\n"
00130         << " it is " << x.indexmin() << "\n";
00131         ad_exit(1);
00132      }
00133      if (x.size() <n)
00134      {
00135        cerr << "Error -- the size of the independent_variables"
00136         " which is " << x.size() << " must be >= " << n << "\n"
00137         << " the number of independent variables in fmin\n";
00138         ad_exit(1);
00139      }
00140      if (g.indexmin() !=1)
00141      {
00142        cerr << "Error -- minimum valid index"
00143          " for the gradient vector in fmin must be 1\n"
00144         << " it is " << g.indexmin() << "\n";
00145         ad_exit(1);
00146      }
00147      if (g.size() <n)
00148      {
00149        cerr << "Error -- the size of the gradient vector"
00150         " which is " << g.size() << " must be >=\n"
00151         << " the number of independent variables in fmin\n";
00152         ad_exit(1);
00153      }
00154      for (i=1; i<=n; i++)
00155            xx.elem(i)=x.elem(i);
00156       itn=0;
00157       icc=0;
00158        for (i=1; i< 11; i++)
00159           funval.elem(i)=0.;
00160       ihang = 0;
00161       llog=1;
00162       np=n+1;
00163       n1=n-1;
00164       nn=n*np/2;
00165       is=n;
00166       iu=n;
00167       iv=n+n;
00168       ib=iv+n;
00169       iexit=0;
00170       dmin=1;
00171       if (dmin <= 0.)
00172          goto label7020;
00173       if(dfn == 0.)
00174          z = 0.0;
00175       for (i=1; i<=n; i++)
00176       {
00177         xsave.elem(i)=x.elem(i);
00178         x.elem(i)=xx.elem(i);
00179       }
00180       ireturn=1;
00181       return;
00182   call1:
00183       for (i=1; i<=n; i++)
00184       {
00185         x.elem(i)=xsave.elem(i);
00186       }
00187       ireturn=3;
00188       fbest = f;
00189       for ( i=1; i<=n; i++)
00190          gbest.elem(i)=g.elem(i);
00191       funval.elem(10) = f;
00192       df=dfn;
00193       if(dfn == 0.0)
00194          df = f - z;
00195       if(dfn < 0.0)
00196          df=fabs(df * f);
00197       if(df <= 0.0)
00198          df=1.;
00199 label20:
00200       ic=0;
00201       iconv = 1;
00202       for ( i=1; i<=9; i++)
00203          funval.elem(i)= funval.elem(i+1);
00204       funval.elem(10) = f;
00205       if ( itn>15 && fabs( funval.elem(1)-funval.elem(10))< min_improve )
00206          ihang = 1;
00207       gmax = 0;
00208       for ( i=1; i<=n; i++)
00209       {
00210         if(fabs(g.elem(i)) > crit) iconv = 2;
00211         if(fabs(g.elem(i)) > fabs(gmax) ) gmax = g.elem(i);
00212       }
00213       if( iconv == 1 || ihang == 1)
00214       {
00215          ireturn=-1;
00216          goto label92;
00217       }
00218       if(iprint == 0)
00219          goto label21;
00220       if(itn == 0)
00221          goto label7000;
00222       if( (itn%iprint) != 0)
00223          goto label21;
00224       if (llog) goto label7010;
00225 #if !defined (_MSC_VER)  && !defined (__GNUC__)
00226         if (!scroll_flag) clrscr();
00227 #endif
00228 label7003:
00229       if (iprint!=0)
00230       {
00231         if (ad_printf)
00232           (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00233               n, itn, ifn);
00234         if (ad_printf)
00235 (*ad_printf)("Function value %12.4le; maximum gradient component mag %12.4le\n",
00236             f, gmax);
00237       }
00238 /*label7002:*/
00239       /* Printing Statistics table */
00240       if(iprint>0)
00241       {
00242         fmmdisp(x, g, n, this->scroll_flag,noprintx);
00243       }
00244 label21 :
00245       itn++;
00246       rrr=update(n,itn-1,xm,g,xstep,xy,x,xold,gold,xrho);
00247       for (i=1; i<=n; i++)
00248          x.elem(i)=xx.elem(i);
00249       for (i=1; i<=n; i++)
00250       {
00251         w.elem(i)=-g.elem(i);
00252       }
00253       for (i=1; i<=n; i++)
00254       {
00255         w(n+i)=rrr(i);
00256       }
00257       gs=0.0;
00258       for (i=1; i<=n; i++)
00259          gs+=w.elem(is+i)*g.elem(i);
00260       iexit=2;
00261       if(gs >= 0.0)
00262       {
00263         double a=1.1*gs/norm2(g);
00264         for (i=1; i<=n; i++)
00265            w.elem(is+i)-=a*g.elem(i);
00266         gs=0.0;
00267         for (i=1; i<=n; i++)
00268           gs+=w.elem(is+i)*g.elem(i);
00269         if(gs >= 0.0)
00270         {
00271           cout << "Kludge didn't work" << endl;
00272           goto label92;
00273         }
00274       }
00275       gso=gs;
00276       alpha=-2.0*df/gs;
00277       if(alpha > 1.0)
00278         alpha=1.0;
00279       df=f;
00280       tot=0.0;
00281 label30:
00282       iexit=3;
00283       if( ifn >= maxfn)
00284       {
00285          maxfn_flag=1;
00286          goto label92;
00287       }
00288       else
00289       {
00290          maxfn_flag=0;
00291          iexit=1;
00292       }
00293       if(quit_flag && use_control_c) goto label92;
00294       for (i=1; i<=n; i++)
00295       {
00296         z=alpha*w.elem(is+i);
00297         xx.elem(i)+=z;
00298       }
00299       for (i=1; i<=n; i++)
00300       {
00301         xsave.elem(i)=x.elem(i);
00302         gsave.elem(i)=g.elem(i);
00303         x.elem(i)=xx.elem(i);
00304         fsave = f;
00305       }
00306       fsave = f;
00307       ireturn=2;
00308       return;
00309   call2:
00310       for (i=1; i<=n; i++)
00311       {
00312         x.elem(i)=xsave.elem(i);
00313         w.elem(i)=g.elem(i);
00314         g.elem(i)=gsave.elem(i);
00315       }
00316       fy = f;
00317       f = fsave;
00318       ireturn=-1;
00319       if (fy < fbest)
00320       {
00321         fbest=fy;
00322         for (i=1; i<=n; i++)
00323         {
00324           x.elem(i)=xx.elem(i);
00325           gbest.elem(i)=w.elem(i);
00326         }
00327       }
00328 #if defined(_MSC_VER)
00329        if (kbhit())
00330 #else
00331        if(ctlc_flag && use_control_c)
00332 #endif
00333        {
00334           #if defined(__DJGPP__)
00335             int c = toupper(getxkey());
00336           #else
00337             int c = toupper(getch());
00338           #endif
00339           if ( c == 'C')
00340           {
00341             for (i=1; i<=n; i++)
00342             {
00343               x.elem(i)=xx.elem(i);
00344             }
00345             ireturn = 3;
00346             derch(f, x , w, n, ireturn);
00347             return;
00348           }
00349           else
00350           {
00351             if ( c == 'Q'|| c == 'N')
00352             {
00353               quit_flag=c;
00354               goto label92;
00355             }
00356             else
00357             {
00358               quit_flag=0;
00359             }
00360           }
00361        }
00362        icc+=1;
00363        if( icc >= 5)
00364           icc=0;
00365       ic++;
00366       if( ic >imax)
00367       {
00368          if (iprint>0)
00369          {
00370            if (ad_printf)
00371              (*ad_printf)("  ic > imax  in fminim is answer attained ?\n" );
00372            fmmdisp(x, g, n, this->scroll_flag,noprintx);
00373          }
00374          ihflag=1;
00375          ihang=1;
00376          goto label92;
00377       }
00378       ifn++;
00379       gys=0.0;
00380       for (i=1; i<= n; i++)
00381          gys+=w.elem(i)*w.elem(is+i);
00382       if(fy>=f)
00383          goto label40;
00384       if(fabs(gys/gso)<=.9)
00385          goto label50;
00386       if(fabs(gys/gso)<=.95 && ic > 4)
00387          goto label50;
00388       if(gys>0.0)
00389          goto  label40;
00390       tot+=alpha;
00391       z=10.0;
00392       if(gs<gys)
00393          z=gys/(gs-gys);
00394       if(z>10.0)
00395          z=10.0;
00396       alpha=alpha*z;
00397       if (alpha == 0.)
00398       {
00399          ialph=1;
00400          #ifdef __ZTC__
00401          if (ireturn <= 0)
00402          {
00403            disp_close();
00404          }
00405          #endif
00406          return;
00407       }
00408       f=fy;
00409       gs=gys;
00410       goto label30;
00411 label40:
00412       for (i=1;i<=n;i++)
00413          xx.elem(i)-=alpha*w.elem(is+i);
00414       z=3.0*(f-fy)/alpha+gys+gs;
00415       zz=dafsqrt(z*z-gs*gys);
00416       z=1.0-(gys+zz-z)/(2.0*zz+gys-gs);
00417       if (fabs(fy-1.e+95) < 1.e-66)
00418       {
00419         alpha*=.001;
00420       }
00421       else
00422       {
00423         alpha=alpha*z;
00424       }
00425       if (alpha == 0.)
00426       {
00427          ialph=1;
00428          #ifdef __ZTC__
00429          if (ireturn <= 0)
00430          {
00431            disp_close();
00432          }
00433          #endif
00434          return;
00435       }
00436       goto label30;
00437 label50:
00438       alpha+=tot;
00439       f=fy;
00440       df-=f;
00441       dgs=gys-gso;
00442       link=1;
00443       if(dgs+alpha*gso>0.0)
00444          goto label52;
00445       for (i=1;i<=n;i++)
00446          w.elem(iu+i)=w.elem(i)-g.elem(i);
00447       sig=1.0/(alpha*dgs);
00448       goto label70;
00449 label52:
00450       zz=alpha/(dgs-alpha*gso);
00451       z=dgs*zz-1.0;
00452       for (i=1;i<=n;i++)
00453          w.elem(iu+i)=z*g.elem(i)+w.elem(i);
00454       sig=1.0/(zz*dgs*dgs);
00455       goto label70;
00456 label60:
00457       link=2;
00458       for (i=1;i<=n;i++)
00459          w.elem(iu+i)=g.elem(i);
00460       if(dgs+alpha*gso>0.0)
00461          goto label62;
00462       sig=1.0/gso;
00463       goto  label70;
00464 label62:
00465       sig=-zz;
00466       goto label70;
00467 label65:
00468       for (i=1;i<=n;i++)
00469          g.elem(i)=w.elem(i);
00470       goto  label20;
00471 label70:
00472       if (link == 1) goto label60;
00473       if (link == 2) goto label65;
00474 /*label90:*/
00475       for (i=1;i<=n;i++)
00476          g.elem(i)=w.elem(i);
00477 label92:
00478       if (iprint>0)
00479       {
00480         if (ihang == 1)
00481           if (ad_printf)
00482             (*ad_printf)(
00483            "Function minimizer not making progress ... is minimum attained?\n");
00484       }
00485       if(iexit == 2)
00486       {
00487         if (iprint>0)
00488         {
00489           if (ad_printf)
00490             (*ad_printf)("*** grad transpose times delta x greater >= 0\n"
00491            " --- convergence critera may be too strict\n");
00492           ireturn=-1;
00493         }
00494       }
00495 #     if defined (_MSC_VER)  && !defined (__WAT32__)
00496         if (scroll_flag == 0) clrscr();
00497 #     endif
00498       if (maxfn_flag == 1)
00499       {
00500         if (iprint>0)
00501         {
00502           if (ad_printf)
00503             (*ad_printf)("Maximum number of function evaluations exceeded");
00504         }
00505       }
00506       if (iprint>0)
00507       {
00508         if (quit_flag == 'Q' && use_control_c)
00509           if (ad_printf) (*ad_printf)("User initiated interrupt");
00510       }
00511       if(iprint == 0) goto label777;
00512       if (ad_printf) (*ad_printf)(" - final statistics:\n");
00513       if (ad_printf)
00514         (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00515           n, itn, ifn);
00516       if (ad_printf)
00517         (*ad_printf)(
00518         "Function value %12.4le; maximum gradient component mag %12.4le\n",
00519         f, gmax);
00520       if (ad_printf)
00521         (*ad_printf)("Exit code = %ld;  converg criter %12.4le\n",iexit,crit);
00522 
00523       fmmdisp(x, g, n, this->scroll_flag,noprintx);
00524 label777:
00525          if (ireturn <= 0)
00526          #ifdef DIAG
00527            if (ad_printf) (*ad_printf)("Final values of h in fmin:\n");
00528            cout << h << "\n";
00529          #endif
00530          #ifdef __ZTC__
00531          {
00532            disp_close();
00533          }
00534          #endif
00535          return;
00536 label7000:
00537       if (iprint>0)
00538       {
00539 #     if defined (_MSC_VER)  && !defined (__WAT32__)
00540         if (!scroll_flag) clrscr();
00541 #endif
00542         if (ad_printf) (*ad_printf)("Initial statistics: ");
00543       }
00544       goto label7003;
00545 label7010:
00546    if (iprint>0)
00547    {
00548 #     if defined (_MSC_VER)  && !defined (__WAT32__)
00549      if (!scroll_flag)  clrscr();
00550 #endif
00551      if (ad_printf) (*ad_printf)("Intermediate statistics: ");
00552    }
00553    llog=0;
00554    goto label7003;
00555 label7020:
00556    if (iprint>0)
00557    {
00558      if (ad_printf) (*ad_printf)("*** hessian not positive definite\n");
00559    }
00560          #ifdef __ZTC__
00561          if (ireturn <= 0)
00562          {
00563            disp_close();
00564          }
00565          #endif
00566          return;
00567    }
00568 
00573 dvector update(int nvar, int iter, int m, const dvector& g, const dmatrix& _s,
00574   dmatrix& y, const dvector& x, const dvector& _xold, const dvector& _gold,
00575   const dvector& _xrho)
00576   {
00577     dvector& xold= (dvector&) _xold;
00578     dmatrix& s= (dmatrix&) _s;
00579     dvector& gold= (dvector&) _gold;
00580     dvector& xrho=(dvector&)_xrho;
00581     dvector beta(1,nvar);
00582     dvector alpha(0,m);
00583     dvector r(1,nvar);
00584     dvector t(1,nvar);
00585     int m1=m+1;
00586     if (iter<1)
00587     {
00588       xold=x;
00589       gold=g;
00590       r=g;
00591     }
00592     else
00593     {
00594       int k=iter-1;
00595       int k1=k%(m1);
00596       y(k1)=g-gold;
00597       s(k1)=x-xold;
00598       xrho(k1)=1./(y(k1)*s(k1));
00599       xold=x;
00600       gold=g;
00601       int i;
00602       int lb=k-m+1;
00603       if (lb <0) lb=0;
00604       t=g;
00605       for (i=k;i>=lb;i--)
00606       {
00607         int i1=i%(m1);
00608         //int i2=(i+1)%(m1);
00609         {
00610           alpha(i-lb)=xrho(i1)*(s(i1)*t);
00611         }
00612         t-=alpha(i-lb)*y(i1);
00613       }
00614       r=t;
00615       for (i=lb;i<=k;i++)
00616       {
00617         int i1=i%(m1);
00618         r+= (alpha(i-lb)-xrho(i1)*(y(i1)*r)) * s(i1);
00619       }
00620     }
00621     return -1.0*r;
00622   }