ADMB Documentation  11.1.2247
 All Classes Files Functions Variables Typedefs Friends Defines
xmodelm3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: xmodelm3.cpp 2244 2014-08-13 19:49:44Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <sstream>
00008 using std::istringstream;
00009 
00010 #include <admodel.h>
00011 
00012 #  include <df1b2fun.h>
00013 #  include <adrndeff.h>
00014 
00015 void check_java_flags(int& start_flag,int& quit_flag,int& der_flag,
00016   int& next_flag);
00017 
00018 void ad_update_function_minimizer_report(int feval,int iter,int phase,
00019   double fval, double gmax,const char * cbuf);
00020 void vm_initialize(void);
00021 
00022 void set_initial_simplex(const tdmatrix& p, const dvector& y,int nvar,
00023   const dvector& x, double delta);
00024 
00025 int get_option_number(const char * option_name,const char * error_message,
00026   int& option_value);
00027 
00028 int get_option_number(const char * option_name,const char * error_message,
00029 #ifdef __BORLANDC__
00030   long int& option_value);
00031 #else
00032   long long int& option_value);
00033 #endif
00034 
00035 extern int traceflag;
00036 
00037 void tracing_message(int traceflag,const char *s);
00038 
00039   int function_minimizer::inner_opt(void)
00040   {
00041     return inner_opt_flag;
00042   }
00043 
00044   int function_minimizer::inner_opt_flag=0;
00045 
00046 
00047   int function_minimizer::bad_step_flag=0;
00048 
00049   void function_minimizer::minimize(void)
00050   {
00051     int nopt=0;
00052     int on=0;
00053     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-shess"))>-1)
00054     {
00055       laplace_approximation_calculator::sparse_hessian_flag=1;
00056     }
00057     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-pis"))>-1)
00058     {
00059      laplace_approximation_calculator::print_importance_sampling_weights_flag=1;
00060     }
00061     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-sp"))>-1)
00062     {
00063       laplace_approximation_calculator::saddlepointflag=1;
00064     }
00065 #    if defined(__MINI_MAX__)
00066         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mm"))>-1)
00067         {
00068           laplace_approximation_calculator::saddlepointflag=2;
00069         }
00070 #    else
00071         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-mm"))>-1)
00072         {
00073            cerr << "option -mm MINI_MAX not defined " << endl;
00074            ad_exit(1);
00075         }
00076 #    endif
00077 
00078     //initial_params::read(); // read in the values for the initial parameters
00079     if (initial_params::restart_phase)
00080     {
00081       initial_params::current_phase = initial_params::restart_phase;
00082       initial_params::restart_phase=0;
00083     }
00084     int allphases=initial_params::max_number_phases;
00085     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxph",nopt))>-1)
00086     {
00087       if (!nopt)
00088       {
00089         cerr << "Usage -maxph option needs integer  -- ignored" << endl;
00090       }
00091       else
00092       {
00093         int jj=atoi(ad_comm::argv[on+1]);
00094         if (jj<=0)
00095         {
00096           cerr << "Usage -maxph option needs positive integer  -- ignored.\n";
00097         }
00098         else
00099         {
00100           if (jj>allphases)
00101           {
00102             allphases=jj;
00103           }
00104         }
00105       }
00106       if (allphases>initial_params::max_number_phases)
00107       {
00108         initial_params::max_number_phases=allphases;
00109       }
00110     }
00111     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-ndv",nopt))>-1)
00112     {
00113       if (!nopt)
00114       {
00115         cerr << "Usage -ndv option needs integer  -- ignored" << endl;
00116       }
00117       else
00118       {
00119         int jj=atoi(ad_comm::argv[on+1]);
00120         if (jj<=0)
00121         {
00122           cerr << "Usage -ndv option needs positive integer  -- ignored.\n";
00123         }
00124         else
00125         {
00126           gradient_structure::NUM_DEPENDENT_VARIABLES=jj;
00127         }
00128       }
00129     }
00130 
00131     // set the maximum number of function evaluations by command line
00132     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-maxfn",nopt))>-1)
00133     {
00134       if (!nopt)
00135       {
00136         cerr << "Usage -maxph option needs integer  -- ignored" << endl;
00137       }
00138       else
00139       {
00140         int _maxfn=atoi(ad_comm::argv[on+1]);
00141         if (_maxfn<0)
00142         {
00143           cerr << "Usage -maxfn option needs positive integer  -- ignored.\n";
00144         }
00145         else
00146         {
00147           maxfn=_maxfn;
00148         }
00149       }
00150     }
00151     double _crit=0;
00152     // set the maximum number of function evaluations by command line
00153     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-crit",nopt))>-1)
00154     {
00155       if (!nopt)
00156       {
00157         cerr << "Usage -crit option needs number  -- ignored" << endl;
00158       }
00159       else
00160       {
00161         istringstream ist(ad_comm::argv[on+1]);
00162         ist >> _crit;
00163 
00164         if (_crit<=0)
00165         {
00166           cerr << "Usage -crit option needs positive number  -- ignored.\n";
00167           _crit=0.0;
00168         }
00169       }
00170     }
00171     int bandwidth=0;
00172     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-bw",nopt))>-1)
00173     {
00174       if (!nopt)
00175       {
00176         cerr << "Usage -bw option needs number  -- ignored" << endl;
00177       }
00178       else
00179       {
00180         istringstream ist(ad_comm::argv[on+1]);
00181         ist >> bandwidth;
00182 
00183         if (bandwidth<=0)
00184         {
00185           cerr << "Usage -bw option needs positive number  -- ignored" << endl;
00186           bandwidth=0.0;
00187         }
00188         else
00189         {
00190           ad_comm::bandwidth=bandwidth;
00191         }
00192       }
00193     }
00194     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-phase"))>-1)
00195     {
00196       int jj=atoi(ad_comm::argv[on+1]);
00197       if (jj <=0)
00198       {
00199         cerr << " Invalid option following command line option -phase -- "
00200           << endl << " phase set equal to 1" << endl;
00201       }
00202       if (jj>allphases)
00203       {
00204         jj=allphases;
00205       }
00206       if (jj<=0)
00207       {
00208         jj=1;
00209       }
00210       initial_params::current_phase = jj;
00211       cout << "Set current phase to " << jj << endl;
00212     }
00213     if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-lapqd"))>-1)
00214     {
00215       ADqd_flag=1;
00216     }
00217 
00218     tracing_message(traceflag,"A2");
00219     while (initial_params::current_phase <= allphases)
00220     {
00221       between_phases_calculations();
00222 
00223       if (random_effects_flag)
00224         initial_params::set_inactive_random_effects();
00225 
00226       int nvar=initial_params::nvarcalc(); // get the number of active
00227              // parameters
00228       if (!nvar)
00229       {
00230         cerr << "Error -- no active parameters. There must be at least 1"
00231              << endl;
00232         exit(1);
00233       }
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       if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-uhess"))>-1)
00239       {
00240         int ierr=0;
00241         ifstream ifs("vector");
00242         if (!ifs)
00243         {
00244           cerr << "couldn't open file vector" << endl;
00245           ierr=1;
00246         }
00247         dvector zz(1,x.indexmax());
00248         if(ierr==0)
00249         {
00250           ifs >> zz;
00251           if (!ifs)
00252           {
00253             cerr << "couldn't read vector" << endl;
00254             ierr=1;
00255           }
00256         }
00257         if (ierr==0)
00258         {
00259           dvector xsave(1,x.indexmax());
00260           do
00261           {
00262             double delta=0;
00263             cout << "enter delta" << endl;
00264             cin >> delta;
00265             xsave=x;
00266             x+=delta*zz;
00267             initial_params::reset(x);    // get the initial values into the
00268             userfunction();
00269             x=xsave;
00270           }
00271           while(1);
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 lmnflag=-1;
00408               int nsteps=5;
00409               if ( (lmnflag=option_match(ad_comm::argc,ad_comm::argv,
00410                  "-lmn2",nopt))>-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       gradient_structure::set_CMPDIF_BUFFER_SIZE(lssz);
00555     }
00556 
00557     nopt=get_option_number("-gbs",
00558       "-gbs option needs positive integer -- ignored",lssz);
00559     if (nopt>-1 && lssz>0) {
00560       gradient_structure::set_GRADSTACK_BUFFER_SIZE
00561         (lssz/sizeof(grad_stack_entry));
00562     }
00563 
00564     if (!sz)
00565     {
00566       pgs = new gradient_structure;
00567     }
00568     else
00569     {
00570       pgs = new gradient_structure(sz);
00571     }
00572   }
00573 
00577 function_minimizer::~function_minimizer()
00578 {
00579   if (multinomial_weights)
00580   {
00581     delete multinomial_weights;
00582     multinomial_weights = 0;
00583   }
00584   if (lapprox)
00585   {
00586     delete lapprox;
00587     lapprox = 0;
00588   }
00589   delete pgs;
00590   pgs = 0;
00591   if (negdirections)
00592   {
00593     delete negdirections;
00594     negdirections = 0;
00595   }
00596 }
00597 
00598 void function_minimizer::set_initial_simplex(const dmatrix& _p,
00599   const dvector& _y, int nvar, const dvector& x, double delta)
00600   {
00601     dvector& y=(dvector&) _y;
00602     dmatrix& p=(dmatrix&) _p;
00603     int i;
00604     p(1)=x;
00605     for (i=2;i<=nvar+1;i++)
00606     {
00607       p(i)=x;
00608       p(i,i-1)+=delta;
00609     }
00610     dvector xx(1,nvar);
00611     double vf=0;
00612     for (i=1;i<=nvar+1;i++)
00613     {
00614       xx=p(i);
00615       vf=value(initial_params::reset(dvar_vector(xx)));
00616       *objective_function_value::pobjfun=0.0;
00617       userfunction();
00618       vf+=value(*objective_function_value::pobjfun);
00619       y(i)=vf;
00620     }
00621   }
00622 
00623 int get_option_number(const char * option_name,const char * error_message,
00624   int& option_value)
00625 {
00626   int on1;
00627   int nopt = 0;
00628   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00629   {
00630     if (!nopt)
00631     {
00632       if (ad_printf)
00633         (*ad_printf)("%s\n",error_message);
00634       else
00635         cerr << error_message << endl;
00636       on1=-1;
00637     }
00638     else
00639     {
00640       option_value=atoi(ad_comm::argv[on1+1]);
00641     }
00642   }
00643   return on1;
00644 }
00645 
00646 int get_option_number(const char * option_name,const char * error_message,
00647 #ifdef __BORLANDC__
00648   long int& option_value)
00649 #else
00650   long long int& option_value)
00651 #endif
00652 {
00653   int on1;
00654   int nopt = 0;
00655   if ( (on1=option_match(ad_comm::argc,ad_comm::argv,option_name,nopt))>-1)
00656   {
00657     if (!nopt)
00658     {
00659       if (ad_printf)
00660         (*ad_printf)("%s\n",error_message);
00661       else
00662         cerr << error_message << endl;
00663       on1=-1;
00664     }
00665     else
00666     {
00667 #if defined(__BORLANDC__) || defined(_MSC_VER)
00668       option_value=atol(ad_comm::argv[on1+1]);
00669 #else
00670       option_value=atoll(ad_comm::argv[on1+1]);
00671 #endif
00672     }
00673   }
00674   return on1;
00675 }
00676 
00677 
00678 void function_minimizer::other_separable_stuff_begin(void)
00679 {
00680   if (lapprox)
00681   {
00682     lapprox->separable_calls_counter++;
00683     /*
00684     lapprox->separable_call_level++;
00685     //lapprox->build_up_nested_shape();
00686     lapprox->nested_separable_calls_counter
00687       (lapprox->separable_call_level)++;
00688     //clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00689     lapprox->nested_tree_position(lapprox->separable_call_level)++;
00690     */
00691   }
00692 }
00693 
00694 void function_minimizer::other_separable_stuff_end(void)
00695 {
00696   /*
00697   if (lapprox)
00698   {
00699     lapprox->build_up_nested_shape();
00700     clean(lapprox->nested_tree_position,lapprox->separable_call_level);
00701     lapprox->separable_call_level--;
00702   }
00703   */
00704 }
00705 
00706 
00707 void function_minimizer::begin_gauss_hermite_stuff(void)
00708 {
00709   int nsc=lapprox->separable_calls_counter;
00710   int is=0;
00711   if (lapprox->gh->mi==0)
00712   {
00713     is=lapprox->gh->is;
00714   }
00715   else
00716   {
00717     is=lapprox->gh->mi->get_offset()+1;
00718   }
00719   lapprox->gh->gauss_hermite_values(nsc,is)=
00720     *objective_function_value::pobjfun;
00721 }
00722 
00723 void function_minimizer::start_get_importance_sampling_comnponent(void)
00724 {
00725   int nsc=lapprox->separable_calls_counter;
00726   int isc=lapprox->importance_sampling_counter;
00727   (*lapprox->importance_sampling_components)(nsc,isc)=
00728      *objective_function_value::pobjfun;
00729 }
00730 
00731 void function_minimizer::end_get_importance_sampling_comnponent(void)
00732 {
00733   int nsc=lapprox->separable_calls_counter;
00734   int is=lapprox->importance_sampling_counter;
00735   if (lapprox->saddlepointflag==2)
00736   {
00737     (*lapprox->importance_sampling_components)(nsc,is)=
00738        (-1)* *objective_function_value::pobjfun-
00739        (*lapprox->importance_sampling_components)(nsc,is);
00740   }
00741   else
00742   {
00743     (*lapprox->importance_sampling_components)(nsc,is)=
00744        *objective_function_value::pobjfun-
00745        (*lapprox->importance_sampling_components)(nsc,is);
00746   }
00747 }
00748 
00749 void function_minimizer::begin_funnel_stuff(void)
00750 {
00751   lapprox->separable_calls_counter++;
00752   if (lapprox->hesstype==2)
00753   {
00754     if (lapprox->in_gauss_hermite_phase)
00755     {
00756        begin_gauss_hermite_stuff();
00757     }
00758     else if (lapprox->num_importance_samples &&
00759       lapprox->importance_sampling_components)
00760     {
00761       if (lapprox->block_diagonal_flag==2)
00762       {
00763         start_get_importance_sampling_comnponent();
00764       }
00765     }
00766   }
00767 }
00768 
00769 void function_minimizer::get_function_difference(void)
00770 {
00771   int nsc=lapprox->separable_calls_counter;
00772   (*(lapprox->separable_function_difference))(nsc)=
00773     value(*objective_function_value::pobjfun);
00774     value(*objective_function_value::pobjfun)=0.0;
00775 }
00776 void function_minimizer::end_df1b2_funnel_stuff(void)
00777 {
00778   if (lapprox->in_gauss_hermite_phase)
00779   {
00780     end_gauss_hermite_stuff();
00781   }
00782   else
00783   {
00784     if (lapprox->hesstype==2)
00785     {
00786       if (lapprox->num_importance_samples &&
00787         lapprox->importance_sampling_components)
00788       {
00789         if (lapprox->block_diagonal_flag==2)
00790         {
00791           end_get_importance_sampling_comnponent();
00792         }
00793       }
00794       if (!lapprox->no_function_component_flag)
00795       {
00796         if (lapprox->saddlepointflag!=2)
00797         {
00798           get_function_difference();
00799         }
00800         else if (inner_opt()!=0)
00801         {
00802           get_function_difference();
00803         }
00804       }
00805     }
00806   }
00807 }
00808 
00809 void function_minimizer::end_gauss_hermite_stuff(void)
00810 {
00811   int nsc=lapprox->separable_calls_counter;
00812   int is=0;
00813   if (lapprox->gh->mi==0)
00814   {
00815     is=lapprox->gh->is;
00816   }
00817   else
00818   {
00819     is=lapprox->gh->mi->get_offset()+1;
00820   }
00821   lapprox->gh->gauss_hermite_values(nsc,is)=
00822     *objective_function_value::pobjfun-
00823     lapprox->gh->gauss_hermite_values(nsc,is);
00824 }
00825 
00826 void print_is_diagnostics(laplace_approximation_calculator *lapprox)
00827 {
00828   if (lapprox->is_diagnostics_flag)
00829   {
00830     if (lapprox->importance_sampling_values)
00831     {
00832       int mmin=lapprox->importance_sampling_values->indexmin();
00833       int mmax=lapprox->importance_sampling_values->indexmax();
00834       double mn= mean(*lapprox->importance_sampling_values);
00835       dmatrix tmp(1,2,mmin,mmax);
00836       tmp(2)=*lapprox->importance_sampling_values-mn;
00837       tmp(1).fill_seqadd(1,1);
00838       tmp=trans(sort(trans(tmp),2));
00839       ofstream ofs("is_diagnostics");
00840       ofs << "number of importance samples "
00841           << lapprox->num_importance_samples << endl;
00842       ofs << "importance_sampling_values" << endl;
00843       ofs << *lapprox->importance_sampling_values << endl<< endl;;
00844       ofs << "normalized importance_sampling_values" << endl;
00845       ofs << *lapprox->importance_sampling_values-mn << endl<< endl;;
00846       ofs << "sorted normalized importance_sampling_values" << endl;
00847       ofs << setw(9) << tmp << endl<< endl;;
00848       ofs << "epsilon(1).indexmax()  "
00849           << lapprox->epsilon(1).indexmax() << endl;
00850       ofs << lapprox->epsilon << endl;
00851       dmatrix plotstuff(1,2,mmin,mmax);
00852       plotstuff(1)=*lapprox->importance_sampling_weights;
00853       plotstuff(2)=*lapprox->importance_sampling_values-mn;
00854       ofs << " weight   value " << endl;
00855       ofs << setw(9) << sort(trans(plotstuff),2) << endl;
00856     }
00857   }
00858 }