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