ADMB Documentation  11.1.1928
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp10.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp10.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if defined(USE_LAPLACE)
00012 #  include <admodel.h>
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015 
00020 void report_calling_set(laplace_approximation_calculator *lapprox)
00021 {
00022   ofstream ofs("callset.rpt");
00023 
00024   imatrix& callset=(*lapprox->calling_set);
00025 
00026   ofs << "Total num_separable calls " <<  callset(0,0)-1 << endl;
00027 
00028   for (int i=1;i<=callset.indexmax();i++)
00029   {
00030     ofs << "Variable " << i << " num calls = " << callset(i)(0) << endl;
00031     ofs << callset(i)(1,callset(i).indexmax())<< endl;
00032   }
00033 }
00034 
00039 void check_order(ivector& v)
00040 {
00041   int mmin=v.indexmin();
00042   int mmax=v.indexmax();
00043   int bad=0;
00044   for (int i=mmin;i<=mmax-1;i++)
00045   {
00046     if (v(i+1)<v(i))
00047     {
00048       bad=1;
00049       break;
00050     }
00051     if (bad)
00052     {
00053       v=sort(v);
00054     }
00055   }
00056 }
00057 
00062 int common(ivector& v,ivector& w)
00063 {
00064   check_order(v);
00065   check_order(w);
00066   //int vmin=v.indexmin();
00067   int wmin=w.indexmin();
00068   int vmax=v.indexmax();
00069   int wmax=w.indexmax();
00070   int common_flag=0;
00071   int i=wmin; int j=wmin;
00072   do
00073   {
00074     if (v(i)==w(j))
00075     {
00076       common_flag=1;
00077       break;
00078     }
00079     else if (v(i)>w(j))
00080     {
00081       if (j<wmax)
00082         j++;
00083       else
00084         break;
00085     }
00086     else if (v(i)<w(j))
00087     {
00088       if (i<vmax)
00089         i++;
00090       else
00091         break;
00092     }
00093   }
00094   while(1);
00095   return common_flag;
00096 }
00097 
00102 void laplace_approximation_calculator::
00103   check_hessian_type2(function_minimizer * pfmin)
00104 {
00105   int ip = 0;
00106   if (quadratic_prior::get_num_quadratic_prior()>0)
00107   {
00108     hesstype=4;
00109     if (allocated(Hess))
00110     {
00111       if (Hess.indexmax()!=usize)
00112       {
00113         Hess.deallocate();
00114         Hess.allocate(1,usize,1,usize);
00115       }
00116     }
00117     else
00118     {
00119        Hess.allocate(1,usize,1,usize);
00120     }
00121     if (allocated(Hessadjoint))
00122     {
00123       if (Hessadjoint.indexmax()!=usize)
00124       {
00125         Hessadjoint.deallocate();
00126         Hessadjoint.allocate(1,usize,1,usize);
00127       }
00128     }
00129     else
00130     {
00131        Hessadjoint.allocate(1,usize,1,usize);
00132     }
00133     return;
00134   }
00135   else
00136   {
00137     int nv=initial_df1b2params::set_index();
00138     if (allocated(used_flags))
00139     {
00140       if (used_flags.indexmax() != nv)
00141       {
00142         used_flags.safe_deallocate();
00143       }
00144     }
00145     if (!allocated(used_flags))
00146     {
00147       used_flags.safe_allocate(1,nv);
00148     }
00149 
00150     //for (ip=1;ip<=num_der_blocks;ip++)
00151     {
00152       used_flags.initialize();
00153       // do we need to reallocate memory for df1b2variables?
00154       check_for_need_to_reallocate(ip);
00155       df1b2_gradlist::set_no_derivatives();
00156       //cout << re_objective_function_value::pobjfun << endl;
00157       //cout << re_objective_function_value::pobjfun->ptr << endl;
00158       (*re_objective_function_value::pobjfun)=0;
00159       df1b2variable pen=0.0;
00160       df1b2variable zz=0.0;
00161 
00162       initial_df1b2params::reset(y,pen);
00163       // call function to do block diagonal newton-raphson
00164       // the step vector from the newton-raphson is in the vector step
00165       df1b2_gradlist::set_no_derivatives();
00166 
00167       funnel_init_var::lapprox=this;
00168       block_diagonal_flag=5;
00169 
00170       quadratic_prior::in_qp_calculations=1;
00171       pfmin->pre_user_function();
00172       quadratic_prior::in_qp_calculations=0;
00173 
00174       int non_block_diagonal=0;
00175       for (int i=xsize+1;i<=xsize+usize;i++)
00176       {
00177         if (used_flags(i)>1)
00178         {
00179           non_block_diagonal=1;
00180           break;
00181         }
00182       }
00183       if (non_block_diagonal)
00184       {
00185         if (bw< usize/2)
00186         {
00187           hesstype=3;  //banded
00188           if (bHess)
00189           {
00190             if (bHess->bandwidth() !=bw)
00191             {
00192               delete bHess;
00193               bHess = new banded_symmetric_dmatrix(1,usize,bw);
00194               if (bHess==0)
00195               {
00196                 cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00197                 ad_exit(1);
00198               }
00199             }
00200           }
00201           else
00202           {
00203             bHess = new banded_symmetric_dmatrix(1,usize,bw);
00204             if (bHess==0)
00205             {
00206               cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00207               ad_exit(1);
00208             }
00209           }
00210           if (bHessadjoint)
00211           {
00212             if (bHessadjoint->bandwidth() !=bw)
00213             {
00214               delete bHessadjoint;
00215               bHessadjoint = new banded_symmetric_dmatrix(1,usize,bw);
00216               if (bHessadjoint==0)
00217               {
00218                 cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00219                 ad_exit(1);
00220               }
00221             }
00222           }
00223           else
00224           {
00225             bHessadjoint = new banded_symmetric_dmatrix(1,usize,bw);
00226             if (bHessadjoint==0)
00227             {
00228               cerr << "Error allocating banded_symmetric_dmatrix" << endl;
00229               ad_exit(1);
00230             }
00231           }
00232         }
00233         else
00234         {
00235           //check_sparse_matrix_structure();
00236           hesstype=4;  // band is so wide so use full matrix
00237           if (bHess)
00238           {
00239             delete bHess;
00240             bHess=0;
00241           }
00242 
00243           if (bHessadjoint)
00244           {
00245             delete bHessadjoint;
00246             bHessadjoint=0;
00247           }
00248 
00249           if (allocated(Hess))
00250           {
00251             if (Hess.indexmax() != usize)
00252             {
00253               Hess.deallocate();
00254               Hess.allocate(1,usize,1,usize);
00255             }
00256           }
00257           else
00258           {
00259             Hess.allocate(1,usize,1,usize);
00260           }
00261           if (allocated(Hessadjoint))
00262           {
00263             if (Hessadjoint.indexmax() != usize)
00264             {
00265               Hessadjoint.deallocate();
00266               Hessadjoint.allocate(1,usize,1,usize);
00267             }
00268           }
00269           else
00270           {
00271             Hessadjoint.allocate(1,usize,1,usize);
00272           }
00273         }
00274       }
00275       else
00276       {
00277         hesstype=2;
00278       }
00279       if (hesstype==2 && num_importance_samples>0)
00280       {
00281         if (importance_sampling_components)
00282         {
00283           delete importance_sampling_components;
00284           importance_sampling_components=0;
00285         }
00286         importance_sampling_components=
00287           new dvar_matrix(1,pmin->lapprox->num_separable_calls,
00288             1,num_importance_samples);
00289       }
00290       // check for containg block diagonal structure
00291       used_flags(1);
00292       if (calling_set)
00293       {
00294         delete calling_set;
00295       }
00296       int mmin=used_flags.indexmin()-1;
00297       int mmax=used_flags.indexmax();
00298       ivector tmp(mmin,mmax);
00299       tmp(mmin)=0;
00300       tmp(mmin+1,mmax)=used_flags;
00301       {
00302         calling_set=new imatrix(mmin,mmax,0,tmp);
00303         calling_set->initialize();
00304         (*calling_set)(0,0)=1;
00305       }
00306       used_flags.initialize();
00307       quadratic_prior::in_qp_calculations=1;
00308       pfmin->pre_user_function();
00309       quadratic_prior::in_qp_calculations=0;
00310       report_calling_set(this);
00311 
00312       if (hesstype==2 && (num_importance_samples>0 || use_gauss_hermite>0))
00313       {
00314         const ivector & itmp=(*num_local_re_array)(1,num_separable_calls);
00315         const ivector & itmpf=(*num_local_fixed_array)(1,num_separable_calls);
00316 
00317         // ****************************************************
00318         // ****************************************************
00319         if (use_gauss_hermite>0)
00320         {
00321           if (gh)
00322           {
00323             delete gh;
00324             gh=0;
00325           }
00326           gh=new gauss_hermite_stuff(this,use_gauss_hermite,
00327             num_separable_calls,itmp);
00328         }
00329 
00330         if (block_diagonal_vch)
00331         {
00332           delete block_diagonal_vch;
00333           block_diagonal_vch=0;
00334         }
00335 
00336         block_diagonal_vch = new dvar3_array(1,num_separable_calls,
00337           1,itmp,1,itmp);
00338         if (block_diagonal_ch)
00339         {
00340           delete block_diagonal_ch;
00341           block_diagonal_ch=0;
00342         }
00343         block_diagonal_ch = new d3_array(1,num_separable_calls,
00344           1,itmp,1,itmp);
00345         if (block_diagonal_hessian)
00346         {
00347           delete block_diagonal_hessian;
00348           block_diagonal_hessian=0;
00349         }
00350         block_diagonal_hessian = new d3_array(1,num_separable_calls,
00351           1,itmp,1,itmp);
00352         if (block_diagonal_hessian ==0)
00353         {
00354           cerr << "error_allocating d3_array" << endl;
00355           ad_exit(1);
00356         }
00357         block_diagonal_re_list = new imatrix(1,num_separable_calls,
00358           1,itmp);
00359         if (block_diagonal_re_list ==0)
00360         {
00361           cerr << "error_allocating imatrix" << endl;
00362           ad_exit(1);
00363         }
00364         block_diagonal_fe_list = new imatrix(1,num_separable_calls,
00365           1,itmpf);
00366         if (block_diagonal_fe_list ==0)
00367         {
00368           cerr << "error_allocating imatrix" << endl;
00369           ad_exit(1);
00370         }
00371         // ****************************************************
00372         if (block_diagonal_Dux)
00373         {
00374           delete block_diagonal_Dux;
00375           block_diagonal_Dux=0;
00376         }
00377         block_diagonal_Dux = new d3_array(1,num_separable_calls,
00378           1,itmp,1,itmpf);
00379         if (block_diagonal_Dux ==0)
00380         {
00381           cerr << "error_allocating d3_array" << endl;
00382           ad_exit(1);
00383         }
00384 
00385         // ****************************************************
00386         // ****************************************************
00387         if (block_diagonal_vhessian)
00388         {
00389           delete block_diagonal_vhessian;
00390           block_diagonal_vhessian=0;
00391         }
00392         block_diagonal_vhessian = new dvar3_array(1,num_separable_calls,
00393           1,itmp,1,itmp);
00394         if (block_diagonal_vhessian ==0)
00395         {
00396           cerr << "error_allocating d3_array" << endl;
00397           ad_exit(1);
00398         }
00399 
00400         if (block_diagonal_vhessianadjoint)
00401         {
00402           delete block_diagonal_vhessianadjoint;
00403           block_diagonal_vhessianadjoint=0;
00404         }
00405         block_diagonal_vhessianadjoint = new d3_array(1,num_separable_calls,
00406           1,itmp,1,itmp);
00407         if (block_diagonal_vhessianadjoint ==0)
00408         {
00409           cerr << "error_allocating d3_array" << endl;
00410           ad_exit(1);
00411         }
00412       }
00413       funnel_init_var::lapprox=0;
00414       block_diagonal_flag=0;
00415       pen.deallocate();
00416     }
00417   }
00418 }
00419 #endif // if defined(USE_LAPLACE)