ADMB Documentation  11.1.1922
 All Classes Files Functions Variables Typedefs Friends Defines
conjprod.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: conjprod.cpp 1893 2014-04-16 20:28:36Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 // this is to get UNIX systems to use getchar
00012 // #define UNIXKLUDGE
00013 
00014 //#define DIAG
00015 #ifdef _MSC_VER
00016   #include <conio.h>
00017 #else
00018   #include <unistd.h>
00019 #endif
00020 
00021 #include <fvar.hpp>
00022 extern int ctlc_flag;
00023 
00024 #ifdef __TURBOC__
00025 #pragma hdrstop
00026 #include <conio.h>
00027 #include <iomanip.h>
00028 #endif
00029 
00030 #include <ctype.h>
00031 #include <string.h>
00032 
00033 //#define CUBIC_INTERPOLATION
00034 
00035 #ifdef __ZTC__
00036 #include <conio.h>
00037 #include <iomanip.hpp>
00038 #include <disp.h>
00039 
00040   // preprocessor defines and function definitions that emulate Borland
00041   // text screen managment functions using the Zortech disp_ package
00042   // note the order of includes is important - this stuff should come
00043   // after stdio.h, disp.h and anything that might also include these
00044 
00045 //#define if (ad_printf) (*ad_printf) disp_if (ad_printf) (*ad_printf)
00046 
00047   void gotoxy(int x, int y);
00048   void clrscr();
00049 
00050   struct text_info
00051   {
00052     unsigned char  winleft,   wintop;
00053     unsigned char  winright,  winbottom;
00054   //unsigned char  attribute, normattr;
00055     unsigned char  currmode;
00056     unsigned char  screenheight;
00057     unsigned char  screenwidth;
00058     unsigned char  curx, cury;
00059   };
00060 
00061   void gettextinfo(struct text_info *r);
00062 #endif
00063 
00064 #if defined(UNIXKLUDGE) || \
00065   ((defined(__SUN__) || defined(__GNU__)) && !defined(__GNUDOS__))
00066 #include <iostream.h>
00067 #include <signal.h>
00068 #define getch getchar
00069   #if defined(__HP__) || defined(linux)
00070   extern "C" void onintr(int k);
00071   #else
00072   extern "C" int onintr(int* k);
00073   #endif
00074 #endif
00075 
00076 #if defined(_MSC_VER)
00077   void __cdecl clrscr(){}
00078 #else
00079   extern "C" void clrscr(){}
00080 #endif
00081 
00082 #ifdef __NDPX__
00083 
00084 #include <iostream.hxx>
00085   extern "C" {
00086 #include <grex.h>
00087     void clrscr();
00088   };
00089 
00090   void gotoxy(int x, int y);
00091 
00092   struct text_info
00093   {
00094     unsigned char  winleft,   wintop;
00095     unsigned char  winright,  winbottom;
00096   //unsigned char  attribute, normattr;
00097   //unsigned char  currmode;
00098     unsigned char  screenheight;
00099     unsigned char  screenwidth;
00100     unsigned char  curx, cury;
00101   };
00102 
00103   void gettextinfo(struct text_info *r);
00104 #endif
00105 
00106 class fmmc;
00107 
00108 double do_interpolate(const double& fret, const double& left_bracket,
00109   double& left_bracket_value, const dvector& left_bracket_gradient,
00110   double& right_bracket, const double& right_bracket_value,
00111   dvector& right_bracket_gradient, const dvector& theta, const dvector& d,
00112   const int& J, long int& ifn, const double& crit1,
00113   int& int_flag, const double& rho_1, const double& Psi_2, const dvector& g1);
00114 
00115 void do_extrapolate(const double& left_bracket,
00116   const double& left_bracket_value, dvector& left_bracket_gradient,
00117   const double& right_bracket, double& right_bracket_value,
00118   const dvector& right_bracket_gradient, const dvector& theta,
00119   dvector& d, const int& J, const double& rho_0, long int& ifn,
00120   const int& ifnex, int& ext_flag, const double& rho_1, const double& rf,
00121   const dvector& g1);
00122 
00123 double mylinmin(const double& fret, const double& Phi_i, const dvector& theta1,
00124   const dvector& q_i, fmmc& cs);
00125 
00126 void bracket_report(const dvector& theta, const double& left_bracket,
00127                     double& right_bracket, const dvector& d);
00128 
00129 double cubic_interpolation(const double& u, const double& v, const double& a,
00130   const double& b, double& ap, const double& bp);
00131 
00132 double Phi(const dvector&);
00133 double min(const double&, const double&);
00134 double max(const double&, const double&);
00135 
00136 #include <math.h>
00137 
00138 #define ITMAX 5000
00139 #define EPS 1.0e-10
00140 #define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n);
00141 
00146 void fmmc::fmin(const double& fret, const dvector& p, const dvector& gg)
00147 {
00148   dfn=0.0;;
00149   maxfn_flag=0;
00150   iexit=0;
00151   ihflag=0;
00152 #ifdef __ZTC__
00153     if (disp_inited == 0)
00154     {
00155        disp_open();
00156        disp_usebios();
00157     }
00158 #endif
00159   int n=p.size();
00160   dvector& xi=*(this->xi);
00161   dvector& h=*(this->h);
00162   dvector& g=*(this->g);
00163   double& fp=this->fp;
00164   //int& its=this->its;
00165   int& J=this->J;
00166   if (this->frp_flag > 0) this->ifn++;
00167 
00168   if (ireturn >= 3)
00169   {
00170     {
00171       derch(fret, p , gg, n, ireturn);
00172     }
00173     return;
00174   }
00175 
00176   if (this->frp_flag>1)
00177   {
00178     if (fret < fbest)
00179     {
00180 #ifdef DIAG
00181         cerr << " Improving fbest  old fbest = " << fbest << "\n" <<
00182                 "                  new fbest = " << fret << "\n";
00183 #endif
00184 
00185       int con_flag;
00186       for (int jj=gg.indexmin();jj<=gg.indexmax();jj++)
00187       {
00188         con_flag=1;
00189         if (fabs(gg(jj)) > crit)
00190         {
00191           con_flag=0;
00192           break;
00193         }
00194       }
00195 
00196       if ( con_flag==1)
00197       {
00198         this->ireturn=-1;
00199         {
00200           if (iprint>0)
00201           {
00202             if (ad_printf)
00203               (*ad_printf)("Gradient magnitude criterion satisfied\n");
00204             if (ad_printf)
00205           (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00206             n, iter, ifn);
00207             if (ad_printf)
00208    (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00209                fret, max(fabs(gg)) );
00210             fmmdisp(p, gg, n, this->scroll_flag); //fmc);
00211           }
00212         }
00213         return;
00214       }
00215       fbest=fret;
00216       *xbest=p;
00217       *gbest=gg;
00218     }
00219   }
00220 
00221   if (this->frp_flag==1) goto label800;
00222   if (this->frp_flag==2) goto label2000;
00223   if (this->frp_flag==3) goto label3000;
00224 
00225 #if defined(__SUN__) && !defined(__GNUDOS__) && !defined(_MSC_VER)
00226   #ifdef __HP__
00227   signal(SIGINT, &onintr);
00228   #else
00229   signal(SIGINT, (SIG_PF)&onintr);
00230   #endif
00231 #endif
00232 #if (defined(__GNU_) && !defined(__GNUDOS__)) || defined(UNIXKLUDGE)
00233   signal(SIGINT, &onintr);
00234 #endif
00235   this->J=1;
00236   this->rho_min=1.e-10;
00237   this->converge_flag=0;
00238 
00239   this->frp_flag=1;
00240   this->ireturn=1;
00241   return;
00242 label800:
00243 
00244   fbest=fret;
00245   *xbest=p;
00246   *gbest=gg;
00247 
00248   this->frp_flag=0;
00249   this->ireturn=0;
00250   xi=gg;
00251   this->fp=fret;
00252   //(*dfunc)(p,xi);
00253 
00254   if (iprint>0)
00255   {
00256     if (!(iter%iprint)&& (iprint > 0) )
00257     {
00258 #if !defined (__WAT32__) && !defined (_MSC_VER)
00259         if (!scroll_flag) clrscr();
00260 #endif
00261       if (ad_printf) (*ad_printf)("Initial statistics: ");
00262       if (ad_printf)
00263  (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00264        n, iter, ifn);
00265       if (ad_printf)
00266  (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00267        fbest, max(fabs(*gbest)) );
00268       fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00269     }
00270   }
00271 
00272   {
00273     for (int j=1;j<=n;j++)
00274     {
00275       g[j] = -xi[j];
00276       xi[j]=h[j]=g[j];
00277     }
00278   }
00279 
00280   this->ifnex=0;
00281 
00282   //this->its=1;
00283 
00284 label1000:
00285     iter++;
00286   label5:
00287 #if ((defined(__SUN__) || defined(__GNU__)) && !defined(__GNUDOS__)) \
00288   || defined(UNIXKLUDGE)
00289     if (ctlc_flag)
00290 #else
00291     if ( kbhit() )
00292 #endif
00293     {
00294       int c = toupper(getch());
00295       if ( c == 'C')
00296       {
00297         ireturn = 3;
00298         {
00299           derch(fret, p , gg, n, ireturn);
00300         }
00301         return;
00302       }
00303       else
00304       {
00305         if ( c == 'Q' || c == 'N' )
00306         {
00307           quit_flag=c;
00308           this->ireturn=-1;
00309           {
00310             if (iprint>0)
00311             {
00312               if (ad_printf) (*ad_printf)("User initiated interrupt\n");
00313               if (ad_printf) (*ad_printf)(" - final statistics:\n");
00314               if (ad_printf)
00315  (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00316                 n, iter, ifn);
00317               if (ad_printf)
00318  (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00319                 fbest, max(fabs(*gbest)) );
00320               fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00321             }
00322           }
00323           p=*xbest;
00324           gg=*gbest;
00325           return;
00326         }
00327         else
00328         {
00329           quit_flag=0;
00330         }
00331       }
00332     }
00333     {
00334       mylinmin(fret,fp,p,-1.*(g),*this);
00335     }
00336     if (this->lin_flag ==0) goto label10;
00337 
00338       p=*(this->extx);
00339       this->frp_flag=2;
00340       this->ireturn=1;
00341       return;
00342 
00343 
00344     label2000:
00345       this->ireturn=0;
00346       this->frp_flag=0;
00347       //this->extf=fcomp( *(this->extx),*(this->extg) );
00348       this->extf=fret;
00349       *(this->extx)=p;
00350       *(this->extg)=gg;
00351       goto label5;
00352   label10:
00353 
00354     {
00355       for ( int i=1; i<=9; i++)
00356       {
00357         (*funval)[i]= (*funval)[i+1];
00358       }
00359     }
00360 
00361     (*funval)[10] = fret;
00362 
00363     if ( iter>15 && fabs( (*funval)[1]-(*funval)[10])< min_improve )
00364     {
00365       ihang = 1;
00366       this->ireturn=-1;
00367       {
00368         if (iprint>0)
00369         {
00370           if (ad_printf)
00371  (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00372             n, iter, ifn);
00373           if (ad_printf)
00374  (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00375              fbest, max(fabs(*gbest)) );
00376           fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00377         }
00378       }
00379       p=*xbest;
00380       gg=*gbest;
00381       return;
00382     }
00383 
00384 
00385   /*  Moved this up to entry point
00386     {
00387       int con_flag;
00388       for (int jj=gg.indexmin();jj<=gg.indexmax();jj++)
00389       {
00390         con_flag=1;
00391         if (fabs(g(jj)) > crit)
00392         {
00393           con_flag=0;
00394           break;
00395         }
00396       }
00397 
00398       if ( con_flag==1)
00399       {
00400         this->ireturn=-1;
00401         {
00402           if (iprint>0)
00403           {
00404             if (ad_printf)
00405  (*ad_printf)("Gradient magnitude criterion satisfied\n");
00406             if (ad_printf)
00407  (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00408               n, iter, ifn);
00409             if (ad_printf)
00410  (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00411              fbest, max(fabs(*gbest)) );
00412             fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00413           }
00414         }
00415         p= *xbest;
00416         gg= *gbest;
00417         return;
00418       }
00419     }
00420  */
00421 
00422 
00423     if ( ifn > maxfn )
00424     {
00425       if (iprint>0)
00426       {
00427         if (ad_printf)
00428         {
00429           (*ad_printf)("Maximum number of function evaluations exceeded\n");
00430           (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00431             n, iter, ifn);
00432         (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00433             fbest, max(fabs(*gbest)) );
00434         }
00435         fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00436       }
00437       p=*xbest;
00438       gg=*gbest;
00439       this->ireturn=-1;
00440       return;
00441     }
00442 
00443     if (iprint>0)
00444     {
00445       if (!(iter%iprint)&&(iprint>0))
00446       {
00447         {
00448 #if !defined (__WAT32__) && !defined (_MSC_VER)
00449             if (!scroll_flag) clrscr();
00450 #endif
00451           if (ad_printf) (*ad_printf)("Intermediate statistics: ");
00452         }
00453         if (ad_printf)
00454  (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00455           n, iter, ifn);
00456         if (ad_printf)
00457  (*ad_printf)("Function value %le; maximum gradient component mag %le\n",
00458           fbest, max(fabs(*gbest)) );
00459         fmmdisp(*xbest, *gbest, n, this->scroll_flag); //fmc);
00460       }
00461     }
00462 
00463     this->frp_flag=3;
00464     this->ireturn=1;
00465     // return;
00466   label3000:
00467     this->frp_flag=0;
00468     this->ireturn=0;
00469     xi=gg;
00470     this->fp=fret;
00471 
00472     //fp=fcomp(p,xi);
00473 
00474     this->dgg=this->gg=0.0;
00475     {
00476       for (int j=1;j<=n;j++)
00477       {
00478         this->gg += g[j]*g[j];
00479 /*      dgg += xi[j]*xi[j];  */
00480         this->dgg += (xi[j]+g[j])*xi[j];
00481       }
00482     }
00483     if (this->gg == 0.0)
00484     {
00485       this->ireturn=-1;
00486       return;
00487     }
00488     this->gam=this->dgg/this->gg;
00489     {
00490       for (int j=1;j<=n;j++)
00491       {
00492         g[j] = -xi[j]; // g seems to hold the negative gradient
00493         xi[j]=h[j]=g[j]+this->gam*h[j];
00494       }
00495     }
00496 //  if (this->iter <= ITMAX) goto label1000;
00497     goto label1000;
00498 
00499   cerr << "Too many iterations in FRPRMN\n";
00500 }
00501 
00502 #undef ITMAX
00503 #undef EPS
00504 #undef FREEALL
00505 
00511 double mylinmin(const double& fret, const double& Phi_i, const dvector& theta1,
00512   const dvector& q_i,
00513   fmmc& cs)
00514 {
00515   // d is the current direction for the 1 dimensional search
00516   // q_i is the gradient at the current best point
00517 
00518   dvector& theta= *(cs.theta);
00519   if (cs.lin_flag==1)goto label1;
00520   if (cs.lin_flag==2)goto label2;
00521   if (cs.lin_flag==3)goto label3;
00522   {
00523     *(cs.theta) = theta1;
00524     *(cs.d)=*(cs.xi);
00525     cs.rho_0=pow(2.,-cs.J)/norm(q_i);
00526     cs.gamma=q_i* *(cs.d);
00527     cs.lin_flag=1;
00528     *(cs.extx)=theta+cs.rho_0* *(cs.d);
00529   }
00530   return 0;
00531 label1:
00532   {
00533     cs.lin_flag=0;
00534     cs.Psi_0=cs.extf;
00535     *(cs.g2)=*(cs.extg);
00536 
00537     //double rho_star=gamma*rho_0*rho_0/(2*(gamma*rho_0+Phi_i-Psi_0));
00538 
00539     cs.ext_flag=0;
00540     cs.int_flag=0;
00541     cs.dir_deriv=*(cs.d)* *(cs.g2);
00542     cs.left_bracket_value=Phi_i;
00543     *(cs.left_bracket_gradient)=q_i;
00544     cs.left_bracket=0;
00545   }
00546 
00547 
00548 //cout << "Check " << cs.Psi_0 << " " << Phi_i << "  " << cs.dir_deriv << "\n";
00549 
00550   if (!(cs.Psi_0 < Phi_i && cs.dir_deriv < 0)) goto label555;
00551   // Extrapolate to get a right bracket
00552 
00553     label100:
00554     {
00555 #ifdef DIAG
00556         cerr << "doing extrapolation\n";
00557 #endif
00558       do_extrapolate(cs.left_bracket,cs.left_bracket_value,
00559         *(cs.left_bracket_gradient), cs.right_bracket,
00560         cs.right_bracket_value,*(cs.right_bracket_gradient),theta,
00561         *(cs.d),cs.J,cs.rho_0,cs.ifn,cs.ifnex,cs.ext_flag,cs.rho_1,
00562         cs.Psi_1,*(cs.grad) );
00563     }
00564     if (cs.ext_flag==0)
00565     {
00566       goto label120;
00567     }
00568     else
00569     {
00570       {
00571        cs.lin_flag=2;
00572        *(cs.extx)=theta+cs.rho_1* *(cs.d);
00573        //extx=theta+rho_1*d;
00574       }
00575       return 0;
00576     label2:
00577       cs.lin_flag=0;
00578       //Psi_1=fcomp(theta+rho_1*d,grad);
00579       cs.Psi_1=cs.extf;
00580       *(cs.grad)= *(cs.extg);
00581     }
00582     goto label100;
00583 
00584   //label120:
00585   //  int itemp=1;
00586 
00587 label555:
00588   //else
00589   // We have a bracket
00590   {
00591 #ifdef DIAG
00592         cerr << "We have a right bracket\n";
00593 #endif
00594      cs.right_bracket=cs.rho_0;
00595      cs.right_bracket_value=cs.Psi_0;
00596      *(cs.right_bracket_gradient)=*(cs.g2);
00597   }
00598  label120:
00599   // cout << "After extrapolation\n";
00600   bracket_report(theta,cs.left_bracket,cs.right_bracket,*(cs.d));
00601 
00602   // We have a bracket [theta,theta+rho_0*d] so interpolate
00603   {
00604   label1100:
00605    {
00606     //cout << "Intrap\n";
00607     cs.rho_i=do_interpolate(fret,cs.left_bracket,cs.left_bracket_value,
00608       *(cs.left_bracket_gradient),cs.right_bracket,cs.right_bracket_value,
00609       *(cs.right_bracket_gradient),theta,*(cs.d),cs.J,cs.ifn,cs.crit1,
00610       cs.int_flag,cs.rho_1,cs.Psi_1,*(cs.grad) );
00611    }
00612     if (cs.int_flag==0)
00613     {
00614       goto label1120;
00615     }
00616     else
00617     {
00618      {
00619       cs.lin_flag=3;
00620       *(cs.extx)=theta+cs.rho_1* *(cs.d) ;
00621      }
00622       return 0;
00623     label3:
00624       cs.lin_flag=0;
00625       //Psi_1=fcomp(theta+rho_1*d,grad);
00626       cs.Psi_1=cs.extf;
00627       *(cs.grad)=*(cs.extg);
00628     }
00629     goto label1100;
00630   label1120:
00631     int itemp=1;;
00632   }
00633   theta=theta+cs.rho_i* *(cs.d) ;
00634   return cs.converge_flag;
00635 }
00636 
00641 double do_interpolate(const double& fret, const double& left_bracket,
00642   double& left_bracket_value, const dvector& left_bracket_gradient,
00643   double& right_bracket, const double& right_bracket_value,
00644   dvector& right_bracket_gradient, const dvector& theta, const dvector& d,
00645   const int& _J, long int& ifn, const double& crit1,
00646   int& int_flag, const double& rho_1, const double& Psi_2, const dvector& g1)
00647 {
00648   double rho_min=1.e-10;
00649   int& J = (int&) _J;
00650   static double rho_star;
00651   static double dir_deriv;
00652   //double Psi_2;
00653   //dvector g1(1,d.size());
00654   static double gamma;
00655   static double gamma1;
00656   static double rho_0;
00657   if (int_flag ==1) goto label200;
00658   J+=1;
00659 label120:
00660   // do
00661   //{
00662     gamma=left_bracket_gradient*d;
00663 #ifdef CUBIC_INTERPOLATION
00664       gamma1=right_bracket_gradient*d;
00665 #endif
00666     rho_0=right_bracket-left_bracket;
00667 
00668 #ifdef DIAG
00669       cout << "f'(x) " << gamma  << "   "
00670          << "f'(y) " << gamma1 << endl;
00671 
00672       cout << "f(x)  " << left_bracket_value  << "   "
00673          << "f(y)  " << right_bracket_value << endl;
00674 
00675       cout << "x     " << left_bracket  << "   "
00676          << "y     " << right_bracket << endl;
00677 #endif
00678 
00679     if (gamma==0.0)
00680     {
00681       rho_star=0.;
00682     }
00683     else
00684     {
00685 #ifdef CUBIC_INTERPOLATION
00686         rho_star=cubic_interpolation(left_bracket,right_bracket,
00687           left_bracket_value,right_bracket_value,gamma,gamma1);
00688 #else
00689         double step=gamma*rho_0*rho_0/
00690           (2*(gamma*rho_0+left_bracket_value-right_bracket_value));
00691          // rho_star=left_bracket+gamma*rho_0*rho_0/
00692          //(2*(gamma*rho_0+left_bracket_value-right_bracket_value));
00693          double width=right_bracket-left_bracket;
00694          if (step<.25*width)step=.25*width;
00695          if (step>.75*width)step=.75*width;
00696          rho_star=left_bracket+step;
00697 
00698 #endif
00699     }
00700 
00701     if (rho_star < left_bracket )
00702     {
00703       cout << " rho_star out of range ..value changed\n";
00704     }
00705 
00706     if (rho_star > right_bracket)
00707     {
00708       cout << " rho_star out of range ..value changed\n";
00709       rho_star=.001*left_bracket+.999*right_bracket;
00710     }
00711 
00712     if (rho_star <= rho_min)
00713     {
00714       return rho_min;
00715     }
00716     int_flag=1;
00717     rho_1=rho_star;
00718     return rho_min;
00719 label200:
00720     int_flag=0;
00721     rho_1=rho_star;
00722     //Psi_2=fcomp(theta+rho_star*d,g1);
00723     //J+=1;
00724 
00725     dir_deriv=d*g1;
00726 
00727     //cout << "Check2 " << Psi_2 << " " << left_bracket_value << "  " <<
00728     //   d*g1 << "\n";
00729 
00730     if (Psi_2 < left_bracket_value && (d*g1) < 0)
00731     {
00732 #ifdef DIAG
00733         cout << " Before interpolation -- new left side\n";
00734         bracket_report( theta,left_bracket,right_bracket,d);
00735 #endif
00736 
00737       left_bracket =rho_star;
00738       left_bracket_value=Psi_2;
00739       left_bracket_gradient=g1;
00740 
00741 #ifdef DIAG
00742         cout << " After interpolation -- new left side\n";
00743         bracket_report( theta,left_bracket,right_bracket,d);
00744 #endif
00745     }
00746     else
00747     {
00748 #ifdef DIAG
00749         cout << " Before interpolation -- new right side\n";
00750         bracket_report( theta,left_bracket,right_bracket,d);
00751 #endif
00752       right_bracket =rho_star;
00753       right_bracket_value=Psi_2;
00754       right_bracket_gradient=g1;
00755 
00756 
00757 #ifdef DIAG
00758         cout << " After interpolation -- new right side\n";
00759         bracket_report( theta,left_bracket,right_bracket,d);
00760 #endif
00761     }
00762 
00763   //while (dir_deriv > crit1
00764   //   && (right_bracket-left_bracket)> 1.e-10);
00765 
00766    double cos_theta=d*g1/(norm(d)*norm(g1));
00767 #ifdef DIAG
00768      cerr << " The cosine of angle between the search direction and \n"
00769        " the gradient is " << cos_theta << "\n";
00770 #endif
00771 
00772   if (cos_theta > crit1
00773      && (right_bracket-left_bracket)> 1.e-10)
00774   {
00775     if (cos_theta > .05)
00776     {
00777       J+=1;
00778     }
00779     goto label120;
00780   }
00781 
00782   if (cos_theta < -crit1 && (right_bracket-left_bracket)> 1.e-10)
00783   {
00784     goto label120;
00785   }
00786 
00787 #ifdef DIAG
00788      if (cos_theta < crit1)
00789      {
00790       cerr << " Leaving interpolation routine because"
00791         " the cosine of angle between the \n search direction and "
00792        " the gradient is < crit1 = " << crit1  << "\n";
00793      }
00794 #endif
00795 
00796 
00797   fret=Psi_2;
00798   return rho_star;
00799 }
00800 
00805 void do_extrapolate(const double& left_bracket,
00806   const double& left_bracket_value, dvector& left_bracket_gradient,
00807   const double& right_bracket, double& right_bracket_value,
00808   const dvector& right_bracket_gradient, const dvector& theta,
00809   dvector& d, const int& J, const double& rho_0, long int& ifn,
00810   const int& ifnex, int& ext_flag, const double& rho_1, const double& rf,
00811   const dvector& g1)
00812 {
00813   if (ext_flag==1) goto label1500;
00814   J=J/2;
00815   rho_1=4.*rho_0;
00816   //dvector g1(1,d.size());
00817 #ifndef __NDPX__
00818   double Psi_1;
00819 #endif
00820 label1000:
00821     ext_flag=1;
00822     return;
00823 label1500:
00824 #ifdef __NDPX__
00825   double Psi_1;
00826 #endif
00827     ext_flag=0;
00828     //Psi_1=fcomp(theta+rho_1*d,g1);
00829     Psi_1=rf;
00830     ifnex++;
00831     if (Psi_1 >= left_bracket_value || d*g1 >0) goto label2000;
00832     left_bracket_value=Psi_1;
00833     left_bracket_gradient=g1;
00834     left_bracket=rho_1;
00835     rho_1=4*rho_1;
00836     goto label1000;
00837 label2000:
00838 
00839   right_bracket=rho_1;
00840   right_bracket_value=Psi_1;
00841   right_bracket_gradient=g1;
00842 }
00843 
00848 void bracket_report(const dvector& theta, const double& left_bracket,
00849                     double& right_bracket, const dvector& d)
00850 {
00851   double f=0;
00852   double fp1=0;
00853   double fp2=0;
00854   dvector g(1,d.size());
00855   dvector u(1,d.size());
00856   ivector ii(1,3);
00857   int one=1;
00858   ii.fill_seqadd(one,one);
00859 
00860 #ifdef DIAG
00861     cout << "lb " <<  setprecision(4) << setw(12) << left_bracket
00862       << "rb " << setprecision(4) << setw(12)  << right_bracket
00863       << setprecision(4) << setw(12) << theta(ii) <<"\n";
00864 #endif
00865 }
00866 
00871 double cubic_interpolation(const double& u, const double& v, const double& aa,
00872   const double& bb, double& ap, const double& bp)
00873 {
00874   //cout <<"Begin cube\n";
00875   dmatrix M(1,4,1,4);
00876 
00877   M(1,1)=1;M(2,1)=0;M(3,1)=1;M(4,1)=0; // First column
00878 
00879   M(1,2)=u;M(2,2)=1;M(3,2)=v;M(4,2)=1; // Second column
00880 
00881   for (int i=3;i<=4;i++)
00882   {
00883     M(1,i)=u*M(1,i-1);
00884     M(2,i)=(i-1)*M(1,i-1);
00885     M(3,i)=v*M(3,i-1);
00886     M(4,i)=(i-1)*M(3,i-1);
00887   }
00888 
00889   dvector d(1,4);
00890   d(1)=aa;
00891   d(2)=ap;
00892   d(3)=bb;
00893   d(4)=bp;
00894 
00895   dvector e = inv(M)*d;
00896   double a,b,c;
00897 
00898   a=3*e(4);
00899   b=2*e(3);
00900   c=e(2);
00901 
00902   double q;
00903   if (a !=0)
00904   {
00905     if (b>0)
00906     {
00907       double y=b*b-4*a*c;
00908       if (y<0) return (u+v)/2.;
00909       q=-.5*(b+sqrt(y));
00910     }
00911     else
00912     {
00913       double y=b*b-4*a*c;
00914       if (y<0) return (u+v)/2.;
00915       q=-.5*(b-sqrt(b*b-4*a*c));
00916     }
00917     double x1,x2;
00918     x1=q/a;       // x1 and x2 are the two roots of the quadratic
00919     x2=c/q;       // equation that is the max andmin of the cubic
00920                   // polynomial
00921     double sgn1,sgn2;
00922     sgn1=b+2*a*x1;
00923     sgn2=b+2*a*x2;
00924     if (sgn1>0)
00925     {
00926       return x1;
00927     }
00928     else
00929     {
00930       return x2;
00931     }
00932   }
00933   else
00934   {
00935     // coeffcient of the quadratic term = 0
00936     return -c/b;
00937   }
00938     //cout <<"end cube\n";
00939 }
00940 
00941 #undef CUBIC_INTERPOLATION
00942 
00947 fmmc::fmmc(const int& n)
00948 {
00949   ctlc_flag = 0;
00950   maxfn=500;
00951   lin_flag=0;
00952   ext_flag=0;
00953   int_flag=0;
00954   ifnex=0;
00955   frp_flag=0;
00956   quit_flag=0;
00957   #ifdef __ZTC__
00958     scroll_flag=0;
00959   #else
00960     scroll_flag=1;
00961   #endif
00962   crit=1.e-4;
00963   crit1=1.e-2;
00964   min_improve=1.e-6;
00965   ireturn=0;
00966   iprint=1;
00967   imax=30;
00968   ihang=0;
00969   iter=0;
00970   ifn=0;
00971   left_bracket_gradient=new dvector(1,n);
00972   right_bracket_gradient=new dvector(1,n);
00973   funval=new dvector(1,10);
00974   g=new dvector(1,n);
00975   h=new dvector(1,n);
00976   xi=new dvector(1,n);
00977   d=new dvector(1,n);
00978   extx=new dvector(1,n);
00979   g2=new dvector(1,n);
00980   grad=new dvector(1,n);
00981   extg=new dvector(1,n);
00982   theta=new dvector(1,n);
00983   xbest=new dvector(1,n);
00984   gbest=new dvector(1,n);
00985 }
00986 
00991 fmmc::~fmmc(void)
00992 {
00993   delete left_bracket_gradient;
00994   delete right_bracket_gradient;
00995   delete funval;
00996   delete g;
00997   delete h;
00998   delete xi;
00999   delete d;
01000   delete extx;
01001   delete g2;
01002   delete grad;
01003   delete extg;
01004   delete theta;
01005   delete xbest;
01006   delete gbest;
01007 }
01008 
01013 void derch(const double& f, const dvector& _x, const dvector& _gg, int n,
01014   const int & _ireturn)
01015 {
01016   int& ireturn=(int&) _ireturn;
01017   dvector& x = (dvector&) _x;
01018   dvector& gg = (dvector&) _gg;
01019   static long int i, n1 ,n2;
01020   static double fsave;
01021   static double s, f1, f2, g2, xsave;
01022   static long int j = 1;
01023   static int si;
01024   si=gg.indexmax();
01025   static dvector g(1,si);
01026 
01027   if (ireturn == 4 ) goto label4;
01028   else if (ireturn == 5) goto label5;
01029   g=gg;
01030   while (j > 0)
01031   {
01032     //cout << "\nEntering derivative checker.\n";
01033     cout << "\n Enter index (1 ... "<< n <<") of derivative to check.";
01034     cout << "  To check all derivatives, enter 0: ";
01035     cin >> j;
01036     if (j <= 0)
01037     {
01038       cout << "\n   Checking all derivatives. Press X to terminate checking.\n";
01039       n1 = 1;
01040       n2 = n;
01041     }
01042     else
01043     {
01044       n1 = j;
01045       n2 = j;
01046     }
01047     cout << "   Enter step size (to quit derivative checker, enter 0): ";
01048     cin >> s;
01049     cout <<
01050        "\n       X           Function     Analytical     Finite Diff;  Index\n";
01051 
01052     if (s <= 0) ad_exit(0);
01053 
01054     for (i=n1; i<=n2; i++)
01055     {
01056       xsave=x(i);
01057       x(i)=xsave+s;
01058       fsave = f;
01059       ireturn = 4; // fgcomp(&f1,x,g1,n, params, vars);
01060       return;
01061 
01062     label4:
01063       f1 = f;
01064       x(i)=xsave-s;
01065       ireturn= 5; //fgcomp(&f2,x,g1,n, params, vars);
01066       return;
01067 
01068     label5:
01069       f2 = f;
01070       f = fsave;
01071       x(i)=xsave;
01072       g2=(f1-f2)/(2.*s);
01073       if (ad_printf) (*ad_printf)("  %12.5e  %12.5e  %12.5e  %12.5e ; %5d \n",
01074               x(i), f, g(i), g2, i);
01075     } // for loop
01076   } // while (j > 0)
01077 //  ireturn = 2;
01078   ad_exit(0);
01079 }
01080 // #undef DIAG