ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
hybmcmc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: hybmcmc.cpp 2695 2014-11-19 20:13:10Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <sstream>
00012 using std::istringstream;
00013 
00014 #  include <df1b2fun.h>
00015 #  include <adrndeff.h>
00016 #include <admodel.h>
00017 
00018 #if defined(_MSC_VER)
00019   #include <conio.h>
00020 #endif
00021 
00022 #ifndef OPT_LIB
00023   #include <cassert>
00024   #include <climits>
00025 #endif
00026 
00027 double better_rand(long int&);
00028 void store_mcmc_values(const ofstream& ofs);
00029 void set_labels_for_mcmc(void);
00030 
00031 void print_hist_data(const dmatrix& hist, const dmatrix& values,
00032   const dvector& h, dvector& m, const dvector& s, const dvector& parsave,
00033   long int iseed, double size_scale);
00034 
00035 int minnz(const dvector& x);
00036 int maxnz(const dvector& xa);
00037 
00038 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS,double s,
00039   int mcmc2_flag);
00040 
00041 int read_hist_data(const dmatrix& hist, const dvector& h,
00042   dvector& m, const dvector& s, const dvector& parsave,long int& iseed,
00043   const double& size_scale);
00044 
00045 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
00046   const dmatrix& values, dmatrix& hist, const dvector& h,int slots,
00047   double total_spread,int probflag=0);
00048 
00049 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& hist,
00050   dvector& mcmc_values,double llc, const dvector& h,int nslots,
00051   double total_spreadd,int probflag=0);
00052 
00053 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00054   const dmatrix& s_covar, adstring& prog_name);
00055 
00056 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00057   const adstring& prog_name);
00058 
00059 void read_hessian_matrix_and_scale(int nvar, const dmatrix& S,
00060   const dvector& pen_vector);
00061 
00062 dvector new_probing_bounded_multivariate_normal(int nvar, const dvector& a1,
00063   const dvector& b1, dmatrix& ch, const double& wght,double pprobe,
00064   random_number_generator& rng);
00065 
00066 void new_probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00067   const dvector& b1, dmatrix& ch, const double& wght, const dvector& _y,
00068   double pprobe, random_number_generator& rng);
00069 
00070 //void newton_raftery_bayes_estimate(double cbf,int ic, const dvector& lk,
00071 //double d);
00072 
00073 void ad_update_mcmc_stats_report
00074   (int feval,int iter,double fval,double gmax);
00075 
00076 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00077   double fval,double gmax,const char * cbuf);
00078 void ad_update_mcmc_report(dmatrix& m,int i,int j,int ff=0);
00079 void ad_update_mcmchist_report(dmatrix& mcmc_values,ivector& number_offsets,
00080   dvector& mean_mcmc_values,dvector& h,int ff=0);
00081 
00086 void function_minimizer::hybrid_mcmc_routine(int nmcmc,int iseed0,double dscale,
00087   int restart_flag)
00088 {
00089   robust_hybrid_flag=0;
00090   uostream * pofs_psave=NULL;
00091   dmatrix mcmc_display_matrix;
00092   //int mcmc_save_index=1;
00093   //int mcmc_wrap_flag=0;
00094 
00095   int on2=-1;
00096   //int nvar1=0;
00097   if ( (on2=option_match(ad_comm::argc,ad_comm::argv,"-nosdmcmc"))>-1)
00098   {
00099     //int no_sd_mcmc = 1;
00100   }
00101   if (mcmc2_flag==1)
00102   {
00103     initial_params::restore_start_phase();
00104     //nvar1=initial_params::nvarcalc(); // get the number of active parameters
00105   }
00106 
00107   if (stddev_params::num_stddev_params==0)
00108   {
00109     cerr << " You must declare at least one object of type sdreport "
00110          << endl << " to do the mcmc calculations" << endl;
00111      return;
00112   }
00113   {
00114     //ofstream of_bf("testbf");
00115 
00116     //if (adjm_ptr) set_labels_for_mcmc();
00117 
00118     ivector number_offsets;
00119     dvector lkvector;
00120     //double current_bf=0;
00121     //double lcurrent_bf=0;
00122     //double size_scale=1.0;
00123     //double total_spread=200;
00124     //double total_spread=2500;
00125     //uostream * pofs_sd = NULL;
00126 
00127     initial_params::set_inactive_random_effects();
00128     //int nvar_x=initial_params::nvarcalc();
00129     initial_params::set_active_random_effects();
00130 
00131     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00132     dmatrix s_covar;
00133     dvector s_mean;
00134     //int ncsim=25000;
00135     //int nslots=800;
00136     //int nslots=3600;
00137     //int initial_nsim=4800;
00138     //int ntmp=0;
00139     //int ncor=0;
00140     //double bfsum=0;
00141     //int ibfcount=0;
00142     //double lbmax;
00143 
00144     s_covar.allocate(1,nvar,1,nvar);
00145     s_mean.allocate(1,nvar);
00146     s_mean.initialize();
00147     s_covar.initialize();
00148 
00149     int ndvar=stddev_params::num_stddev_calc();
00150     /*int numdvar=*/stddev_params::num_stddev_number_calc();
00151 
00152     if (mcmc2_flag==0)
00153     {
00154       initial_params::set_inactive_random_effects();
00155       nvar=initial_params::nvarcalc(); // get the number of active parameters
00156     }
00157 
00158     independent_variables parsave(1,nvar);
00159     initial_params::restore_start_phase();
00160 
00161     dvector x(1,nvar);
00162     //dvector scale(1,nvar);
00163     dmatrix values;
00164     //int have_hist_flag=0;
00165     initial_params::xinit(x);
00166     dvector pen_vector(1,nvar);
00167     {
00168       initial_params::reset(dvar_vector(x),pen_vector);
00169       //cout << pen_vector << endl << endl;
00170     }
00171 
00172     initial_params::mc_phase=0;
00173     //initial_params::stddev_scale(scale,x);
00174     initial_params::mc_phase=1;
00175     dvector bmn(1,nvar);
00176     dvector mean_mcmc_values(1,ndvar);
00177     dvector s(1,ndvar);
00178     dvector h(1,ndvar);
00179     //dvector h;
00180     dvector square_mcmc_values(1,ndvar);
00181     square_mcmc_values.initialize();
00182     mean_mcmc_values.initialize();
00183     bmn.initialize();
00184     int use_empirical_flag=0;
00185     int diag_option=0;
00186     //int topt=0;
00187     int hybnstep=10;
00188     double hybeps=0.1;
00189     double _hybeps=-1.0;
00190     int old_Hybrid_bounded_flag=-1;
00191 
00192     int on,nopt = 0;
00193     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hyeps",nopt))>-1)
00194     {
00195       if (!nopt)
00196       {
00197         cerr << "Usage -hyeps option needs number  -- ignored" << endl;
00198       }
00199       else
00200       {
00201         istringstream ist(ad_comm::argv[on+1]);
00202         ist >> _hybeps;
00203 
00204         if (_hybeps<=0)
00205         {
00206           cerr << "Usage -hyeps option needs positive number  -- ignored\n";
00207           _hybeps=-1.0;
00208         }
00209       }
00210     }
00211     if (_hybeps>0.0)
00212       hybeps=_hybeps;
00213 
00214     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hynstep",nopt))>-1)
00215     {
00216       if (nopt)
00217       {
00218         int iii=atoi(ad_comm::argv[on+1]);
00219         if (iii < 1 )
00220         {
00221           cerr << " -hynstep argument must be integer between > 0 --"
00222                   " using default of 10" << endl;
00223         }
00224         else
00225         {
00226           hybnstep=iii;
00227         }
00228       }
00229     }
00230 
00231     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcdiag"))>-1)
00232     {
00233       diag_option=1;
00234       cout << " Setting covariance matrix to diagonal with entries " << dscale
00235            << endl;
00236     }
00237     int mcrestart_flag=option_match(ad_comm::argc,ad_comm::argv,"-mcr");
00238     dmatrix S(1,nvar,1,nvar);
00239     dvector old_scale(1,nvar);
00240     int old_nvar;
00241     if (!diag_option)
00242     {
00243       int rescale_bounded_flag=0;
00244       double rescale_bounded_power=0.5;
00245       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcrb",nopt))>-1)
00246       {
00247         if (nopt)
00248         {
00249           int iii=atoi(ad_comm::argv[on+1]);
00250           if (iii < 1 || iii > 9)
00251           {
00252             cerr << " -mcrb argument must be integer between 1 and 9 --"
00253                     " using default of 5" << endl;
00254             rescale_bounded_power=0.5;
00255           }
00256           else
00257             rescale_bounded_power=iii/10.0;
00258         }
00259         else
00260         {
00261           rescale_bounded_power=0.5;
00262         }
00263         rescale_bounded_flag=1;
00264       }
00265       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcec"))>-1)
00266       {
00267         use_empirical_flag=1;
00268       }
00269       if (use_empirical_flag)
00270       {
00271         read_empirical_covariance_matrix(nvar,S,ad_comm::adprogram_name);
00272       }
00273       else if (!rescale_bounded_flag)
00274       {
00275         if (mcmc2_flag==0)
00276         {
00277           read_covariance_matrix(S,nvar,old_Hybrid_bounded_flag,old_scale);
00278         }
00279         else
00280         {
00281           int tmp_nvar = 0;
00282           adstring tmpstring = ad_comm::adprogram_name + ".bgs";
00283           uistream uis((char*)(tmpstring));
00284           if (!uis)
00285           {
00286             cerr << "error opening file " << tmpstring << endl;
00287             ad_exit(1);
00288           }
00289           uis >> tmp_nvar;
00290           if (!uis)
00291           {
00292             cerr << "error reading from file " << tmpstring << endl;
00293             ad_exit(1);
00294           }
00295           if (tmp_nvar != nvar)
00296           {
00297             cerr << "size error reading from " <<  tmpstring << endl;
00298             ad_exit(1);
00299           }
00300           uis >> S;
00301           if (!uis)
00302           {
00303             cerr << "error reading from file " << tmpstring << endl;
00304             ad_exit(1);
00305           }
00306           dvector tmp=read_old_scale(old_nvar);
00307           old_scale=1.0;
00308           old_scale(1,old_nvar)=tmp;
00309         }
00310       }
00311       else
00312       {
00313         read_hessian_matrix_and_scale1(nvar,S,rescale_bounded_power,
00314           mcmc2_flag);
00315         //read_hessian_matrix_and_scale(nvar,S,pen_vector);
00316       }
00317 
00318       // maybe we dont want to do this?
00319       /*
00320       {  // scale covariance matrix for model space
00321         dmatrix tmp(1,nvar,1,nvar);
00322         for (int i=1;i<=nvar;i++)
00323         {
00324           tmp(i,i)=S(i,i)*(scale(i)*scale(i));
00325           for (int j=1;j<i;j++)
00326           {
00327             tmp(i,j)=S(i,j)*(scale(i)*scale(j));
00328             tmp(j,i)=tmp(i,j);
00329           }
00330         }
00331         S=tmp;
00332       }
00333       */
00334     }
00335     else
00336     {
00337       S.initialize();
00338       for (int i=1;i<=nvar;i++)
00339       {
00340         S(i,i)=dscale;
00341       }
00342     }
00343 
00344     // for hybrid mcmc option always save output
00345     //if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcsave"))>-1)
00346     if ( mcrestart_flag>-1)
00347     {
00348       // check that nvar is correct
00349       uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00350       if (!uis)
00351       {
00352         cerr << "Error trying to open file" <<
00353           ad_comm::adprogram_name + adstring(".psv") <<
00354           " for mcrestart" <<   endl;
00355         ad_exit(1);
00356       } else {
00357         int nv1 = 0;
00358         uis >> nv1;
00359         if (nv1 !=nvar) {
00360           cerr << "wrong number of independent variables in" <<
00361             ad_comm::adprogram_name + adstring(".psv") << endl
00362            << " I found " << nv1 << " instead of " << nvar  << endl;
00363           ad_exit(1);
00364         }
00365         // get last x vector from file
00366 #ifndef OPT_LIB
00367         assert(parsave.size() >= 0);
00368 #endif
00369         std::streamoff sz = (unsigned int)parsave.size() * sizeof(double);
00370         // backup from end of file
00371         uis.seekg(-sz, ios::end);
00372         uis >> parsave;
00373        cout << " restored "  << parsave(parsave.indexmin()) << " "
00374             << parsave(parsave.indexmax()) << endl;
00375         if (!uis)
00376         {
00377           cerr << "error resotring last mcmc poistion from file "
00378             << ad_comm::adprogram_name + adstring(".psv") << endl;
00379           ad_exit(1);
00380         }
00381         int ii=1;
00382         dvector x0(1,nvar);
00383         initial_params::restore_all_values(parsave,ii);
00384         //x0.initialize();
00385         //dvector pen_vector(1,nvar);
00386         //initial_params::reset(dvar_vector(parsave),pen_vector);
00387         //initial_params::xinit(x0);
00388         //cout << " x0 " << x0(x0.indexmin()) << endl;
00389       }
00390     }
00391 
00392 
00393     if ( mcrestart_flag>-1)
00394     {
00395       pofs_psave= new uostream( (char*)(ad_comm::adprogram_name
00396           + adstring(".psv")),ios::app);
00397      /*
00398       pofs_psave->seekp(0,std::ios::end);
00399       *pofs_psave << 123;
00400       delete pofs_psave;
00401       uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00402       if (!uis)
00403       {
00404         cerr << "Error trying to open file" <<
00405           ad_comm::adprogram_name + adstring(".psv") <<
00406           " for mcrestart" <<   endl;
00407         ad_exit(1);
00408       } else {
00409         int nv1;
00410         uis >> nv1;
00411         if (nv1 !=nvar) {
00412           cerr << "wrong number of independent variables in" <<
00413             ad_comm::adprogram_name + adstring(".psv") << endl
00414            << " I found " << nv1 << " instead of " << nvar  << endl;
00415           ad_exit(1);
00416         }
00417       }
00418      */
00419     }
00420     else
00421     {
00422       pofs_psave=
00423         new uostream((char*)(ad_comm::adprogram_name + adstring(".psv")));
00424     }
00425 
00426     if (!pofs_psave|| !(*pofs_psave))
00427     {
00428       cerr << "Error trying to open file" <<
00429         ad_comm::adprogram_name + adstring(".psv") << endl;
00430       ad_exit(1);
00431     }
00432     if (mcrestart_flag == -1 )
00433     {
00434         (*pofs_psave) << nvar;
00435     }
00436 
00437       // need to rescale the hessian
00438       // get the current scale
00439       dvector x0(1,nvar);
00440       dvector current_scale(1,nvar);
00441       initial_params::xinit(x0);
00442       // cout << "starting with " << x0(x0.indexmin()) << " "
00443         //    << x0(x0.indexmax()) << endl;
00444       int mctmp=initial_params::mc_phase;
00445       initial_params::mc_phase=0;
00446       initial_params::stddev_scale(current_scale,x);
00447       initial_params::mc_phase=mctmp;
00448       int i;
00449       for (i=1;i<=nvar;i++)
00450       {
00451         for (int j=1;j<=nvar;j++)
00452         {
00453           S(i,j)*=old_scale(i)*old_scale(j);
00454         }
00455       }
00456       for (i=1;i<=nvar;i++)
00457       {
00458         for (int j=1;j<=nvar;j++)
00459         {
00460           S(i,j)/=current_scale(i)*current_scale(j);
00461         }
00462       }
00463 
00464       //cout << sort(eigenvalues(S)) << endl;
00465       dmatrix chd = choleski_decomp(S);
00466       //dmatrix tchd = trans(chd);
00467       //dmatrix chdinv=inv(chd);
00468       //int sgn;
00469       // ***************************************************************
00470       // ***************************************************************
00471       // NEW HYBRID MCMC
00472       // ***************************************************************
00473       // ***************************************************************
00474       independent_variables z(1,nvar);
00475       z=x0;
00476       gradient_structure::set_NO_DERIVATIVES();
00477       dvector g(1,nvar);
00478       //cout << "initial ll value " << get_hybrid_monte_carlo_value(nvar,z,g)
00479       //     << endl;
00480       dvector y(1,nvar);
00481       y.initialize();
00482 
00483       if (mcmc2_flag==0)
00484       {
00485         initial_params::set_inactive_random_effects();
00486       }
00487 
00488       dvector p(1,nvar);  // momentum
00489       int iseed=2197;
00490       if (iseed0) iseed=iseed0;
00491       if ( mcrestart_flag>-1)
00492       {
00493         ifstream ifs("hybrid_seed");
00494         ifs >> iseed;
00495         if (!ifs)
00496         {
00497            cerr << "error reading random number seed" << endl;
00498         }
00499       }
00500       random_number_generator rng(iseed);
00501       gradient_structure::set_YES_DERIVATIVES();
00502       ofstream ogs("sims");
00503       initial_params::xinit(x0);
00504 
00505       z=x0+chd*y;
00506       /*double llc=*/get_hybrid_monte_carlo_value(nvar,z,g);
00507       /*double llbest=*/get_hybrid_monte_carlo_value(nvar,z,g);
00508       //lbmax=llbest;
00509 
00510 
00511        int number_sims;
00512        if (nmcmc<=0)
00513        {
00514          number_sims=  100000;
00515        }
00516        else
00517        {
00518          number_sims=  nmcmc;
00519        }
00520        //double hybeps2=0.5*hybeps;
00521        double beginprior=get_hybrid_monte_carlo_value(nvar,z,g);
00522        dvector Fbegin=g*chd;
00523        // use trand(chd) ?
00524        //dvector Fbegin=tchd*g;
00525        dvector F(1,nvar);
00526        F=Fbegin;
00527        p.fill_randn(rng);
00528        if (robust_hybrid_flag)
00529        {
00530          double choice=randu(rng);
00531          if (choice<0.05)
00532          {
00533            p*=3.0;
00534          }
00535        }
00536        dmatrix xvalues(1,number_sims,1,nvar);
00537        dvector yold(1,nvar);
00538        yold=y;
00539        double pprob;
00540        if (robust_hybrid_flag==0)
00541        {
00542          pprob=0.5*norm2(p);
00543        }
00544        else
00545        {
00546          double r2=0.5*norm2(p);
00547          pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00548        }
00549        double Hbegin=beginprior+pprob;
00550        double tmpprior = 0;
00551        int ii=1;
00552        initial_params::copy_all_values(parsave,ii);
00553        // detmine whether to go forward or backward
00554 
00555        double iaccept=0.0;
00556        for (int is=1;is<=number_sims;is++)
00557        {
00558          int forflag=1;
00559          //double rnd=randu(rng);
00560          //if (rnd<0.5) forflag=0;
00561          double hstep,hstep2;
00562          //if (forflag)
00563            hstep=hybeps;
00564          //else
00565          // hstep=-hybeps;
00566          hstep2=0.5*hstep;
00567          // randomize the number of steps
00568          double rnd2=randn(rng);
00569 #ifdef OPT_LIB
00570          int hnsteps = (int)(exp(0.2 * rnd2) * hybnstep);
00571 #else
00572          double _hnsteps=exp(0.2 * rnd2) * hybnstep;
00573          assert(_hnsteps > 0 && _hnsteps <= (double)INT_MAX);
00574          int hnsteps = (int)_hnsteps;
00575 #endif
00576          for (int i=1;i<=hnsteps;i++)
00577          {
00578            if (forflag==1)
00579            {
00580              dvector phalf=p-hstep2*F;
00581              if (robust_hybrid_flag==0)
00582              {
00583                y+=hstep*phalf;
00584              }
00585              else
00586              {
00587                //pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00588                double r2=0.5*norm2(phalf);
00589                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00590                double xx=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z;
00591                dvector zz=xx*phalf;
00592                y+=hstep*zz;
00593              }
00594              z=x0+chd*y;
00595              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00596              F=g*chd;
00597              //F=tchd*g;
00598              p=phalf-hstep2*F;
00599            }
00600            else
00601            {
00602              dvector phalf=p+hstep2*F;
00603              if (robust_hybrid_flag==0)
00604              {
00605                y-=hstep*phalf;
00606              }
00607              else
00608              {
00609                double r2=0.5*norm2(phalf);
00610                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00611                dvector zz=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z*phalf;
00612                y-=hstep*phalf;
00613              }
00614              z=x0+chd*y;
00615              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00616              F=g*chd;
00617              //F=tchd*g;
00618              p=phalf+hstep2*F;
00619            }
00620          }
00621          if (robust_hybrid_flag==0)
00622          {
00623            pprob=0.5*norm2(p);
00624          }
00625          else
00626          {
00627            double r2=0.5*norm2(p);
00628            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00629          }
00630          double Ham=tmpprior+pprob;
00631          double rr=randu(rng);
00632          double pp=exp(Hbegin-Ham);
00633          p.fill_randn(rng);
00634          //p*=1.2;
00635          if (robust_hybrid_flag)
00636          {
00637            double choice=randu(rng);
00638            if (choice<0.05)
00639            {
00640              p*=3.0;
00641            }
00642          }
00643          if (robust_hybrid_flag==0)
00644          {
00645            pprob=0.5*norm2(p);
00646          }
00647          else
00648          {
00649            double r2=0.5*norm2(p);
00650            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00651          }
00652          if ((is%50)==1)
00653            //  cout << iaccept/is << " " << Hbegin-Ham << " " << Ham << endl;
00654            cout << " hybrid sim " << is <<  "  accept rate " << iaccept/is
00655                 << "  Hbegin-Ham " << Hbegin-Ham << "  Ham " << Ham << endl;
00656          if (rr<pp)
00657          {
00658            iaccept++;
00659            yold=y;
00660            beginprior=tmpprior;
00661            Hbegin=beginprior+pprob;
00662            Fbegin=F;
00663            ii=1;
00664            initial_params::copy_all_values(parsave,ii);
00665          }
00666          else
00667          {
00668            y=yold;
00669            z=x0+chd*y;
00670            Hbegin=beginprior+pprob;
00671            F=Fbegin;
00672          }
00673          (*pofs_psave) << parsave;
00674        }
00675       // cout << " saved  " << parsave(parsave.indexmin()) << " "
00676         //    << parsave(parsave.indexmax()) << endl;
00677        //double ll=get_hybrid_monte_carlo_value(nvar,parsave,g);
00678        //cout << "ll  " << ll << endl;
00679     // ***************************************************************
00680     // ***************************************************************
00681     // ***************************************************************
00682     ofstream ofs("hybrid_seed");
00683     int seed=(int) (10000*randu(rng));
00684     ofs << seed;
00685     if (pofs_psave)
00686     {
00687       delete pofs_psave;
00688       pofs_psave=NULL;
00689     }
00690   }
00691 }
00692 
00697 double function_minimizer::get_hybrid_monte_carlo_value(int nvar,
00698   const independent_variables& x,dvector& g)
00699 {
00700   //initial_params::xinit(x);
00701   double f=0.0;
00702   if (mcmc2_flag==0 && lapprox)
00703   {
00704     cerr << "error not implemented" << endl;
00705     ad_exit(1);
00706     g=(*lapprox)(x,f,this);
00707   }
00708   else
00709   {
00710     dvariable vf=0.0;
00711     dvar_vector vx=dvar_vector(x);
00712     vf=initial_params::reset(vx);
00713     *objective_function_value::pobjfun=0.0;
00714     userfunction();
00715     dvar_vector d(1,nvar);
00716     initial_params::stddev_vscale(d,vx);
00717     vf-=sum(log(d));
00718     vf+=*objective_function_value::pobjfun;
00719     f=value(vf);
00720     gradcalc(nvar,g);
00721   }
00722   return f;
00723 }