ADMB Documentation  11.1.2491
 All Classes Files Functions Variables Typedefs Friends Defines
xmodelm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xmodelm3.cpp 2470 2014-10-03 23:21:25Z 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 #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     int bandwidth=0;
00171     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-bw",nopt))>-1)
00172     {
00173       if (!nopt)
00174       {
00175         cerr << "Usage -bw option needs number  -- ignored" << endl;
00176       }
00177       else
00178       {
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           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       size_t _nvar=initial_params::nvarcalc(); // get the number of active
00226              // parameters
00227       if (_nvar < 1)
00228       {
00229         cerr << "Error -- no active parameters. There must be at least 1"
00230              << endl;
00231         exit(1);
00232       }
00233       const int nvar = (int)_nvar;
00234       dvector g(1,nvar);
00235       independent_variables x(1,nvar);
00236       tracing_message(traceflag,"B2");
00237       initial_params::xinit(x);    // get the initial values into the
00238 
00240       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-uhess"))>-1)
00241       {
00242         int ierr=0;
00243         ifstream ifs("vector");
00244         if (!ifs)
00245         {
00246           cerr << "couldn't open file vector" << endl;
00247           ierr=1;
00248         }
00249         dvector zz(1,x.indexmax());
00250         if(ierr==0)
00251         {
00252           ifs >> zz;
00253           if (!ifs)
00254           {
00255             cerr << "couldn't read vector" << endl;
00256             ierr=1;
00257           }
00258         }
00259         if (ierr==0)
00260         {
00261           dvector xsave(1,x.indexmax());
00262           for (;;)
00263           {
00264             double delta=0;
00265             cout << "enter delta" << endl;
00266             cin >> delta;
00267             xsave=x;
00268             x+=delta*zz;
00269             initial_params::reset(x);    // get the initial values into the
00270             userfunction();
00271             x=xsave;
00272           }
00273         }
00274       }
00275 
00276       double f=0.0;
00277 
00278       int simpflag = -1;
00279       if ( (simpflag=option_match(ad_comm::argc,ad_comm::argv,"-simplex"))>-1)
00280       {
00281         gradient_structure::set_NO_DERIVATIVES();
00282         double delta=1.e-4;
00283         double ftol=1.e-16;
00284         dmatrix p(1,nvar+1,1,nvar);
00285         dvector y(1,nvar+1);
00286         set_initial_simplex(p,y,nvar,x,delta);
00287         adamoeba(p,y,nvar,ftol,maxfn);
00288         double ymin=min(y);
00289         for (int i=1;i<=nvar+1;i++)
00290         if (ymin==y(i))
00291         {
00292           x=p(i);
00293           break;
00294         }
00295         cerr << "The -simplex option is deprecated. The user should port "
00296              << "to the -neldmead option." << endl;
00297       }
00298       if ( (simpflag=option_match(ad_comm::argc,ad_comm::argv,"-neldmead"))>-1)
00299       {
00300         gradient_structure::set_NO_DERIVATIVES();
00301         double delta=1.e-4;
00302         double ftol=1.e-16;
00303         dvector mincords = x;
00304         double ynewlo;
00305         double* pynewlo = &ynewlo;
00306         int icount, numres, ifault;
00307         int* picount = &icount;
00308         int* pnumres = &numres;
00309         int* pifault = &ifault;
00310         neldmead(nvar,mincords,mincords,pynewlo,ftol,delta,picount,pnumres,
00311           pifault);
00312         x = mincords;
00313       }
00314       int lmnflag = -1;
00315       int lmnsteps=10;
00316       if ( (lmnflag=option_match(ad_comm::argc,ad_comm::argv,"-lmn",nopt))>-1)
00317       {
00318         if (random_effects_flag)
00319         {
00320           cerr << "At present you can not use the -lmn option for the outer"
00321                << endl << " optimiation in a random-effects model" << endl;
00322           ad_exit(1);
00323         }
00324         if (!nopt)
00325         {
00326           cerr << "Usage -lmn option needs integer  -- set to default 10.\n";
00327         }
00328         else
00329         {
00330           int jj=atoi(ad_comm::argv[lmnflag+1]);
00331           if (jj<=0)
00332           {
00333             cerr << "Usage -lmn option needs positive integer"
00334                     " -- set to default 10.\n";
00335           }
00336           else
00337           {
00338             lmnsteps=jj;
00339           }
00340         }
00341       }
00342       if (lmnflag==-1)
00343       {
00344         // *********************************************************
00345         // block for quasi newton minimization
00346         if (negative_eigenvalue_flag)
00347         {
00348           trust_region_update(nvar,_crit,x,g,f);
00349         }
00350 #if defined(USE_ADPVM)
00351         if (ad_comm::pvm_manager)
00352         {
00353           if (random_effects_flag)
00354           {
00355             if (maxfn>0)
00356             {
00357               switch (ad_comm::pvm_manager->mode)
00358               {
00359               case 1: // master
00360                 quasi_newton_block_pvm_master_random_effects(nvar,_crit,x,g,f);
00361                 break;
00362               case 2: // slave
00363               // these don't exist yet
00364                 function_evaluation_block_pvm_slave_random_effects
00365                   (nvar,_crit,x,g,f);
00366                 break;
00367               default:
00368                 cerr << "error illega value for pvm_manager->mode" << endl;
00369                 exit(1);
00370               }
00371             }
00372           }
00373           else
00374           {
00375             if (maxfn>0)
00376             {
00377               switch (ad_comm::pvm_manager->mode)
00378               {
00379               case 1: // master
00380                 quasi_newton_block_pvm_master(nvar,_crit,x,g,f);
00381                 break;
00382               case 2: // slave
00383                 function_evaluation_block_pvm_slave();
00384                 break;
00385               default:
00386                 cerr << "error illega value for pvm_manager->mode" << endl;
00387                 exit(1);
00388               }
00389             }
00390           }
00391         }
00392         else
00393 #endif
00394         {
00395           do
00396           {
00397             if (spminflag)
00398             {
00399               repeatminflag=1;
00400               spminflag=0;
00401             }
00402             else
00403             {
00404               repeatminflag=0;
00405             }
00406             if (maxfn>0)
00407             {
00408               int lmnflag=-1;
00409               int nsteps=5;
00410               if ( (lmnflag=option_match(ad_comm::argc,ad_comm::argv,
00411                  "-lmn2",nopt))>-1)
00412               {
00413                 if (!nopt)
00414                 {
00415                   cerr << "Usage -lmn option needs integer"
00416                      "  -- set to default 5" << endl;
00417                 }
00418                 else
00419                 {
00420                   int jj=atoi(ad_comm::argv[lmnflag+1]);
00421                   if (jj<=0)
00422                   {
00423                     cerr << "Usage -lmn option needs positive integer "
00424                      " -- set to default 5" << endl;
00425                   }
00426                   else
00427                   {
00428                     nsteps=jj;
00429                   }
00430                 }
00431               }
00432               if (lmnflag<0)
00433               {
00434                 quasi_newton_block(nvar,_crit,x,g,f);
00435               }
00436               else
00437               {
00438                 limited_memory_quasi_newton_block(nvar,_crit,x,g,f,nsteps);
00439               }
00440             }
00441           }
00442           while(repeatminflag);
00443         }
00444       } // end block for quasi newton minimization
00445       else
00446       {  // block for limited memory quasi newton minimization
00447         if (maxfn>0)
00448         {
00449           function_minimizer::limited_memory_quasi_newton(x,lmnsteps);
00450         }
00451       }
00452       // end block for limited memory quasi newton minimization
00453       // *********************************************************
00454       tracing_message(traceflag,"M2");
00455 
00456       gradient_structure::set_NO_DERIVATIVES();
00457       initial_params::reset(dvar_vector(x));
00458       *objective_function_value::pobjfun=0.0;
00459       if (!random_effects_flag || !lapprox)
00460       {
00461 #if defined(USE_ADPVM)
00462         if (ad_comm::pvm_manager)
00463         {
00464           switch (ad_comm::pvm_manager->mode)
00465           {
00466           case 1:
00467             pvm_master_function_evaluation_no_derivatives(f,x,nvar);
00468             *objective_function_value::pobjfun=f;
00469             break;
00470           case 2:
00471             pvm_slave_function_evaluation_no_derivatives();
00472             break;
00473           default:
00474             cerr << "Illegal value for ad_comm::pvm_manager->mode" << endl;
00475             ad_exit(1);
00476           }
00477         }
00478         else
00479         {
00480 #endif  //#if defined(USE_ADPVM)
00481           userfunction();
00482 #if defined(USE_ADPVM)
00483         }
00484 #endif  //#if defined(USE_ADPVM)
00485       }
00486       else
00487       {
00488         (*lapprox)(x,f,this);
00489         *objective_function_value::pobjfun=f;
00490         initial_params::set_inactive_only_random_effects();
00491         print_is_diagnostics(lapprox);
00492       }
00493       initial_params::save();
00494       report(g);
00495       // in case the user changes some initial_params in the report section
00496       // call reset again
00497       initial_params::reset(dvar_vector(x));
00498       report_function_minimizer_stats();
00499       if (quit_flag=='Q') break;
00500       if (!quit_flag || quit_flag == 'N')
00501       {
00502         initial_params::current_phase++;
00503       }
00504     }
00505     if (initial_params::current_phase >
00506       initial_params::max_number_phases)
00507     {
00508       initial_params::current_phase =
00509         initial_params::max_number_phases;
00510     }
00511     tracing_message(traceflag,"N2");
00512   }
00513   void function_minimizer::set_multinomial_weights(dvector& d)
00514   {
00515     multinomial_weights=new dvector(d);
00516   }
00517 
00518   function_minimizer::function_minimizer(long int sz)
00519   {
00520     lapprox=0;
00521     multinomial_weights=0;
00522     //cout << lapprox << endl;
00523     maxfn  = 1000;
00524     iprint = 1;
00525     crit   = 0.0001;
00526     imax   = 30;
00527     dfn    = 0.01;
00528     iexit  = 0;
00529     ihflag = 0;
00530     ihang  = 0;
00531     scroll_flag = 1;
00532     maxfn_flag=0;
00533     quit_flag=0;
00534     min_improve=.0;
00535     negdirections=0;
00536     spminflag=0;
00537     repeatminflag=0;
00538 
00539     int ssz;
00540 
00541     int nopt=get_option_number("-ams",
00542       "-ams option needs positive integer -- ignored",ssz);
00543     if (nopt>-1 && ssz>0) {
00544       sz=ssz;
00545     }
00546 
00547 #ifdef __BORLANDC__
00548     long int lssz;
00549 #else
00550     long long int lssz;
00551 #endif
00552     nopt=get_option_number("-cbs",
00553       "-cbs option needs positive integer -- ignored",lssz);
00554     if (nopt>-1 && lssz>0) {
00555       const size_t size = (size_t)lssz;
00556       gradient_structure::set_CMPDIF_BUFFER_SIZE(size);
00557     }
00558 
00559     nopt=get_option_number("-gbs",
00560       "-gbs option needs positive integer -- ignored",lssz);
00561     if (nopt>-1 && lssz>0) {
00562       const size_t size = (size_t)lssz / sizeof(grad_stack_entry);
00563       gradient_structure::set_GRADSTACK_BUFFER_SIZE(size);
00564     }
00565 
00566     if (!sz)
00567     {
00568       pgs = new gradient_structure;
00569     }
00570     else
00571     {
00572       pgs = new gradient_structure(sz);
00573     }
00574   }
00575 
00579 function_minimizer::~function_minimizer()
00580 {
00581   if (multinomial_weights)
00582   {
00583     delete multinomial_weights;
00584     multinomial_weights = 0;
00585   }
00586   if (lapprox)
00587   {
00588     delete lapprox;
00589     lapprox = 0;
00590   }
00591   delete pgs;
00592   pgs = 0;
00593   if (negdirections)
00594   {
00595     delete negdirections;
00596     negdirections = 0;
00597   }
00598 }
00599 
00600 void function_minimizer::set_initial_simplex(const dmatrix& _p,
00601   const dvector& _y, int nvar, const dvector& x, double delta)
00602   {
00603     dvector& y=(dvector&) _y;
00604     dmatrix& p=(dmatrix&) _p;
00605     int i;
00606     p(1)=x;
00607     for (i=2;i<=nvar+1;i++)
00608     {
00609       p(i)=x;
00610       p(i,i-1)+=delta;
00611     }
00612     dvector xx(1,nvar);
00613     double vf=0;
00614     for (i=1;i<=nvar+1;i++)
00615     {
00616       xx=p(i);
00617       vf=value(initial_params::reset(dvar_vector(xx)));
00618       *objective_function_value::pobjfun=0.0;
00619       userfunction();
00620       vf+=value(*objective_function_value::pobjfun);
00621       y(i)=vf;
00622     }
00623   }
00624 
00625 int get_option_number(const char * option_name,const char * error_message,
00626   int& option_value)
00627 {
00628   int on1;
00629   int nopt = 0;
00630   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00631   {
00632     if (!nopt)
00633     {
00634       if (ad_printf)
00635         (*ad_printf)("%s\n",error_message);
00636       else
00637         cerr << error_message << endl;
00638       on1=-1;
00639     }
00640     else
00641     {
00642       option_value=atoi(ad_comm::argv[on1+1]);
00643     }
00644   }
00645   return on1;
00646 }
00647 
00648 int get_option_number(const char * option_name,const char * error_message,
00649 #ifdef __BORLANDC__
00650   long int& option_value)
00651 #else
00652   long long int& option_value)
00653 #endif
00654 {
00655   int on1;
00656   int nopt = 0;
00657   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00658   {
00659     if (!nopt)
00660     {
00661       if (ad_printf)
00662         (*ad_printf)("%s\n",error_message);
00663       else
00664         cerr << error_message << endl;
00665       on1=-1;
00666     }
00667     else
00668     {
00669 #if defined(__BORLANDC__) || defined(_MSC_VER)
00670       option_value=atol(ad_comm::argv[on1+1]);
00671 #else
00672       option_value=atoll(ad_comm::argv[on1+1]);
00673 #endif
00674     }
00675   }
00676   return on1;
00677 }
00678 
00679 
00680 void function_minimizer::other_separable_stuff_begin(void)
00681 {
00682   if (lapprox)
00683   {
00684     lapprox->separable_calls_counter++;
00685     /*
00686     lapprox->separable_call_level++;
00687     //lapprox->build_up_nested_shape();
00688     lapprox->nested_separable_calls_counter
00689       (lapprox->separable_call_level)++;
00690     //clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00691     lapprox->nested_tree_position(lapprox->separable_call_level)++;
00692     */
00693   }
00694 }
00695 
00696 void function_minimizer::other_separable_stuff_end(void)
00697 {
00698   /*
00699   if (lapprox)
00700   {
00701     lapprox->build_up_nested_shape();
00702     clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00703     lapprox->separable_call_level--;
00704   }
00705   */
00706 }
00707 
00708 
00709 void function_minimizer::begin_gauss_hermite_stuff(void)
00710 {
00711   int nsc=lapprox->separable_calls_counter;
00712   int is=0;
00713   if (lapprox->gh->mi==0)
00714   {
00715     is=lapprox->gh->is;
00716   }
00717   else
00718   {
00719     is=lapprox->gh->mi->get_offset()+1;
00720   }
00721   lapprox->gh->gauss_hermite_values(nsc,is)=
00722     *objective_function_value::pobjfun;
00723 }
00724 
00725 void function_minimizer::start_get_importance_sampling_comnponent(void)
00726 {
00727   int nsc=lapprox->separable_calls_counter;
00728   int isc=lapprox->importance_sampling_counter;
00729   (*lapprox->importance_sampling_components)(nsc,isc)=
00730      *objective_function_value::pobjfun;
00731 }
00732 
00733 void function_minimizer::end_get_importance_sampling_comnponent(void)
00734 {
00735   int nsc=lapprox->separable_calls_counter;
00736   int is=lapprox->importance_sampling_counter;
00737   if (lapprox->saddlepointflag==2)
00738   {
00739     (*lapprox->importance_sampling_components)(nsc,is)=
00740        (-1)* *objective_function_value::pobjfun-
00741        (*lapprox->importance_sampling_components)(nsc,is);
00742   }
00743   else
00744   {
00745     (*lapprox->importance_sampling_components)(nsc,is)=
00746        *objective_function_value::pobjfun-
00747        (*lapprox->importance_sampling_components)(nsc,is);
00748   }
00749 }
00750 
00751 void function_minimizer::begin_funnel_stuff(void)
00752 {
00753   lapprox->separable_calls_counter++;
00754   if (lapprox->hesstype==2)
00755   {
00756     if (lapprox->in_gauss_hermite_phase)
00757     {
00758        begin_gauss_hermite_stuff();
00759     }
00760     else if (lapprox->num_importance_samples &&
00761       lapprox->importance_sampling_components)
00762     {
00763       if (lapprox->block_diagonal_flag==2)
00764       {
00765         start_get_importance_sampling_comnponent();
00766       }
00767     }
00768   }
00769 }
00770 
00771 void function_minimizer::get_function_difference(void)
00772 {
00773   int nsc=lapprox->separable_calls_counter;
00774   (*(lapprox->separable_function_difference))(nsc)=
00775     value(*objective_function_value::pobjfun);
00776     value(*objective_function_value::pobjfun)=0.0;
00777 }
00778 void function_minimizer::end_df1b2_funnel_stuff(void)
00779 {
00780   if (lapprox->in_gauss_hermite_phase)
00781   {
00782     end_gauss_hermite_stuff();
00783   }
00784   else
00785   {
00786     if (lapprox->hesstype==2)
00787     {
00788       if (lapprox->num_importance_samples &&
00789         lapprox->importance_sampling_components)
00790       {
00791         if (lapprox->block_diagonal_flag==2)
00792         {
00793           end_get_importance_sampling_comnponent();
00794         }
00795       }
00796       if (!lapprox->no_function_component_flag)
00797       {
00798         if (lapprox->saddlepointflag!=2)
00799         {
00800           get_function_difference();
00801         }
00802         else if (inner_opt()!=0)
00803         {
00804           get_function_difference();
00805         }
00806       }
00807     }
00808   }
00809 }
00810 
00811 void function_minimizer::end_gauss_hermite_stuff(void)
00812 {
00813   int nsc=lapprox->separable_calls_counter;
00814   int is=0;
00815   if (lapprox->gh->mi==0)
00816   {
00817     is=lapprox->gh->is;
00818   }
00819   else
00820   {
00821     is=lapprox->gh->mi->get_offset()+1;
00822   }
00823   lapprox->gh->gauss_hermite_values(nsc,is)=
00824     *objective_function_value::pobjfun-
00825     lapprox->gh->gauss_hermite_values(nsc,is);
00826 }
00827 
00828 void print_is_diagnostics(laplace_approximation_calculator *lapprox)
00829 {
00830   if (lapprox->is_diagnostics_flag)
00831   {
00832     if (lapprox->importance_sampling_values)
00833     {
00834       int mmin=lapprox->importance_sampling_values->indexmin();
00835       int mmax=lapprox->importance_sampling_values->indexmax();
00836       double mn= mean(*lapprox->importance_sampling_values);
00837       dmatrix tmp(1,2,mmin,mmax);
00838       tmp(2)=*lapprox->importance_sampling_values-mn;
00839       tmp(1).fill_seqadd(1,1);
00840       tmp=trans(sort(trans(tmp),2));
00841       ofstream ofs("is_diagnostics");
00842       ofs << "number of importance samples "
00843           << lapprox->num_importance_samples << endl;
00844       ofs << "importance_sampling_values" << endl;
00845       ofs << *lapprox->importance_sampling_values << endl<< endl;;
00846       ofs << "normalized importance_sampling_values" << endl;
00847       ofs << *lapprox->importance_sampling_values-mn << endl<< endl;;
00848       ofs << "sorted normalized importance_sampling_values" << endl;
00849       ofs << setw(9) << tmp << endl<< endl;;
00850       ofs << "epsilon(1).indexmax()  "
00851           << lapprox->epsilon(1).indexmax() << endl;
00852       ofs << lapprox->epsilon << endl;
00853       dmatrix plotstuff(1,2,mmin,mmax);
00854       plotstuff(1)=*lapprox->importance_sampling_weights;
00855       plotstuff(2)=*lapprox->importance_sampling_values-mn;
00856       ofs << " weight   value " << endl;
00857       ofs << setw(9) << sort(trans(plotstuff),2) << endl;
00858     }
00859   }
00860 }