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