ADMB Documentation  11.5.3150
 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   ifn(0)
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 }