ADMB Documentation  11.1.2397
 All Classes Files Functions Variables Typedefs Friends Defines
modspmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: modspmin.cpp 2299 2014-09-05 18:18:07Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 #  include <df1b2fun.h>
00009 #  include <adrndeff.h>
00010 
00011 //#define NO_MCMC
00012 #if ( (defined(_WINDOWS) || defined(_Windows)) && !defined(BORBUGS))
00013 #  include <windows.h>
00014 #endif
00015 
00016 void ADSleep(unsigned int x);
00017 
00018   void check_java_flags(int& start_flag,int& quit_flag,int& der_flag,
00019     int& next_flag);
00020   void vm_initialize(void);
00021   void test_mcmc_options_window(void);
00022   void ad_open_mcmc_options_window(void);
00023   void ad_open_mcmchist_window(void);
00024   void ad_update_mcmc_report(double * v,int l);
00025   void write_banner_stuff(void);
00026   void calljunk(double x){;}
00027 //#if defined (AD_DEMO)
00028   void adwait(double sec);
00029 //#endif
00030   int function_minimizer::have_constraints=0;
00031   int function_minimizer::first_hessian_flag=0;
00032   //int function_minimizer::in_likeprof_flag=0;
00033 
00034 class admb_javapointers;
00035 extern admb_javapointers * adjm_ptr;
00036 
00037   void function_minimizer::computations(int argc,char * argv[])
00038   {
00039     //traceflag=1;
00040     tracing_message(traceflag,"A1");
00041     //if (option_match(argc,argv,"-gui")>-1)
00042     //{
00043     //  vm_initialize();
00044     //  cout << " called vm_initialize() " << endl;
00045     //}
00046 #if defined (AD_DEMO)
00047      write_banner_stuff();
00048 #endif
00049     if (option_match(argc,argv,"-mceval") == -1)
00050     {
00051         computations1(argc,argv);
00052     }
00053     else
00054     {
00055       initial_params::mceval_phase=1;
00056       mcmc_eval();
00057       initial_params::mceval_phase=0;
00058     }
00059     other_calculations();
00060 
00061     final_calcs();
00062     // clean up if have random effects
00063      // cleanup_laplace_stuff(lapprox);
00064   }
00065 
00066   void function_minimizer::computations1(int argc,char * argv[])
00067   {
00068     tracing_message(traceflag,"B1");
00069 
00070     int on=-1;
00071     int nopt=-1;
00072 #if defined(USE_ADPVM)
00073     if (ad_comm::pvm_manager)
00074     {
00075       switch (ad_comm::pvm_manager->mode)
00076       {
00077       case 1: //master
00078         pvm_params::send_all_to_slaves();
00079         break;
00080       case 2: //slave
00081         pvm_params::get_all_from_master();
00082         break;
00083       default:
00084         cerr << "Illegal value for ad_comm::pvm_manager->mode"
00085          " value was " << ad_comm::pvm_manager->mode << endl;
00086         ad_exit(1);
00087       }
00088     }
00089 #endif  // #if defined(USE_ADPVM)
00090 
00091     set_runtime();
00092 
00093     if ( (on=option_match(argc,argv,"-hbf",nopt))>-1)
00094     {
00095       gradient_structure::Hybrid_bounded_flag=1;
00096     }
00097 
00098     // Sets the maximum number of function evaluation as determined from the
00099     // command line
00100     if ( (on=option_match(argc,argv,"-maxfn",nopt))>-1)
00101     {
00102       if (nopt ==1)
00103       {
00104         set_runtime_maxfn(argv[on+1]);
00105       }
00106       else
00107       {
00108         cerr << "Wrong number of options to -mafxn -- must be 1"
00109           " you have " << nopt << endl;
00110       }
00111     }
00112 
00113     if ( (on=option_match(argc,argv,"-ttr",nopt))>-1)
00114     {
00115       test_trust_flag=1;
00116     }
00117 
00118     if ( (on=option_match(argc,argv,"-crit",nopt))>-1)
00119     {
00120       if (nopt ==1)
00121       {
00122         set_runtime_crit(argv[on+1]);
00123       }
00124       else
00125       {
00126         cerr << "Wrong number of options to -crit -- must be 1"
00127           " you have " << nopt << endl;
00128       }
00129     }
00130 
00131     stddev_params::get_stddev_number_offset();
00132 
00133     tracing_message(traceflag,"C1");
00134 
00135     repeatminflag=0;
00136     do
00137     {
00138      /*
00139       if (spminflag)
00140       {
00141         repeatminflag=1;
00142         spminflag=0;
00143       }
00144       else
00145       {
00146         repeatminflag=0;
00147       }
00148       */
00149 
00150       if (option_match(argc,argv,"-noest") == -1)
00151       {
00152         if (!function_minimizer::have_constraints)
00153         {
00154           minimize();
00155         }
00156         else
00157         {
00158           constraints_minimize();
00159         }
00160       }
00161       else
00162       {
00163         initial_params::current_phase=initial_params::max_number_phases;
00164       }
00165       tracing_message(traceflag,"D1");
00166 
00167       //double ratio=100.*gradient_structure::max_last_offset/12000.0;
00168       tracing_message(traceflag,"E1");
00169       if (option_match(argc,argv,"-est") == -1)
00170       {
00171         if (!quit_flag)
00172         {
00173           // save the sparse Hessian for the random effects
00174           if (lapprox && lapprox->sparse_hessian_flag)
00175           {
00176             if (lapprox->sparse_triplet2)
00177             {
00178               dcompressed_triplet& dct=*(lapprox->sparse_triplet2);
00179               adstring tmpstring = ad_comm::adprogram_name + ".shess";
00180               uostream uos((char*)(tmpstring));
00181               uos << dct.get_n()  << dct.indexmin() << dct.indexmax()
00182                   << dct.get_coords() << dct.get_x();
00183             }
00184           }
00185 
00186           int on=option_match(argc,argv,"-nohess");
00187           int on1=option_match(argc,argv,"-noest");
00188           if (on==-1 && on1==-1)
00189           {
00190             if (option_match(argc,argv,"-sdonly")==-1)
00191             {
00192               hess_routine();
00193             }
00194             // set this flag so that variables only needed for their std devs
00195             // will be calculated
00196             initial_params::sd_phase=1;
00197 #if defined(USE_ADPVM)
00198             if (ad_comm::pvm_manager)
00199             {
00200               if (ad_comm::pvm_manager->mode==1)  //master
00201               {
00202                 depvars_routine();
00203                 hess_inv();
00204                 if (spminflag==0)
00205                 {
00206                   sd_routine();
00207                 }
00208               }
00209             }
00210             else
00211 #endif
00212             {
00213               depvars_routine();
00214               hess_inv();
00215               if (spminflag==0)
00216               {
00217                 sd_routine();
00218               }
00219             }
00220           }
00221           else
00222           {
00223             initial_params::sd_phase=1;
00224           }
00225           if (spminflag==0)
00226           {
00227             if ( (on=option_match(argc,argv,"-lprof"))>-1)
00228             {
00229               if (likeprof_params::num_likeprof_params)
00230               {
00231     #if defined(USE_ADPVM)
00232                 if (ad_comm::pvm_manager)
00233                 {
00234                   switch (ad_comm::pvm_manager->mode)
00235                   {
00236                   case 1: // master
00237                     likeprof_routine(ffbest);
00238                     break;
00239                   case 2: // slave
00240                     pvm_slave_likeprof_routine();
00241                     break;
00242                   default:
00243                     cerr << "error illega value for pvm_manager->mode" << endl;
00244                     exit(1);
00245                   }
00246                 }
00247                 else
00248     #endif
00249                 {
00250                   const double f = value(*objective_function_value::pobjfun);
00251                   likeprof_routine(f);
00252                 }
00253               }
00254             }
00255             int nopt=0;
00256             int on2=-1;
00257             int nopt2=-1;
00258 
00259             // stuff for mcmc
00260             //cout << "checking for mcmc" << endl;
00261             if ( (on=option_match(argc,argv,"-mcmc",nopt))>-1 ||
00262                  (on=option_match(argc,argv,"-mcmc2",nopt))>-1)
00263             {
00264               if ( (on2=option_match(argc,argv,"-mcmc2",nopt2))>-1)
00265                 mcmc2_flag=1;
00266               else
00267                 mcmc2_flag=0;
00268 
00269  #if defined(USE_ADPVM)
00270               if (ad_comm::pvm_manager)
00271               {
00272                 switch (ad_comm::pvm_manager->mode)
00273                 {
00274                 case 1: // master
00275                   pvm_master_mcmc_computations();
00276                 break;
00277                 case 2: // slave
00278                   pvm_slave_mcmc_routine();
00279                   break;
00280                 default:
00281                   cerr << "error illega value for pvm_manager->mode" << endl;
00282                   exit(1);
00283                 }
00284               }
00285               else
00286  #endif
00287               {
00288                 mcmc_computations();
00289               }
00290             }
00291             if ( (on=option_match(argc,argv,"-sob",nopt))>-1)
00292             {
00293               int nsob=0;
00294               //int iseed0=0;
00295               //double dscale=1.0;
00296               if (nopt)
00297               {
00298                 nsob=atoi(argv[on+1]);
00299                 if (nsob <=0)
00300                 {
00301                   cerr << " Invalid option following command line option -sob"
00302                           " -- "
00303                     << endl << " ignored" << endl;
00304                 }
00305               }
00306               if ( (on=option_match(argc,argv,"-mcr",nopt))>-1)
00307               {
00308     #if defined(NO_MCMC)
00309                 cerr << "mcmc and sob option not supported you must purchase"
00310                         " this as an add on" << endl;
00311                 exit(1);
00312     #else
00313                 //sob_routine(nsob,dscale,1);
00314                 //sobol_importance_routine(nsob,iseed0,dscale,1);
00315     #endif
00316               }
00317               else
00318               {
00319     #if defined(NO_MCMC)
00320                 cerr << "mcmc and sob option not supported you must purchase"
00321                         " this as an add on" << endl;
00322                 exit(1);
00323     #else
00324                 //sobol_importance_routine(nsob,iseed0,dscale,0);
00325     #endif
00326               }
00327             }
00328             initial_params::sd_phase=0;
00329           }
00330         }
00331       }
00332     }
00333     while(spminflag || repeatminflag);
00334   }
00335 
00336   void function_minimizer::computations(void)
00337   {
00338     // for now just do parameter estimates
00339     //function_minimizer::minimize();
00340     minimize();
00341   }
00342 
00343 void write_banner_stuff(void)
00344 {
00345   if (ad_printf)
00346   {
00347     char banner0[56]={"*****************************************************"};
00348     char banner1[56]={"This is the open source version of AD Model Builder"};
00349     char banner1a[58]={"You can freely use AD Model Builder"};
00350     char banner2[30]={"http://www.admb-project.org/"};
00351     char banner3[55]={"http://www.admb-project.org/"};
00352     char banner4[60]={"users@admb-project.org   http://www.admb-project.org/"};
00353     (*ad_printf)("\n%s\n", banner0);
00354     (*ad_printf)("%s\n", banner1);
00355     (*ad_printf)("%s\n", banner1a);
00356     (*ad_printf)("%s\n", banner2);
00357     (*ad_printf)("%s\n", banner3);
00358     (*ad_printf)("%s\n", banner4);
00359     //(*ad_printf)("%s\n", adcopy);
00360     (*ad_printf)("%s\n", banner0);
00361     (*ad_printf)("%s\n\n", banner0);
00362   }
00363   adwait(2.5);
00364 }
00365 
00366   void test_mcmc_options_window(void)
00367   {
00368     dvector v(1,1000);
00369     random_number_generator rng(908);
00370 
00371     for (int i=5;i<=1000;i++)
00372     {
00373       rng.reinitialize(908);
00374       v(1,i).fill_randn(rng);
00375       for (int j=2;j<=i;j++)
00376       {
00377         v(j)=0.9*v(j-1)+0.435889*v(j);
00378       }
00379 
00380       //ad_update_mcmc_report(&(v(1)),i);
00381       ADSleep(500);
00382     }
00383   }
00384 
00385   void function_minimizer::set_runtime(void){;}
00386 
00387   void function_minimizer::set_runtime_maxfn(const char * s)
00388   {
00389     adstring opt="{" + adstring(s) + "}";
00390     dvector temp1((char*)(opt));
00391     if (allocated(maximum_function_evaluations))
00392       maximum_function_evaluations.deallocate();
00393     maximum_function_evaluations.allocate(temp1.indexmin(),temp1.indexmax());
00394     maximum_function_evaluations=temp1;
00395   }
00396 
00397   void function_minimizer::set_runtime_crit(const char * s)
00398   {
00399     adstring opt="{" + adstring(s) + "}";
00400     dvector temp1((char*)(opt));
00401     if (allocated(convergence_criteria)) convergence_criteria.deallocate();
00402     convergence_criteria.allocate(temp1.indexmin(),temp1.indexmax());
00403     convergence_criteria=temp1;
00404   }
00405 
00406   void function_minimizer::mcmc_computations(void)
00407   {
00408     int ton,tnopt = 0;
00409     ton=option_match(ad_comm::argc,ad_comm::argv,"-mcmc",tnopt);
00410     if (ton<0)
00411     {
00412       ton=option_match(ad_comm::argc,ad_comm::argv,"-mcmc2",tnopt);
00413     }
00414     int on=ton;
00415     int nopt=tnopt;
00416 
00417     if (on>-1)
00418     {
00419       /*
00420       if (adjm_ptr)
00421       {
00422         ad_open_mcmc_options_window();
00423         ad_open_mcmchist_window();
00424         //test_mcmc_options_window();
00425       }
00426       */
00427       int nmcmc=0;
00428       int iseed0=0;
00429       double dscale=1.0;
00430       if (nopt)
00431       {
00432         nmcmc=(int)atof(ad_comm::argv[on+1]);
00433         if (nmcmc <=0)
00434         {
00435           cerr << " Invalid option following command line option -mcmc -- "
00436             << endl << " ignored" << endl;
00437         }
00438       }
00439       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmult",nopt))>-1)
00440       {
00441         if (nopt)
00442         {
00443           char * end;
00444           dscale=strtod(ad_comm::argv[on+1],&end);
00445           if (!dscale)
00446           {
00447             cerr << "Invalid argument to option -mcmult" << endl;
00448             dscale=1.0;
00449           }
00450         }
00451       }
00452       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcseed",nopt))>-1)
00453       {
00454         if (nopt)
00455         {
00456           iseed0=atoi(ad_comm::argv[on+1]);
00457           if (iseed0 <=0)
00458           {
00459             cerr << " Invalid option following command line option -mcseed -- "
00460               << endl << " ignored" << endl;
00461             iseed0=0;
00462           }
00463         }
00464       }
00465       int hybrid_flag=0;
00466       if (option_match(ad_comm::argc,ad_comm::argv,"-hybrid") > -1)
00467       {
00468         hybrid_flag=1;
00469         gradient_structure::Hybrid_bounded_flag=1;
00470       }
00471       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcr",nopt))>-1)
00472       {
00473         if (hybrid_flag==0)
00474         {
00475           mcmc_routine(nmcmc,iseed0,dscale,1);
00476         }
00477         else
00478         {
00479           hybrid_mcmc_routine(nmcmc,iseed0,dscale,1);
00480         }
00481       }
00482       else
00483       {
00484         if (hybrid_flag==0)
00485         {
00486           mcmc_routine(nmcmc,iseed0,dscale,0);
00487         }
00488         else
00489         {
00490           hybrid_mcmc_routine(nmcmc,iseed0,dscale,0);
00491         }
00492       }
00493     }
00494   }
00495 
00496   void function_minimizer::pvm_master_mcmc_computations(void)
00497   {
00498     int on,nopt = 0;
00499     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmc",nopt))>-1)
00500     {
00501       /*
00502       if (adjm_ptr)
00503       {
00504         ad_open_mcmc_options_window();
00505         ad_open_mcmchist_window();
00506         //test_mcmc_options_window();
00507       }
00508       */
00509       int nmcmc=0;
00510       int iseed0=0;
00511       double dscale=1.0;
00512       if (nopt)
00513       {
00514         nmcmc=(int)atof(ad_comm::argv[on+1]);
00515         if (nmcmc <=0)
00516         {
00517           cerr << " Invalid option following command line option -mcmc -- "
00518             << endl << " ignored" << endl;
00519         }
00520       }
00521       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmult",nopt))>-1)
00522       {
00523         if (nopt)
00524         {
00525           char * end;
00526           dscale=strtod(ad_comm::argv[on+1],&end);
00527           if (!dscale)
00528           {
00529             cerr << "Invalid argument to option -mcmult" << endl;
00530             dscale=1.0;
00531           }
00532         }
00533       }
00534       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcseed",nopt))>-1)
00535       {
00536         if (nopt)
00537         {
00538           iseed0=atoi(ad_comm::argv[on+1]);
00539           if (iseed0 <=0)
00540           {
00541             cerr << " Invalid option following command line option -mcseed -- "
00542               << endl << " ignored" << endl;
00543             iseed0=0;
00544           }
00545         }
00546       }
00547       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcr",nopt))>-1)
00548       {
00549 #if defined(NO_MCMC)
00550         cerr << "mcmc option not supported you must purchase"
00551                 " this as an add on" << endl;
00552         exit(1);
00553 #else
00554         //mcmc_routine(nmcmc,iseed0,dscale,1);
00555         #if defined(USE_ADPVM)
00556         pvm_master_mcmc_routine(nmcmc,iseed0,dscale,1);
00557         #endif
00558 #endif
00559       }
00560       else
00561       {
00562 #if defined(NO_MCMC)
00563         cerr << "mcmc option not supported you must purchase"
00564                 " this as an add on" << endl;
00565         exit(1);
00566 #else
00567         //mcmc_routine(nmcmc,iseed0,dscale,0);
00568         #if defined(USE_ADPVM)
00569         pvm_master_mcmc_routine(nmcmc,iseed0,dscale,0);
00570         #endif
00571 #endif
00572       }
00573     }
00574   }
00575   void function_minimizer::pvm_slave_mcmc_computations(void)
00576   {
00577     pvm_slave_mcmc_routine();
00578   }