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