ADMB Documentation  11.1.1810
 All Classes Files Functions Variables Typedefs Friends Defines
xxmcmc2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xxmcmc2.cpp 1807 2014-03-27 21:40:03Z jianelli $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 //ofstream tmpof("testmc");
00010 
00011 #if defined(__GNU__) || defined(UNIXKLUDGE) || defined(__SUN__)
00012   #define getch getchar
00013 #endif
00014 
00015 #ifdef __GNUDOS__
00016   #include <gccmanip.h>
00017 #endif
00018 
00019 #if defined (__ZTC__) || defined(__TURBOC__) || defined(__WAT32__) \
00020   || defined (_MSC_VER)
00021 #  if  !defined(__linux__)
00022 #    include <conio.h>
00023 #  endif
00024 #endif
00025 
00026 //double better_rand(long int&);
00027 void store_mcmc_values(const ofstream& ofs);
00028 void set_labels_for_mcmc(void);
00029 void save_mcmc_for_gui(const dvector& mcmc_values,dmatrix &mdm,int& ids);
00030 void save_mcmc_for_gui1(const dvector& mcmc_values,
00031   dmatrix &mdm,int& ids,int& iwrap,ivector& no);
00032 
00033 void check_java_flags(int& start_flag,int& quit_flag,int& der_flag,
00034   int& next_flag);
00035 void print_hist_data(const dmatrix& hist, const dmatrix& values,
00036                      const dvector& h, dvector& m, const dvector& s,
00037                      const dvector& parsave, long int iseed, double size_scale);
00038 
00039 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS,double s,
00040   int flag=0);
00041 
00042 int minnz(const dvector& x);
00043 int maxnz(const dvector& xa);
00044 
00045 int read_hist_data(const dmatrix& hist, const dvector& h,
00046   dvector& m, const dvector& s, const dvector& parsave,long int& iseed,
00047   const double& size_scale);
00048 
00049 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
00050                            const dmatrix& values, dmatrix& hist, 
00051                            const dvector& h,
00052                            int slots,double total_spread,int probflag=0);
00053 
00054 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& hist,
00055   dvector& mcmc_values,double llc, const dvector& h,int nslots,
00056   double total_spreadd,int probflag=0);
00057 
00058 void add_guihist_values(const dvector& s, const dvector& m,
00059   const dmatrix& _hist, dvector& mcmcnumber_values, double llc,
00060   const dvector& h,int nslots,double total_spread);
00061 
00062 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00063                                        const dmatrix& s_covar,
00064                                        adstring& prog_name);
00065 
00066 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00067                                       const adstring& prog_name);
00068 
00069 
00070 void read_hessian_matrix_and_scale(int nvar, const dmatrix& S,
00071                                    const dvector& pen_vector);
00072 
00073 int user_stop(void);
00074 
00075 extern int ctlc_flag;
00076 class admb_javapointers;
00077 extern admb_javapointers * adjm_ptr;
00078 
00079 dvector new_probing_bounded_multivariate_normal(int nvar, const dvector& a1,
00080                                                 const dvector& b1,
00081   dmatrix& ch, const double& wght,double pprobe, random_number_generator& rng);
00082 
00083 void new_probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00084                                                   const dvector& b1,
00085   dmatrix& ch, const double& wght, const dvector& _y,double pprobe,
00086   random_number_generator& rng);
00087 
00088 //void  newton_raftery_bayes_estimate(double cbf,int ic, 
00089 // const dvector& lk,double d);
00090 void  newton_raftery_bayes_estimate_new(double cbf,int ic, const dvector& lk,
00091                                         double d);
00092 
00093 void ad_update_mcmc_stats_report
00094   (int feval,int iter,double fval,double gmax);
00095 
00096 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00097   double fval,double gmax,const char * cbuf);
00098 void ad_update_mcmc_report(dmatrix& m,int i,int j,int ff=0);
00099 void ad_update_mcmchist_report(dmatrix& mcmc_values,ivector& number_offsets,
00100   dvector& mean_mcmc_values,dvector& h,int ff=0);
00101 
00102 void ADSleep(int);
00103 
00104 #if !defined(USE_ADPVM)
00105 void function_minimizer::pvm_master_mcmc_routine(int nmcmc,int iseed0,
00106   double dscale,int restart_flag){;}
00107 #else
00108 void function_minimizer::pvm_master_mcmc_routine(int nmcmc,int iseed0,
00109                                               double dscale, int restart_flag)
00110 {
00111   uostream * pofs_psave=NULL;
00112   dmatrix mcmc_display_matrix;
00113   int mcmc_save_index=1;
00114   int mcmc_wrap_flag=0;
00115   int mcmc_gui_length=10000;
00116   int no_sd_mcmc=1; // Ianelli changed to make default behavior have sdmcmc off
00117 
00118   int on2=-1;
00119   if ( (on2=option_match(ad_comm::argc,ad_comm::argv,"-sdmcmc"))>-1)
00120     no_sd_mcmc=0;
00121 
00122   if (stddev_params::num_stddev_params==0)
00123   {
00124     cerr << " You must declare at least one object of type sdreport "
00125          << endl << " to do the mcmc calculations" << endl;
00126      return;
00127   }
00128   {
00129     //ofstream of_bf("testbf");
00130 
00131     //if (adjm_ptr) set_labels_for_mcmc();
00132 
00133     ivector number_offsets;
00134     dvector lkvector;
00135     //double current_bf=0;
00136     double lcurrent_bf=0;
00137     double size_scale=1.0;
00138     double total_spread=200;
00139     //double total_spread=2500;
00140     uostream * pofs_sd = NULL;
00141     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00142     int scov_option=0;
00143     dmatrix s_covar;
00144     dvector s_mean;
00145     int on=-1;
00146     int ncsim=25000;
00147     int nslots=800;
00148     //int nslots=3600;
00149     int initial_nsim=4800;
00150     int ntmp=0;
00151     int ncor=0;
00152     double bfsum=0;
00153     int ibfcount=0;
00154     double llbest;
00155     double lbmax;
00156 
00157     //if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcscov",ntmp))>-1)
00158     //{
00159     scov_option=1;
00160     s_covar.allocate(1,nvar,1,nvar);
00161     s_mean.allocate(1,nvar);
00162     s_mean.initialize();
00163     s_covar.initialize();
00164 
00165     int ndvar=stddev_params::num_stddev_calc();
00166     int numdvar=stddev_params::num_stddev_number_calc();
00167     /*
00168     if (adjm_ptr)
00169     {
00170       mcmc_display_matrix.allocate(1,numdvar,1,mcmc_gui_length);
00171       number_offsets.allocate(1,numdvar);
00172       number_offsets=stddev_params::copy_all_number_offsets();
00173     }
00174     */
00175     dvector x(1,nvar);
00176     dvector scale(1,nvar);
00177     dmatrix values;
00178     int have_hist_flag=0;
00179     initial_params::xinit(x);
00180     dvector pen_vector(1,nvar);
00181     {
00182       initial_params::reset(dvar_vector(x),pen_vector);
00183       cout << pen_vector << endl << endl;
00184     }
00185 
00186     initial_params::mc_phase=0;
00187     initial_params::stddev_scale(scale,x);
00188     initial_params::mc_phase=1;
00189     dvector bmn(1,nvar);
00190     dvector mean_mcmc_values(1,ndvar);
00191     dvector s(1,ndvar);
00192     dvector h(1,ndvar);
00193     //dvector h;
00194     dvector square_mcmc_values(1,ndvar);
00195     square_mcmc_values.initialize();
00196     mean_mcmc_values.initialize();
00197     bmn.initialize();
00198     int use_empirical_flag=0;
00199     int diag_option=0;
00200     int topt=0;
00201     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcdiag"))>-1)
00202     {
00203       diag_option=1;
00204       cout << " Setting covariance matrix to diagonal with entries " << dscale
00205            << endl;
00206     }
00207     dmatrix S(1,nvar,1,nvar);
00208     dvector sscale(1,nvar);
00209     if (!diag_option)
00210     {
00211       int on,nopt;
00212       int rescale_bounded_flag=0;
00213       double rescale_bounded_power=0.5;
00214       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcrb",nopt))>-1)
00215       {
00216         if (nopt)
00217         {
00218           int iii=atoi(ad_comm::argv[on+1]);
00219           if (iii < 1 || iii > 9)
00220           {
00221             cerr << " -mcrb argument must be integer between 1 and 9 --"
00222                     " using default of 5" << endl;
00223             rescale_bounded_power=0.5;
00224           }
00225           else
00226             rescale_bounded_power=iii/10.0;
00227         }
00228         else
00229         {
00230           rescale_bounded_power=0.5;
00231         }
00232         rescale_bounded_flag=1;
00233       }
00234       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcec"))>-1)
00235       {
00236         use_empirical_flag=1;
00237       }
00238       if (use_empirical_flag)
00239       {
00240         read_empirical_covariance_matrix(nvar,S,ad_comm::adprogram_name);
00241       }
00242       else if (!rescale_bounded_flag)
00243       {
00244         int tmp;
00245         read_covariance_matrix(S,nvar,tmp,sscale);
00246       }
00247       else
00248       {
00249         read_hessian_matrix_and_scale1(nvar,S,rescale_bounded_power,
00250           mcmc2_flag);
00251         //read_hessian_matrix_and_scale(nvar,S,pen_vector);
00252       }
00253 
00254       {  // scale covariance matrix for model space
00255         dmatrix tmp(1,nvar,1,nvar);
00256         for (int i=1;i<=nvar;i++)
00257         {
00258           tmp(i,i)=S(i,i)*(scale(i)*scale(i));
00259           for (int j=1;j<i;j++)
00260           {
00261             tmp(i,j)=S(i,j)*(scale(i)*scale(j));
00262             tmp(j,i)=tmp(i,j);
00263           }
00264         }
00265         S=tmp;
00266       }
00267     }
00268     else
00269     {
00270       S.initialize();
00271       for (int i=1;i<=nvar;i++)
00272       {
00273         S(i,i)=dscale;
00274       }
00275     }
00276 
00277     cout << sort(eigenvalues(S)) << endl;
00278     dmatrix chd = choleski_decomp( (dscale*2.4/sqrt(double(nvar))) * S);
00279     dmatrix chdinv=inv(chd);
00280     int sgn;
00281 
00282     dmatrix symbds(1,2,1,nvar);
00283     initial_params::set_all_simulation_bounds(symbds);
00284     ofstream ofs_sd1((char*)(ad_comm::adprogram_name + adstring(".mc2")));
00285 
00286     {
00287       long int iseed=0;
00288       int number_sims;
00289       if (nmcmc<=0)
00290       {
00291         number_sims=  100000;
00292       }
00293       else
00294       {
00295         number_sims=  nmcmc;
00296       }
00297       //cin >> iseed;
00298       if (iseed0<=0)
00299       {
00300         iseed=-36519;
00301       }
00302       else
00303       {
00304         iseed=-iseed0;
00305       }
00306       if (iseed>0)
00307       {
00308         iseed=-iseed;
00309       }
00310       cout << "Initial seed value " << iseed << endl;
00311       random_number_generator rng(iseed);
00312       rng.better_rand();
00313       //better_rand(iseed);
00314       double lprob=0.0;
00315       double lpinv=0.0;
00316       double lprob3=0.0;
00317       // get lower and upper bounds
00318 
00319       independent_variables y(1,nvar);
00320       independent_variables parsave(1,nvar);
00321 
00322       // read in the mcmc values to date
00323       int ii=1;
00324       dmatrix hist;
00325       if (restart_flag)
00326       {
00327         int tmp=0;
00328         if (!no_sd_mcmc) {
00329           hist.allocate(1,ndvar,-nslots,nslots);
00330           tmp=read_hist_data(hist,h,mean_mcmc_values,s,parsave,iseed,
00331             size_scale);
00332           values.allocate(1,ndvar,-nslots,nslots);
00333           for (int i=1;i<=ndvar;i++)
00334           {
00335             values(i).fill_seqadd(mean_mcmc_values(i)-0.5*total_spread*s(i)
00336               +.5*h(i),h(i));
00337           }
00338         }
00339         if (iseed>0)
00340         {
00341           iseed=-iseed;
00342         }
00343         double br=rng.better_rand();
00344         if (tmp) have_hist_flag=1;
00345         chd=size_scale*chd;
00346         chdinv=chdinv/size_scale;
00347       }
00348       else
00349       {
00350         int on=-1;
00351         int nopt=0;
00352         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcpin",nopt))>-1)
00353         {
00354           if (nopt)
00355           {
00356             cifstream cif((char *)ad_comm::argv[on+1]);
00357             if (!cif)
00358             {
00359               cerr << "Error trying to open mcmc par input file "
00360                    << ad_comm::argv[on+1] << endl;
00361               exit(1);
00362             }
00363             cif >> parsave;
00364             if (!cif)
00365             {
00366               cerr << "Error reading from mcmc par input file "
00367                    << ad_comm::argv[on+1] << endl;
00368               exit(1);
00369             }
00370           }
00371           else
00372           {
00373             cerr << "Illegal option with -mcpin" << endl;
00374           }
00375         }
00376         else
00377         {
00378           ii=1;
00379           initial_params::copy_all_values(parsave,ii);
00380         }
00381       }
00382 
00383       ii=1;
00384       initial_params::restore_all_values(parsave,ii);
00385 
00386       gradient_structure::set_NO_DERIVATIVES();
00387       ofstream ogs("sims");
00388       ogs << nvar << " " << number_sims << endl;
00389       initial_params::xinit(y);
00390 
00391       send_int_to_slaves(1);
00392       double llc=-pvm_master_get_monte_carlo_value(nvar,y);
00393       send_int_to_slaves(1);
00394       llbest=-pvm_master_get_monte_carlo_value(nvar,y);
00395 
00396 
00397 
00398       lbmax=llbest;
00399       // store current mcmc variable values in param_values
00400       //dmatrix store_mcmc_values(1,number_sims,1,ndvar);
00401 #if defined(USE_BAYES_FACTORS)
00402       lkvector.allocate(1,number_sims);
00403 #endif
00404       dvector mcmc_values(1,ndvar);
00405       dvector mcmc_number_values;
00406       //if (adjm_ptr) mcmc_number_values.allocate(1,numdvar);
00407       int offs=1;
00408       stddev_params::copy_all_values(mcmc_values,offs);
00409 
00410       /*
00411       if (adjm_ptr)
00412       {
00413         offs=1;
00414         stddev_params::copy_all_number_values(mcmc_number_values,offs);
00415       }
00416       */
00417       int change_ball=2500;
00418       int nopt;
00419       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcscale",nopt))>-1)
00420       {
00421         if (nopt)
00422         {
00423           int iii=atoi(ad_comm::argv[on+1]);
00424           if (iii <=0)
00425           {
00426             cerr << " Invalid option following command line option -mcball -- "
00427               << endl << " ignored" << endl;
00428           }
00429           else
00430             change_ball=iii;
00431         }
00432       }
00433       int iac=0;
00434       int liac=0;
00435       int isim=0;
00436       int itmp=0;
00437       double logr;
00438       int u_option=0;
00439       double ll;
00440       int s_option=1;
00441       int psvflag=0;
00442       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcu"))>-1)
00443       {
00444         u_option=1;
00445       }
00446       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcnoscale"))>-1)
00447       {
00448         s_option=0;
00449       }
00450       //cout << llc << " " << llc << endl;
00451       int iac_old=0;
00452       int i_old=0;
00453 
00454       {
00455        if (!restart_flag)
00456        {
00457          pofs_sd =
00458            new uostream((char*)(ad_comm::adprogram_name + adstring(".mcm")));
00459        }
00460 
00461       int mcsave_flag=0;
00462       int mcrestart_flag=option_match(ad_comm::argc,ad_comm::argv,"-mcr");
00463 
00464       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcsave"))>-1)
00465       {
00466         int jj=(int)atof(ad_comm::argv[on+1]);
00467         if (jj <=0)
00468         {
00469           cerr << " Invalid option following command line option -mcsave -- "
00470             << endl;
00471         }
00472         else
00473         {
00474           mcsave_flag=jj;
00475           if ( mcrestart_flag>-1)
00476           {
00477             // check that nvar is correct
00478             {
00479               uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00480               if (!uis)
00481               {
00482                 cerr << "Error trying to open file" <<
00483                   ad_comm::adprogram_name + adstring(".psv") <<
00484                   " for mcrestart" <<   endl;
00485                 cerr << " I am starting a new file " << endl;
00486                 psvflag=1;
00487               }
00488               else
00489               {
00490                 int nv1;
00491                 uis >> nv1;
00492                 if (nv1 !=nvar)
00493                 {
00494                   cerr << "wrong number of independent variables in" <<
00495                     ad_comm::adprogram_name + adstring(".psv") <<
00496                   cerr << " I am starting a new file " << endl;
00497                   psvflag=1;
00498                 }
00499               }
00500             }
00501 
00502             if (!psvflag) {
00503               pofs_psave=
00504                 new uostream(
00505                   (char*)(ad_comm::adprogram_name + adstring(".psv")),ios::app);
00506             } else {
00507               pofs_psave=
00508                 new uostream((char*)(ad_comm::adprogram_name +
00509                              adstring(".psv")));
00510             }
00511 
00512           } else {
00513             pofs_psave=
00514               new uostream((char*)(ad_comm::adprogram_name + adstring(".psv")));
00515           }
00516           if (!pofs_psave|| !(*pofs_psave))
00517           {
00518             cerr << "Error trying to open file" <<
00519               ad_comm::adprogram_name + adstring(".psv") << endl;
00520             mcsave_flag=0;
00521             if (pofs_psave)
00522             {
00523               delete pofs_psave;
00524               pofs_psave=NULL;
00525             }
00526           }
00527           else
00528           {
00529             if (psvflag || (mcrestart_flag == -1) )
00530               (*pofs_psave) << nvar;
00531           }
00532         }
00533       }
00534 
00535       double pprobe=0.05;
00536       int probe_flag=0;
00537       int nopt=0;
00538       on = option_match(ad_comm::argc, ad_comm::argv, "-mcprobe", nopt);
00539       if (on == -1)
00540       {
00541         on = option_match(ad_comm::argc,ad_comm::argv,"-mcgrope",nopt);
00542       }
00543       if (on > -1)
00544       {
00545         probe_flag=1;
00546         if (nopt)
00547         {
00548           char * end;
00549           pprobe=strtod(ad_comm::argv[on+1],&end);
00550           if (pprobe<=0.00001 || pprobe > .499)
00551           {
00552             cerr << "Invalid argument to option -mcgrope" << endl;
00553             pprobe=-1.0;
00554             probe_flag=0;
00555           }
00556         }
00557       }
00558 
00559        int start_flag;
00560        int java_quit_flag=0;
00561        int der_flag,next_flag;
00562        for (int i=1;i<=number_sims;i++)
00563        {
00564          if (user_stop()) break;
00565          if (java_quit_flag) break;
00566 
00567         if (!((i-1)%200))
00568         {
00569           double ratio = double(iac)/i;
00570           iac_old=iac-iac_old;
00571           i_old=i-i_old;
00572           cout << llc << " " << llc << endl;
00573           double tratio=double(liac)/200;
00574           liac=0;
00575           cout << " mcmc sim " << i <<  "  acceptance rate "
00576                << ratio << " " << tratio << endl;
00577 
00578          /*
00579           if (adjm_ptr && i>1)
00580           {
00581             ad_update_mcmc_report(mcmc_display_matrix,mcmc_save_index,
00582               mcmc_wrap_flag);
00583 
00584             ad_update_mcmc_stats_report
00585               (i,number_sims,100.*tratio,100.*ratio);
00586 
00587             if (allocated(hist)) ad_update_mcmchist_report(hist,
00588               number_offsets,mean_mcmc_values,h);
00589             check_java_flags(start_flag,java_quit_flag,der_flag,next_flag);
00590             ADSleep(50);
00591           }
00592           */
00593 
00594           if (i>50 && s_option && i<change_ball && !restart_flag)
00595           {
00596             if (tratio < .15)
00597             {
00598               chd=.2*chd;
00599               size_scale*=0.2;
00600               chdinv=chdinv/0.2;
00601               cout << "decreasing step size " << ln_det(chd,itmp) << endl;
00602             }
00603             if (tratio > .6)
00604             {
00605               chd=2.*chd;
00606               size_scale*=2.0;
00607               chdinv=chdinv/2.;
00608               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00609             }
00610             else if (tratio > .5)
00611             {
00612               chd=1.5*chd;
00613               size_scale*=1.5;
00614               chdinv=chdinv/1.5;
00615               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00616             }
00617             else if (tratio > .4)
00618             {
00619               chd=1.2*chd;
00620               size_scale*=1.2;
00621               chdinv=chdinv/1.2;
00622               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00623             }
00624           }
00625         }
00626         ii=1;
00627         initial_params::restore_all_values(parsave,ii);
00628         initial_params::set_all_simulation_bounds(symbds);
00629 
00630         // option of generating uniform or normal random variables
00631         dvector bmn1(1,nvar);
00632         if (!u_option)
00633         {
00634           if (!probe_flag)
00635             bmn1=bounded_multivariate_normal(nvar,symbds(1),symbds(2),
00636               chd,lprob,rng);
00637           else
00638             bmn1=new_probing_bounded_multivariate_normal(
00639               nvar,symbds(1),symbds(2),chd,lprob,pprobe,rng);
00640 
00641           initial_params::add_random_vector(bmn1);
00642           initial_params::xinit(y);
00643           // get the simulation bounds for the inverse transition
00644           initial_params::set_all_simulation_bounds(symbds);
00645           if (!probe_flag)
00646             bounded_multivariate_normal_mcmc(nvar,symbds(1),symbds(2),chd,
00647               lpinv,-1*(chdinv*bmn1),rng);
00648           else
00649             new_probing_bounded_multivariate_normal_mcmc(nvar,symbds(1),
00650               symbds(2), chd,lpinv,-1*(chdinv*bmn1),pprobe,rng);
00651 
00652           send_int_to_slaves(1);
00653           ll=-pvm_master_get_monte_carlo_value(nvar,y);
00654           //cout << ll << " " << llc << endl;
00655           double ldiff=lprob-lpinv;
00656           logr= ll - ldiff - llc;
00657         }
00658         else
00659         {
00660           dvector bmn1=bounded_multivariate_uniform(nvar,symbds(1),symbds(2),
00661                                                     chd, lprob,rng);
00662           initial_params::add_random_vector(bmn1);
00663           initial_params::xinit(y);
00664           // get the simulation bounds for the inverse transition
00665           initial_params::set_all_simulation_bounds(symbds);
00666           bounded_multivariate_uniform_mcmc(nvar,symbds(1),symbds(2),chd,
00667                                             lpinv,-1*(chdinv*bmn1),rng);
00668           send_int_to_slaves(1);
00669           ll=-pvm_master_get_monte_carlo_value(nvar,y);
00670           double ldiff=lprob-lpinv;
00671           logr= ll - ldiff - llc;
00672         }
00673         //cout << logr << endl;
00674         // decide whether to accept the new point
00675         isim++;
00676         double br=rng.better_rand();
00677         if (logr>=0 || br< exp(logr) )
00678         {
00679           ii=1;
00680           // save current parameter values
00681           initial_params::copy_all_values(parsave,ii);
00682           ii=1;
00683           // save current mcmc values
00684           stddev_params::copy_all_values(mcmc_values,ii);
00685           /*
00686           if (adjm_ptr)
00687           {
00688             ii=1;
00689             stddev_params::copy_all_number_values(mcmc_number_values,ii);
00690           }
00691           */
00692           // update prob density for current point
00693           llc =ll;
00694           liac++;
00695           iac++;
00696         }
00697         //tmpof << mcmc_values(1) << endl;
00698         //tofs << mcmc_values<< " "  << llc << endl;
00699         int mmin=mcmc_values.indexmin();
00700         int mmax=mcmc_values.indexmax();
00701         double lbf=llbest-llc;
00702         if (lbf>lbmax) lbmax=lbf;
00703         //of_bf << lbf << endl;
00704         bfsum+=exp(lbf);
00705         ibfcount+=1;
00706 #if defined(USE_BAYES_FACTORS)
00707         lkvector(ibfcount)=llc;
00708         //current_bf=1.0/(bfsum/double(ibfcount))*exp(llbest);
00709         lcurrent_bf=-1.*log(bfsum/double(ibfcount))+llbest;
00710 #endif
00711         if (mcsave_flag && (!((i-1)%mcsave_flag)))
00712         {
00713           (*pofs_psave) << parsave;
00714         }
00715        /*
00716         if (adjm_ptr)
00717         {
00718           save_mcmc_for_gui1(mcmc_values,mcmc_display_matrix,
00719             mcmc_save_index,mcmc_wrap_flag,number_offsets);
00720         }
00721         */
00722         if (!no_sd_mcmc)
00723         {
00724           if (!have_hist_flag)
00725           {
00726             for (ii=mmin;ii<=mmax;ii++)
00727             {
00728               (*pofs_sd) << float(mcmc_values(ii));
00729             }
00730             (*pofs_sd) << (float)(llc);
00731             mean_mcmc_values+=mcmc_values;
00732             square_mcmc_values+=square(mcmc_values);
00733           }
00734           else
00735           {
00736             add_hist_values(s,mean_mcmc_values,hist,mcmc_values,llc,
00737               h,nslots,total_spread);
00738           }
00739         }
00740         //if (scov_option)
00741         {
00742           ncor++;
00743           s_mean+=parsave;
00744           s_covar+=outer_prod(parsave,parsave);
00745         }
00746         if (!no_sd_mcmc && iac>5 && isim>initial_nsim)
00747         {
00748           if (!have_hist_flag)
00749           {
00750             have_hist_flag=1;
00751             delete pofs_sd;
00752             pofs_sd=NULL;
00753             mean_mcmc_values/=double(isim);
00754             square_mcmc_values/=double(isim);
00755             s=2.*sqrt(1.e-20+square_mcmc_values-square(mean_mcmc_values));
00756             make_preliminary_hist(s,mean_mcmc_values,isim,values,hist,h,
00757               nslots,total_spread);
00758           }
00759         }
00760        }
00761       }
00762       if (!no_sd_mcmc && !have_hist_flag)
00763       {
00764         have_hist_flag=1;
00765         delete pofs_sd;
00766         pofs_sd=NULL;
00767         mean_mcmc_values/=double(isim);
00768         square_mcmc_values/=double(isim);
00769         s=2.*sqrt(1.e-20+square_mcmc_values-square(mean_mcmc_values));
00770         make_preliminary_hist(s,mean_mcmc_values,isim,values,hist,h,nslots,
00771           total_spread);
00772       }
00773       if (!no_sd_mcmc)
00774         if (iac>5)
00775           print_hist_data(hist,values,h,mean_mcmc_values,s,parsave,iseed,
00776             size_scale);
00777       cout << iac/double(isim) << endl;
00778       initial_params::mc_phase=0;
00779      /*
00780       if (adjm_ptr)
00781       {
00782         ad_update_mcmc_report(mcmc_display_matrix,mcmc_save_index,
00783           mcmc_wrap_flag,1);
00784 
00785         if (allocated(hist)) ad_update_mcmchist_report(hist,
00786           number_offsets,mean_mcmc_values,h,1);
00787       }
00788       */
00789     }
00790 
00791     write_empirical_covariance_matrix(ncor,s_mean,s_covar,
00792       ad_comm::adprogram_name);
00793     //newton_raftery_bayes_estimate(current_bf,ibfcount,exp(lkvector),.01);
00794 #if defined(USE_BAYES_FACTORS)
00795     cout << "log current bayes factor " << lcurrent_bf << endl;
00796     newton_raftery_bayes_estimate_new(lcurrent_bf,ibfcount,lkvector,.01);
00797 #endif
00798     if (pofs_psave)
00799     {
00800       delete pofs_psave;
00801       pofs_psave=NULL;
00802     }
00803   }
00804   send_int_to_slaves(0);
00805 }
00806 #endif