ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
model7.cpp
Go to the documentation of this file.
00001 
00007 #  include <df1b2fun.h>
00008 #  include <admodel.h>
00009 //#include <parallel.h>
00010 #include <signal.h>
00011 #ifdef __MINGW64__
00012   #include <cassert>
00013 #endif
00014 
00015 #if defined(_MSC_VER)
00016 void strip_full_path(const adstring& _s)
00017 {
00018   adstring& s = (adstring&)_s;
00019   size_t n = s.size();
00020   size_t i = n - 1;
00021   for (; i >= 1; i--)
00022   {
00023     if ( s(i) == '\\' || s(i) == '/' || s(i) == ':') break;
00024   }
00025   s = s(i + 1, n);
00026 }
00027 #endif
00028 
00029 void set_signal_handlers(void)
00030 {
00031   signal(SIGFPE,exit_handler);
00032   signal(SIGSEGV,exit_handler);
00033   signal(SIGILL,exit_handler);
00034   signal(SIGINT,exit_handler);
00035 }
00036 
00037 ad_comm::ad_comm(int _argc,char * _argv[])
00038 {
00039   if (option_match(_argc,_argv,"-version") > -1
00040    || option_match(_argc,_argv,"--version") > -1)
00041   {
00042     void banner(const adstring& program_name);
00043     banner(_argv[0]);
00044 
00045     exit(0);
00046   }
00047 
00048   ad_comm::argc=_argc;
00049   ad_comm::argv=_argv;
00050   if (option_match(_argc,_argv,"-time")>-1)
00051   {
00052     time_flag=1;
00053   }
00054   else
00055   {
00056     time_flag=0;
00057   }
00058   if (time_flag)
00059   {
00060     if (!ptm)
00061     {
00062       ptm=new adtimer();
00063     }
00064     if (!ptm1)
00065     {
00066       ptm1=new adtimer();
00067     }
00068   }
00069   no_atlas_flag=0;
00070   if (option_match(_argc,_argv,"-noatlas")>-1) no_atlas_flag=1;
00071 
00072 #if defined(USE_ADPVM)
00073   int pvm_flag=0;
00074   if (option_match(_argc,_argv,"-slave")>-1)  pvm_flag=2;
00075   if (option_match(_argc,_argv,"-master")>-1) pvm_flag=1;
00076 
00077   if (pvm_flag)
00078     pvm_manager = new adpvm_manager(pvm_flag);
00079   else
00080     pvm_manager = NULL;
00081 
00082   if (pvm_manager)
00083   {
00084     if (pvm_manager->mode==2)  //slave
00085     {
00086       int on=0; int nopt=0;
00087       if ( (on=option_match(_argc,_argv,"-slave",nopt))>-1)
00088       {
00089         if (nopt ==1)
00090         {
00091           pvm_manager->slave_number=atoi(ad_comm::argv[on+1]);
00092         }
00093         else
00094         {
00095           cerr << "Wrong number of options to -slave -- must be 1"
00096                   " you have " << nopt << endl;
00097           ad_exit(1);
00098         }
00099       }
00100       if ( (on=option_match(_argc,_argv,"-slavedir",nopt))>-1)
00101       {
00102         if (nopt ==1)
00103         {
00104           ad_chdir(_argv[on+1]);
00105         }
00106         else
00107         {
00108           cerr << "Wrong number of options to -slavedir -- must be 1"
00109                   " you have " << nopt << endl;
00110         }
00111       }
00112     }
00113   }
00114 #endif
00115 
00116   /*
00117     if (option_match(_argc,_argv,"-gui")>-1)
00118     {
00119       void vm_initialize(void);
00120       vm_initialize();
00121     }
00122   */
00123   set_signal_handlers();
00124   adprogram_name=_argv[0];
00125   //int len=strlen(_argv[0]);
00126   //for (int i=1;i<=len;i++) adprogram_name[i]=tolower(adprogram_name[i]);
00127 #if defined(_MSC_VER)
00128   strip_full_path(adprogram_name);
00129 #endif
00130   adstring workdir;
00131   ad_getcd(workdir);
00132 //#define STR(x) STR2(x)
00133 //#define STR2(x) #x
00134   if (_argc>1)
00135   {
00136     if (option_match(_argc,_argv,"-?")>-1
00137      || option_match(_argc,_argv,"-help")>-1
00138      || option_match(_argc,_argv,"--help")>-1)
00139     {
00140       // remove path (if user runs -help)
00141       for (size_t i = adprogram_name.size(); i >= 1; i--)
00142       {
00143 #ifdef _WIN32
00144         if (adprogram_name(i) == '\\')
00145 #else
00146         if (adprogram_name(i) == '/')
00147 #endif
00148         {
00149           adprogram_name=adprogram_name(i+1,adprogram_name.size());
00150           break;
00151         }
00152       }
00153 
00154       //(*ad_printf)(" %s", (char*)admb_banner);
00155       (*ad_printf)( "Usage: %s [options]\n\n",(char*)(adprogram_name));
00156 
00157       (*ad_printf)( "Options:\n");
00158       (*ad_printf)( " -ainp FILE      change default ascii input parameter "
00159       "filename to FILE\n");
00160       (*ad_printf)( " -binp FILE      change default binary input parameter "
00161       "filename to FILE\n");
00162       (*ad_printf)( " -est            only do the parameter estimation\n");
00163       (*ad_printf)( " -noest          do not do the parameter estimation "
00164       "(optimization) \n");
00165       (*ad_printf)( " -ind FILE       change default input data filename to "
00166       "FILE\n");
00167       (*ad_printf)( " -lmn N          use limited memory quasi newton -- keep "
00168       "N steps\n");
00169       (*ad_printf)( " -lmn2 N         use other limited memory quasi newton -- "
00170       "keep N steps\n");
00171       (*ad_printf)( " -ilmn N         use other limited memory quasi newton "
00172       "for random effects models - keep N steps\n");
00173       (*ad_printf)( " -dd N           check derivatives after N function "
00174       "evaluations\n");
00175       (*ad_printf)( " -lprof          perform profile likelihood "
00176       "calculations\n");
00177       (*ad_printf)( " -maxph N        increase the maximum phase number to "
00178       "N\n");
00179       (*ad_printf)( " -mcdiag         use diagonal covariance matrix for mcmc "
00180       "with diagonal values 1\n");
00181       (*ad_printf)( " -mcmc [N]       perform markov chain monte carlo with N "
00182       "simulations\n");
00183       (*ad_printf)( " -mcmult N       multiplier N for mcmc default\n");
00184       (*ad_printf)( " -mcr            resume previous mcmc\n");
00185       (*ad_printf)( " -mcrb  N        reduce amount of correlation in the "
00186       "covariance matrix 1<=N<=9\n");
00187       (*ad_printf)( " -mcnoscale      don't rescale step size for mcmc "
00188       "depending on acceptance rate\n");
00189       (*ad_printf)( " -nosdmcmc       turn off mcmc histogram calcs to make "
00190       "mcsave run faster\n");
00191       (*ad_printf)( " -mcprobe N      use probing strategy for mcmc with "
00192       "factor N\n");
00193       (*ad_printf)( " -mcgrope N      Deprecated, same as -mcprobe\n");
00194       (*ad_printf)( " -mcseed N       seed for random number generator for "
00195       "markov chain monte carlo\n");
00196       (*ad_printf)( " -mcscale N      rescale step size for first N "
00197       "evaluations\n");
00198       (*ad_printf)( " -mcsave N       save the parameters for every Nth "
00199       "simulation\n");
00200       (*ad_printf)( " -mceval         go through the saved mcmc values from a "
00201       "previous mcsave\n");
00202       (*ad_printf)( " -mcu            use uniformaly distributed steps for "
00203       "mcmc instead of random normal\n");
00204       (*ad_printf)( " -crit N1,N2,... set gradient magnitude convergence "
00205       "criterion to N\n");
00206       (*ad_printf)( " -iprint N       print out function minimizer report "
00207       "every N iterations\n");
00208       (*ad_printf)( " -maxfn N1,N2,.. set maximum number opf function eval's "
00209       "to N\n");
00210       (*ad_printf)( " -rs             if function minimizer can't make "
00211       "progress rescale and try again\n");
00212       //(*ad_printf)( " -sp             for DLL running from splus write to "
00213       //"command window\n");
00214       (*ad_printf)( " -nox            suppress vector and gradient values in "
00215       "minimizer screen report\n");
00216       (*ad_printf)( " -phase N        start minimization in phase N\n");
00217       (*ad_printf)( " -simplex        use simplex for minimization -- "
00218       "deprecated, use -neldmead\n");
00219       (*ad_printf)( " -neldmead       use Nelder-Mead simplex algorithm for "
00220       "minimization\n");
00221       (*ad_printf)( " -nohess         don't do hessian or delta method for std "
00222       "dev\n");
00223       (*ad_printf)( " -eigvec         calculate eigenvectors of the Hessian\n");
00224       (*ad_printf)( " -sdonly         do delta method for std dev estimates "
00225       "without redoing hessian\n");
00226       (*ad_printf)( " -ams N          set arrmblsize to N "
00227       "(ARRAY_MEMBLOCK_SIZE)\n");
00228       (*ad_printf)( " -cbs N          set CMPDIF_BUFFER_SIZE to N "
00229       "(ARRAY_MEMBLOCK_SIZE)\n");
00230       (*ad_printf)( " -mno N          set the maximum number of independent "
00231       "variables to N\n");
00232       (*ad_printf)( " -mdl N          set the maximum number of dvariables to "
00233       "N\n");
00234       (*ad_printf)( " -gbs N          set GRADSTACK_BUFFER_SIZE to N "
00235       "(ARRAY_MEMBLOCK_SIZE)\n");
00236 #if defined(USE_ADPVM)
00237       (*ad_printf)( " -master         run as PVM master program\n");
00238       (*ad_printf)( " -slave          run as PVM slave program\n");
00239       (*ad_printf)( " -pvmtime        record timing information for PVM "
00240       "performance analysis\n");
00241 #endif
00242       (*ad_printf)( " -info           show how to cite ADMB, license, and "
00243       "acknowledgements\n");
00244       (*ad_printf)( " -version        show version information\n");
00245       (*ad_printf)( " -help           show this message\n\n");
00246     //if (function_minimizer::random_effects_flag)
00247     //{
00248       (*ad_printf)( "Random effects options if applicable\n");
00249       (*ad_printf)( " -nr N           maximum number of Newton-Raphson "
00250       "steps\n");
00251       (*ad_printf)( " -imaxfn N       maximum number of evals in quasi-Newton "
00252       "inner optimization\n");
00253       (*ad_printf)( " -is N           set importance sampling size to N\n");
00254       (*ad_printf)( " -isf N          set importance sampling size funnel "
00255       "blocks to N\n");
00256       (*ad_printf)( " -isdiag         print importance sampling diagnostics\n");
00257       (*ad_printf)( " -hybrid         do hybrid Monte Carlo version of MCMC\n");
00258       (*ad_printf)( " -hbf            set the hybrid bounded flag for bounded "
00259       "parameters\n");
00260       (*ad_printf)( " -hyeps          mean step size for hybrid Monte Carlo\n");
00261       (*ad_printf)( " -hynstep        number of steps for hybrid Monte "
00262       "Carlo\n");
00263       (*ad_printf)( " -noinit         do not initialize RE before inner "
00264       "optimization\n");
00265       (*ad_printf)( " -ndi N          set maximum number of separable calls\n");
00266       (*ad_printf)( " -ndb N          set number of blocks for RE derivatives "
00267       "(reduce temp file size)\n");
00268       (*ad_printf)( " -ddnr           use high precision Newton-Raphson, for "
00269       "banded Hessian case only\n");
00270       (*ad_printf)( " -nrdbg          verbose reporting for debugging "
00271       "newton-raphson\n");
00272 #  if defined(__MINI_MAX__)
00273       (*ad_printf)( " -mm N           do minimax optimization\n");
00274 #  endif
00275       (*ad_printf)( " -shess          use sparse Hessian structure inner "
00276       "optimzation\n\n");
00277 
00278       (*ad_printf)("Read online documentation at http://admb-project.org\n");
00279       (*ad_printf)("Contact <users@admb-project.org> for help.\n");
00280     //}
00281       ad_exit(0);
00282     }
00283     else if (option_match(_argc,_argv,"-info") > -1)
00284     {
00285       (*ad_printf)("ADMB Information\n");
00286       (*ad_printf)("================\n\n");
00287 
00288       (*ad_printf)("How to Cite ADMB\n");
00289       (*ad_printf)("----------------\n\n");
00290 
00291       (*ad_printf)("Fournier, D.A., H.J. Skaug, J. Ancheta, J. Ianelli, "
00292       "A. Magnusson, M.N. Maunder,\n");
00293       (*ad_printf)("A. Nielsen, and J. Sibert. 2012. AD Model Builder: using "
00294       "automatic\n");
00295       (*ad_printf)("differentiation for statistical inference of highly "
00296       "parameterized complex\n");
00297       (*ad_printf)("nonlinear models. Optim. Methods Softw. 27:233-249.\n\n");
00298 
00299       //(*ad_printf)(" %s", (char*)admb_banner);
00300       (*ad_printf)("License\n");
00301       (*ad_printf)("-------\n\n");
00302 
00303       (*ad_printf)("Copyright (c) 2008-2013\n");
00304       (*ad_printf)("Regents of the University of California and ADMB "
00305       "Foundation\n\n");
00306       (*ad_printf)("ADMB is free software and comes with ABSOLUTELY NO "
00307       "WARRANTY.\n");
00308       (*ad_printf)("You are welcome to redistribute it under certain "
00309       "conditions.\n\n");
00310       (*ad_printf)("AD Model Builder, or ADMB, was developed by David Fournier "
00311       "of Otter Research\n");
00312       (*ad_printf)("Ltd, Sidney, BC, Canada. In 2007, scientists from the "
00313       "University of Hawai'i at\n");
00314       (*ad_printf)("Manoa Pelagic Fisheries Research Program (John Sibert and "
00315       "Anders Nielsen) and\n");
00316       (*ad_printf)("the Inter-American Tropical Tuna Commission (Mark "
00317       "Maunder), in consultation with\n");
00318       (*ad_printf)("scientists from NOAA Fisheries (Richard Methot), created "
00319       "the non-profit ADMB\n");
00320       (*ad_printf)("Foundation (admb-foundation.org) with the goal of "
00321       "increasing the number of ADMB\n");
00322       (*ad_printf)("users by making the software free and open source. In "
00323       "partnership with NOAA\n");
00324       (*ad_printf)("Fisheries and the National Center for Ecological Analysis "
00325       "and Synthesis (NCEAS,\n");
00326       (*ad_printf)("www.nceas.ucsb.edu), the ADMB Foundation obtained funding "
00327       "from the Gordon and\n");
00328       (*ad_printf)("Betty Moore Foundation (www.moore.org) to acquire the "
00329       "copyright to the ADMB\n");
00330       (*ad_printf)("software suite, in order to make it broadly and freely "
00331       "available to the research\n");
00332       (*ad_printf)("community. In 2008 the copyright was transferred from "
00333       "Otter Research Ltd to the\n");
00334       (*ad_printf)("University of California. The binary files were released "
00335       "in November 2008 and\n");
00336       (*ad_printf)("the source code was released in December 2009. More "
00337       "information about the ADMB\n");
00338       (*ad_printf)("Project can be found at admb-project.org.\n\n");
00339       (*ad_printf)("ADMB was originally developed by David Fournier of Otter "
00340       "Research Ltd.\n\n");
00341       (*ad_printf)("It is now maintained by the ADMB Core Team, whose members "
00342       "are listed on\n");
00343       (*ad_printf)("http://admb-project.org/developers/core-team.\n");
00344 
00345       ad_exit(0);
00346     }
00347   }
00348   allocate();
00349 }
00350 
00351 ad_comm::ad_comm(void)
00352 {
00353   allocate();
00354 }
00355 
00356 void ad_comm::allocate(void)
00357 {
00358 #if defined (_WIN32)
00359   directory_prefix='\\';
00360 #else
00361   directory_prefix='/';
00362 #endif
00363   adstring tmpstring;
00364 
00365 #if defined(_MSC_VER)
00366   //remove path
00367   for (int i = (int)adprogram_name.size(); i >= 1; i--)
00368   {
00369     if (adprogram_name(i)==directory_prefix)
00370     {
00371       adprogram_name=adprogram_name(i+1,adprogram_name.size());
00372       break;
00373     }
00374   }
00375 #endif
00376 
00377 #if defined(_WIN32)
00378   // strip off the .exe
00379   #ifdef __MINGW64__
00380   size_t _n = adprogram_name.size();
00381   assert(_n <= INT_MAX);
00382   int n = (int)_n;
00383   #else
00384   int n = (int)adprogram_name.size();
00385   #endif
00386   if (n > 4)
00387   {
00388     if (adprogram_name(n - 3) == '.'
00389         && tolower(adprogram_name(n - 2)) == 'e'
00390         && tolower(adprogram_name(n - 1)) == 'x'
00391         && tolower(adprogram_name(n)) == 'e')
00392     {
00393       n -= 4;
00394     }
00395   }
00396   adprogram_name=adprogram_name(1,n);
00397 #endif
00398 
00399   // change the working directory name
00400   if (argc > 1)
00401   {
00402     int on=0;
00403     if ( (on=option_match(argc,argv,"-wd"))>-1)
00404     {
00405       if (on>argc-2 || argv[on+1][0] == '-')
00406       {
00407         cerr << "Invalid input data command line option"
00408                 " -- ignored" << endl;
00409       }
00410       else
00411       {
00412         tmpstring = adstring(argv[on+1]);
00413         wd_flag=1;
00414       }
00415     }
00416   }
00417   if (length(tmpstring))
00418   {
00419     if (tmpstring(length(tmpstring)) == directory_prefix)
00420     {
00421       adprogram_name=tmpstring + adprogram_name;
00422       working_directory_path = tmpstring;
00423     }
00424     else
00425     {
00426       adprogram_name=tmpstring + directory_prefix + adprogram_name;
00427       working_directory_path = tmpstring + directory_prefix;
00428     }
00429   }
00430 
00431   tmpstring=adprogram_name + adstring(".dat");
00432   if (argc > 1)
00433   {
00434     int on=0;
00435     if ( (on=option_match(argc,argv,"-ind"))>-1)
00436     {
00437       if (on>argc-2 || argv[on+1][0] == '-')
00438       {
00439         cerr << "Invalid input data command line option"
00440                 " -- ignored" << endl;
00441       }
00442       else
00443       {
00444         tmpstring = adstring(argv[on+1]);
00445       }
00446     }
00447   }
00448   global_datafile= new cifstream(tmpstring);
00449   if (!global_datafile)
00450   {
00451     cerr << "Error trying to open data input file "
00452          << tmpstring << endl;
00453   }
00454   else
00455   {
00456     if (!(*global_datafile))
00457     {
00458       cerr << "Error trying to open data input file "
00459            << tmpstring << endl;
00460       delete global_datafile;
00461       global_datafile=NULL;
00462     }
00463   }
00464   adstring ts=adprogram_name + adstring(".log");
00465   global_logfile= new ofstream( (char*)ts);
00466 
00467   int biopt=-1;
00468   int aiopt=-1;
00469   biopt=option_match(argc,argv,"-binp");
00470   aiopt=option_match(argc,argv,"-ainp");
00471 
00472   tmpstring=adprogram_name + adstring(".bin");
00473   if (!global_bparfile && aiopt == -1)
00474   {
00475     if (biopt>-1)
00476     {
00477       if (biopt>argc-2 || argv[biopt+1][0] == '-')
00478       {
00479         cerr << "Invalid input parameter file command line option"
00480                 " -- ignored" << endl;
00481       }
00482       else
00483       {
00484         tmpstring = adstring(argv[biopt+1]);
00485       }
00486     }
00487     global_bparfile= new uistream(tmpstring);
00488     if (global_bparfile)
00489     {
00490       if (!(*global_bparfile))
00491       {
00492         if (biopt>-1)
00493         {
00494           cerr << "Error trying to open binary inoput par file "
00495                << tmpstring << endl;
00496           exit(1);
00497         }
00498         delete global_bparfile;
00499         global_bparfile=NULL;
00500       }
00501     }
00502   }
00503   tmpstring=adprogram_name + adstring(".pin");
00504   if (!global_parfile)
00505   {
00506     if (aiopt>-1)
00507     {
00508       if (aiopt>argc-2 || argv[aiopt+1][0] == '-')
00509       {
00510         cerr << "Invalid input parameter file command line option"
00511                 " -- ignored" << endl;
00512       }
00513       else
00514       {
00515         tmpstring = adstring(argv[aiopt+1]);
00516       }
00517     }
00518     global_parfile= new cifstream(tmpstring);
00519     if (global_parfile)
00520     {
00521       if (!(*global_parfile))
00522       {
00523         if (aiopt>-1)
00524         {
00525           cerr << "Error trying to open ascii inoput par file "
00526                << tmpstring << endl;
00527           exit(1);
00528         }
00529         delete global_parfile;
00530         global_parfile=NULL;
00531       }
00532     }
00533   }
00534 }
00535 
00539 ad_comm::~ad_comm()
00540 {
00541   if (ptm)
00542   {
00543     delete ptm;
00544     ptm=0;
00545   }
00546   if (ptm1)
00547   {
00548     delete ptm1;
00549     ptm1=0;
00550   }
00551   if (global_datafile)
00552   {
00553     delete global_datafile;
00554     global_datafile=NULL;
00555   }
00556   if (global_parfile)
00557   {
00558     delete global_parfile;
00559     global_parfile=NULL;
00560   }
00561   if (global_logfile)
00562   {
00563     delete global_logfile;
00564     global_logfile=NULL;
00565   }
00566 }
00567 
00568 void function_minimizer::pre_userfunction(void)
00569 {
00570   if (lapprox)
00571   {
00572     if (lapprox->hesstype==2)
00573     {
00574       //lapprox->num_separable_calls=0;
00575       lapprox->separable_calls_counter=0;
00576     }
00577   }
00578   userfunction();
00579   if (lapprox)
00580   {
00581     if (lapprox->hesstype==2)
00582     {
00583       lapprox->num_separable_calls=lapprox->separable_calls_counter;
00584 
00585       double tmp=0.0;
00586       int inner_opt_value=inner_opt();
00587       if (lapprox->saddlepointflag==2)
00588       {
00589         if (inner_opt_value !=0 )
00590         {
00591           for (int i = 1; i <= lapprox->num_separable_calls; i++)
00592           {
00593             tmp-=(*lapprox->separable_function_difference)(i);
00594           }
00595           value(*objective_function_value::pobjfun)=tmp;
00596         }
00597       }
00598       else
00599       {
00600         for (int i = 1; i <= lapprox->num_separable_calls; i++)
00601         {
00602           tmp+=(*lapprox->separable_function_difference)(i);
00603         }
00604         value(*objective_function_value::pobjfun)=tmp;
00605       }
00606     }
00607   }
00608 }