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