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