ADMB Documentation  11.1.2500
 All Classes Files Functions Variables Typedefs Friends Defines
hybmcmc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: hybmcmc.cpp 2463 2014-10-03 00:20:26Z 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 on,nopt = 0;
00244       int rescale_bounded_flag=0;
00245       double rescale_bounded_power=0.5;
00246       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcrb",nopt))>-1)
00247       {
00248         if (nopt)
00249         {
00250           int iii=atoi(ad_comm::argv[on+1]);
00251           if (iii < 1 || iii > 9)
00252           {
00253             cerr << " -mcrb argument must be integer between 1 and 9 --"
00254                     " using default of 5" << endl;
00255             rescale_bounded_power=0.5;
00256           }
00257           else
00258             rescale_bounded_power=iii/10.0;
00259         }
00260         else
00261         {
00262           rescale_bounded_power=0.5;
00263         }
00264         rescale_bounded_flag=1;
00265       }
00266       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcec"))>-1)
00267       {
00268         use_empirical_flag=1;
00269       }
00270       if (use_empirical_flag)
00271       {
00272         read_empirical_covariance_matrix(nvar,S,ad_comm::adprogram_name);
00273       }
00274       else if (!rescale_bounded_flag)
00275       {
00276         if (mcmc2_flag==0)
00277         {
00278           read_covariance_matrix(S,nvar,old_Hybrid_bounded_flag,old_scale);
00279         }
00280         else
00281         {
00282           int tmp_nvar = 0;
00283           adstring tmpstring = ad_comm::adprogram_name + ".bgs";
00284           uistream uis((char*)(tmpstring));
00285           if (!uis)
00286           {
00287             cerr << "error opening file " << tmpstring << endl;
00288             ad_exit(1);
00289           }
00290           uis >> tmp_nvar;
00291           if (!uis)
00292           {
00293             cerr << "error reading from file " << tmpstring << endl;
00294             ad_exit(1);
00295           }
00296           if (tmp_nvar != nvar)
00297           {
00298             cerr << "size error reading from " <<  tmpstring << endl;
00299             ad_exit(1);
00300           }
00301           uis >> S;
00302           if (!uis)
00303           {
00304             cerr << "error reading from file " << tmpstring << endl;
00305             ad_exit(1);
00306           }
00307           dvector tmp=read_old_scale(old_nvar);
00308           old_scale=1.0;
00309           old_scale(1,old_nvar)=tmp;
00310         }
00311       }
00312       else
00313       {
00314         read_hessian_matrix_and_scale1(nvar,S,rescale_bounded_power,
00315           mcmc2_flag);
00316         //read_hessian_matrix_and_scale(nvar,S,pen_vector);
00317       }
00318 
00319       // maybe we dont want to do this?
00320       /*
00321       {  // scale covariance matrix for model space
00322         dmatrix tmp(1,nvar,1,nvar);
00323         for (int i=1;i<=nvar;i++)
00324         {
00325           tmp(i,i)=S(i,i)*(scale(i)*scale(i));
00326           for (int j=1;j<i;j++)
00327           {
00328             tmp(i,j)=S(i,j)*(scale(i)*scale(j));
00329             tmp(j,i)=tmp(i,j);
00330           }
00331         }
00332         S=tmp;
00333       }
00334       */
00335     }
00336     else
00337     {
00338       S.initialize();
00339       for (int i=1;i<=nvar;i++)
00340       {
00341         S(i,i)=dscale;
00342       }
00343     }
00344 
00345     // for hybrid mcmc option always save output
00346     //if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcsave"))>-1)
00347     if ( mcrestart_flag>-1)
00348     {
00349       // check that nvar is correct
00350       uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00351       if (!uis)
00352       {
00353         cerr << "Error trying to open file" <<
00354           ad_comm::adprogram_name + adstring(".psv") <<
00355           " for mcrestart" <<   endl;
00356         ad_exit(1);
00357       } else {
00358         int nv1 = 0;
00359         uis >> nv1;
00360         if (nv1 !=nvar) {
00361           cerr << "wrong number of independent variables in" <<
00362             ad_comm::adprogram_name + adstring(".psv") << endl
00363            << " I found " << nv1 << " instead of " << nvar  << endl;
00364           ad_exit(1);
00365         }
00366         // get last x vector from file
00367         int sz=parsave.size()*sizeof(double);
00368         // backup from end of file
00369         uis.seekg(-sz,ios::end);
00370         uis >> parsave;
00371        cout << " restored "  << parsave(parsave.indexmin()) << " "
00372             << parsave(parsave.indexmax()) << endl;
00373         if (!uis)
00374         {
00375           cerr << "error resotring last mcmc poistion from file "
00376             << ad_comm::adprogram_name + adstring(".psv") << endl;
00377           ad_exit(1);
00378         }
00379         int ii=1;
00380         dvector x0(1,nvar);
00381         initial_params::restore_all_values(parsave,ii);
00382         //x0.initialize();
00383         //dvector pen_vector(1,nvar);
00384         //initial_params::reset(dvar_vector(parsave),pen_vector);
00385         //initial_params::xinit(x0);
00386         //cout << " x0 " << x0(x0.indexmin()) << endl;
00387       }
00388     }
00389 
00390 
00391     if ( mcrestart_flag>-1)
00392     {
00393       pofs_psave= new uostream( (char*)(ad_comm::adprogram_name
00394           + adstring(".psv")),ios::app);
00395      /*
00396       pofs_psave->seekp(0,std::ios::end);
00397       *pofs_psave << 123;
00398       delete pofs_psave;
00399       uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00400       if (!uis)
00401       {
00402         cerr << "Error trying to open file" <<
00403           ad_comm::adprogram_name + adstring(".psv") <<
00404           " for mcrestart" <<   endl;
00405         ad_exit(1);
00406       } else {
00407         int nv1;
00408         uis >> nv1;
00409         if (nv1 !=nvar) {
00410           cerr << "wrong number of independent variables in" <<
00411             ad_comm::adprogram_name + adstring(".psv") << endl
00412            << " I found " << nv1 << " instead of " << nvar  << endl;
00413           ad_exit(1);
00414         }
00415       }
00416      */
00417     }
00418     else
00419     {
00420       pofs_psave=
00421         new uostream((char*)(ad_comm::adprogram_name + adstring(".psv")));
00422     }
00423 
00424     if (!pofs_psave|| !(*pofs_psave))
00425     {
00426       cerr << "Error trying to open file" <<
00427         ad_comm::adprogram_name + adstring(".psv") << endl;
00428       ad_exit(1);
00429     }
00430     if (mcrestart_flag == -1 )
00431     {
00432         (*pofs_psave) << nvar;
00433     }
00434 
00435       // need to rescale the hessian
00436       // get the current scale
00437       dvector x0(1,nvar);
00438       dvector current_scale(1,nvar);
00439       initial_params::xinit(x0);
00440       // cout << "starting with " << x0(x0.indexmin()) << " "
00441         //    << x0(x0.indexmax()) << endl;
00442       int mctmp=initial_params::mc_phase;
00443       initial_params::mc_phase=0;
00444       initial_params::stddev_scale(current_scale,x);
00445       initial_params::mc_phase=mctmp;
00446       int i;
00447       for (i=1;i<=nvar;i++)
00448       {
00449         for (int j=1;j<=nvar;j++)
00450         {
00451           S(i,j)*=old_scale(i)*old_scale(j);
00452         }
00453       }
00454       for (i=1;i<=nvar;i++)
00455       {
00456         for (int j=1;j<=nvar;j++)
00457         {
00458           S(i,j)/=current_scale(i)*current_scale(j);
00459         }
00460       }
00461 
00462       //cout << sort(eigenvalues(S)) << endl;
00463       dmatrix chd = choleski_decomp(S);
00464       //dmatrix tchd = trans(chd);
00465       //dmatrix chdinv=inv(chd);
00466       //int sgn;
00467       // ***************************************************************
00468       // ***************************************************************
00469       // NEW HYBRID MCMC
00470       // ***************************************************************
00471       // ***************************************************************
00472       independent_variables z(1,nvar);
00473       z=x0;
00474       gradient_structure::set_NO_DERIVATIVES();
00475       dvector g(1,nvar);
00476       //cout << "initial ll value " << get_hybrid_monte_carlo_value(nvar,z,g)
00477       //     << endl;
00478       dvector y(1,nvar);
00479       y.initialize();
00480 
00481       if (mcmc2_flag==0)
00482       {
00483         initial_params::set_inactive_random_effects();
00484       }
00485 
00486       dvector p(1,nvar);  // momentum
00487       int iseed=2197;
00488       if (iseed0) iseed=iseed0;
00489       if ( mcrestart_flag>-1)
00490       {
00491         ifstream ifs("hybrid_seed");
00492         ifs >> iseed;
00493         if (!ifs)
00494         {
00495            cerr << "error reading random number seed" << endl;
00496         }
00497       }
00498       random_number_generator rng(iseed);
00499       gradient_structure::set_YES_DERIVATIVES();
00500       ofstream ogs("sims");
00501       initial_params::xinit(x0);
00502 
00503       z=x0+chd*y;
00504       /*double llc=*/get_hybrid_monte_carlo_value(nvar,z,g);
00505       /*double llbest=*/get_hybrid_monte_carlo_value(nvar,z,g);
00506       //lbmax=llbest;
00507 
00508 
00509        int number_sims;
00510        if (nmcmc<=0)
00511        {
00512          number_sims=  100000;
00513        }
00514        else
00515        {
00516          number_sims=  nmcmc;
00517        }
00518        //double hybeps2=0.5*hybeps;
00519        double beginprior=get_hybrid_monte_carlo_value(nvar,z,g);
00520        dvector Fbegin=g*chd;
00521        // use trand(chd) ?
00522        //dvector Fbegin=tchd*g;
00523        dvector F(1,nvar);
00524        F=Fbegin;
00525        p.fill_randn(rng);
00526        if (robust_hybrid_flag)
00527        {
00528          double choice=randu(rng);
00529          if (choice<0.05)
00530          {
00531            p*=3.0;
00532          }
00533        }
00534        dmatrix xvalues(1,number_sims,1,nvar);
00535        dvector yold(1,nvar);
00536        yold=y;
00537        double pprob;
00538        if (robust_hybrid_flag==0)
00539        {
00540          pprob=0.5*norm2(p);
00541        }
00542        else
00543        {
00544          double r2=0.5*norm2(p);
00545          pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00546        }
00547        double Hbegin=beginprior+pprob;
00548        double tmpprior = 0;
00549        int ii=1;
00550        initial_params::copy_all_values(parsave,ii);
00551        // detmine whether to go forward or backward
00552 
00553        double iaccept=0.0;
00554        for (int is=1;is<=number_sims;is++)
00555        {
00556          int forflag=1;
00557          //double rnd=randu(rng);
00558          //if (rnd<0.5) forflag=0;
00559          double hstep,hstep2;
00560          //if (forflag)
00561            hstep=hybeps;
00562          //else
00563          // hstep=-hybeps;
00564          hstep2=0.5*hstep;
00565          // randomize the number of steps
00566          double rnd2=randn(rng);
00567 #ifdef OPT_LIB
00568          int hnsteps = (int)(exp(0.2 * rnd2) * hybnstep);
00569 #else
00570          double _hnsteps=exp(0.2 * rnd2) * hybnstep;
00571          assert(_hnsteps > 0 && _hnsteps <= (double)INT_MAX);
00572          int hnsteps = (int)_hnsteps;
00573 #endif
00574          for (int i=1;i<=hnsteps;i++)
00575          {
00576            if (forflag==1)
00577            {
00578              dvector phalf=p-hstep2*F;
00579              if (robust_hybrid_flag==0)
00580              {
00581                y+=hstep*phalf;
00582              }
00583              else
00584              {
00585                //pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00586                double r2=0.5*norm2(phalf);
00587                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00588                double xx=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z;
00589                dvector zz=xx*phalf;
00590                y+=hstep*zz;
00591              }
00592              z=x0+chd*y;
00593              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00594              F=g*chd;
00595              //F=tchd*g;
00596              p=phalf-hstep2*F;
00597            }
00598            else
00599            {
00600              dvector phalf=p+hstep2*F;
00601              if (robust_hybrid_flag==0)
00602              {
00603                y-=hstep*phalf;
00604              }
00605              else
00606              {
00607                double r2=0.5*norm2(phalf);
00608                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00609                dvector zz=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z*phalf;
00610                y-=hstep*phalf;
00611              }
00612              z=x0+chd*y;
00613              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00614              F=g*chd;
00615              //F=tchd*g;
00616              p=phalf+hstep2*F;
00617            }
00618          }
00619          if (robust_hybrid_flag==0)
00620          {
00621            pprob=0.5*norm2(p);
00622          }
00623          else
00624          {
00625            double r2=0.5*norm2(p);
00626            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00627          }
00628          double Ham=tmpprior+pprob;
00629          double rr=randu(rng);
00630          double pp=exp(Hbegin-Ham);
00631          p.fill_randn(rng);
00632          //p*=1.2;
00633          if (robust_hybrid_flag)
00634          {
00635            double choice=randu(rng);
00636            if (choice<0.05)
00637            {
00638              p*=3.0;
00639            }
00640          }
00641          if (robust_hybrid_flag==0)
00642          {
00643            pprob=0.5*norm2(p);
00644          }
00645          else
00646          {
00647            double r2=0.5*norm2(p);
00648            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00649          }
00650          if ((is%50)==1)
00651            //  cout << iaccept/is << " " << Hbegin-Ham << " " << Ham << endl;
00652            cout << " hybrid sim " << is <<  "  accept rate " << iaccept/is
00653                 << "  Hbegin-Ham " << Hbegin-Ham << "  Ham " << Ham << endl;
00654          if (rr<pp)
00655          {
00656            iaccept++;
00657            yold=y;
00658            beginprior=tmpprior;
00659            Hbegin=beginprior+pprob;
00660            Fbegin=F;
00661            ii=1;
00662            initial_params::copy_all_values(parsave,ii);
00663          }
00664          else
00665          {
00666            y=yold;
00667            z=x0+chd*y;
00668            Hbegin=beginprior+pprob;
00669            F=Fbegin;
00670          }
00671          (*pofs_psave) << parsave;
00672        }
00673       // cout << " saved  " << parsave(parsave.indexmin()) << " "
00674         //    << parsave(parsave.indexmax()) << endl;
00675        //double ll=get_hybrid_monte_carlo_value(nvar,parsave,g);
00676        //cout << "ll  " << ll << endl;
00677     // ***************************************************************
00678     // ***************************************************************
00679     // ***************************************************************
00680     ofstream ofs("hybrid_seed");
00681     int seed=(int) (10000*randu(rng));
00682     ofs << seed;
00683     if (pofs_psave)
00684     {
00685       delete pofs_psave;
00686       pofs_psave=NULL;
00687     }
00688   }
00689 }
00690 
00695 double function_minimizer::get_hybrid_monte_carlo_value(int nvar,
00696   const independent_variables& x,dvector& g)
00697 {
00698   //initial_params::xinit(x);
00699   double f=0.0;
00700   if (mcmc2_flag==0 && lapprox)
00701   {
00702     cerr << "error not implemented" << endl;
00703     ad_exit(1);
00704     g=(*lapprox)(x,f,this);
00705   }
00706   else
00707   {
00708     dvariable vf=0.0;
00709     dvar_vector vx=dvar_vector(x);
00710     vf=initial_params::reset(vx);
00711     *objective_function_value::pobjfun=0.0;
00712     userfunction();
00713     dvar_vector d(1,nvar);
00714     initial_params::stddev_vscale(d,vx);
00715     vf-=sum(log(d));
00716     vf+=*objective_function_value::pobjfun;
00717     f=value(vf);
00718     gradcalc(nvar,g);
00719   }
00720   return f;
00721 }