ADMB Documentation  11.1.1920
 All Classes Files Functions Variables Typedefs Friends Defines
xxmcmc.cpp
Go to the documentation of this file.
00001 
00007 #include <admodel.h>
00008 
00009 #if defined(_MSC_VER)
00010   #include <conio.h>
00011 #else
00012   #define getch getchar
00013 #endif
00014 
00015 double better_rand(long int&);
00016 void store_mcmc_values(const ofstream& ofs);
00017 void set_labels_for_mcmc(void);
00018 void save_mcmc_for_gui(const dvector& mcmc_values,dmatrix &mdm,int& ids);
00019 void save_mcmc_for_gui1(const dvector& mcmc_values,
00020   dmatrix &mdm,int& ids,int& iwrap,ivector& no);
00021 
00022 void check_java_flags(int& start_flag,int& quit_flag,int& der_flag,
00023   int& next_flag);
00024 void print_hist_data(const dmatrix& hist, const dmatrix& values,
00025   const dvector& h, dvector& m, const dvector& s, const dvector& parsave,
00026   long int iseed, double size_scale);
00027 
00028 int minnz(const dvector& x);
00029 int maxnz(const dvector& xa);
00030 
00031 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS, double s,
00032   int mcmc2_flag);
00033 
00034 int read_hist_data(const dmatrix& hist, const dvector& h, dvector& m,
00035   const dvector& s, const dvector& parsave, long int& iseed,
00036   const double& size_scale);
00037 
00038 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
00039   const dmatrix& values, dmatrix& hist, const dvector& h,int slots,
00040   double total_spread,int probflag=0);
00041 
00042 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& hist,
00043   dvector& mcmc_values,double llc, const dvector& h,int nslots,
00044   double total_spreadd,int probflag=0);
00045 
00046 void add_guihist_values(const dvector& s, const dvector& m,
00047   const dmatrix& _hist, dvector& mcmcnumber_values, double llc,
00048   const dvector& h, int nslots, double total_spread);
00049 
00050 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00051   const dmatrix& s_covar, adstring& prog_name);
00052 
00053 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00054   const adstring& prog_name);
00055 
00056 void read_hessian_matrix_and_scale(int nvar, const dmatrix& S,
00057   const dvector& pen_vector);
00058 
00059 int user_stop(void);
00060 
00061 extern int ctlc_flag;
00062 
00063 dvector new_probing_bounded_multivariate_normal(int nvar, const dvector& a1,
00064   const dvector& b1, dmatrix& ch, const double& wght,double pprobe,
00065    random_number_generator& rng);
00066  // const random_number_generator& rng);
00067 
00068 void new_probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00069   const dvector& b1, dmatrix& ch, const double& wght, const dvector& _y,
00070   double pprobe, random_number_generator& rng);
00071 
00072 //void newton_raftery_bayes_estimate(double cbf,int ic, const dvector& lk,
00073 //  double d);
00074 void newton_raftery_bayes_estimate_new(double cbf,int ic, const dvector& lk,
00075   double d);
00076 
00077 void ad_update_mcmc_stats_report
00078   (int feval,int iter,double fval,double gmax);
00079 
00080 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00081   double fval,double gmax,const char * cbuf);
00082 void ad_update_mcmc_report(dmatrix& m,int i,int j,int ff=0);
00083 void ad_update_mcmchist_report(dmatrix& mcmc_values,ivector& number_offsets,
00084   dvector& mean_mcmc_values,dvector& h,int ff=0);
00085 
00086 void ADSleep(int);
00087 
00123 void function_minimizer::mcmc_routine(int nmcmc,int iseed0, double dscale,
00124   int restart_flag)
00125  {
00126   uostream * pofs_psave=NULL;
00127   dmatrix mcmc_display_matrix;
00128   //int mcmc_save_index=1;
00129   //int mcmc_wrap_flag=0;
00130   //int mcmc_gui_length=10000;
00131   int no_sd_mcmc=0;
00132 
00133   int on2=-1;
00134   if ( (on2=option_match(ad_comm::argc,ad_comm::argv,"-nosdmcmc"))>-1)
00135   {
00136     no_sd_mcmc=1;
00137   }
00138   if (mcmc2_flag==1)
00139   {
00140     initial_params::restore_start_phase();
00141     //get the number of active parameters
00142     //int nvar1=initial_params::nvarcalc();
00143   }
00144 
00145   if (stddev_params::num_stddev_params==0)
00146   {
00147     cerr << " You must declare at least one object of type sdreport "
00148          << endl << " to do the mcmc calculations" << endl;
00149      return;
00150   }
00151   {
00152     ivector number_offsets;
00153     dvector lkvector;
00154     //double current_bf=0;
00155 #if defined(USE_BAYES_FACTORS)
00156     double lcurrent_bf=0;
00157 #endif
00158     double size_scale=1.0;
00159     double total_spread=200;
00160     //double total_spread=2500;
00161     uostream * pofs_sd = NULL;
00162 
00163 #if defined(USE_LAPLACE)
00164     initial_params::set_inactive_random_effects();
00165     //int nvar_x=initial_params::nvarcalc();
00166     initial_params::set_active_random_effects();
00167     int nvar_re=initial_params::nvarcalc();
00168 #endif
00169 
00170     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00171     //int scov_option=0;
00172     dmatrix s_covar;
00173     dvector s_mean;
00174     int on=-1;
00175     int nslots=800;
00176     //int nslots=3600;
00177     int initial_nsim=4800;
00178     int ncor=0;
00179     double bfsum=0;
00180     int ibfcount=0;
00181     double llbest;
00182     double lbmax;
00183 
00184     //int ntmp=0;
00185     //if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcscov",ntmp))>-1)
00186     //{
00187     //scov_option=1;
00188     s_covar.allocate(1,nvar,1,nvar);
00189     s_mean.allocate(1,nvar);
00190     s_mean.initialize();
00191     s_covar.initialize();
00192 
00193     int ndvar=stddev_params::num_stddev_calc();
00194     /*
00195     int numdvar=stddev_params::num_stddev_number_calc();
00196     if (adjm_ptr)
00197     {
00198       mcmc_display_matrix.allocate(1,numdvar,1,mcmc_gui_length);
00199       number_offsets.allocate(1,numdvar);
00200       number_offsets=stddev_params::copy_all_number_offsets();
00201     }
00202     */
00203 #if defined(USE_LAPLACE)
00204     if (mcmc2_flag==0)
00205     {
00206       initial_params::set_inactive_random_effects();
00207       nvar=initial_params::nvarcalc(); // get the number of active parameters
00208     }
00209 #endif
00210     dvector x(1,nvar);
00211     dvector scale(1,nvar);
00212     dmatrix values;
00213     int have_hist_flag=0;
00214     initial_params::xinit(x);
00215     dvector pen_vector(1,nvar);
00216     {
00217       initial_params::reset(dvar_vector(x),pen_vector);
00218       //cout << pen_vector << endl << endl;
00219     }
00220 
00221     initial_params::mc_phase=0;
00222     initial_params::stddev_scale(scale,x);
00223     initial_params::mc_phase=1;
00224     dvector bmn(1,nvar);
00225     dvector mean_mcmc_values(1,ndvar);
00226     dvector s(1,ndvar);
00227     dvector h(1,ndvar);
00228     //dvector h;
00229     dvector square_mcmc_values(1,ndvar);
00230     square_mcmc_values.initialize();
00231     mean_mcmc_values.initialize();
00232     bmn.initialize();
00233     int use_empirical_flag=0;
00234     int diag_option=0;
00235     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcdiag"))>-1)
00236     {
00237       diag_option=1;
00238       cout << " Setting covariance matrix to diagonal with entries " << dscale
00239            << endl;
00240     }
00241     dmatrix S(1,nvar,1,nvar);
00242     dvector sscale(1,nvar);
00243     if (!diag_option)
00244     {
00245       int on,nopt = 0;
00246       int rescale_bounded_flag=0;
00247       double rescale_bounded_power=0.5;
00248       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcrb",nopt))>-1)
00249       {
00250         if (nopt)
00251         {
00252           int iii=atoi(ad_comm::argv[on+1]);
00253           if (iii < 1 || iii > 9)
00254           {
00255             cerr << " -mcrb argument must be integer between 1 and 9 --"
00256                     " using default of 5" << endl;
00257             rescale_bounded_power=0.5;
00258           }
00259           else
00260             rescale_bounded_power=iii/10.0;
00261         }
00262         else
00263         {
00264           rescale_bounded_power=0.5;
00265         }
00266         rescale_bounded_flag=1;
00267       }
00268       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcec"))>-1)
00269       {
00270         use_empirical_flag=1;
00271       }
00272       if (use_empirical_flag)
00273       {
00274         read_empirical_covariance_matrix(nvar,S,ad_comm::adprogram_name);
00275       }
00276       else if (!rescale_bounded_flag)
00277       {
00278         int hbf;
00279         if (mcmc2_flag==0)
00280         {
00281           read_covariance_matrix(S,nvar,hbf,sscale);
00282         }
00283         else
00284         {
00285           int tmp_nvar = 0;
00286           adstring tmpstring = ad_comm::adprogram_name + ".bgs";
00287           uistream uis((char*)(tmpstring));
00288           if (!uis)
00289           {
00290             cerr << "error opening file " << tmpstring << endl;
00291             ad_exit(1);
00292           }
00293           uis >> tmp_nvar;
00294           if (!uis)
00295           {
00296             cerr << "error reading from file " << tmpstring << endl;
00297             ad_exit(1);
00298           }
00299           if (tmp_nvar != nvar)
00300           {
00301             cerr << "size error reading from " <<  tmpstring << endl;
00302             ad_exit(1);
00303           }
00304           uis >> S;
00305           if (!uis)
00306           {
00307             cerr << "error reading from file " << tmpstring << endl;
00308             ad_exit(1);
00309           }
00310         }
00311       }
00312       else
00313       {
00314         read_hessian_matrix_and_scale1(nvar,S,rescale_bounded_power,
00315           mcmc2_flag);
00316         //read_hessian_matrix_and_scale(nvar,S,pen_vector);
00317       }
00318 
00319       {  // scale covariance matrix for model space
00320         dmatrix tmp(1,nvar,1,nvar);
00321         for (int i=1;i<=nvar;i++)
00322         {
00323           tmp(i,i)=S(i,i)*(scale(i)*scale(i));
00324           for (int j=1;j<i;j++)
00325           {
00326             tmp(i,j)=S(i,j)*(scale(i)*scale(j));
00327             tmp(j,i)=tmp(i,j);
00328           }
00329         }
00330         S=tmp;
00331       }
00332     }
00333     else
00334     {
00335       S.initialize();
00336       for (int i=1;i<=nvar;i++)
00337       {
00338         S(i,i)=dscale;
00339       }
00340     }
00341 
00342     cout << sort(eigenvalues(S)) << endl;
00343     dmatrix chd = choleski_decomp( (dscale*2.4/sqrt(double(nvar))) * S);
00344     dmatrix chdinv=inv(chd);
00345 
00346     dmatrix symbds(1,2,1,nvar);
00347     initial_params::set_all_simulation_bounds(symbds);
00348     ofstream ofs_sd1((char*)(ad_comm::adprogram_name + adstring(".mc2")));
00349 
00350     {
00351       long int iseed=0;
00352       int number_sims;
00353       if (nmcmc<=0)
00354       {
00355         number_sims=  100000;
00356       }
00357       else
00358       {
00359         number_sims=  nmcmc;
00360       }
00361       //cin >> iseed;
00362       if (iseed0<=0)
00363       {
00364         iseed=-36519;
00365       }
00366       else
00367       {
00368         iseed=-iseed0;
00369       }
00370       if (iseed>0)
00371       {
00372         iseed=-iseed;
00373       }
00374       cout << "Initial seed value " << iseed << endl;
00375       random_number_generator rng(iseed);
00376       rng.better_rand();
00377       double lprob=0.0;
00378       double lpinv=0.0;
00379       // get lower and upper bounds
00380 
00381       independent_variables y(1,nvar);
00382 #if defined(USE_LAPLACE)
00383       independent_variables parsave(1,nvar_re);
00384       initial_params::restore_start_phase();
00385 #else
00386       independent_variables parsave(1,nvar);
00387 #endif
00388 
00389       // read in the mcmc values to date
00390       int ii=1;
00391       dmatrix hist;
00392       if (restart_flag)
00393       {
00394         int tmp=0;
00395         if (!no_sd_mcmc) {
00396           hist.allocate(1,ndvar,-nslots,nslots);
00397           tmp=read_hist_data(hist,h,mean_mcmc_values,s,parsave,iseed,
00398             size_scale);
00399           values.allocate(1,ndvar,-nslots,nslots);
00400           for (int i=1;i<=ndvar;i++)
00401           {
00402             values(i).fill_seqadd(mean_mcmc_values(i)-0.5*total_spread*s(i)
00403               +.5*h(i),h(i));
00404           }
00405         }
00406         if (iseed>0)
00407         {
00408           iseed=-iseed;
00409         }
00410         /*double br=*/rng.better_rand();
00411         if (tmp) have_hist_flag=1;
00412         chd=size_scale*chd;
00413         chdinv=chdinv/size_scale;
00414       }
00415       else
00416       {
00417         int on=-1;
00418         int nopt=0;
00419         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcpin",nopt))>-1)
00420         {
00421           if (nopt)
00422           {
00423             cifstream cif((char *)ad_comm::argv[on+1]);
00424             if (!cif)
00425             {
00426               cerr << "Error trying to open mcmc par input file "
00427                    << ad_comm::argv[on+1] << endl;
00428               exit(1);
00429             }
00430             cif >> parsave;
00431             if (!cif)
00432             {
00433               cerr << "Error reading from mcmc par input file "
00434                    << ad_comm::argv[on+1] << endl;
00435               exit(1);
00436             }
00437           }
00438           else
00439           {
00440             cerr << "Illegal option with -mcpin" << endl;
00441           }
00442         }
00443         else
00444         {
00445           ii=1;
00446           initial_params::copy_all_values(parsave,ii);
00447         }
00448       }
00449 
00450       ii=1;
00451       initial_params::restore_all_values(parsave,ii);
00452 
00453 #if defined(USE_LAPLACE)
00454       if (mcmc2_flag==0)
00455       {
00456         initial_params::set_inactive_random_effects();
00457       }
00458 #endif
00459 
00460       gradient_structure::set_NO_DERIVATIVES();
00461       ofstream ogs("sims");
00462       ogs << nvar << " " << number_sims << endl;
00463       initial_params::xinit(y);
00464       double llc=-get_monte_carlo_value(nvar,y);
00465       llbest=-get_monte_carlo_value(nvar,y);
00466       lbmax=llbest;
00467       // store current mcmc variable values in param_values
00468       //dmatrix store_mcmc_values(1,number_sims,1,ndvar);
00469 #if defined(USE_BAYES_FACTORS)
00470       lkvector.allocate(1,number_sims);
00471 #endif
00472       dvector mcmc_values(1,ndvar);
00473       dvector mcmc_number_values;
00474       //if (adjm_ptr) mcmc_number_values.allocate(1,numdvar);
00475       int offs=1;
00476       stddev_params::copy_all_values(mcmc_values,offs);
00477 
00478       /*
00479       if (adjm_ptr)
00480       {
00481         offs=1;
00482         stddev_params::copy_all_number_values(mcmc_number_values,offs);
00483       }
00484       */
00485       int change_ball=2500;
00486       int nopt = 0;
00487       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcscale",nopt))>-1)
00488       {
00489         if (nopt)
00490         {
00491           int iii=atoi(ad_comm::argv[on+1]);
00492           if (iii <=0)
00493           {
00494             cerr << " Invalid option following command line option -mcscale -- "
00495               << endl << " ignored" << endl;
00496           }
00497           else
00498             change_ball=iii;
00499         }
00500       }
00501       int iac=0;
00502       int liac=0;
00503       int isim=0;
00504       int itmp=0;
00505       double logr;
00506       int u_option=0;
00507       double ll;
00508       int s_option=1;
00509       int psvflag=0;
00510       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcu"))>-1)
00511       {
00512         u_option=1;
00513       }
00514       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcnoscale"))>-1)
00515       {
00516         s_option=0;
00517       }
00518       //cout << llc << " " << llc << endl;
00519       int iac_old=0;
00520       int i_old=0;
00521 
00522       {
00523        if (!restart_flag)
00524        {
00525          pofs_sd =
00526            new uostream((char*)(ad_comm::adprogram_name + adstring(".mcm")));
00527        }
00528 
00529       int mcsave_flag=0;
00530       int mcrestart_flag=option_match(ad_comm::argc,ad_comm::argv,"-mcr");
00531 
00532       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcsave"))>-1)
00533       {
00534         int jj=(int)atof(ad_comm::argv[on+1]);
00535         if (jj <=0)
00536         {
00537           cerr << " Invalid option following command line option -mcsave -- "
00538             << endl;
00539         }
00540         else
00541         {
00542           mcsave_flag=jj;
00543           if ( mcrestart_flag>-1)
00544           {
00545             // check that nvar is correct
00546             {
00547               uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00548               if (!uis)
00549               {
00550                 cerr << "Error trying to open file" <<
00551                   ad_comm::adprogram_name + adstring(".psv") <<
00552                   " for mcrestart" <<   endl;
00553                 cerr << " I am starting a new file " << endl;
00554                 psvflag=1;
00555               }
00556               else
00557               {
00558                 int nv1 = 0;
00559                 uis >> nv1;
00560                 if (nv1 !=nvar)
00561                 {
00562                   cerr << "wrong number of independent variables in" <<
00563                     ad_comm::adprogram_name + adstring(".psv") <<
00564                   cerr << " I am starting a new file " << endl;
00565                   psvflag=1;
00566                 }
00567               }
00568             }
00569 
00570             if (!psvflag) {
00571               pofs_psave=
00572                 new uostream(
00573                   (char*)(ad_comm::adprogram_name + adstring(".psv")),ios::app);
00574             } else {
00575               pofs_psave= new uostream((char*)(ad_comm::adprogram_name
00576                 + adstring(".psv")));
00577             }
00578           } else {
00579             pofs_psave=
00580               new uostream((char*)(ad_comm::adprogram_name + adstring(".psv")));
00581           }
00582           if (!pofs_psave|| !(*pofs_psave))
00583           {
00584             cerr << "Error trying to open file" <<
00585               ad_comm::adprogram_name + adstring(".psv") << endl;
00586             mcsave_flag=0;
00587             if (pofs_psave)
00588             {
00589               delete pofs_psave;
00590               pofs_psave=NULL;
00591             }
00592           }
00593           else
00594           {
00595             if (psvflag || (mcrestart_flag == -1) )
00596             {
00597 #           if defined(USE_LAPLACE)
00598               (*pofs_psave) << nvar_re;
00599 #           else
00600               (*pofs_psave) << nvar;
00601 #           endif
00602             }
00603           }
00604         }
00605       }
00606 
00607       double pprobe=0.05;
00608       int probe_flag=0;
00609       int nopt=0;
00610       on = option_match(ad_comm::argc, ad_comm::argv, "-mcprobe", nopt);
00611       if (on == -1)
00612       {
00613         on = option_match(ad_comm::argc,ad_comm::argv,"-mcgrope",nopt);
00614       }
00615       if (on > -1)
00616       {
00617         probe_flag=1;
00618         if (nopt)
00619         {
00620           char* end = 0;
00621           pprobe=strtod(ad_comm::argv[on+1],&end);
00622           if (pprobe<=0.00001 || pprobe > .499)
00623           {
00624             cerr << "Invalid argument to option -mcprobe" << endl;
00625             pprobe=-1.0;
00626             probe_flag=0;
00627           }
00628         }
00629       }
00630 
00631        int java_quit_flag=0;
00632        for (int i=1;i<=number_sims;i++)
00633        {
00634          if (user_stop()) break;
00635          if (java_quit_flag) break;
00636 
00637         if (!((i-1)%200))
00638         {
00639           double ratio = double(iac)/i;
00640           iac_old=iac-iac_old;
00641           i_old=i-i_old;
00642           cout << llc << " " << llc << endl;
00643           double tratio=double(liac)/200;
00644           liac=0;
00645           cout << " mcmc sim " << i <<  "  acceptance rate "
00646                << ratio << " " << tratio << endl;
00647 
00648          /*
00649           int start_flag;
00650           int der_flag,next_flag;
00651           if (adjm_ptr && i>1)
00652           {
00653             ad_update_mcmc_report(mcmc_display_matrix,mcmc_save_index,
00654               mcmc_wrap_flag);
00655 
00656             ad_update_mcmc_stats_report
00657               (i,number_sims,100.*tratio,100.*ratio);
00658 
00659             if (allocated(hist)) ad_update_mcmchist_report(hist,
00660               number_offsets,mean_mcmc_values,h);
00661             check_java_flags(start_flag,java_quit_flag,der_flag,next_flag);
00662             ADSleep(50);
00663           }
00664           */
00665 
00666           if (i>50 && s_option && i<change_ball && !restart_flag)
00667           {
00668             if (tratio < .15)
00669             {
00670               chd=.2*chd;
00671               size_scale*=0.2;
00672               chdinv=chdinv/0.2;
00673               cout << "decreasing step size " << ln_det(chd,itmp) << endl;
00674             }
00675             if (tratio > .6)
00676             {
00677               chd=2.*chd;
00678               size_scale*=2.0;
00679               chdinv=chdinv/2.;
00680               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00681             }
00682             else if (tratio > .5)
00683             {
00684               chd=1.5*chd;
00685               size_scale*=1.5;
00686               chdinv=chdinv/1.5;
00687               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00688             }
00689             else if (tratio > .4)
00690             {
00691               chd=1.2*chd;
00692               size_scale*=1.2;
00693               chdinv=chdinv/1.2;
00694               cout << "increasing step size " << ln_det(chd,itmp) << endl;
00695             }
00696           }
00697         }
00698         ii=1;
00699 #if defined(USE_LAPLACE)
00700         if (random_effects_flag)
00701         {
00702           initial_params::restore_start_phase();
00703           initial_params::restore_all_values(parsave,ii);
00704           if (mcmc2_flag==0)
00705           {
00706             initial_params::set_inactive_random_effects();
00707           }
00708         }
00709         else
00710         {
00711           initial_params::restore_all_values(parsave,ii);
00712         }
00713 #else
00714         initial_params::restore_all_values(parsave,ii);
00715 #endif
00716         initial_params::set_all_simulation_bounds(symbds);
00717 
00718         // option of generating uniform or normal random variables
00719         dvector bmn1(1,nvar);
00720         if (!u_option)
00721         {
00722           if (!probe_flag)
00723             bmn1=bounded_multivariate_normal(nvar,symbds(1),symbds(2),
00724               chd,lprob,rng);
00725           else
00726             bmn1=new_probing_bounded_multivariate_normal(
00727               nvar,symbds(1),symbds(2),chd,lprob,pprobe,rng);
00728 
00729           initial_params::add_random_vector(bmn1);
00730           initial_params::xinit(y);
00731           // get the simulation bounds for the inverse transition
00732           initial_params::set_all_simulation_bounds(symbds);
00733           if (!probe_flag)
00734             bounded_multivariate_normal_mcmc(nvar,symbds(1),symbds(2),chd,
00735               lpinv,-1*(chdinv*bmn1),rng);
00736           else
00737             new_probing_bounded_multivariate_normal_mcmc(nvar,symbds(1),
00738               symbds(2), chd,lpinv,-1*(chdinv*bmn1),pprobe,rng);
00739 
00740           ll=-get_monte_carlo_value(nvar,y);
00741 #if defined(USE_LAPLACE)
00742           if (random_effects_flag)
00743           {
00744             initial_params::restore_start_phase();
00745           }
00746 #endif
00747           //cout << ll << " " << llc << endl;
00748           double ldiff=lprob-lpinv;
00749           logr= ll - ldiff - llc;
00750         }
00751         else
00752         {
00753           dvector bmn1=bounded_multivariate_uniform(nvar,symbds(1),symbds(2),
00754             chd, lprob,rng);
00755           initial_params::add_random_vector(bmn1);
00756           initial_params::xinit(y);
00757           // get the simulation bounds for the inverse transition
00758           initial_params::set_all_simulation_bounds(symbds);
00759           bounded_multivariate_uniform_mcmc(nvar,symbds(1),symbds(2),chd,
00760                                             lpinv,-1*(chdinv*bmn1),rng);
00761           ll=-get_monte_carlo_value(nvar,y);
00762           double ldiff=lprob-lpinv;
00763           logr= ll - ldiff - llc;
00764         }
00765         //cout << logr << endl;
00766         // decide whether to accept the new point
00767         isim++;
00768         double br=rng.better_rand();
00769         if (logr>=0 || br< exp(logr) )
00770         {
00771           ii=1;
00772           // save current parameter values
00773           initial_params::copy_all_values(parsave,ii);
00774           ii=1;
00775           // save current mcmc values
00776           stddev_params::copy_all_values(mcmc_values,ii);
00777 #         if defined(USE_LAPLACE)
00778             if (random_effects_flag)
00779             {
00780               if (mcmc2_flag==0)
00781               {
00782                 initial_params::set_inactive_random_effects();
00783               }
00784             }
00785 #         endif
00786           // update prob density for current point
00787           llc =ll;
00788           liac++;
00789           iac++;
00790         }
00791         int mmin=mcmc_values.indexmin();
00792         int mmax=mcmc_values.indexmax();
00793         double lbf=llbest-llc;
00794         if (lbf>lbmax) lbmax=lbf;
00795         //of_bf << lbf << endl;
00796         bfsum+=exp(lbf);
00797         ibfcount+=1;
00798 #if defined(USE_BAYES_FACTORS)
00799         lkvector(ibfcount)=llc;
00800         //current_bf=1.0/(bfsum/double(ibfcount))*exp(llbest);
00801         lcurrent_bf=-1.*log(bfsum/double(ibfcount))+llbest;
00802 #endif
00803         if (mcsave_flag && (!((i-1)%mcsave_flag)))
00804         {
00805           (*pofs_psave) << parsave;
00806         }
00807        /*
00808         if (adjm_ptr)
00809         {
00810           save_mcmc_for_gui1(mcmc_values,mcmc_display_matrix,
00811             mcmc_save_index,mcmc_wrap_flag,number_offsets);
00812         }
00813         */
00814         if (!no_sd_mcmc)
00815         {
00816           if (!have_hist_flag)
00817           {
00818             for (ii=mmin;ii<=mmax;ii++)
00819             {
00820               (*pofs_sd) << float(mcmc_values(ii));
00821             }
00822             (*pofs_sd) << (float)(llc);
00823             mean_mcmc_values+=mcmc_values;
00824             square_mcmc_values+=square(mcmc_values);
00825           }
00826           else
00827           {
00828             add_hist_values(s,mean_mcmc_values,hist,mcmc_values,llc,
00829               h,nslots,total_spread);
00830           }
00831         }
00832         //if (scov_option)
00833         {
00834           ncor++;
00835           s_mean+=parsave;
00836           s_covar+=outer_prod(parsave,parsave);
00837         }
00838         if (!no_sd_mcmc && iac>5 && isim>initial_nsim)
00839         {
00840           if (!have_hist_flag)
00841           {
00842             have_hist_flag=1;
00843             delete pofs_sd;
00844             pofs_sd=NULL;
00845             mean_mcmc_values/=double(isim);
00846             square_mcmc_values/=double(isim);
00847 #if defined(USE_DDOUBLE)
00848             s=2.*sqrt(double(1.e-20)+square_mcmc_values
00849               -square(mean_mcmc_values));
00850 #else
00851             s=2.*sqrt(1.e-20+square_mcmc_values-square(mean_mcmc_values));
00852 #endif
00853             make_preliminary_hist(s,mean_mcmc_values,isim,values,hist,h,
00854               nslots,total_spread);
00855           }
00856         }
00857        }
00858       }
00859       if (!no_sd_mcmc && !have_hist_flag)
00860       {
00861         have_hist_flag=1;
00862         delete pofs_sd;
00863         pofs_sd=NULL;
00864         mean_mcmc_values/=double(isim);
00865         square_mcmc_values/=double(isim);
00866 #if defined(USE_DDOUBLE)
00867         s=2.*sqrt(double(1.e-20)+square_mcmc_values-square(mean_mcmc_values));
00868 #else
00869         s=2.*sqrt(1.e-20+square_mcmc_values-square(mean_mcmc_values));
00870 #endif
00871         make_preliminary_hist(s,mean_mcmc_values,isim,values,hist,h,nslots,
00872           total_spread);
00873       }
00874       if (!no_sd_mcmc)
00875         if (iac>5)
00876           print_hist_data(hist,values,h,mean_mcmc_values,s,parsave,iseed,
00877             size_scale);
00878       cout << iac/double(isim) << endl;
00879       initial_params::mc_phase=0;
00880      /*
00881       if (adjm_ptr)
00882       {
00883         ad_update_mcmc_report(mcmc_display_matrix,mcmc_save_index,
00884           mcmc_wrap_flag,1);
00885 
00886         if (allocated(hist)) ad_update_mcmchist_report(hist,
00887           number_offsets,mean_mcmc_values,h,1);
00888       }
00889       */
00890     }
00891 
00892     write_empirical_covariance_matrix(ncor,s_mean,s_covar,
00893       ad_comm::adprogram_name);
00894     //newton_raftery_bayes_estimate(current_bf,ibfcount,exp(lkvector),.01);
00895 #if defined(USE_BAYES_FACTORS)
00896     cout << "log current bayes factor " << lcurrent_bf << endl;
00897     newton_raftery_bayes_estimate_new(lcurrent_bf,ibfcount,lkvector,.01);
00898 #endif
00899     if (pofs_psave)
00900     {
00901       delete pofs_psave;
00902       pofs_psave=NULL;
00903     }
00904   }
00905 }
00906 
00917 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00918   const dmatrix& s_covar,adstring& prog_name)
00919 {
00920   uostream ofs((char*)(ad_comm::adprogram_name + adstring(".ecm")));
00921   dvector tmp=s_mean/ncor;
00922   int nvar=s_mean.indexmax();
00923   //ofs << ncor << " " << nvar << endl;
00924   dmatrix sigma=s_covar/ncor -outer_prod(tmp,tmp);
00925   cout << "In write empirical covariance matrix" << endl;
00926   cout << sort(eigenvalues(sigma)) << endl;
00927   dvector std(1,nvar);
00928   ofs << sigma;
00929   /*
00930   int i;
00931   for (i=1;i<=nvar;i++)
00932   {
00933     std(i)=sigma(i,i);
00934     for (int j=1;j<=i;j++)
00935     {
00936       sigma(i,j)/=sqrt(std(i)*std(j));
00937       sigma(j,i)=sigma(i,j);
00938     }
00939   }
00940   for (i=1;i<=nvar;i++)
00941   {
00942     ofs << setw(10) << setscientific() << std(i);
00943     for (int j=1;j<=i;j++)
00944     {
00945       ofs << " " << setfixed() << setw(7) << setprecision(4)
00946            << sigma(i,j);
00947     }
00948     ofs << endl;
00949   }
00950  */
00951 }
00952 
00961 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00962   const adstring& prog_name)
00963 {
00964   adstring infile_name=ad_comm::adprogram_name + adstring(".ecm");
00965   uistream ifs((char*)(infile_name));
00966   if (!ifs)
00967   {
00968     cerr << "Error opening file " << infile_name << endl;
00969   }
00970   ifs  >> S;
00971   /*
00972   int ncor=0;
00973   int _nvar=0;
00974   ifs >> ncor >> _nvar;
00975   if (nvar != _nvar)
00976   {
00977     cerr << "wromng number of indepdendent variables in routine"
00978       " read_empirical_covariance_matrix" << endl;
00979   }
00980   dvector std(1,nvar);
00981   int i;
00982   for (i=1;i<=nvar;i++)
00983   {
00984     ifs >> std(i);
00985     for (int j=1;j<=i;j++)
00986     {
00987       ifs  >> S(i,j);
00988       S(j,i)=S(i,j);
00989     }
00990   }
00991   if (!ifs)
00992   {
00993     cerr << "Error reading from file " << infile_name << endl;
00994   }
00995   for (i=1;i<=nvar;i++)
00996   {
00997     for (int j=1;j<=i;j++)
00998     {
00999       S(i,j)*=sqrt(std(i)*std(j));
01000       S(j,i)=S(i,j);
01001     }
01002   }
01003  */
01004 }
01005 
01006 void print_hist_data(const dmatrix& hist, const dmatrix& values,
01007   const dvector& h, dvector& m, const dvector& s, const dvector& parsave,
01008   long int iseed, double size_scale)
01009 {
01010   ofstream ofs((char*)(ad_comm::adprogram_name + adstring(".hst")));
01011   int nsdvars=stddev_params::num_stddev_calc();
01012   adstring_array param_labels(1,nsdvars);
01013   ivector param_size(1,nsdvars);
01014   int ii=1;
01015   //int max_name_length=0;
01016   int i;
01017   for (i=0;i< stddev_params::num_stddev_params;i++)
01018   {
01019     param_labels[ii]=
01020       stddev_params::stddevptr[i]->label();
01021     param_size[ii]=
01022       stddev_params::stddevptr[i]->size_count();
01023 /*
01024     if (max_name_length<param_labels[ii].size())
01025     {
01026       max_name_length=param_labels[ii].size();
01027     }
01028 */
01029     ii++;
01030   }
01031   //int end_stdlabels=ii-1;
01032 
01033   int lc=1;
01034   int ic=1;
01035   ivector mmin(1,nsdvars);
01036   ivector mmax(1,nsdvars);
01037 
01038   int nsim = sum(hist(1));
01039   for (i=1;i<=nsdvars;i++)
01040   {
01041     mmin(i)=minnz(hist(i));
01042     mmax(i)=maxnz(hist(i));
01043   }
01044   ofs << "# samples sizes" << endl;
01045   ofs << nsim << endl;
01046   ofs << "# step size scaling factor" << endl;
01047   ofs << size_scale << endl;
01048   ofs << "# step sizes" << endl;
01049   ofs << h << endl;
01050   ofs << "# means" << endl;
01051   ofs << m << endl;
01052   ofs << "# standard devs" << endl;
01053   ofs << s << endl;
01054   ofs << "# lower bounds" << endl;
01055   ofs << mmin << endl;
01056   ofs << "# upper bounds" << endl;
01057   ofs << mmax<< endl;
01058   ofs << "#number of parameters" << endl;
01059   ofs << parsave.indexmax() << endl;
01060   ofs << "#current parameter values for mcmc restart" << endl;
01061   ofs << parsave << endl;
01062   ofs << "#random number seed" << endl;
01063   ofs << iseed << endl;
01064   for (i=1;i<=nsdvars;i++)
01065   {
01066     ofs << "#" << param_labels[lc];
01067     if (param_size[lc]>1)
01068     {
01069       ofs << "[" << ic << "]";
01070     }
01071     ofs << endl;
01072 
01073     if (++ic> param_size[lc])
01074     {
01075       lc++;
01076       ic=1;
01077     }
01078     for (ii=mmin(i);ii<=mmax(i);ii++)
01079     {
01080       ofs << values(i,ii) << " " << hist(i,ii)/(nsim*h(i)) << endl;
01081     }
01082     if (i<nsdvars) ofs << endl;
01083   }
01084 }
01085 
01086 int read_hist_data(const dmatrix& _hist, const dvector& h, dvector& m,
01087   const dvector& s, const dvector& parsave, long int& iseed,
01088   const double& size_scale)
01089 {
01090   dmatrix& hist= (dmatrix&) _hist;
01091   cifstream cif((char*)(ad_comm::adprogram_name + adstring(".hst")));
01092   int nsdvars=stddev_params::num_stddev_calc();
01093   int nsim=0;
01094   int ii=1;
01095   int i;
01096   ivector mmin(1,nsdvars);
01097   ivector mmax(1,nsdvars);
01098   //ofs << # samples sizes << endl;
01099   cif >> nsim;
01100   cif >> size_scale;
01101   //ofs << # step sizes << endl;
01102   cif >> h;
01103   //ofs << # means << endl;
01104   cif >> m;
01105   //ofs << # standard devs << endl;
01106   cif >> s;
01107   //ofs << # lower bounds << endl;
01108   cif >> mmin;
01109   //ofs << # upper bounds << endl;
01110   cif >> mmax;
01111   int num_vars=0;
01112   cif >> num_vars;
01113   int flag=1;
01114   if (num_vars!=parsave.indexmax())
01115   {
01116     cerr << "wrong number of independent variables in file"
01117        << ad_comm::adprogram_name + adstring(".hst") << endl;
01118     flag=0;
01119   }
01120   if (flag)
01121   {
01122   cif >> parsave;
01123     cif >> iseed;
01124     double tmp=0.0;
01125     hist.initialize();
01126     for (i=1;i<=nsdvars;i++)
01127     {
01128       for (ii=mmin(i);ii<=mmax(i);ii++)
01129       {
01130         cif >> tmp >> hist(i,ii);
01131       }
01132       hist(i)*=nsim*h(i);
01133     }
01134   }
01135   return flag;
01136 }
01137 
01138 int minnz(const dvector& x)
01139 {
01140   int mmin=x.indexmin();
01141   int mmax=x.indexmax();
01142   int m=mmin;
01143   for (int ii=mmin;ii<=mmax;ii++)
01144   {
01145     if (x(ii))
01146     {
01147       m=ii;
01148       if (m>mmin) m--;
01149       break;
01150     }
01151   }
01152   return m;
01153 }
01154 
01155 int maxnz(const dvector& x)
01156 {
01157   int mmin=x.indexmin();
01158   int mmax=x.indexmax();
01159   int m=mmax;
01160   for (int ii=mmax;ii>=mmin;ii--)
01161   {
01162     if (x(ii))
01163     {
01164       m=ii;
01165       if (m<mmax) m++;
01166       break;
01167     }
01168   }
01169   return m;
01170 }
01171 
01172 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& _hist,
01173   dvector& mcmc_values, double llc, const dvector& h, int nslots,
01174   double total_spread,int probflag)
01175 {
01176   dmatrix& hist= (dmatrix&) _hist;
01177   int nsdvars=stddev_params::num_stddev_calc();
01178   for (int ii=1;ii<=nsdvars;ii++)
01179   {
01180     int x;
01181     double xx=(mcmc_values(ii)-m(ii))/h(ii);
01182     if (xx>0.0)
01183       x=int (xx+0.5);
01184     else
01185       x=int(xx-0.5);
01186 
01187     if (x<-nslots)
01188     {
01189       x=-nslots;
01190     }
01191     if (x>nslots)
01192     {
01193       x=nslots;
01194     }
01195     if (!probflag)
01196     {
01197       hist(ii,x)+=1;
01198     }
01199     else
01200     {
01201       hist(ii,x)+=exp(llc);
01202     }
01203   }
01204 }
01205 
01206 /*
01207 void add_guihist_values(const dvector& s, const dvector& m,
01208   const dmatrix& _hist,dvector& mcmcnumber_values,double llc,
01209   const dvector& h,int nslots,double total_spread)
01210 {
01211   dmatrix& hist= (dmatrix&) _hist;
01212   int nsdvars=stddev_params::num_stddev_number_calc();
01213   for (int ii=1;ii<=nsdvars;ii++)
01214   {
01215     int x=int((mcmcnumber_values(ii)-m(ii))/h(ii));
01216     //cout << "xxx = " << xxx << endl;
01217     char ch;
01218     //cin >> ch;
01219 
01220     if (x<1)
01221     {
01222       x=-nslots;
01223     }
01224     if (x>nslots)
01225     {
01226       x=nslots;
01227     }
01228     {
01229       hist(ii,x)+=1;
01230     }
01231   }
01232 }
01233 */
01234 
01235 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
01236   const dmatrix& _values, dmatrix& hist, const dvector& _h, int nslots,
01237   double total_spread, int probflag)
01238 {
01239   ADUNCONST(dmatrix,values)
01240   dvector& h = (dvector&) _h;
01241   int nsdvars=stddev_params::num_stddev_calc();
01242   dvector mcmc_values(1,nsdvars);
01243   values.allocate(1,nsdvars,-nslots,nslots);
01244   h=total_spread/(2*nslots+1)*s;
01245   hist.allocate(1,nsdvars,-nslots,nslots);
01246   hist.initialize();
01247   uistream ifs((char*)(ad_comm::adprogram_name + adstring(".mcm")));
01248   int i;
01249   double llc;
01250   for (i=1;i<=nsdvars;i++)
01251   {
01252     values(i).fill_seqadd(m(i)-.5*total_spread*s(i)+.5*h(i),h(i));
01253   }
01254 
01255   if (!ifs)
01256   {
01257     cerr << "Error trying to open file "
01258          << ad_comm::adprogram_name + adstring(".mcm");
01259     return;
01260   }
01261   if (!ifs)
01262   {
01263     cerr << "Error trying to read number of simulations from file "
01264          << ad_comm::adprogram_name + adstring(".mcm");
01265     return;
01266   }
01267   for (i=1;i<=nsim;i++)
01268   {
01269     float ftmp = 0.0;
01270     int ii;
01271     int mmin=mcmc_values.indexmin();
01272     int mmax=mcmc_values.indexmax();
01273     for (ii=mmin;ii<=mmax;ii++)
01274     {
01275       ifs >>  ftmp;
01276       mcmc_values(ii)=double(ftmp);
01277     }
01278     ifs >>  ftmp;
01279     llc=double(ftmp);
01280     for (ii=1;ii<=nsdvars;ii++)
01281     {
01282       int x;
01283       double xx=(mcmc_values(ii)-m(ii))/h(ii);
01284       if (xx>0.0)
01285         x=int (xx+0.5);
01286       else
01287         x=int(xx-0.5);
01288       if (x<-nslots)
01289       {
01290         x=-nslots;
01291       }
01292       if (x>nslots)
01293       {
01294         x=nslots;
01295       }
01296       if (!probflag)
01297       {
01298         hist(ii,x)+=1;
01299       }
01300       else
01301       {
01302         hist(ii,x)+=exp(llc);
01303       }
01304     }
01305     if (!ifs)
01306     {
01307       cerr << "Error trying to read data from file "
01308          << ad_comm::adprogram_name + adstring(".mcm");
01309       return;
01310     }
01311   }
01312 }
01313 
01314 void read_covariance_matrix(const dmatrix& S,int nvar,int& oldHbf,
01315   dvector & sscale)
01316 {
01317   adstring tmpstring="admodel.cov";
01318   if (ad_comm::wd_flag)
01319     tmpstring = ad_comm::adprogram_name + ".cov";
01320   uistream cif((char*)tmpstring);
01321   if (!cif)
01322   {
01323     cerr << "Error trying to open file " << tmpstring
01324          << "  for reading" << endl;
01325   }
01326   int tmp_nvar = 0;
01327   cif >> tmp_nvar;
01328   if (nvar !=tmp_nvar)
01329   {
01330     cerr << "Incorrect number of independent variables in file"
01331       " model.cov" << endl;
01332     exit(1);
01333   }
01334   cif >> S;
01335   if (!cif)
01336   {
01337     cerr << "error reading covariance matrix from " << tmpstring
01338          << endl;
01339     exit(1);
01340   }
01341   cif >> oldHbf;
01342   if (!cif)
01343   {
01344     cerr << "error reading old_hybrid_bounded_flag from " << tmpstring
01345          << endl;
01346     exit(1);
01347   }
01348   cif >> sscale;
01349   if (!cif)
01350   {
01351     cerr << "error reading sscale from " << tmpstring
01352          << endl;
01353     exit(1);
01354   }
01355 }
01356 void read_hessian_matrix_and_scale(int nvar, const dmatrix& _SS,
01357   const dvector& pen_vector)
01358 {
01359   dmatrix& S= (dmatrix&) _SS;
01360   adstring tmpstring="admodel.hes";
01361   if (ad_comm::wd_flag)
01362      tmpstring = ad_comm::adprogram_name + ".hes";
01363   uistream cif((char*)tmpstring);
01364 
01365   if (!cif)
01366   {
01367     cerr << "Error trying to open file " << tmpstring
01368          << "  for reading" << endl;
01369   }
01370   int tmp_nvar = 0;
01371   cif >> tmp_nvar;
01372   if (nvar !=tmp_nvar)
01373   {
01374     cerr << "Incorrect number of independent variables in file"
01375       " admodel.hes" << endl;
01376     exit(1);
01377   }
01378   cif >> S;
01379   if (!cif)
01380   {
01381     cerr << "error reading covariance matrix from model.cov" << endl;
01382     exit(1);
01383   }
01384   cifstream cifs((char*)(ad_comm::adprogram_name + adstring(".bvs")));
01385   dvector tmp(1,nvar);
01386   cifs >> tmp;
01387   dvector wts=pen_vector/.16;
01388   dvector diag_save(1,nvar);
01389   //int neg_flag;
01390   //double base=5.0;
01391   double dmin=min(eigenvalues(S));
01392   cout << "Smallest eigenvalue = " << dmin << endl;
01393   for (int i=1;i<=nvar;i++)
01394   {
01395     if (tmp(i)>0)
01396     {
01397 #if defined(USE_DDOUBLE)
01398       S(i,i)/=pow(double(10.0),tmp(i));
01399 #else
01400       S(i,i)/=pow(10.0,tmp(i));
01401 #endif
01402     }
01403   }
01404   dmin=min(eigenvalues(S));
01405   if (dmin<0)
01406   {
01407     cout << "Smallest eigenvalue = " << dmin << endl;
01408     exit(1);
01409   }
01410   /*
01411   do
01412   {
01413     cout << wts << endl << endl;
01414     for (int i=1;i<=nvar;i++)
01415     {
01416       diag_save(i)=S(i,i);
01417       if (wts(i)>0)
01418       {
01419         S(i,i)/=pow(base,wts(i));
01420         cout << "  wts(" << i << ") = " << wts(i) << endl;
01421       }
01422     }
01423     dmin=min(eigenvalues(S));
01424     if (dmin<0)
01425     {
01426       cout << "Smallest eigenvalue = " << dmin << endl;
01427       neg_flag=1;
01428       base=sqrt(base);
01429       cout << "base = " << base << endl;
01430       for (int i=1;i<=nvar;i++)
01431       {
01432         S(i,i)=diag_save(i);
01433       }
01434       dmin=min(eigenvalues(S));
01435       cout << "XX Smallest eigenvalue = " << dmin << endl;
01436     }
01437     else
01438     {
01439       neg_flag=0;
01440     }
01441   }
01442   while (neg_flag);
01443  */
01444   S=inv(S);
01445 }
01446 
01447 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS,
01448   double rbp,int mcmc2_flag)
01449 {
01450   dmatrix& S= (dmatrix&) _SS;
01451   adstring tmpstring="admodel.hes";
01452   if (mcmc2_flag)
01453   {
01454     tmpstring = ad_comm::adprogram_name + ".bgs";
01455   }
01456   else
01457   {
01458     if (ad_comm::wd_flag)
01459        tmpstring = ad_comm::adprogram_name + ".hes";
01460   }
01461   uistream cif((char*)tmpstring);
01462 
01463   if (!cif)
01464   {
01465     cerr << "Error trying to open file " << tmpstring
01466          << "  for reading" << endl;
01467   }
01468   int tmp_nvar = 0;
01469   cif >> tmp_nvar;
01470   if (nvar !=tmp_nvar)
01471   {
01472     cerr << "Incorrect number of independent variables in file"
01473       " admodel.hes" << endl;
01474     exit(1);
01475   }
01476   cif >> S;
01477   if (!cif)
01478   {
01479     cerr << "error reading covariance matrix from model.cov" << endl;
01480     exit(1);
01481   }
01482   dmatrix Save=1*S;
01483   // for mcmc2 option Hessian is already inverted.
01484   if (mcmc2_flag==0)
01485   {
01486     S=inv(S);
01487   }
01488   int mmin=S.indexmin();
01489   int mmax=S.indexmax();
01490   dvector diag(mmin,mmax);
01491   int i,j;
01492   for (i=mmin;i<=mmax;i++)
01493   {
01494     diag(i)=sqrt(S(i,i));
01495   }
01496   for (i=mmin;i<=mmax;i++)
01497     for (j=mmin;j<=mmax;j++)
01498       S(i,j)/=diag(i)*diag(j);
01499 
01500   dmatrix ch=choleski_decomp(S);
01501   ofstream ofs("corrtest");
01502   ofs << "corr matrix" << endl;
01503   ofs << S << endl;
01504   ofs << "choleski decomp of corr matrix" << endl;
01505   ofs << ch << endl;
01506   dmatrix tmp(mmin,mmax,mmin,mmax);
01507 
01508   for (i=mmin;i<=mmax;i++)
01509     tmp(i)=ch(i)/norm(ch(i));
01510   ofs << "tmp" << endl;
01511   ofs << tmp << endl;
01512 
01513   for (i=mmin;i<=mmax;i++)
01514     tmp(i)/=tmp(i,i);
01515 
01516 
01517   if (rbp<=0.0 || rbp >= 1.0)
01518     rbp=0.5;
01519   for (i=mmin;i<=mmax;i++)
01520   {
01521     for (j=mmin;j<=mmax;j++)
01522     {
01523       if (tmp(i,j)>=1)
01524          tmp(i,j)=pow(tmp(i,j),rbp);
01525       else if (tmp(i,j)<-1)
01526          tmp(i,j)=-pow(-tmp(i,j),rbp);
01527     }
01528   }
01529 
01530   for (i=mmin;i<=mmax;i++)
01531     tmp(i)/=norm(tmp(i));
01532 
01533   dmatrix T=tmp*trans(tmp);
01534 
01535   ofs << "T-S" << endl;
01536   ofs << T-S << endl;
01537 
01538   S=T;
01539   ofs << "modified corr matrix" << endl;
01540   ofs << S << endl;
01541   for (i=mmin;i<=mmax;i++)
01542     for (j=mmin;j<=mmax;j++)
01543       S(i,j)*=diag(i)*diag(j);
01544 
01545   ofs << "modified S" << endl;
01546   ofs << S << endl;
01547 
01548   ofs << "S* modified S" << endl;
01549   ofs << S*Save << endl;
01550 }
01551 
01552 int user_stop(void)
01553 {
01554   int quit_flag=0;
01555 #if defined(_MSC_VER)
01556   if (kbhit())
01557 #else
01558   if(ctlc_flag)
01559 #endif
01560   {
01561 #if defined(__DJGPP__)
01562     int c = toupper(getxkey());
01563 #else
01564     int c = toupper(getch());
01565 #endif
01566     if (c == 'Q')
01567     {
01568       quit_flag=1;
01569     }
01570   }
01571   return quit_flag;
01572 }
01573 
01574 /*
01575 void newton_raftery_bayes_estimate(double cbf, int ic, const dvector& lk,
01576   double d)
01577 {
01578   double d1=1.0-d;
01579   double cbold=cbf;
01580   do
01581   {
01582     cout << "initial bayes factor" << cbf << endl;
01583     double sum1=0;
01584     double sum2=0;
01585 
01586     for (int i=1;i<=ic;i++)
01587     {
01588       sum1+=1./(d1/cbf+d/lk(i));
01589 
01590       sum2+=1./(d1/cbf*lk(i)+d);
01591     }
01592     sum1+=d*ic*cbf;
01593     sum2+=d*ic;
01594 
01595     cbf=sum1/sum2;
01596     double diff=log(cbf)-log(cbold);
01597     if (fabs(diff)<1.e-5) break;
01598     cbold=cbf;
01599  }
01600  while(1);
01601 }
01602 */
01603 
01604 void newton_raftery_bayes_estimate_new(double lcbf, int ic, const dvector& lk,
01605   double d)
01606 {
01607   double d1=1.0-d;
01608   double lcbold=lcbf;
01609   do
01610   {
01611     cout << "initial log bayes factor" << lcbf << endl;
01612     double sum1=0;
01613     double sum2=0;
01614 
01615     for (int i=1;i<=ic;i++)
01616     {
01617       double dtmp=exp(lcbf-lk(i));
01618       sum1+=1./(d1+d*dtmp);
01619       sum2+=1./(d1/dtmp+d);
01620     }
01621     sum1+=d*ic;
01622     sum2+=d*ic;
01623     lcbf=lcbf+log(sum1)-log(sum2);
01624     double diff=lcbf-lcbold;
01625     if (fabs(diff)<1.e-5) break;
01626     lcbold=lcbf;
01627   }
01628   while(1);
01629 }
01630 
01631 /*
01632 void save_mcmc_for_gui(const dvector& mcmc_number_values,
01633   dmatrix &mdm,int& ids)
01634 {
01635   int imax=mdm.colmax();
01636   if (ids>imax) ids=1;
01637   int rmax=mcmc_number_values.indexmax();
01638   for (int i=1;i<=rmax;i++)
01639     mdm(i,ids)=mcmc_number_values(i);
01640   ids++;
01641 }
01642 */
01643 
01644 void save_mcmc_for_gui1(const dvector& mcmc_values,
01645   dmatrix &mdm,int& ids,int& iwrap,ivector& no)
01646 {
01647   int rmax=no.indexmax();
01648   int imax=mdm.colmax();
01649   if (ids>imax)
01650   {
01651     ids=1;
01652     iwrap=1;
01653   }
01654   for (int i=1;i<=rmax;i++)
01655     mdm(i,ids)=mcmc_values(no(i));
01656   ids++;
01657 }
01658 
01659 dvector read_old_scale(int & old_nvar)
01660 {
01661   adstring tmpstring="admodel.cov";
01662   if (ad_comm::wd_flag)
01663     tmpstring = ad_comm::adprogram_name + ".cov";
01664   uistream cif((char*)tmpstring);
01665   if (!cif)
01666   {
01667     cerr << "Error trying to open file " << tmpstring
01668          << "  for reading" << endl;
01669   }
01670   cif >> old_nvar;
01671   dmatrix S(1,old_nvar,1,old_nvar);
01672   cif >> S;
01673   if (!cif)
01674   {
01675     cerr << "error reading covariance matrix from " << tmpstring
01676          << endl;
01677     exit(1);
01678   }
01679   int oldHbf;
01680   cif >> oldHbf;
01681   if (!cif)
01682   {
01683     cerr << "error reading old_hybrid_bounded_flag from " << tmpstring
01684          << endl;
01685     exit(1);
01686   }
01687   dvector sscale(1,old_nvar);
01688   cif >> sscale;
01689   if (!cif)
01690   {
01691     cerr << "error reading sscale from " << tmpstring
01692          << endl;
01693     exit(1);
01694   }
01695   return sscale;
01696 }