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