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