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