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