ADMB Documentation  11.1x.2735
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2fnl2.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2fnl2.cpp 2615 2014-11-10 22:35:40Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fnl.h>
00012 #include <adrndeff.h>
00013 #ifndef OPT_LIB
00014   #include <cassert>
00015   #include <climits>
00016 #endif
00017 
00018 int pool_check_flag=0;
00019 
00020 extern int noboundepen_flag;
00021 
00026 void laplace_approximation_calculator::do_separable_stuff(void)
00027 {
00028   df1b2variable& ff=  *re_objective_function_value::pobjfun;
00029   if (block_diagonal_flag==0)
00030   {
00031     f1b2gradlist->reset();
00032     f1b2gradlist->list.initialize();
00033     f1b2gradlist->list2.initialize();
00034     f1b2gradlist->list3.initialize();
00035     f1b2gradlist->nlist.initialize();
00036     f1b2gradlist->nlist2.initialize();
00037     f1b2gradlist->nlist3.initialize();
00038     funnel_init_var::num_vars=0;
00039     funnel_init_var::num_active_parameters=0;
00040     funnel_init_var::num_inactive_vars=0;
00041     if (funnel_init_var::funnel_constraints_penalty)
00042     {
00043       delete funnel_init_var::funnel_constraints_penalty;
00044       funnel_init_var::funnel_constraints_penalty=0;
00045     }
00046     return;
00047   }
00048   if (funnel_init_var::funnel_constraints_penalty)
00049   {
00050     if (noboundepen_flag==0)
00051     {
00052       ff+=*funnel_init_var::funnel_constraints_penalty;
00053     }
00054     delete funnel_init_var::funnel_constraints_penalty;
00055     funnel_init_var::funnel_constraints_penalty=0;
00056   }
00057 
00058   // this should not be called a block diagopnal flag but it is for now.
00059   //if (pool_check_flag)
00060    // check_pool_depths();
00061   switch (block_diagonal_flag)
00062   {
00063   case 1:
00064     switch(hesstype)
00065     {
00066     case 2:
00067       do_separable_stuff_newton_raphson_block_diagonal(ff);
00068       break;
00069     case 3: //banded
00070     case 4: // full matrix
00071       do_separable_stuff_newton_raphson_banded(ff);
00072   //if (pool_check_flag)
00073    // check_pool_depths();
00074       break;
00075     default:
00076        cerr << "Illegal value for hesstype" << endl;
00077        ad_exit(1);
00078     }
00079     break;
00080   case 2:
00081     switch(hesstype)
00082     {
00083     case 2:
00084       ++separable_calls_counter;
00085       if (saddlepointflag==2)
00086       {
00087         do_separable_stuff_x_u_block_diagonal(ff);
00088       }
00089       else
00090       {
00091         do_separable_stuff_laplace_approximation_block_diagonal(ff);
00092       }
00093       break;
00094     case 3:
00095     case 4: // full matrix
00096       do_separable_stuff_laplace_approximation_banded(ff);
00097       break;
00098     default:
00099        cerr << "Illegal value for hesstype" << endl;
00100        ad_exit(1);
00101     }
00102     break;
00103   case 3:
00104     switch(hesstype)
00105     {
00106     case 2:
00107       do_separable_stuff_laplace_approximation_importance_sampling_adjoint(ff);
00108       break;
00109     case 3:
00110     case 4: // full matrix
00111       do_separable_stuff_laplace_approximation_banded_adjoint(ff);
00112       break;
00113     default:
00114        cerr << "Illegal value for hesstype" << endl;
00115        ad_exit(1);
00116     }
00117     break;
00118   case 5:   // get hessian type information
00119       do_separable_stuff_hessian_type_information();
00120       break;
00121   case 6:   // get hessian type information
00122       get_block_diagonal_hessian(ff);
00123       break;
00124   default:
00125     cerr << "illegal value for block_diagonal_flag = "
00126       << block_diagonal_flag << endl;
00127     ad_exit(1);
00128   }
00129   df1b2variable::pool=df1b2variable::adpool_vector[0];
00130   df1b2variable::nvar=df1b2variable::pool->nvar;
00131   if  (df1b2variable::pool->size<=0)
00132   {
00133     cerr << "this can't happen" << endl;
00134     ad_exit(1);
00135   }
00136 }
00137 
00142 void laplace_approximation_calculator::
00143   do_separable_stuff_newton_raphson_block_diagonal(df1b2variable& ff)
00144 {
00145   //***********************************************************
00146   //***********************************************************
00147   set_dependent_variable(ff);
00148   df1b2_gradlist::set_no_derivatives();
00149   df1b2variable::passnumber=1;
00150   df1b2_gradcalc1();
00151 
00152   init_df1b2vector & locy= *funnel_init_var::py;
00153   imatrix& list=*funnel_init_var::plist;
00154   int num_local_re=0;
00155   int num_fixed_effects=0;
00156 
00157   //cout << list << endl;
00158 #ifndef OPT_LIB
00159   assert(funnel_init_var::num_active_parameters <= INT_MAX);
00160 #endif
00161   ivector lre_index(1, (int)funnel_init_var::num_active_parameters);
00162   ivector lfe_index(1, (int)funnel_init_var::num_active_parameters);
00163 
00164   for (int i=1;i<=(int)funnel_init_var::num_active_parameters;i++)
00165   {
00166     if (list(i,1)>xsize)
00167     {
00168       lre_index(++num_local_re)=i;
00169     }
00170     else if (list(i,1)>0)
00171     {
00172       lfe_index(++num_fixed_effects)=i;
00173     }
00174   }
00175 
00176   if (num_local_re > 0)
00177   {
00178     //cout << " num_local_re = " << num_local_re << endl;
00179     //cout << " num_fixed_effects = " << num_fixed_effects << endl;
00180     //cout << " list = " << endl << list << endl;
00181     //cout << " lre_index= " << endl <<  lre_index << endl;
00182     //cout << " lfe_index= " << endl <<  lfe_index << endl;
00183     //cout << "the number of local res is " << num_local_re << endl;
00184     dmatrix local_Hess(1,num_local_re,1,num_local_re);
00185     dvector local_grad(1,num_local_re);
00186     //dmatrix local_Dux(1,num_local_re,1,
00187      // funnel_init_var::num_active_parameters-num_local_re);
00188     //dmatrix Hess1(1,funnel_init_var::num_vars,1,funnel_init_var::num_vars);
00189     local_Hess.initialize();
00190 
00191     for (int i=1;i<=num_local_re;i++)
00192     {
00193       int lrei=lre_index(i);
00194       for (int j=1;j<=num_local_re;j++)
00195       {
00196         int lrej=lre_index(j);
00197         int i2=list(lrei,2);
00198         int j2=list(lrej,2);
00199         //Hess(i1-xsize,j1-xsize)+=locy(i2).u_bar[j2-1];
00200         local_Hess(i,j)+=locy(i2).u_bar[j2-1];
00201       }
00202     }
00203      // i<=funnel_init_var::num_vars;i++)
00204     for (int i=1;i<=num_local_re;i++)
00205     {
00206       int lrei=lre_index(i);
00207       //int i1=list(lrei,1);
00208       int i2=list(lrei,2);
00209       //grad(i1-xsize)= re_objective_function_value::pobjfun->u_dot[i2-1];
00210       //grad(i1-xsize)= ff.u_dot[i2-1];
00211       local_grad(i)= ff.u_dot[i2-1];
00212     }
00213 
00214     have_bounded_random_effects=0;
00215     if (have_bounded_random_effects)
00216     {
00217       for (int i=1;i<=num_local_re;i++)
00218       {
00219         int lrei=lre_index(i);
00220         int i1=list(lrei,1);
00221         for (int j=1;j<=num_local_re;j++)
00222         {
00223           int lrej=lre_index(j);
00224           int j1=list(lrej,1);
00225           local_Hess(i,j)*=scale(i1-xsize)*scale(j1-xsize);
00226         }
00227       }
00228 
00229       for (int i=1;i<=num_local_re;i++)
00230       {
00231         int lrei=lre_index(i);
00232         int i1=list(lrei,1);
00233         local_Hess(i,i)+=local_grad(i)*curv(i1-xsize);
00234       }
00235 
00236       for (int i=1;i<=num_local_re;i++)
00237       {
00238         int lrei=lre_index(i);
00239         int i1=list(lrei,1);
00240         local_grad(i)*=scale(i1-xsize);
00241       }
00242     }
00243 
00244     double mg=max(fabs(local_grad));
00245     if (max_separable_g< mg) max_separable_g=mg;
00246     dvector local_step=-solve(local_Hess,local_grad);
00247 
00248     for (int i=1;i<=num_local_re;i++)
00249     {
00250       int lrei=lre_index(i);
00251       int i1=list(lrei,1);
00252       //int i2=list(lrei,2);
00253       step(i1-xsize)=local_step(i);
00254     }
00255   }
00256 
00257   f1b2gradlist->reset();
00258   f1b2gradlist->list.initialize();
00259   f1b2gradlist->list2.initialize();
00260   f1b2gradlist->list3.initialize();
00261   f1b2gradlist->nlist.initialize();
00262   f1b2gradlist->nlist2.initialize();
00263   f1b2gradlist->nlist3.initialize();
00264   funnel_init_var::num_vars=0;
00265   funnel_init_var::num_active_parameters=0;
00266   funnel_init_var::num_inactive_vars=0;
00267 }