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