ADMB Documentation  11.1x.2723
 All Classes Files Functions Variables Typedefs Friends Defines
xmodelm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xmodelm3.cpp 2688 2014-11-18 21:39:48Z 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     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   {
00519     lapprox=0;
00520     multinomial_weights=0;
00521     //cout << lapprox << endl;
00522     maxfn  = 1000;
00523     iprint = 1;
00524     crit   = 0.0001;
00525     imax   = 30;
00526     dfn    = 0.01;
00527     iexit  = 0;
00528     ihflag = 0;
00529     ihang  = 0;
00530     scroll_flag = 1;
00531     maxfn_flag=0;
00532     quit_flag=0;
00533     min_improve=.0;
00534     negdirections=0;
00535     spminflag=0;
00536     repeatminflag=0;
00537 
00538     int ssz;
00539 
00540     int nopt=get_option_number("-ams",
00541       "-ams option needs positive integer -- ignored",ssz);
00542     if (nopt>-1 && ssz>0) {
00543       sz=ssz;
00544     }
00545 
00546 #ifdef __BORLANDC__
00547     long int lssz;
00548 #else
00549     long long int lssz;
00550 #endif
00551     nopt=get_option_number("-cbs",
00552       "-cbs option needs positive integer -- ignored",lssz);
00553     if (nopt>-1 && lssz>0) {
00554       const size_t size = (size_t)lssz;
00555       gradient_structure::set_CMPDIF_BUFFER_SIZE(size);
00556     }
00557 
00558     nopt=get_option_number("-gbs",
00559       "-gbs option needs positive integer -- ignored",lssz);
00560     if (nopt>-1 && lssz>0) {
00561       const size_t size = (size_t)lssz / sizeof(grad_stack_entry);
00562       gradient_structure::set_GRADSTACK_BUFFER_SIZE(size);
00563     }
00564 
00565     if (!sz)
00566     {
00567       pgs = new gradient_structure;
00568     }
00569     else
00570     {
00571       pgs = new gradient_structure(sz);
00572     }
00573   }
00574 
00578 function_minimizer::~function_minimizer()
00579 {
00580   if (multinomial_weights)
00581   {
00582     delete multinomial_weights;
00583     multinomial_weights = 0;
00584   }
00585   if (lapprox)
00586   {
00587     delete lapprox;
00588     lapprox = 0;
00589   }
00590   delete pgs;
00591   pgs = 0;
00592   if (negdirections)
00593   {
00594     delete negdirections;
00595     negdirections = 0;
00596   }
00597 }
00598 
00599 void function_minimizer::set_initial_simplex(const dmatrix& _p,
00600   const dvector& _y, int nvar, const dvector& x, double delta)
00601   {
00602     dvector& y=(dvector&) _y;
00603     dmatrix& p=(dmatrix&) _p;
00604     int i;
00605     p(1)=x;
00606     for (i=2;i<=nvar+1;i++)
00607     {
00608       p(i)=x;
00609       p(i,i-1)+=delta;
00610     }
00611     dvector xx(1,nvar);
00612     double vf=0;
00613     for (i=1;i<=nvar+1;i++)
00614     {
00615       xx=p(i);
00616       vf=value(initial_params::reset(dvar_vector(xx)));
00617       *objective_function_value::pobjfun=0.0;
00618       userfunction();
00619       vf+=value(*objective_function_value::pobjfun);
00620       y(i)=vf;
00621     }
00622   }
00623 
00624 int get_option_number(const char * option_name,const char * error_message,
00625   int& option_value)
00626 {
00627   int on1;
00628   int nopt = 0;
00629   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00630   {
00631     if (!nopt)
00632     {
00633       if (ad_printf)
00634         (*ad_printf)("%s\n",error_message);
00635       else
00636         cerr << error_message << endl;
00637       on1=-1;
00638     }
00639     else
00640     {
00641       option_value=atoi(ad_comm::argv[on1+1]);
00642     }
00643   }
00644   return on1;
00645 }
00646 
00647 int get_option_number(const char * option_name,const char * error_message,
00648 #ifdef __BORLANDC__
00649   long int& option_value)
00650 #else
00651   long long int& option_value)
00652 #endif
00653 {
00654   int on1;
00655   int nopt = 0;
00656   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00657   {
00658     if (!nopt)
00659     {
00660       if (ad_printf)
00661         (*ad_printf)("%s\n",error_message);
00662       else
00663         cerr << error_message << endl;
00664       on1=-1;
00665     }
00666     else
00667     {
00668 #if defined(__BORLANDC__) || defined(_MSC_VER)
00669       option_value=atol(ad_comm::argv[on1+1]);
00670 #else
00671       option_value=atoll(ad_comm::argv[on1+1]);
00672 #endif
00673     }
00674   }
00675   return on1;
00676 }
00677 
00678 
00679 void function_minimizer::other_separable_stuff_begin(void)
00680 {
00681   if (lapprox)
00682   {
00683     lapprox->separable_calls_counter++;
00684     /*
00685     lapprox->separable_call_level++;
00686     //lapprox->build_up_nested_shape();
00687     lapprox->nested_separable_calls_counter
00688       (lapprox->separable_call_level)++;
00689     //clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00690     lapprox->nested_tree_position(lapprox->separable_call_level)++;
00691     */
00692   }
00693 }
00694 
00695 void function_minimizer::other_separable_stuff_end(void)
00696 {
00697   /*
00698   if (lapprox)
00699   {
00700     lapprox->build_up_nested_shape();
00701     clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00702     lapprox->separable_call_level--;
00703   }
00704   */
00705 }
00706 
00707 
00708 void function_minimizer::begin_gauss_hermite_stuff(void)
00709 {
00710   int nsc=lapprox->separable_calls_counter;
00711   int is=0;
00712   if (lapprox->gh->mi==0)
00713   {
00714     is=lapprox->gh->is;
00715   }
00716   else
00717   {
00718     is=lapprox->gh->mi->get_offset()+1;
00719   }
00720   lapprox->gh->gauss_hermite_values(nsc,is)=
00721     *objective_function_value::pobjfun;
00722 }
00723 
00724 void function_minimizer::start_get_importance_sampling_comnponent(void)
00725 {
00726   int nsc=lapprox->separable_calls_counter;
00727   int isc=lapprox->importance_sampling_counter;
00728   (*lapprox->importance_sampling_components)(nsc,isc)=
00729      *objective_function_value::pobjfun;
00730 }
00731 
00732 void function_minimizer::end_get_importance_sampling_comnponent(void)
00733 {
00734   int nsc=lapprox->separable_calls_counter;
00735   int is=lapprox->importance_sampling_counter;
00736   if (lapprox->saddlepointflag==2)
00737   {
00738     (*lapprox->importance_sampling_components)(nsc,is)=
00739        (-1)* *objective_function_value::pobjfun-
00740        (*lapprox->importance_sampling_components)(nsc,is);
00741   }
00742   else
00743   {
00744     (*lapprox->importance_sampling_components)(nsc,is)=
00745        *objective_function_value::pobjfun-
00746        (*lapprox->importance_sampling_components)(nsc,is);
00747   }
00748 }
00749 
00750 void function_minimizer::begin_funnel_stuff(void)
00751 {
00752   lapprox->separable_calls_counter++;
00753   if (lapprox->hesstype==2)
00754   {
00755     if (lapprox->in_gauss_hermite_phase)
00756     {
00757        begin_gauss_hermite_stuff();
00758     }
00759     else if (lapprox->num_importance_samples &&
00760       lapprox->importance_sampling_components)
00761     {
00762       if (lapprox->block_diagonal_flag==2)
00763       {
00764         start_get_importance_sampling_comnponent();
00765       }
00766     }
00767   }
00768 }
00769 
00770 void function_minimizer::get_function_difference(void)
00771 {
00772   int nsc=lapprox->separable_calls_counter;
00773   (*(lapprox->separable_function_difference))(nsc)=
00774     value(*objective_function_value::pobjfun);
00775     value(*objective_function_value::pobjfun)=0.0;
00776 }
00777 void function_minimizer::end_df1b2_funnel_stuff(void)
00778 {
00779   if (lapprox->in_gauss_hermite_phase)
00780   {
00781     end_gauss_hermite_stuff();
00782   }
00783   else
00784   {
00785     if (lapprox->hesstype==2)
00786     {
00787       if (lapprox->num_importance_samples &&
00788         lapprox->importance_sampling_components)
00789       {
00790         if (lapprox->block_diagonal_flag==2)
00791         {
00792           end_get_importance_sampling_comnponent();
00793         }
00794       }
00795       if (!lapprox->no_function_component_flag)
00796       {
00797         if (lapprox->saddlepointflag!=2)
00798         {
00799           get_function_difference();
00800         }
00801         else if (inner_opt()!=0)
00802         {
00803           get_function_difference();
00804         }
00805       }
00806     }
00807   }
00808 }
00809 
00810 void function_minimizer::end_gauss_hermite_stuff(void)
00811 {
00812   int nsc=lapprox->separable_calls_counter;
00813   int is=0;
00814   if (lapprox->gh->mi==0)
00815   {
00816     is=lapprox->gh->is;
00817   }
00818   else
00819   {
00820     is=lapprox->gh->mi->get_offset()+1;
00821   }
00822   lapprox->gh->gauss_hermite_values(nsc,is)=
00823     *objective_function_value::pobjfun-
00824     lapprox->gh->gauss_hermite_values(nsc,is);
00825 }
00826 
00827 void print_is_diagnostics(laplace_approximation_calculator *lapprox)
00828 {
00829   if (lapprox->is_diagnostics_flag)
00830   {
00831     if (lapprox->importance_sampling_values)
00832     {
00833       int mmin=lapprox->importance_sampling_values->indexmin();
00834       int mmax=lapprox->importance_sampling_values->indexmax();
00835       double mn= mean(*lapprox->importance_sampling_values);
00836       dmatrix tmp(1,2,mmin,mmax);
00837       tmp(2)=*lapprox->importance_sampling_values-mn;
00838       tmp(1).fill_seqadd(1,1);
00839       tmp=trans(sort(trans(tmp),2));
00840       ofstream ofs("is_diagnostics");
00841       ofs << "number of importance samples "
00842           << lapprox->num_importance_samples << endl;
00843       ofs << "importance_sampling_values" << endl;
00844       ofs << *lapprox->importance_sampling_values << endl<< endl;;
00845       ofs << "normalized importance_sampling_values" << endl;
00846       ofs << *lapprox->importance_sampling_values-mn << endl<< endl;;
00847       ofs << "sorted normalized importance_sampling_values" << endl;
00848       ofs << setw(9) << tmp << endl<< endl;;
00849       ofs << "epsilon(1).indexmax()  "
00850           << lapprox->epsilon(1).indexmax() << endl;
00851       ofs << lapprox->epsilon << endl;
00852       dmatrix plotstuff(1,2,mmin,mmax);
00853       plotstuff(1)=*lapprox->importance_sampling_weights;
00854       plotstuff(2)=*lapprox->importance_sampling_values-mn;
00855       ofs << " weight   value " << endl;
00856       ofs << setw(9) << sort(trans(plotstuff),2) << endl;
00857     }
00858   }
00859 }