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