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