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