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