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