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