ADMB Documentation  11.5.3260
 All Classes Files Functions Variables Typedefs Friends Defines
hybmcmc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
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     int nvar_re=initial_params::nvarcalc();
00131 
00132     int nvar=initial_params::nvarcalc(); // get the number of active parameters
00133     dmatrix s_covar;
00134     dvector s_mean;
00135     //int ncsim=25000;
00136     //int nslots=800;
00137     //int nslots=3600;
00138     //int initial_nsim=4800;
00139     //int ntmp=0;
00140     //int ncor=0;
00141     //double bfsum=0;
00142     //int ibfcount=0;
00143     //double lbmax;
00144 
00145     s_covar.allocate(1,nvar,1,nvar);
00146     s_mean.allocate(1,nvar);
00147     s_mean.initialize();
00148     s_covar.initialize();
00149 
00150     int ndvar=stddev_params::num_stddev_calc();
00151     //int numdvar=stddev_params::num_stddev_number_calc();
00152 
00153     if (mcmc2_flag==0)
00154     {
00155       initial_params::set_inactive_random_effects();
00156       nvar=initial_params::nvarcalc(); // get the number of active parameters
00157     }
00158 
00159     independent_variables parsave(1,nvar_re);
00160     initial_params::restore_start_phase();
00161 
00162     dvector x(1,nvar);
00163     //dvector scale(1,nvar);
00164     dmatrix values;
00165     //int have_hist_flag=0;
00166     initial_params::xinit(x);
00167     dvector pen_vector(1,nvar);
00168     {
00169       initial_params::reset(dvar_vector(x),pen_vector);
00170       //cout << pen_vector << endl << endl;
00171     }
00172 
00173     initial_params::mc_phase=0;
00174     //initial_params::stddev_scale(scale,x);
00175     initial_params::mc_phase=1;
00176     dvector bmn(1,nvar);
00177     dvector mean_mcmc_values(1,ndvar);
00178     dvector s(1,ndvar);
00179     dvector h(1,ndvar);
00180     //dvector h;
00181     dvector square_mcmc_values(1,ndvar);
00182     square_mcmc_values.initialize();
00183     mean_mcmc_values.initialize();
00184     bmn.initialize();
00185     int use_empirical_flag=0;
00186     int diag_option=0;
00187     //int topt=0;
00188     int hybnstep=10;
00189     double hybeps=0.1;
00190     double _hybeps=-1.0;
00191     int old_Hybrid_bounded_flag=-1;
00192 
00193     int on,nopt = 0;
00194     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hyeps",nopt))>-1)
00195     {
00196       if (!nopt)
00197       {
00198         cerr << "Usage -hyeps option needs number  -- ignored" << endl;
00199       }
00200       else
00201       {
00202         istringstream ist(ad_comm::argv[on+1]);
00203         ist >> _hybeps;
00204 
00205         if (_hybeps<=0)
00206         {
00207           cerr << "Usage -hyeps option needs positive number  -- ignored\n";
00208           _hybeps=-1.0;
00209         }
00210       }
00211     }
00212     if (_hybeps>0.0)
00213       hybeps=_hybeps;
00214 
00215     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-hynstep",nopt))>-1)
00216     {
00217       if (nopt)
00218       {
00219         int iii=atoi(ad_comm::argv[on+1]);
00220         if (iii < 1 )
00221         {
00222           cerr << " -hynstep argument must be integer between > 0 --"
00223                   " using default of 10" << endl;
00224         }
00225         else
00226         {
00227           hybnstep=iii;
00228         }
00229       }
00230     }
00231 
00232     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mcdiag"))>-1)
00233     {
00234       diag_option=1;
00235       cout << " Setting covariance matrix to diagonal with entries " << dscale
00236            << endl;
00237     }
00238     int mcrestart_flag=option_match(ad_comm::argc,ad_comm::argv,"-mcr");
00239     dmatrix S(1,nvar,1,nvar);
00240     dvector old_scale(1,nvar);
00241     int old_nvar;
00242     if (!diag_option)
00243     {
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 #ifndef OPT_LIB
00368         assert(parsave.size() >= 0);
00369 #endif
00370         std::streamoff sz = (unsigned int)parsave.size() * sizeof(double);
00371         // backup from end of file
00372         uis.seekg(-sz, ios::end);
00373         uis >> parsave;
00374        cout << " restored "  << parsave(parsave.indexmin()) << " "
00375             << parsave(parsave.indexmax()) << endl;
00376         if (!uis)
00377         {
00378           cerr << "error resotring last mcmc poistion from file "
00379             << ad_comm::adprogram_name + adstring(".psv") << endl;
00380           ad_exit(1);
00381         }
00382         int ii=1;
00383         dvector x0(1,nvar);
00384         initial_params::restore_all_values(parsave,ii);
00385         //x0.initialize();
00386         //dvector pen_vector(1,nvar);
00387         //initial_params::reset(dvar_vector(parsave),pen_vector);
00388         //initial_params::xinit(x0);
00389         //cout << " x0 " << x0(x0.indexmin()) << endl;
00390       }
00391     }
00392 
00393 
00394     if ( mcrestart_flag>-1)
00395     {
00396       pofs_psave= new uostream( (char*)(ad_comm::adprogram_name
00397           + adstring(".psv")),ios::app);
00398      /*
00399       pofs_psave->seekp(0,std::ios::end);
00400       *pofs_psave << 123;
00401       delete pofs_psave;
00402       uistream uis((char*)(ad_comm::adprogram_name + adstring(".psv")));
00403       if (!uis)
00404       {
00405         cerr << "Error trying to open file" <<
00406           ad_comm::adprogram_name + adstring(".psv") <<
00407           " for mcrestart" <<   endl;
00408         ad_exit(1);
00409       } else {
00410         int nv1;
00411         uis >> nv1;
00412         if (nv1 !=nvar) {
00413           cerr << "wrong number of independent variables in" <<
00414             ad_comm::adprogram_name + adstring(".psv") << endl
00415            << " I found " << nv1 << " instead of " << nvar  << endl;
00416           ad_exit(1);
00417         }
00418       }
00419      */
00420     }
00421     else
00422     {
00423       pofs_psave=
00424         new uostream((char*)(ad_comm::adprogram_name + adstring(".psv")));
00425     }
00426 
00427     if (!pofs_psave|| !(*pofs_psave))
00428     {
00429       cerr << "Error trying to open file" <<
00430         ad_comm::adprogram_name + adstring(".psv") << endl;
00431       ad_exit(1);
00432     }
00433     if (mcrestart_flag == -1 )
00434     {
00435         (*pofs_psave) << nvar;
00436     }
00437 
00438       // need to rescale the hessian
00439       // get the current scale
00440       dvector x0(1,nvar);
00441       dvector current_scale(1,nvar);
00442       initial_params::xinit(x0);
00443       // cout << "starting with " << x0(x0.indexmin()) << " "
00444         //    << x0(x0.indexmax()) << endl;
00445       int mctmp=initial_params::mc_phase;
00446       initial_params::mc_phase=0;
00447       initial_params::stddev_scale(current_scale,x);
00448       initial_params::mc_phase=mctmp;
00449       for (int 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 (int 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 }