ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
modspmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: modspmin.cpp 1971 2014-05-02 18:45:14Z 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(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 #if defined(ADMB_REDEMO)
00136   if (return_int(*pflag5) != 1912)
00137   {
00138     for (int i=0;i<initial_params::num_initial_params;i++)
00139     {
00140       initial_params::varsptr[i]=0;
00141     }
00142   }
00143 
00144 #endif // #if defined(ADMB_REDEMO)
00145 
00146     repeatminflag=0;
00147     do
00148     {
00149      /*
00150       if (spminflag)
00151       {
00152         repeatminflag=1;
00153         spminflag=0;
00154       }
00155       else
00156       {
00157         repeatminflag=0;
00158       }
00159       */
00160 
00161       if (option_match(argc,argv,"-noest") == -1)
00162       {
00163         if (!function_minimizer::have_constraints)
00164         {
00165           minimize();
00166         }
00167         else
00168         {
00169           constraints_minimize();
00170         }
00171       }
00172       else
00173       {
00174         initial_params::current_phase=initial_params::max_number_phases;
00175       }
00176       tracing_message(traceflag,"D1");
00177 
00178       //double ratio=100.*gradient_structure::max_last_offset/12000.0;
00179       tracing_message(traceflag,"E1");
00180       if (option_match(argc,argv,"-est") == -1)
00181       {
00182         if (!quit_flag)
00183         {
00184           // save the sparse Hessian for the random effects
00185           if (lapprox && lapprox->sparse_hessian_flag)
00186           {
00187             if (lapprox->sparse_triplet2)
00188             {
00189               dcompressed_triplet& dct=*(lapprox->sparse_triplet2);
00190               adstring tmpstring = ad_comm::adprogram_name + ".shess";
00191               uostream uos((char*)(tmpstring));
00192               uos << dct.get_n()  << dct.indexmin() << dct.indexmax()
00193                   << dct.get_coords() << dct.get_x();
00194             }
00195           }
00196 
00197           int on=-1;
00198           int on1=-1;
00199           on=option_match(argc,argv,"-nohess");
00200           on1=option_match(argc,argv,"-noest");
00201           if (on==-1 && on1==-1)
00202           {
00203             if (option_match(argc,argv,"-sdonly")==-1)
00204             {
00205               hess_routine();
00206             }
00207             // set this flag so that variables only needed for their std devs
00208             // will be calculated
00209             initial_params::sd_phase=1;
00210 #if defined(USE_ADPVM)
00211             if (ad_comm::pvm_manager)
00212             {
00213               if (ad_comm::pvm_manager->mode==1)  //master
00214               {
00215                 depvars_routine();
00216                 hess_inv();
00217                 if (spminflag==0)
00218                 {
00219                   sd_routine();
00220                 }
00221               }
00222             }
00223             else
00224 #endif
00225             {
00226               depvars_routine();
00227               hess_inv();
00228               if (spminflag==0)
00229               {
00230                 sd_routine();
00231               }
00232             }
00233           }
00234           else
00235           {
00236             initial_params::sd_phase=1;
00237           }
00238           if (spminflag==0)
00239           {
00240             if ( (on=option_match(argc,argv,"-lprof"))>-1)
00241             {
00242               if (likeprof_params::num_likeprof_params)
00243               {
00244     #if defined(USE_ADPVM)
00245                 if (ad_comm::pvm_manager)
00246                 {
00247                   switch (ad_comm::pvm_manager->mode)
00248                   {
00249                   case 1: // master
00250                     likeprof_routine(ffbest);
00251                     break;
00252                   case 2: // slave
00253                     pvm_slave_likeprof_routine();
00254                     break;
00255                   default:
00256                     cerr << "error illega value for pvm_manager->mode" << endl;
00257                     exit(1);
00258                   }
00259                 }
00260                 else
00261     #endif
00262                 {
00263                   likeprof_routine(ffbest);
00264                 }
00265               }
00266             }
00267             int nopt=0;
00268             int on2=-1;
00269             int nopt2=-1;
00270 
00271             // stuff for mcmc
00272             //cout << "checking for mcmc" << endl;
00273             if ( (on=option_match(argc,argv,"-mcmc",nopt))>-1 ||
00274                  (on=option_match(argc,argv,"-mcmc2",nopt))>-1)
00275             {
00276               if ( (on2=option_match(argc,argv,"-mcmc2",nopt2))>-1)
00277                 mcmc2_flag=1;
00278               else
00279                 mcmc2_flag=0;
00280 
00281  #if defined(USE_ADPVM)
00282               if (ad_comm::pvm_manager)
00283               {
00284                 switch (ad_comm::pvm_manager->mode)
00285                 {
00286                 case 1: // master
00287                   pvm_master_mcmc_computations();
00288                 break;
00289                 case 2: // slave
00290                   pvm_slave_mcmc_routine();
00291                   break;
00292                 default:
00293                   cerr << "error illega value for pvm_manager->mode" << endl;
00294                   exit(1);
00295                 }
00296               }
00297               else
00298  #endif
00299               {
00300                 mcmc_computations();
00301               }
00302             }
00303             if ( (on=option_match(argc,argv,"-sob",nopt))>-1)
00304             {
00305               int nsob=0;
00306               //int iseed0=0;
00307               //double dscale=1.0;
00308               if (nopt)
00309               {
00310                 nsob=atoi(argv[on+1]);
00311                 if (nsob <=0)
00312                 {
00313                   cerr << " Invalid option following command line option -sob"
00314                           " -- "
00315                     << endl << " ignored" << endl;
00316                 }
00317               }
00318               if ( (on=option_match(argc,argv,"-mcr",nopt))>-1)
00319               {
00320     #if defined(NO_MCMC)
00321                 cerr << "mcmc and sob option not supported you must purchase"
00322                         " this as an add on" << endl;
00323                 exit(1);
00324     #else
00325                 //sob_routine(nsob,dscale,1);
00326                 //sobol_importance_routine(nsob,iseed0,dscale,1);
00327     #endif
00328               }
00329               else
00330               {
00331     #if defined(NO_MCMC)
00332                 cerr << "mcmc and sob option not supported you must purchase"
00333                         " this as an add on" << endl;
00334                 exit(1);
00335     #else
00336                 //sobol_importance_routine(nsob,iseed0,dscale,0);
00337     #endif
00338               }
00339             }
00340             initial_params::sd_phase=0;
00341           }
00342         }
00343       }
00344     }
00345     while(spminflag || repeatminflag);
00346   }
00347 
00348   void function_minimizer::computations(void)
00349   {
00350     // for now just do parameter estimates
00351     //function_minimizer::minimize();
00352     minimize();
00353   }
00354 
00355 void write_banner_stuff(void)
00356 {
00357   if (ad_printf)
00358   {
00359     char banner0[56]={"*****************************************************"};
00360     char banner1[56]={"This is the open source version of AD Model Builder"};
00361     char banner1a[58]={"You can freely use AD Model Builder"};
00362     char banner2[30]={"http://www.admb-project.org/"};
00363     char banner3[55]={"http://www.admb-project.org/"};
00364     char banner4[60]={"users@admb-project.org   http://www.admb-project.org/"};
00365     (*ad_printf)("\n%s\n", banner0);
00366     (*ad_printf)("%s\n", banner1);
00367     (*ad_printf)("%s\n", banner1a);
00368     (*ad_printf)("%s\n", banner2);
00369     (*ad_printf)("%s\n", banner3);
00370     (*ad_printf)("%s\n", banner4);
00371     //(*ad_printf)("%s\n", adcopy);
00372     (*ad_printf)("%s\n", banner0);
00373     (*ad_printf)("%s\n\n", banner0);
00374   }
00375   adwait(2.5);
00376 }
00377 
00378   void test_mcmc_options_window(void)
00379   {
00380     dvector v(1,1000);
00381     random_number_generator rng(908);
00382 
00383     for (int i=5;i<=1000;i++)
00384     {
00385       rng.reinitialize(908);
00386       v(1,i).fill_randn(rng);
00387       for (int j=2;j<=i;j++)
00388       {
00389         v(j)=0.9*v(j-1)+0.435889*v(j);
00390       }
00391 
00392       //ad_update_mcmc_report(&(v(1)),i);
00393       ADSleep(500);
00394     }
00395   }
00396 
00397   void function_minimizer::set_runtime(void){;}
00398 
00399   void function_minimizer::set_runtime_maxfn(const char * s)
00400   {
00401     adstring opt="{" + adstring(s) + "}";
00402     dvector temp1((char*)(opt));
00403     if (allocated(maximum_function_evaluations))
00404       maximum_function_evaluations.deallocate();
00405     maximum_function_evaluations.allocate(temp1.indexmin(),temp1.indexmax());
00406     maximum_function_evaluations=temp1;
00407   }
00408 
00409   void function_minimizer::set_runtime_crit(const char * s)
00410   {
00411     adstring opt="{" + adstring(s) + "}";
00412     dvector temp1((char*)(opt));
00413     if (allocated(convergence_criteria)) convergence_criteria.deallocate();
00414     convergence_criteria.allocate(temp1.indexmin(),temp1.indexmax());
00415     convergence_criteria=temp1;
00416   }
00417 
00418   void function_minimizer::mcmc_computations(void)
00419   {
00420     int ton,tnopt = 0;
00421     ton=option_match(ad_comm::argc,ad_comm::argv,"-mcmc",tnopt);
00422     if (ton<0)
00423     {
00424       ton=option_match(ad_comm::argc,ad_comm::argv,"-mcmc2",tnopt);
00425     }
00426     int on=ton;
00427     int nopt=tnopt;
00428 
00429     if (on>-1)
00430     {
00431       /*
00432       if (adjm_ptr)
00433       {
00434         ad_open_mcmc_options_window();
00435         ad_open_mcmchist_window();
00436         //test_mcmc_options_window();
00437       }
00438       */
00439       int nmcmc=0;
00440       int iseed0=0;
00441       double dscale=1.0;
00442       if (nopt)
00443       {
00444         nmcmc=(int)atof(ad_comm::argv[on+1]);
00445         if (nmcmc <=0)
00446         {
00447           cerr << " Invalid option following command line option -mcmc -- "
00448             << endl << " ignored" << endl;
00449         }
00450       }
00451       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmult",nopt))>-1)
00452       {
00453         if (nopt)
00454         {
00455           char * end;
00456           dscale=strtod(ad_comm::argv[on+1],&end);
00457           if (!dscale)
00458           {
00459             cerr << "Invalid argument to option -mcmult" << endl;
00460             dscale=1.0;
00461           }
00462         }
00463       }
00464       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcseed",nopt))>-1)
00465       {
00466         if (nopt)
00467         {
00468           iseed0=atoi(ad_comm::argv[on+1]);
00469           if (iseed0 <=0)
00470           {
00471             cerr << " Invalid option following command line option -mcseed -- "
00472               << endl << " ignored" << endl;
00473             iseed0=0;
00474           }
00475         }
00476       }
00477       int hybrid_flag=0;
00478       if (option_match(ad_comm::argc,ad_comm::argv,"-hybrid") > -1)
00479       {
00480         hybrid_flag=1;
00481         gradient_structure::Hybrid_bounded_flag=1;
00482       }
00483       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcr",nopt))>-1)
00484       {
00485         if (hybrid_flag==0)
00486         {
00487           mcmc_routine(nmcmc,iseed0,dscale,1);
00488         }
00489         else
00490         {
00491           hybrid_mcmc_routine(nmcmc,iseed0,dscale,1);
00492         }
00493       }
00494       else
00495       {
00496         if (hybrid_flag==0)
00497         {
00498           mcmc_routine(nmcmc,iseed0,dscale,0);
00499         }
00500         else
00501         {
00502           hybrid_mcmc_routine(nmcmc,iseed0,dscale,0);
00503         }
00504       }
00505     }
00506   }
00507 
00508   void function_minimizer::pvm_master_mcmc_computations(void)
00509   {
00510     int on,nopt = 0;
00511     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmc",nopt))>-1)
00512     {
00513       /*
00514       if (adjm_ptr)
00515       {
00516         ad_open_mcmc_options_window();
00517         ad_open_mcmchist_window();
00518         //test_mcmc_options_window();
00519       }
00520       */
00521       int nmcmc=0;
00522       int iseed0=0;
00523       double dscale=1.0;
00524       if (nopt)
00525       {
00526         nmcmc=(int)atof(ad_comm::argv[on+1]);
00527         if (nmcmc <=0)
00528         {
00529           cerr << " Invalid option following command line option -mcmc -- "
00530             << endl << " ignored" << endl;
00531         }
00532       }
00533       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcmult",nopt))>-1)
00534       {
00535         if (nopt)
00536         {
00537           char * end;
00538           dscale=strtod(ad_comm::argv[on+1],&end);
00539           if (!dscale)
00540           {
00541             cerr << "Invalid argument to option -mcmult" << endl;
00542             dscale=1.0;
00543           }
00544         }
00545       }
00546       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcseed",nopt))>-1)
00547       {
00548         if (nopt)
00549         {
00550           iseed0=atoi(ad_comm::argv[on+1]);
00551           if (iseed0 <=0)
00552           {
00553             cerr << " Invalid option following command line option -mcseed -- "
00554               << endl << " ignored" << endl;
00555             iseed0=0;
00556           }
00557         }
00558       }
00559       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcr",nopt))>-1)
00560       {
00561 #if defined(NO_MCMC)
00562         cerr << "mcmc option not supported you must purchase"
00563                 " this as an add on" << endl;
00564         exit(1);
00565 #else
00566         //mcmc_routine(nmcmc,iseed0,dscale,1);
00567         #if defined(USE_ADPVM)
00568         pvm_master_mcmc_routine(nmcmc,iseed0,dscale,1);
00569         #endif
00570 #endif
00571       }
00572       else
00573       {
00574 #if defined(NO_MCMC)
00575         cerr << "mcmc option not supported you must purchase"
00576                 " this as an add on" << endl;
00577         exit(1);
00578 #else
00579         //mcmc_routine(nmcmc,iseed0,dscale,0);
00580         #if defined(USE_ADPVM)
00581         pvm_master_mcmc_routine(nmcmc,iseed0,dscale,0);
00582         #endif
00583 #endif
00584       }
00585     }
00586   }
00587   void function_minimizer::pvm_slave_mcmc_computations(void)
00588   {
00589     pvm_slave_mcmc_routine();
00590   }