ADMB Documentation  11.1.1897
 All Classes Files Functions Variables Typedefs Friends Defines
xmodelm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xmodelm3.cpp 1819 2014-04-01 01:31:22Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <sstream>
00008 using std::istringstream;
00009 
00010 #include <admodel.h>
00011 
00012 #if defined(USE_LAPLACE)
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015 #endif
00016 
00017 void check_java_flags(int& start_flag,int& quit_flag,int& der_flag,
00018   int& next_flag);
00019 
00020 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00021   double fval, double gmax,const char * cbuf);
00022 void vm_initialize(void);
00023 
00024 void set_initial_simplex(const tdmatrix& p, const dvector& y,int nvar,
00025   const dvector& x, double delta);
00026 
00027 int get_option_number(const char * option_name,const char * error_message,
00028   int& option_value);
00029 
00030 int get_option_number(const char * option_name,const char * error_message,
00031 #ifdef __BORLANDC__
00032   long int& option_value);
00033 #else
00034   long long int& option_value);
00035 #endif
00036 
00037 extern int traceflag;
00038 
00039 void tracing_message(int traceflag,const char *s);
00040 
00041 #if defined(USE_LAPLACE)
00042   int function_minimizer::inner_opt(void)
00043   {
00044     return inner_opt_flag;
00045   }
00046 
00047   int function_minimizer::inner_opt_flag=0;
00048 #endif
00049 
00050 
00051   int function_minimizer::bad_step_flag=0;
00052 
00053   void function_minimizer::minimize(void)
00054   {
00055     int nopt=0;
00056     int on=0;
00057 #   if defined(USE_LAPLACE)
00058     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-shess"))>-1)
00059     {
00060       laplace_approximation_calculator::sparse_hessian_flag=1;
00061     }
00062     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-pis"))>-1)
00063     {
00064      laplace_approximation_calculator::print_importance_sampling_weights_flag=1;
00065     }
00066     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-sp"))>-1)
00067     {
00068       laplace_approximation_calculator::saddlepointflag=1;
00069     }
00070 #    if defined(__MINI_MAX__)
00071         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mm"))>-1)
00072         {
00073           laplace_approximation_calculator::saddlepointflag=2;
00074         }
00075 #    else
00076         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mm"))>-1)
00077         {
00078            cerr << "option -mm MINI_MAX not defined " << endl;
00079            ad_exit(1);
00080         }
00081 #    endif
00082 
00083 #   endif
00084     //initial_params::read(); // read in the values for the initial parameters
00085     if (initial_params::restart_phase)
00086     {
00087       initial_params::current_phase = initial_params::restart_phase;
00088       initial_params::restart_phase=0;
00089     }
00090     int allphases=initial_params::max_number_phases;
00091     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxph",nopt))>-1)
00092     {
00093       if (!nopt)
00094       {
00095         cerr << "Usage -maxph option needs integer  -- ignored" << endl;
00096       }
00097       else
00098       {
00099         int jj=atoi(ad_comm::argv[on+1]);
00100         if (jj<=0)
00101         {
00102           cerr << "Usage -maxph option needs positive integer  -- ignored.\n";
00103         }
00104         else
00105         {
00106           if (jj>allphases)
00107           {
00108             allphases=jj;
00109           }
00110         }
00111       }
00112       if (allphases>initial_params::max_number_phases)
00113       {
00114         initial_params::max_number_phases=allphases;
00115       }
00116     }
00117     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ndv",nopt))>-1)
00118     {
00119       if (!nopt)
00120       {
00121         cerr << "Usage -ndv option needs integer  -- ignored" << endl;
00122       }
00123       else
00124       {
00125         int jj=atoi(ad_comm::argv[on+1]);
00126         if (jj<=0)
00127         {
00128           cerr << "Usage -ndv option needs positive integer  -- ignored.\n";
00129         }
00130         else
00131         {
00132           gradient_structure::NUM_DEPENDENT_VARIABLES=jj;
00133         }
00134       }
00135     }
00136 
00137     // set the maximum number of function evaluations by command line
00138     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxfn",nopt))>-1)
00139     {
00140       if (!nopt)
00141       {
00142         cerr << "Usage -maxph option needs integer  -- ignored" << endl;
00143       }
00144       else
00145       {
00146         int _maxfn=atoi(ad_comm::argv[on+1]);
00147         if (_maxfn<0)
00148         {
00149           cerr << "Usage -maxfn option needs positive integer  -- ignored.\n";
00150         }
00151         else
00152         {
00153           maxfn=_maxfn;
00154         }
00155       }
00156     }
00157     double _crit=0;
00158     // set the maximum number of function evaluations by command line
00159     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-crit",nopt))>-1)
00160     {
00161       if (!nopt)
00162       {
00163         cerr << "Usage -crit option needs number  -- ignored" << endl;
00164       }
00165       else
00166       {
00167         istringstream ist(ad_comm::argv[on+1]);
00168         ist >> _crit;
00169 
00170         if (_crit<=0)
00171         {
00172           cerr << "Usage -crit option needs positive number  -- ignored.\n";
00173           _crit=0.0;
00174         }
00175       }
00176     }
00177     int bandwidth=0;
00178     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-bw",nopt))>-1)
00179     {
00180       if (!nopt)
00181       {
00182         cerr << "Usage -bw option needs number  -- ignored" << endl;
00183       }
00184       else
00185       {
00186         istringstream ist(ad_comm::argv[on+1]);
00187         ist >> bandwidth;
00188 
00189         if (bandwidth<=0)
00190         {
00191           cerr << "Usage -bw option needs positive number  -- ignored" << endl;
00192           bandwidth=0.0;
00193         }
00194         else
00195         {
00196           ad_comm::bandwidth=bandwidth;
00197         }
00198       }
00199     }
00200     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-phase"))>-1)
00201     {
00202       int jj=atoi(ad_comm::argv[on+1]);
00203       if (jj <=0)
00204       {
00205         cerr << " Invalid option following command line option -phase -- "
00206           << endl << " phase set equal to 1" << endl;
00207       }
00208       if (jj>allphases)
00209       {
00210         jj=allphases;
00211       }
00212       if (jj<=0)
00213       {
00214         jj=1;
00215       }
00216       initial_params::current_phase = jj;
00217       cout << "Set current phase to " << jj << endl;
00218     }
00219     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-lapqd"))>-1)
00220     {
00221       ADqd_flag=1;
00222     }
00223 
00224     tracing_message(traceflag,"A2");
00225     while (initial_params::current_phase <= allphases)
00226     {
00227       between_phases_calculations();
00228 
00229 #if defined(USE_LAPLACE)
00230       if (random_effects_flag)
00231         initial_params::set_inactive_random_effects();
00232 #endif
00233 
00234       int nvar=initial_params::nvarcalc(); // get the number of active
00235              // parameters
00236       if (!nvar)
00237       {
00238         cerr << "Error -- no active parameters. There must be at least 1"
00239              << endl;
00240         exit(1);
00241       }
00242       dvector g(1,nvar);
00243       independent_variables x(1,nvar);
00244       tracing_message(traceflag,"B2");
00245       initial_params::xinit(x);    // get the initial values into the
00246       dvector zz(1,x.indexmax());
00247       dvector xsave(1,x.indexmax());
00248       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-uhess"))>-1)
00249       {
00250         int ierr=0;
00251         ifstream ifs("vector");
00252         if (!ifs)
00253         {
00254           cerr << "couldn't open file vector" << endl;
00255           ierr=1;
00256         }
00257         if(ierr==0)
00258         {
00259           ifs >> zz;
00260           if (!ifs)
00261           {
00262             cerr << "couldn't read vector" << endl;
00263             ierr=1;
00264           }
00265         }
00266         if (ierr==0)
00267         {
00268           do
00269           {
00270             double delta=0;
00271             cout << "enter delta" << endl;
00272             cin >> delta;
00273             xsave=x;
00274             x+=delta*zz;
00275             initial_params::reset(x);    // get the initial values into the
00276             userfunction();
00277             x=xsave;
00278           }
00279           while(1);
00280         }
00281       }
00282 
00283       double f=0.0;
00284 
00285       int lmnflag = -1;
00286       int simpflag = -1;
00287       int lmnsteps=10;
00288       if ( (simpflag=option_match(ad_comm::argc,ad_comm::argv,"-simplex"))>-1)
00289       {
00290         gradient_structure::set_NO_DERIVATIVES();
00291         double delta=1.e-4;
00292         double ftol=1.e-16;
00293         dmatrix p(1,nvar+1,1,nvar);
00294         dvector y(1,nvar+1);
00295         set_initial_simplex(p,y,nvar,x,delta);
00296         adamoeba(p,y,nvar,ftol,maxfn);
00297         double ymin=min(y);
00298         for (int i=1;i<=nvar+1;i++)
00299         if (ymin==y(i))
00300         {
00301           x=p(i);
00302           break;
00303         }
00304         cerr << "The -simplex option is deprecated. The user should port "
00305              << "to the -neldmead option." << endl;
00306       }
00307       if ( (simpflag=option_match(ad_comm::argc,ad_comm::argv,"-neldmead"))>-1)
00308       {
00309         gradient_structure::set_NO_DERIVATIVES();
00310         double delta=1.e-4;
00311         double ftol=1.e-16;
00312         dvector mincords = x;
00313         double ynewlo;
00314         double* pynewlo = &ynewlo;
00315         int icount, numres, ifault;
00316         int* picount = &icount;
00317         int* pnumres = &numres;
00318         int* pifault = &ifault;
00319         neldmead(nvar,mincords,mincords,pynewlo,ftol,delta,picount,pnumres,
00320           pifault);
00321         x = mincords;
00322       }
00323       if ( (lmnflag=option_match(ad_comm::argc,ad_comm::argv,"-lmn",nopt))>-1)
00324       {
00325 #if   defined(USE_LAPLACE)
00326         if (random_effects_flag)
00327         {
00328           cerr << "At present you can not use the -lmn option for the outer"
00329                << endl << " optimiation in a random-effects model" << endl;
00330           ad_exit(1);
00331         }
00332 #endif
00333         if (!nopt)
00334         {
00335           cerr << "Usage -lmn option needs integer  -- set to default 10.\n";
00336         }
00337         else
00338         {
00339           int jj=atoi(ad_comm::argv[lmnflag+1]);
00340           if (jj<=0)
00341           {
00342             cerr << "Usage -lmn option needs positive integer"
00343                     " -- set to default 10.\n";
00344           }
00345           else
00346           {
00347             lmnsteps=jj;
00348           }
00349         }
00350       }
00351       if (lmnflag==-1)
00352       {
00353         // *********************************************************
00354         // block for quasi newton minimization
00355 #if   defined(USE_LAPLACE)
00356         if (negative_eigenvalue_flag)
00357         {
00358           trust_region_update(nvar,_crit,x,g,f);
00359         }
00360 #endif   //defined(USE_LAPLACE)
00361         if (!ad_comm::pvm_manager)
00362         {
00363           do
00364           {
00365             if (spminflag)
00366             {
00367               repeatminflag=1;
00368               spminflag=0;
00369             }
00370             else
00371             {
00372               repeatminflag=0;
00373             }
00374             if (maxfn>0)
00375             {
00376               int lmnflag=-1;
00377               int nsteps=5;
00378               if ( (lmnflag=option_match(ad_comm::argc,ad_comm::argv,
00379                  "-lmn2",nopt))>-1)
00380               {
00381                 if (!nopt)
00382                 {
00383                   cerr << "Usage -lmn option needs integer"
00384                      "  -- set to default 5" << endl;
00385                 }
00386                 else
00387                 {
00388                   int jj=atoi(ad_comm::argv[lmnflag+1]);
00389                   if (jj<=0)
00390                   {
00391                     cerr << "Usage -lmn option needs positive integer "
00392                      " -- set to defalt 5" << endl;
00393                   }
00394                   else
00395                   {
00396                     nsteps=jj;
00397                   }
00398                 }
00399               }
00400               if (lmnflag<0)
00401               {
00402                 quasi_newton_block(nvar,_crit,x,g,f);
00403               }
00404               else
00405               {
00406                 limited_memory_quasi_newton_block(nvar,_crit,x,g,f,nsteps);
00407               }
00408             }
00409           }
00410           while(repeatminflag);
00411         }
00412         else
00413         {
00414 #if defined(USE_ADPVM)
00415 #if   defined(USE_LAPLACE)
00416       if (random_effects_flag)
00417       {
00418         if (maxfn>0)
00419         {
00420           switch (ad_comm::pvm_manager->mode)
00421           {
00422           case 1: // master
00423             quasi_newton_block_pvm_master_random_effects(nvar,_crit,x,g,f);
00424             break;
00425           case 2: // slave
00426           // these don't exist yet
00427             function_evaluation_block_pvm_slave_random_effects
00428               (nvar,_crit,x,g,f);
00429             break;
00430           default:
00431             cerr << "error illega value for pvm_manager->mode" << endl;
00432             exit(1);
00433           }
00434         }
00435       }
00436       else
00437 #  endif // #if   defined(USE_LAPLACE)
00438       {
00439         if (maxfn>0)
00440         {
00441           switch (ad_comm::pvm_manager->mode)
00442           {
00443           case 1: // master
00444             quasi_newton_block_pvm_master(nvar,_crit,x,g,f);
00445             break;
00446           case 2: // slave
00447             function_evaluation_block_pvm_slave();
00448             break;
00449           default:
00450             cerr << "error illega value for pvm_manager->mode" << endl;
00451             exit(1);
00452           }
00453         }
00454       }
00455 #else
00456           {
00457             cerr << "PVM not included with this distribution" << endl;
00458             ad_exit(1);
00459           }
00460 #endif
00461         }
00462       } // end block for quasi newton minimization
00463       else
00464       {  // block for limited memory quasi newton minimization
00465         if (maxfn>0)
00466         {
00467           function_minimizer::limited_memory_quasi_newton(x,lmnsteps);
00468         }
00469       }
00470       // end block for limited memory quasi newton minimization
00471       // *********************************************************
00472       tracing_message(traceflag,"M2");
00473 
00474       gradient_structure::set_NO_DERIVATIVES();
00475       initial_params::reset(dvar_vector(x));
00476       *objective_function_value::pobjfun=0.0;
00477 #if defined(USE_LAPLACE)
00478       if (!random_effects_flag || !lapprox)
00479       {
00480 #endif
00481 
00482 #if defined(USE_ADPVM)
00483         if (ad_comm::pvm_manager)
00484         {
00485           switch (ad_comm::pvm_manager->mode)
00486           {
00487           case 1:
00488             pvm_master_function_evaluation_no_derivatives(f,x,nvar);
00489             *objective_function_value::pobjfun=f;
00490             break;
00491           case 2:
00492             pvm_slave_function_evaluation_no_derivatives();
00493             break;
00494           default:
00495             cerr << "Illegal value for ad_comm::pvm_manager->mode" << endl;
00496             ad_exit(1);
00497           }
00498         }
00499         else
00500         {
00501 #endif  //#if defined(USE_ADPVM)
00502           userfunction();
00503 #if defined(USE_ADPVM)
00504         }
00505 #endif  //#if defined(USE_ADPVM)
00506 
00507 #if defined(USE_LAPLACE)
00508       }
00509       else
00510       {
00511         (*lapprox)(x,f,this);
00512         *objective_function_value::pobjfun=f;
00513         initial_params::set_inactive_only_random_effects();
00514         print_is_diagnostics(lapprox);
00515       }
00516 #endif
00517       initial_params::save();
00518       report(g);
00519       // in case the user changes some initial_params in the report section
00520       // call reset again
00521       initial_params::reset(dvar_vector(x));
00522       report_function_minimizer_stats();
00523       if (quit_flag=='Q') break;
00524       if (!quit_flag || quit_flag == 'N')
00525       {
00526         initial_params::current_phase++;
00527       }
00528     }
00529     if (initial_params::current_phase >
00530       initial_params::max_number_phases)
00531     {
00532       initial_params::current_phase =
00533         initial_params::max_number_phases;
00534     }
00535     tracing_message(traceflag,"N2");
00536   }
00537 #  if defined(USE_LAPLACE)
00538   void function_minimizer::set_multinomial_weights(dvector& d)
00539   {
00540     multinomial_weights=new dvector(d);
00541   }
00542 #endif
00543 
00544   function_minimizer::function_minimizer(long int sz)
00545   {
00546 #  if defined(USE_LAPLACE)
00547     lapprox=0;
00548     multinomial_weights=0;
00549     //cout << lapprox << endl;
00550 #  endif
00551     maxfn  = 1000;
00552     iprint = 1;
00553     crit   = 0.0001;
00554     imax   = 30;
00555     dfn    = 0.01;
00556     iexit  = 0;
00557     ihflag = 0;
00558     ihang  = 0;
00559     scroll_flag = 1;
00560     maxfn_flag=0;
00561     quit_flag=0;
00562     min_improve=.0;
00563     negdirections=0;
00564     spminflag=0;
00565     repeatminflag=0;
00566 
00567     int ssz;
00568 
00569     int nopt=get_option_number("-ams",
00570       "-ams option needs positive integer -- ignored",ssz);
00571     if (nopt>-1 && ssz>0) {
00572       sz=ssz;
00573     }
00574 
00575 #ifdef __BORLANDC__
00576     long int lssz;
00577 #else
00578     long long int lssz;
00579 #endif
00580     nopt=get_option_number("-cbs",
00581       "-cbs option needs positive integer -- ignored",lssz);
00582     if (nopt>-1 && lssz>0) {
00583       gradient_structure::set_CMPDIF_BUFFER_SIZE(lssz);
00584     }
00585 
00586     nopt=get_option_number("-gbs",
00587       "-gbs option needs positive integer -- ignored",lssz);
00588     if (nopt>-1 && lssz>0) {
00589       gradient_structure::set_GRADSTACK_BUFFER_SIZE
00590         (lssz/sizeof(grad_stack_entry));
00591     }
00592 
00593     if (!sz)
00594     {
00595       pgs = new gradient_structure;
00596     }
00597     else
00598     {
00599       pgs = new gradient_structure(sz);
00600     }
00601   }
00602 
00606 function_minimizer::~function_minimizer()
00607 {
00608 #if defined(USE_LAPLACE)
00609   if (multinomial_weights)
00610   {
00611     delete multinomial_weights;
00612     multinomial_weights = 0;
00613   }
00614   if (lapprox)
00615   {
00616     delete lapprox;
00617     lapprox = 0;
00618   }
00619 #endif
00620   delete pgs;
00621   pgs = 0;
00622   if (negdirections)
00623   {
00624     delete negdirections;
00625     negdirections = 0;
00626   }
00627 }
00628 
00629 void function_minimizer::set_initial_simplex(const dmatrix& _p,
00630   const dvector& _y, int nvar, const dvector& x, double delta)
00631   {
00632     dvector& y=(dvector&) _y;
00633     dmatrix& p=(dmatrix&) _p;
00634     int i;
00635     p(1)=x;
00636     for (i=2;i<=nvar+1;i++)
00637     {
00638       p(i)=x;
00639       p(i,i-1)+=delta;
00640     }
00641     dvector xx(1,nvar);
00642     double vf=0;
00643     for (i=1;i<=nvar+1;i++)
00644     {
00645       xx=p(i);
00646       vf=value(initial_params::reset(dvar_vector(xx)));
00647       *objective_function_value::pobjfun=0.0;
00648       userfunction();
00649       vf+=value(*objective_function_value::pobjfun);
00650       y(i)=vf;
00651     }
00652   }
00653 
00654 int get_option_number(const char * option_name,const char * error_message,
00655   int& option_value)
00656 {
00657   int on1;
00658   int nopt = 0;
00659   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00660   {
00661     if (!nopt)
00662     {
00663       if (ad_printf)
00664         (*ad_printf)("%s\n",error_message);
00665       else
00666         cerr << error_message << endl;
00667       on1=-1;
00668     }
00669     else
00670     {
00671       option_value=atoi(ad_comm::argv[on1+1]);
00672     }
00673   }
00674   return on1;
00675 }
00676 
00677 int get_option_number(const char * option_name,const char * error_message,
00678 #ifdef __BORLANDC__
00679   long int& option_value)
00680 #else
00681   long long int& option_value)
00682 #endif
00683 {
00684   int on1;
00685   int nopt = 0;
00686   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00687   {
00688     if (!nopt)
00689     {
00690       if (ad_printf)
00691         (*ad_printf)("%s\n",error_message);
00692       else
00693         cerr << error_message << endl;
00694       on1=-1;
00695     }
00696     else
00697     {
00698 #if defined(__BORLANDC__) || defined(_MSC_VER)
00699       option_value=atol(ad_comm::argv[on1+1]);
00700 #else
00701       option_value=atoll(ad_comm::argv[on1+1]);
00702 #endif
00703     }
00704   }
00705   return on1;
00706 }
00707 
00708 #if defined(USE_LAPLACE)
00709 
00710 void function_minimizer::other_separable_stuff_begin(void)
00711 {
00712   if (lapprox)
00713   {
00714     lapprox->separable_calls_counter++;
00715     /*
00716     lapprox->separable_call_level++;
00717     //lapprox->build_up_nested_shape();
00718     lapprox->nested_separable_calls_counter
00719       (lapprox->separable_call_level)++;
00720     //clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00721     lapprox->nested_tree_position(lapprox->separable_call_level)++;
00722     */
00723   }
00724 }
00725 
00726 void function_minimizer::other_separable_stuff_end(void)
00727 {
00728   /*
00729   if (lapprox)
00730   {
00731     lapprox->build_up_nested_shape();
00732     clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00733     lapprox->separable_call_level--;
00734   }
00735   */
00736 }
00737 
00738 
00739 void function_minimizer::begin_gauss_hermite_stuff(void)
00740 {
00741   int nsc=lapprox->separable_calls_counter;
00742   int is=0;
00743   if (lapprox->gh->mi==0)
00744   {
00745     is=lapprox->gh->is;
00746   }
00747   else
00748   {
00749     is=lapprox->gh->mi->get_offset()+1;
00750   }
00751   lapprox->gh->gauss_hermite_values(nsc,is)=
00752     *objective_function_value::pobjfun;
00753 }
00754 
00755 void function_minimizer::start_get_importance_sampling_comnponent(void)
00756 {
00757   int nsc=lapprox->separable_calls_counter;
00758   int isc=lapprox->importance_sampling_counter;
00759   (*lapprox->importance_sampling_components)(nsc,isc)=
00760      *objective_function_value::pobjfun;
00761 }
00762 
00763 void function_minimizer::end_get_importance_sampling_comnponent(void)
00764 {
00765   int nsc=lapprox->separable_calls_counter;
00766   int is=lapprox->importance_sampling_counter;
00767   if (lapprox->saddlepointflag==2)
00768   {
00769     (*lapprox->importance_sampling_components)(nsc,is)=
00770        (-1)* *objective_function_value::pobjfun-
00771        (*lapprox->importance_sampling_components)(nsc,is);
00772   }
00773   else
00774   {
00775     (*lapprox->importance_sampling_components)(nsc,is)=
00776        *objective_function_value::pobjfun-
00777        (*lapprox->importance_sampling_components)(nsc,is);
00778   }
00779 }
00780 
00781 void function_minimizer::begin_funnel_stuff(void)
00782 {
00783   lapprox->separable_calls_counter++;
00784   if (lapprox->hesstype==2)
00785   {
00786     if (lapprox->in_gauss_hermite_phase)
00787     {
00788        begin_gauss_hermite_stuff();
00789     }
00790     else if (lapprox->num_importance_samples &&
00791       lapprox->importance_sampling_components)
00792     {
00793       if (lapprox->block_diagonal_flag==2)
00794       {
00795         start_get_importance_sampling_comnponent();
00796       }
00797     }
00798   }
00799 }
00800 
00801 void function_minimizer::get_function_difference(void)
00802 {
00803   int nsc=lapprox->separable_calls_counter;
00804   (*(lapprox->separable_function_difference))(nsc)=
00805     value(*objective_function_value::pobjfun);
00806     value(*objective_function_value::pobjfun)=0.0;
00807 }
00808 void function_minimizer::end_df1b2_funnel_stuff(void)
00809 {
00810   if (lapprox->in_gauss_hermite_phase)
00811   {
00812     end_gauss_hermite_stuff();
00813   }
00814   else
00815   {
00816     if (lapprox->hesstype==2)
00817     {
00818       if (lapprox->num_importance_samples &&
00819         lapprox->importance_sampling_components)
00820       {
00821         if (lapprox->block_diagonal_flag==2)
00822         {
00823           end_get_importance_sampling_comnponent();
00824         }
00825       }
00826       if (!lapprox->no_function_component_flag)
00827       {
00828         if (lapprox->saddlepointflag!=2)
00829         {
00830           get_function_difference();
00831         }
00832         else if (inner_opt()!=0)
00833         {
00834           get_function_difference();
00835         }
00836       }
00837     }
00838   }
00839 }
00840 
00841 void function_minimizer::end_gauss_hermite_stuff(void)
00842 {
00843   int nsc=lapprox->separable_calls_counter;
00844   int is=0;
00845   if (lapprox->gh->mi==0)
00846   {
00847     is=lapprox->gh->is;
00848   }
00849   else
00850   {
00851     is=lapprox->gh->mi->get_offset()+1;
00852   }
00853   lapprox->gh->gauss_hermite_values(nsc,is)=
00854     *objective_function_value::pobjfun-
00855     lapprox->gh->gauss_hermite_values(nsc,is);
00856 }
00857 
00858 void print_is_diagnostics(laplace_approximation_calculator *lapprox)
00859 {
00860   if (lapprox->is_diagnostics_flag)
00861   {
00862     if (lapprox->importance_sampling_values)
00863     {
00864       int mmin=lapprox->importance_sampling_values->indexmin();
00865       int mmax=lapprox->importance_sampling_values->indexmax();
00866       double mn= mean(*lapprox->importance_sampling_values);
00867       dmatrix tmp(1,2,mmin,mmax);
00868       tmp(2)=*lapprox->importance_sampling_values-mn;
00869       tmp(1).fill_seqadd(1,1);
00870       tmp=trans(sort(trans(tmp),2));
00871       ofstream ofs("is_diagnostics");
00872       ofs << "number of importance samples "
00873           << lapprox->num_importance_samples << endl;
00874       ofs << "importance_sampling_values" << endl;
00875       ofs << *lapprox->importance_sampling_values << endl<< endl;;
00876       ofs << "normalized importance_sampling_values" << endl;
00877       ofs << *lapprox->importance_sampling_values-mn << endl<< endl;;
00878       ofs << "sorted normalized importance_sampling_values" << endl;
00879       ofs << setw(9) << tmp << endl<< endl;;
00880       ofs << "epsilon(1).indexmax()  "
00881           << lapprox->epsilon(1).indexmax() << endl;
00882       ofs << lapprox->epsilon << endl;
00883       dmatrix plotstuff(1,2,mmin,mmax);
00884       plotstuff(1)=*lapprox->importance_sampling_weights;
00885       plotstuff(2)=*lapprox->importance_sampling_values-mn;
00886       ofs << " weight   value " << endl;
00887       ofs << setw(9) << sort(trans(plotstuff),2) << endl;
00888     }
00889   }
00890 }
00891 #endif // #if defined(USE_LAPLACE)