ADMB Documentation  11.5.3287
 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     int nvar_re=initial_params::nvarcalc();
00169 
00170     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00171     //int scov_option=0;
00172     dmatrix s_covar;
00173     dvector s_mean;
00174     int on=-1;
00175     int nslots=800;
00176     //int nslots=3600;
00177     int initial_nsim=4800;
00178     int ncor=0;
00179     double bfsum=0;
00180     int ibfcount=0;
00181     double llbest;
00182     double lbmax;
00183 
00184     //int ntmp=0;
00185     //if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcscov",ntmp))>-1)
00186     //{
00187     //scov_option=1;
00188     s_covar.allocate(1,nvar,1,nvar);
00189     s_mean.allocate(1,nvar);
00190     s_mean.initialize();
00191     s_covar.initialize();
00192 
00193     int ndvar=stddev_params::num_stddev_calc();
00194     /*
00195     int numdvar=stddev_params::num_stddev_number_calc();
00196     if (adjm_ptr)
00197     {
00198       mcmc_display_matrix.allocate(1,numdvar,1,mcmc_gui_length);
00199       number_offsets.allocate(1,numdvar);
00200       number_offsets=stddev_params::copy_all_number_offsets();
00201     }
00202     */
00203     if (mcmc2_flag==0)
00204     {
00205       initial_params::set_inactive_random_effects();
00206       nvar=initial_params::nvarcalc(); // get the number of active parameters
00207     }
00208     dvector x(1,nvar);
00209     dvector scale(1,nvar);
00210     dmatrix values;
00211     int have_hist_flag=0;
00212     initial_params::xinit(x);
00213     dvector pen_vector(1,nvar);
00214     {
00215       initial_params::reset(dvar_vector(x),pen_vector);
00216       //cout << pen_vector << endl << endl;
00217     }
00218 
00219     initial_params::mc_phase=0;
00220     initial_params::stddev_scale(scale,x);
00221     initial_params::mc_phase=1;
00222     dvector bmn(1,nvar);
00223     dvector mean_mcmc_values(1,ndvar);
00224     dvector s(1,ndvar);
00225     dvector h(1,ndvar);
00226     //dvector h;
00227     dvector square_mcmc_values(1,ndvar);
00228     square_mcmc_values.initialize();
00229     mean_mcmc_values.initialize();
00230     bmn.initialize();
00231     int use_empirical_flag=0;
00232     int diag_option=0;
00233     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcdiag"))>-1)
00234     {
00235       diag_option=1;
00236       cout << " Setting covariance matrix to diagonal with entries " << dscale
00237            << endl;
00238     }
00239     dmatrix S(1,nvar,1,nvar);
00240     dvector sscale(1,nvar);
00241     if (!diag_option)
00242     {
00243       int nopt = 0;
00244       int rescale_bounded_flag=0;
00245       double rescale_bounded_power=0.5;
00246       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcrb",nopt))>-1)
00247       {
00248         if (nopt)
00249         {
00250           int iii=atoi(ad_comm::argv[on+1]);
00251           if (iii < 1 || iii > 9)
00252           {
00253             cerr << " -mcrb argument must be integer between 1 and 9 --"
00254                     " using default of 5" << endl;
00255             rescale_bounded_power=0.5;
00256           }
00257           else
00258             rescale_bounded_power=iii/10.0;
00259         }
00260         else
00261         {
00262           rescale_bounded_power=0.5;
00263         }
00264         rescale_bounded_flag=1;
00265       }
00266       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcec"))>-1)
00267       {
00268         use_empirical_flag=1;
00269       }
00270       if (use_empirical_flag)
00271       {
00272         read_empirical_covariance_matrix(nvar,S,ad_comm::adprogram_name);
00273       }
00274       else if (!rescale_bounded_flag)
00275       {
00276         int hbf;
00277         if (mcmc2_flag==0)
00278         {
00279           read_covariance_matrix(S,nvar,hbf,sscale);
00280         }
00281         else
00282         {
00283           int tmp_nvar = 0;
00284           adstring tmpstring = ad_comm::adprogram_name + ".bgs";
00285           uistream uis((char*)(tmpstring));
00286           if (!uis)
00287           {
00288             cerr << "error opening file " << tmpstring << endl;
00289             ad_exit(1);
00290           }
00291           uis >> tmp_nvar;
00292           if (!uis)
00293           {
00294             cerr << "error reading from file " << tmpstring << endl;
00295             ad_exit(1);
00296           }
00297           if (tmp_nvar != nvar)
00298           {
00299             cerr << "size error reading from " <<  tmpstring << endl;
00300             ad_exit(1);
00301           }
00302           uis >> S;
00303           if (!uis)
00304           {
00305             cerr << "error reading from file " << tmpstring << endl;
00306             ad_exit(1);
00307           }
00308         }
00309       }
00310       else
00311       {
00312         read_hessian_matrix_and_scale1(nvar,S,rescale_bounded_power,
00313           mcmc2_flag);
00314         //read_hessian_matrix_and_scale(nvar,S,pen_vector);
00315       }
00316 
00317       {  // scale covariance matrix for model space
00318         dmatrix tmp(1,nvar,1,nvar);
00319         for (int i=1;i<=nvar;i++)
00320         {
00321           tmp(i,i)=S(i,i)*(scale(i)*scale(i));
00322           for (int j=1;j<i;j++)
00323           {
00324             tmp(i,j)=S(i,j)*(scale(i)*scale(j));
00325             tmp(j,i)=tmp(i,j);
00326           }
00327         }
00328         S=tmp;
00329       }
00330     }
00331     else
00332     {
00333       S.initialize();
00334       for (int i=1;i<=nvar;i++)
00335       {
00336         S(i,i)=dscale;
00337       }
00338     }
00339 
00340     cout << sort(eigenvalues(S)) << endl;
00341     dmatrix chd = choleski_decomp( (dscale*2.4/sqrt(double(nvar))) * S);
00342     dmatrix chdinv=inv(chd);
00343 
00344     dmatrix symbds(1,2,1,nvar);
00345     initial_params::set_all_simulation_bounds(symbds);
00346     ofstream ofs_sd1((char*)(ad_comm::adprogram_name + adstring(".mc2")));
00347 
00348     {
00349       int number_sims = 100000;
00350       if (nmcmc > 0)
00351       {
00352         number_sims = nmcmc;
00353       }
00354       int iseed=0;
00355       //cin >> iseed;
00356       if (iseed0<=0)
00357       {
00358         iseed=-36519;
00359       }
00360       else
00361       {
00362         iseed=-iseed0;
00363       }
00364       if (iseed>0)
00365       {
00366         iseed=-iseed;
00367       }
00368       cout << "Initial seed value " << iseed << endl;
00369       random_number_generator rng(iseed);
00370       rng.better_rand();
00371       double lprob=0.0;
00372       double lpinv=0.0;
00373       // get lower and upper bounds
00374 
00375       independent_variables y(1,nvar);
00376       independent_variables parsave(1,nvar_re);
00377       initial_params::restore_start_phase();
00378 
00379       // read in the mcmc values to date
00380       int ii=1;
00381       dmatrix hist;
00382       if (restart_flag)
00383       {
00384         int tmp=0;
00385         if (!no_sd_mcmc) {
00386           hist.allocate(1,ndvar,-nslots,nslots);
00387           tmp=read_hist_data(hist,h,mean_mcmc_values,s,parsave,iseed,
00388             size_scale);
00389           values.allocate(1,ndvar,-nslots,nslots);
00390           for (int i=1;i<=ndvar;i++)
00391           {
00392             values(i).fill_seqadd(mean_mcmc_values(i)-0.5*total_spread*s(i)
00393               +.5*h(i),h(i));
00394           }
00395         }
00396         if (iseed>0)
00397         {
00398           iseed=-iseed;
00399         }
00400         /*double br=*/rng.better_rand();
00401         if (tmp) have_hist_flag=1;
00402         chd=size_scale*chd;
00403         chdinv=chdinv/size_scale;
00404       }
00405       else
00406       {
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_re;
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 #ifndef OPT_LIB
00858       assert(isim != 0);
00859 #endif
00860       cout << iac/double(isim) << endl;
00861       initial_params::mc_phase=0;
00862      /*
00863       if (adjm_ptr)
00864       {
00865         ad_update_mcmc_report(mcmc_display_matrix,mcmc_save_index,
00866           mcmc_wrap_flag,1);
00867 
00868         if (allocated(hist)) ad_update_mcmchist_report(hist,
00869           number_offsets,mean_mcmc_values,h,1);
00870       }
00871       */
00872     }
00873 
00874     write_empirical_covariance_matrix(ncor,s_mean,s_covar,
00875       ad_comm::adprogram_name);
00876     //newton_raftery_bayes_estimate(current_bf,ibfcount,exp(lkvector),.01);
00877 #if defined(USE_BAYES_FACTORS)
00878     cout << "log current bayes factor " << lcurrent_bf << endl;
00879     newton_raftery_bayes_estimate_new(lcurrent_bf,ibfcount,lkvector,.01);
00880 #endif
00881     if (pofs_psave)
00882     {
00883       delete pofs_psave;
00884       pofs_psave=NULL;
00885     }
00886   }
00887 }
00888 
00899 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00900   const dmatrix& s_covar,adstring& prog_name)
00901 {
00902   uostream ofs((char*)(ad_comm::adprogram_name + adstring(".ecm")));
00903   dvector tmp=s_mean/ncor;
00904   int nvar=s_mean.indexmax();
00905   //ofs << ncor << " " << nvar << endl;
00906   dmatrix sigma=s_covar/ncor -outer_prod(tmp,tmp);
00907   cout << "In write empirical covariance matrix" << endl;
00908   cout << sort(eigenvalues(sigma)) << endl;
00909   dvector std(1,nvar);
00910   ofs << sigma;
00911   /*
00912   int i;
00913   for (i=1;i<=nvar;i++)
00914   {
00915     std(i)=sigma(i,i);
00916     for (int j=1;j<=i;j++)
00917     {
00918       sigma(i,j)/=sqrt(std(i)*std(j));
00919       sigma(j,i)=sigma(i,j);
00920     }
00921   }
00922   for (i=1;i<=nvar;i++)
00923   {
00924     ofs << setw(10) << setscientific() << std(i);
00925     for (int j=1;j<=i;j++)
00926     {
00927       ofs << " " << setfixed() << setw(7) << setprecision(4)
00928            << sigma(i,j);
00929     }
00930     ofs << endl;
00931   }
00932  */
00933 }
00934 
00943 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00944   const adstring& prog_name)
00945 {
00946   adstring infile_name=ad_comm::adprogram_name + adstring(".ecm");
00947   uistream ifs((char*)(infile_name));
00948   if (!ifs)
00949   {
00950     cerr << "Error opening file " << infile_name << endl;
00951   }
00952   ifs  >> S;
00953   /*
00954   int ncor=0;
00955   int _nvar=0;
00956   ifs >> ncor >> _nvar;
00957   if (nvar != _nvar)
00958   {
00959     cerr << "wromng number of indepdendent variables in routine"
00960       " read_empirical_covariance_matrix" << endl;
00961   }
00962   dvector std(1,nvar);
00963   int i;
00964   for (i=1;i<=nvar;i++)
00965   {
00966     ifs >> std(i);
00967     for (int j=1;j<=i;j++)
00968     {
00969       ifs  >> S(i,j);
00970       S(j,i)=S(i,j);
00971     }
00972   }
00973   if (!ifs)
00974   {
00975     cerr << "Error reading from file " << infile_name << endl;
00976   }
00977   for (i=1;i<=nvar;i++)
00978   {
00979     for (int j=1;j<=i;j++)
00980     {
00981       S(i,j)*=sqrt(std(i)*std(j));
00982       S(j,i)=S(i,j);
00983     }
00984   }
00985  */
00986 }
00987 
00988 void print_hist_data(const dmatrix& hist, const dmatrix& values,
00989   const dvector& h, dvector& m, const dvector& s, const dvector& parsave,
00990   int iseed, double size_scale)
00991 {
00992   ofstream ofs((char*)(ad_comm::adprogram_name + adstring(".hst")));
00993   int nsdvars=stddev_params::num_stddev_calc();
00994   adstring_array param_labels(1,nsdvars);
00995   ivector param_size(1,nsdvars);
00996   int ii=1;
00997   //int max_name_length=0;
00998   int i;
00999   for (i=0;i< stddev_params::num_stddev_params;i++)
01000   {
01001     param_labels[ii]=
01002       stddev_params::stddevptr[i]->label();
01003     param_size[ii]=
01004       stddev_params::stddevptr[i]->size_count();
01005 /*
01006     if (max_name_length<param_labels[ii].size())
01007     {
01008       max_name_length=param_labels[ii].size();
01009     }
01010 */
01011     ii++;
01012   }
01013   //int end_stdlabels=ii-1;
01014 
01015   int lc=1;
01016   int ic=1;
01017   ivector mmin(1,nsdvars);
01018   ivector mmax(1,nsdvars);
01019 
01020   for (i=1;i<=nsdvars;i++)
01021   {
01022     mmin(i)=minnz(hist(i));
01023     mmax(i)=maxnz(hist(i));
01024   }
01025 #ifdef OPT_LIB
01026   int nsim = (int)sum(hist(1));
01027 #else
01028   double _nsim = sum(hist(1));
01029   assert(_nsim <= (double)INT_MAX);
01030   int nsim = (int)_nsim;
01031 #endif
01032   ofs << "# samples sizes" << endl;
01033   ofs << nsim << endl;
01034   ofs << "# step size scaling factor" << endl;
01035   ofs << size_scale << endl;
01036   ofs << "# step sizes" << endl;
01037   ofs << h << endl;
01038   ofs << "# means" << endl;
01039   ofs << m << endl;
01040   ofs << "# standard devs" << endl;
01041   ofs << s << endl;
01042   ofs << "# lower bounds" << endl;
01043   ofs << mmin << endl;
01044   ofs << "# upper bounds" << endl;
01045   ofs << mmax<< endl;
01046   ofs << "#number of parameters" << endl;
01047   ofs << parsave.indexmax() << endl;
01048   ofs << "#current parameter values for mcmc restart" << endl;
01049   ofs << parsave << endl;
01050   ofs << "#random number seed" << endl;
01051   ofs << iseed << endl;
01052   for (i=1;i<=nsdvars;i++)
01053   {
01054     ofs << "#" << param_labels[lc];
01055     if (param_size[lc]>1)
01056     {
01057       ofs << "[" << ic << "]";
01058     }
01059     ofs << endl;
01060 
01061     if (++ic> param_size[lc])
01062     {
01063       lc++;
01064       ic=1;
01065     }
01066     for (ii=mmin(i);ii<=mmax(i);ii++)
01067     {
01068       ofs << values(i,ii) << " " << hist(i,ii)/(nsim*h(i)) << endl;
01069     }
01070     if (i<nsdvars) ofs << endl;
01071   }
01072 }
01073 
01074 int read_hist_data(const dmatrix& _hist, const dvector& h, dvector& m,
01075   const dvector& s, const dvector& parsave, int& iseed,
01076   const double& size_scale)
01077 {
01078   dmatrix& hist= (dmatrix&) _hist;
01079   cifstream cif((char*)(ad_comm::adprogram_name + adstring(".hst")));
01080   int nsdvars=stddev_params::num_stddev_calc();
01081   int nsim=0;
01082   int ii=1;
01083   int i;
01084   ivector mmin(1,nsdvars);
01085   ivector mmax(1,nsdvars);
01086   //ofs << # samples sizes << endl;
01087   cif >> nsim;
01088   cif >> size_scale;
01089   //ofs << # step sizes << endl;
01090   cif >> h;
01091   //ofs << # means << endl;
01092   cif >> m;
01093   //ofs << # standard devs << endl;
01094   cif >> s;
01095   //ofs << # lower bounds << endl;
01096   cif >> mmin;
01097   //ofs << # upper bounds << endl;
01098   cif >> mmax;
01099   int num_vars=0;
01100   cif >> num_vars;
01101   int flag=1;
01102   if (num_vars!=parsave.indexmax())
01103   {
01104     cerr << "wrong number of independent variables in file"
01105        << ad_comm::adprogram_name + adstring(".hst") << endl;
01106     flag=0;
01107   }
01108   if (flag)
01109   {
01110   cif >> parsave;
01111     cif >> iseed;
01112     double tmp=0.0;
01113     hist.initialize();
01114     for (i=1;i<=nsdvars;i++)
01115     {
01116       for (ii=mmin(i);ii<=mmax(i);ii++)
01117       {
01118         cif >> tmp >> hist(i,ii);
01119       }
01120       hist(i)*=nsim*h(i);
01121     }
01122   }
01123   return flag;
01124 }
01125 
01126 int minnz(const dvector& x)
01127 {
01128   int mmin=x.indexmin();
01129   int mmax=x.indexmax();
01130   int m=mmin;
01131   for (int ii=mmin;ii<=mmax;ii++)
01132   {
01133     if (!ISZERO(x(ii)))
01134     {
01135       m=ii;
01136       if (m>mmin) m--;
01137       break;
01138     }
01139   }
01140   return m;
01141 }
01142 
01143 int maxnz(const dvector& x)
01144 {
01145   int mmin=x.indexmin();
01146   int mmax=x.indexmax();
01147   int m=mmax;
01148   for (int ii=mmax;ii>=mmin;ii--)
01149   {
01150     if (!ISZERO(x(ii)))
01151     {
01152       m=ii;
01153       if (m<mmax) m++;
01154       break;
01155     }
01156   }
01157   return m;
01158 }
01159 
01160 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& _hist,
01161   dvector& mcmc_values, double llc, const dvector& h, int nslots,
01162   double total_spread,int probflag)
01163 {
01164   dmatrix& hist= (dmatrix&) _hist;
01165   int nsdvars=stddev_params::num_stddev_calc();
01166   for (int ii=1;ii<=nsdvars;ii++)
01167   {
01168     int x;
01169     double xx=(mcmc_values(ii)-m(ii))/h(ii);
01170     if (xx>0.0)
01171       x=int (xx+0.5);
01172     else
01173       x=int(xx-0.5);
01174 
01175     if (x<-nslots)
01176     {
01177       x=-nslots;
01178     }
01179     if (x>nslots)
01180     {
01181       x=nslots;
01182     }
01183     if (!probflag)
01184     {
01185       hist(ii,x)+=1;
01186     }
01187     else
01188     {
01189       hist(ii,x)+=exp(llc);
01190     }
01191   }
01192 }
01193 
01194 /*
01195 void add_guihist_values(const dvector& s, const dvector& m,
01196   const dmatrix& _hist,dvector& mcmcnumber_values,double llc,
01197   const dvector& h,int nslots,double total_spread)
01198 {
01199   dmatrix& hist= (dmatrix&) _hist;
01200   int nsdvars=stddev_params::num_stddev_number_calc();
01201   for (int ii=1;ii<=nsdvars;ii++)
01202   {
01203     int x=int((mcmcnumber_values(ii)-m(ii))/h(ii));
01204     //cout << "xxx = " << xxx << endl;
01205     char ch;
01206     //cin >> ch;
01207 
01208     if (x<1)
01209     {
01210       x=-nslots;
01211     }
01212     if (x>nslots)
01213     {
01214       x=nslots;
01215     }
01216     {
01217       hist(ii,x)+=1;
01218     }
01219   }
01220 }
01221 */
01222 
01223 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
01224   const dmatrix& _values, dmatrix& hist, const dvector& _h, int nslots,
01225   double total_spread, int probflag)
01226 {
01227   ADUNCONST(dmatrix,values)
01228   dvector& h = (dvector&) _h;
01229   int nsdvars=stddev_params::num_stddev_calc();
01230   dvector mcmc_values(1,nsdvars);
01231   values.allocate(1,nsdvars,-nslots,nslots);
01232   h=total_spread/(2*nslots+1)*s;
01233   hist.allocate(1,nsdvars,-nslots,nslots);
01234   hist.initialize();
01235   uistream ifs((char*)(ad_comm::adprogram_name + adstring(".mcm")));
01236   int i;
01237   double llc;
01238   for (i=1;i<=nsdvars;i++)
01239   {
01240     values(i).fill_seqadd(m(i)-.5*total_spread*s(i)+.5*h(i),h(i));
01241   }
01242 
01243   if (!ifs)
01244   {
01245     cerr << "Error trying to open file "
01246          << ad_comm::adprogram_name + adstring(".mcm");
01247     return;
01248   }
01249   if (!ifs)
01250   {
01251     cerr << "Error trying to read number of simulations from file "
01252          << ad_comm::adprogram_name + adstring(".mcm");
01253     return;
01254   }
01255   for (i=1;i<=nsim;i++)
01256   {
01257     float ftmp = 0.0;
01258     int ii;
01259     int mmin=mcmc_values.indexmin();
01260     int mmax=mcmc_values.indexmax();
01261     for (ii=mmin;ii<=mmax;ii++)
01262     {
01263       ifs >>  ftmp;
01264       mcmc_values(ii)=double(ftmp);
01265     }
01266     ifs >>  ftmp;
01267     llc=double(ftmp);
01268     for (ii=1;ii<=nsdvars;ii++)
01269     {
01270       int x;
01271       double xx=(mcmc_values(ii)-m(ii))/h(ii);
01272       if (xx>0.0)
01273         x=int (xx+0.5);
01274       else
01275         x=int(xx-0.5);
01276       if (x<-nslots)
01277       {
01278         x=-nslots;
01279       }
01280       if (x>nslots)
01281       {
01282         x=nslots;
01283       }
01284       if (!probflag)
01285       {
01286         hist(ii,x)+=1;
01287       }
01288       else
01289       {
01290         hist(ii,x)+=exp(llc);
01291       }
01292     }
01293     if (!ifs)
01294     {
01295       cerr << "Error trying to read data from file "
01296          << ad_comm::adprogram_name + adstring(".mcm");
01297       return;
01298     }
01299   }
01300 }
01301 
01302 void read_covariance_matrix(const dmatrix& S,int nvar,int& oldHbf,
01303   dvector & sscale)
01304 {
01305   adstring tmpstring="admodel.cov";
01306   if (ad_comm::wd_flag)
01307     tmpstring = ad_comm::adprogram_name + ".cov";
01308   uistream cif((char*)tmpstring);
01309   if (!cif)
01310   {
01311     cerr << "Error trying to open file " << tmpstring
01312          << "  for reading" << endl;
01313   }
01314   int tmp_nvar = 0;
01315   cif >> tmp_nvar;
01316   if (nvar !=tmp_nvar)
01317   {
01318     cerr << "Incorrect number of independent variables in file"
01319       " model.cov" << endl;
01320     exit(1);
01321   }
01322   cif >> S;
01323   if (!cif)
01324   {
01325     cerr << "error reading covariance matrix from " << tmpstring
01326          << endl;
01327     exit(1);
01328   }
01329   cif >> oldHbf;
01330   if (!cif)
01331   {
01332     cerr << "error reading old_hybrid_bounded_flag from " << tmpstring
01333          << endl;
01334     exit(1);
01335   }
01336   cif >> sscale;
01337   if (!cif)
01338   {
01339     cerr << "error reading sscale from " << tmpstring
01340          << endl;
01341     exit(1);
01342   }
01343 }
01344 void read_hessian_matrix_and_scale(int nvar, const dmatrix& _SS,
01345   const dvector& pen_vector)
01346 {
01347   dmatrix& S= (dmatrix&) _SS;
01348   adstring tmpstring="admodel.hes";
01349   if (ad_comm::wd_flag)
01350      tmpstring = ad_comm::adprogram_name + ".hes";
01351   uistream cif((char*)tmpstring);
01352 
01353   if (!cif)
01354   {
01355     cerr << "Error trying to open file " << tmpstring
01356          << "  for reading" << endl;
01357   }
01358   int tmp_nvar = 0;
01359   cif >> tmp_nvar;
01360   if (nvar !=tmp_nvar)
01361   {
01362     cerr << "Incorrect number of independent variables in file"
01363       " admodel.hes" << endl;
01364     exit(1);
01365   }
01366   cif >> S;
01367   if (!cif)
01368   {
01369     cerr << "error reading covariance matrix from model.cov" << endl;
01370     exit(1);
01371   }
01372   cifstream cifs((char*)(ad_comm::adprogram_name + adstring(".bvs")));
01373   dvector tmp(1,nvar);
01374   cifs >> tmp;
01375   dvector wts=pen_vector/.16;
01376   dvector diag_save(1,nvar);
01377   //int neg_flag;
01378   //double base=5.0;
01379   double dmin=min(eigenvalues(S));
01380   cout << "Smallest eigenvalue = " << dmin << endl;
01381   for (int i=1;i<=nvar;i++)
01382   {
01383     if (tmp(i)>0)
01384     {
01385 #if defined(USE_DDOUBLE)
01386       S(i,i)/=pow(double(10.0),tmp(i));
01387 #else
01388       S(i,i)/=pow(10.0,tmp(i));
01389 #endif
01390     }
01391   }
01392   dmin=min(eigenvalues(S));
01393   if (dmin<0)
01394   {
01395     cout << "Smallest eigenvalue = " << dmin << endl;
01396     exit(1);
01397   }
01398   /*
01399   do
01400   {
01401     cout << wts << endl << endl;
01402     for (int i=1;i<=nvar;i++)
01403     {
01404       diag_save(i)=S(i,i);
01405       if (wts(i)>0)
01406       {
01407         S(i,i)/=pow(base,wts(i));
01408         cout << "  wts(" << i << ") = " << wts(i) << endl;
01409       }
01410     }
01411     dmin=min(eigenvalues(S));
01412     if (dmin<0)
01413     {
01414       cout << "Smallest eigenvalue = " << dmin << endl;
01415       neg_flag=1;
01416       base=sqrt(base);
01417       cout << "base = " << base << endl;
01418       for (int i=1;i<=nvar;i++)
01419       {
01420         S(i,i)=diag_save(i);
01421       }
01422       dmin=min(eigenvalues(S));
01423       cout << "XX Smallest eigenvalue = " << dmin << endl;
01424     }
01425     else
01426     {
01427       neg_flag=0;
01428     }
01429   }
01430   while (neg_flag);
01431  */
01432   S=inv(S);
01433 }
01434 
01435 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS,
01436   double rbp,int mcmc2_flag)
01437 {
01438   dmatrix& S= (dmatrix&) _SS;
01439   adstring tmpstring="admodel.hes";
01440   if (mcmc2_flag)
01441   {
01442     tmpstring = ad_comm::adprogram_name + ".bgs";
01443   }
01444   else
01445   {
01446     if (ad_comm::wd_flag)
01447        tmpstring = ad_comm::adprogram_name + ".hes";
01448   }
01449   uistream cif((char*)tmpstring);
01450 
01451   if (!cif)
01452   {
01453     cerr << "Error trying to open file " << tmpstring
01454          << "  for reading" << endl;
01455   }
01456   int tmp_nvar = 0;
01457   cif >> tmp_nvar;
01458   if (nvar !=tmp_nvar)
01459   {
01460     cerr << "Incorrect number of independent variables in file"
01461       " admodel.hes" << endl;
01462     exit(1);
01463   }
01464   cif >> S;
01465   if (!cif)
01466   {
01467     cerr << "error reading covariance matrix from model.cov" << endl;
01468     exit(1);
01469   }
01470   dmatrix Save=1*S;
01471   // for mcmc2 option Hessian is already inverted.
01472   if (mcmc2_flag==0)
01473   {
01474     S=inv(S);
01475   }
01476   int mmin=S.indexmin();
01477   int mmax=S.indexmax();
01478   dvector diag(mmin,mmax);
01479   int i,j;
01480   for (i=mmin;i<=mmax;i++)
01481   {
01482     diag(i)=sqrt(S(i,i));
01483   }
01484   for (i=mmin;i<=mmax;i++)
01485     for (j=mmin;j<=mmax;j++)
01486       S(i,j)/=diag(i)*diag(j);
01487 
01488   dmatrix ch=choleski_decomp(S);
01489   ofstream ofs("corrtest");
01490   ofs << "corr matrix" << endl;
01491   ofs << S << endl;
01492   ofs << "choleski decomp of corr matrix" << endl;
01493   ofs << ch << endl;
01494   dmatrix tmp(mmin,mmax,mmin,mmax);
01495 
01496   for (i=mmin;i<=mmax;i++)
01497     tmp(i)=ch(i)/norm(ch(i));
01498   ofs << "tmp" << endl;
01499   ofs << tmp << endl;
01500 
01501   for (i=mmin;i<=mmax;i++)
01502     tmp(i)/=tmp(i,i);
01503 
01504 
01505   if (rbp<=0.0 || rbp >= 1.0)
01506     rbp=0.5;
01507   for (i=mmin;i<=mmax;i++)
01508   {
01509     for (j=mmin;j<=mmax;j++)
01510     {
01511       if (tmp(i,j)>=1)
01512          tmp(i,j)=pow(tmp(i,j),rbp);
01513       else if (tmp(i,j)<-1)
01514          tmp(i,j)=-pow(-tmp(i,j),rbp);
01515     }
01516   }
01517 
01518   for (i=mmin;i<=mmax;i++)
01519     tmp(i)/=norm(tmp(i));
01520 
01521   dmatrix T=tmp*trans(tmp);
01522 
01523   ofs << "T-S" << endl;
01524   ofs << T-S << endl;
01525 
01526   S=T;
01527   ofs << "modified corr matrix" << endl;
01528   ofs << S << endl;
01529   for (i=mmin;i<=mmax;i++)
01530     for (j=mmin;j<=mmax;j++)
01531       S(i,j)*=diag(i)*diag(j);
01532 
01533   ofs << "modified S" << endl;
01534   ofs << S << endl;
01535 
01536   ofs << "S* modified S" << endl;
01537   ofs << S*Save << endl;
01538 }
01539 
01540 int user_stop(void)
01541 {
01542   int quit_flag=0;
01543 #if defined(_MSC_VER)
01544   if (kbhit())
01545 #else
01546   if(ctlc_flag)
01547 #endif
01548   {
01549 #if defined(__DJGPP__)
01550     int c = toupper(getxkey());
01551 #else
01552     int c = toupper(getch());
01553 #endif
01554     if (c == 'Q')
01555     {
01556       quit_flag=1;
01557     }
01558   }
01559   return quit_flag;
01560 }
01561 
01562 /*
01563 void newton_raftery_bayes_estimate(double cbf, int ic, const dvector& lk,
01564   double d)
01565 {
01566   double d1=1.0-d;
01567   double cbold=cbf;
01568   do
01569   {
01570     cout << "initial bayes factor" << cbf << endl;
01571     double sum1=0;
01572     double sum2=0;
01573 
01574     for (int i=1;i<=ic;i++)
01575     {
01576       sum1+=1./(d1/cbf+d/lk(i));
01577 
01578       sum2+=1./(d1/cbf*lk(i)+d);
01579     }
01580     sum1+=d*ic*cbf;
01581     sum2+=d*ic;
01582 
01583     cbf=sum1/sum2;
01584     double diff=log(cbf)-log(cbold);
01585     if (fabs(diff)<1.e-5) break;
01586     cbold=cbf;
01587  }
01588  while(1);
01589 }
01590 */
01591 
01592 #if defined(USE_BAYES_FACTORS)
01593 void newton_raftery_bayes_estimate_new(double lcbf, int ic, const dvector& lk,
01594   double d)
01595 {
01596   double d1=1.0-d;
01597   double lcbold=lcbf;
01598   do
01599   {
01600     cout << "initial log bayes factor" << lcbf << endl;
01601     double sum1=0;
01602     double sum2=0;
01603 
01604     for (int i=1;i<=ic;i++)
01605     {
01606       double dtmp=exp(lcbf-lk(i));
01607       sum1+=1./(d1+d*dtmp);
01608       sum2+=1./(d1/dtmp+d);
01609     }
01610     sum1+=d*ic;
01611     sum2+=d*ic;
01612     lcbf=lcbf+log(sum1)-log(sum2);
01613     double diff=lcbf-lcbold;
01614     if (fabs(diff)<1.e-5) break;
01615     lcbold=lcbf;
01616   }
01617   while(1);
01618 }
01619 #endif
01620 
01621 /*
01622 void save_mcmc_for_gui(const dvector& mcmc_number_values,
01623   dmatrix &mdm,int& ids)
01624 {
01625   int imax=mdm.colmax();
01626   if (ids>imax) ids=1;
01627   int rmax=mcmc_number_values.indexmax();
01628   for (int i=1;i<=rmax;i++)
01629     mdm(i,ids)=mcmc_number_values(i);
01630   ids++;
01631 }
01632 
01633 void save_mcmc_for_gui1(const dvector& mcmc_values,
01634   dmatrix &mdm,int& ids,int& iwrap,ivector& no)
01635 {
01636   int rmax=no.indexmax();
01637   int imax=mdm.colmax();
01638   if (ids>imax)
01639   {
01640     ids=1;
01641     iwrap=1;
01642   }
01643   for (int i=1;i<=rmax;i++)
01644     mdm(i,ids)=mcmc_values(no(i));
01645   ids++;
01646 }
01647 */
01648 
01649 dvector read_old_scale(int & old_nvar)
01650 {
01651   adstring tmpstring="admodel.cov";
01652   if (ad_comm::wd_flag)
01653     tmpstring = ad_comm::adprogram_name + ".cov";
01654   uistream cif((char*)tmpstring);
01655   if (!cif)
01656   {
01657     cerr << "Error trying to open file " << tmpstring
01658          << "  for reading" << endl;
01659   }
01660   cif >> old_nvar;
01661   dmatrix S(1,old_nvar,1,old_nvar);
01662   cif >> S;
01663   if (!cif)
01664   {
01665     cerr << "error reading covariance matrix from " << tmpstring
01666          << endl;
01667     exit(1);
01668   }
01669   int oldHbf;
01670   cif >> oldHbf;
01671   if (!cif)
01672   {
01673     cerr << "error reading old_hybrid_bounded_flag from " << tmpstring
01674          << endl;
01675     exit(1);
01676   }
01677   dvector sscale(1,old_nvar);
01678   cif >> sscale;
01679   if (!cif)
01680   {
01681     cerr << "error reading sscale from " << tmpstring
01682          << endl;
01683     exit(1);
01684   }
01685   return sscale;
01686 }