ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
hybmcmc.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: hybmcmc.cpp 1964 2014-04-30 22:09:22Z 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 double better_rand(long int&);
00023 void store_mcmc_values(const ofstream& ofs);
00024 void set_labels_for_mcmc(void);
00025 
00026 void print_hist_data(const dmatrix& hist, const dmatrix& values,
00027   const dvector& h, dvector& m, const dvector& s, const dvector& parsave,
00028   long int iseed, double size_scale);
00029 
00030 int minnz(const dvector& x);
00031 int maxnz(const dvector& xa);
00032 
00033 void read_hessian_matrix_and_scale1(int nvar, const dmatrix& _SS,double s,
00034   int mcmc2_flag);
00035 
00036 int read_hist_data(const dmatrix& hist, const dvector& h,
00037   dvector& m, const dvector& s, const dvector& parsave,long int& iseed,
00038   const double& size_scale);
00039 
00040 void make_preliminary_hist(const dvector& s, const dvector& m,int nsim,
00041   const dmatrix& values, dmatrix& hist, const dvector& h,int slots,
00042   double total_spread,int probflag=0);
00043 
00044 void add_hist_values(const dvector& s, const dvector& m, const dmatrix& hist,
00045   dvector& mcmc_values,double llc, const dvector& h,int nslots,
00046   double total_spreadd,int probflag=0);
00047 
00048 void write_empirical_covariance_matrix(int ncor, const dvector& s_mean,
00049   const dmatrix& s_covar, adstring& prog_name);
00050 
00051 void read_empirical_covariance_matrix(int nvar, const dmatrix& S,
00052   const adstring& prog_name);
00053 
00054 void read_hessian_matrix_and_scale(int nvar, const dmatrix& S,
00055   const dvector& pen_vector);
00056 
00057 dvector new_probing_bounded_multivariate_normal(int nvar, const dvector& a1,
00058   const dvector& b1, dmatrix& ch, const double& wght,double pprobe,
00059   random_number_generator& rng);
00060 
00061 void new_probing_bounded_multivariate_normal_mcmc(int nvar, const dvector& a1,
00062   const dvector& b1, dmatrix& ch, const double& wght, const dvector& _y,
00063   double pprobe, random_number_generator& rng);
00064 
00065 //void newton_raftery_bayes_estimate(double cbf,int ic, const dvector& lk,
00066 //double d);
00067 
00068 void newton_raftery_bayes_estimate_new(double cbf, int ic, const dvector& lk,
00069   double d);
00070 
00071 void ad_update_mcmc_stats_report
00072   (int feval,int iter,double fval,double gmax);
00073 
00074 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00075   double fval,double gmax,const char * cbuf);
00076 void ad_update_mcmc_report(dmatrix& m,int i,int j,int ff=0);
00077 void ad_update_mcmchist_report(dmatrix& mcmc_values,ivector& number_offsets,
00078   dvector& mean_mcmc_values,dvector& h,int ff=0);
00079 
00080 void ADSleep(int);
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          int hnsteps=hybnstep*exp(.2*rnd2);
00568 
00569          for (int i=1;i<=hnsteps;i++)
00570          {
00571            if (forflag==1)
00572            {
00573              dvector phalf=p-hstep2*F;
00574              if (robust_hybrid_flag==0)
00575              {
00576                y+=hstep*phalf;
00577              }
00578              else
00579              {
00580                //pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00581                double r2=0.5*norm2(phalf);
00582                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00583                double xx=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z;
00584                dvector zz=xx*phalf;
00585                y+=hstep*zz;
00586              }
00587              z=x0+chd*y;
00588              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00589              F=g*chd;
00590              //F=tchd*g;
00591              p=phalf-hstep2*F;
00592            }
00593            else
00594            {
00595              dvector phalf=p+hstep2*F;
00596              if (robust_hybrid_flag==0)
00597              {
00598                y-=hstep*phalf;
00599              }
00600              else
00601              {
00602                double r2=0.5*norm2(phalf);
00603                double z=0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0);
00604                dvector zz=(0.95*exp(-r2)+0.05/27.0*exp(-r2/9.0))/z*phalf;
00605                y-=hstep*phalf;
00606              }
00607              z=x0+chd*y;
00608              tmpprior=get_hybrid_monte_carlo_value(nvar,z,g);
00609              F=g*chd;
00610              //F=tchd*g;
00611              p=phalf+hstep2*F;
00612            }
00613          }
00614          if (robust_hybrid_flag==0)
00615          {
00616            pprob=0.5*norm2(p);
00617          }
00618          else
00619          {
00620            double r2=0.5*norm2(p);
00621            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00622          }
00623          double Ham=tmpprior+pprob;
00624          double rr=randu(rng);
00625          double pp=exp(Hbegin-Ham);
00626          p.fill_randn(rng);
00627          //p*=1.2;
00628          if (robust_hybrid_flag)
00629          {
00630            double choice=randu(rng);
00631            if (choice<0.05)
00632            {
00633              p*=3.0;
00634            }
00635          }
00636          if (robust_hybrid_flag==0)
00637          {
00638            pprob=0.5*norm2(p);
00639          }
00640          else
00641          {
00642            double r2=0.5*norm2(p);
00643            pprob=-log(0.95*exp(-r2)+0.05/3.0*exp(-r2/9.0));
00644          }
00645          if ((is%50)==1)
00646            //  cout << iaccept/is << " " << Hbegin-Ham << " " << Ham << endl;
00647            cout << " hybrid sim " << is <<  "  accept rate " << iaccept/is
00648                 << "  Hbegin-Ham " << Hbegin-Ham << "  Ham " << Ham << endl;
00649          if (rr<pp)
00650          {
00651            iaccept++;
00652            yold=y;
00653            beginprior=tmpprior;
00654            Hbegin=beginprior+pprob;
00655            Fbegin=F;
00656            ii=1;
00657            initial_params::copy_all_values(parsave,ii);
00658          }
00659          else
00660          {
00661            y=yold;
00662            z=x0+chd*y;
00663            Hbegin=beginprior+pprob;
00664            F=Fbegin;
00665          }
00666          (*pofs_psave) << parsave;
00667        }
00668       // cout << " saved  " << parsave(parsave.indexmin()) << " "
00669         //    << parsave(parsave.indexmax()) << endl;
00670        //double ll=get_hybrid_monte_carlo_value(nvar,parsave,g);
00671        //cout << "ll  " << ll << endl;
00672     // ***************************************************************
00673     // ***************************************************************
00674     // ***************************************************************
00675     ofstream ofs("hybrid_seed");
00676     int seed=(int) (10000*randu(rng));
00677     ofs << seed;
00678     if (pofs_psave)
00679     {
00680       delete pofs_psave;
00681       pofs_psave=NULL;
00682     }
00683   }
00684 }
00685 
00690 double function_minimizer::get_hybrid_monte_carlo_value(int nvar,
00691   const independent_variables& x,dvector& g)
00692 {
00693   //initial_params::xinit(x);
00694   double f=0.0;
00695   if (mcmc2_flag==0 && lapprox)
00696   {
00697     cerr << "error not implemented" << endl;
00698     ad_exit(1);
00699     g=(*lapprox)(x,f,this);
00700   }
00701   else
00702   {
00703     dvariable vf=0.0;
00704     dvar_vector vx=dvar_vector(x);
00705     vf=initial_params::reset(vx);
00706     *objective_function_value::pobjfun=0.0;
00707     userfunction();
00708     dvar_vector d(1,nvar);
00709     initial_params::stddev_vscale(d,vx);
00710     vf-=sum(log(d));
00711     vf+=*objective_function_value::pobjfun;
00712     f=value(vf);
00713     gradcalc(nvar,g);
00714   }
00715   return f;
00716 }