ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp7.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp7.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 <df1b2fun.h>
00013 
00018 void laplace_approximation_calculator::
00019    do_separable_stuff_laplace_approximation_banded_adjoint
00020    (const df1b2variable& ff)
00021 {
00022   //***********************************************************
00023   //***********************************************************
00024   set_dependent_variable(ff);
00025   df1b2_gradlist::set_no_derivatives();
00026   df1b2variable::passnumber=1;
00027   df1b2_gradcalc1();
00028 
00029   init_df1b2vector & locy= *funnel_init_var::py;
00030   imatrix& list=*funnel_init_var::plist;
00031   int num_local_re=0;
00032   int num_local_fe=0;
00033 
00034   int i;
00035   ivector lre_index(1,funnel_init_var::num_active_parameters);
00036   ivector lfe_index(1,funnel_init_var::num_active_parameters);
00037 
00038   for (i=1;i<=funnel_init_var::num_active_parameters;i++)
00039   {
00040     if (list(i,1)>xsize)
00041     {
00042       lre_index(++num_local_re)=i;
00043     }
00044     else if (list(i,1)>0)
00045     {
00046       lfe_index(++num_local_fe)=i;
00047     }
00048   }
00049 
00050   if (num_local_re > 0)
00051   {
00052     int j;
00053     if (hesstype==3)
00054     {
00055       for (i=1;i<=num_local_re;i++)
00056       {
00057         int lrei=lre_index(i);
00058         for (j=1;j<=num_local_re;j++)
00059         {
00060           int lrej=lre_index(j);
00061           int i1=list(lrei,1)-xsize;
00062           int i2=list(lrei,2);
00063           int j1=list(lrej,1)-xsize;
00064           int j2=list(lrej,2);
00065           if (i1>=j1) (*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00066         }
00067       }
00068     }
00069     else if (hesstype==4)
00070     {
00071       if (sparse_hessian_flag==0)
00072       {
00073         for (i=1;i<=num_local_re;i++)
00074         {
00075           int lrei=lre_index(i);
00076           for (j=1;j<=num_local_re;j++)
00077           {
00078             int lrej=lre_index(j);
00079             int i1=list(lrei,1)-xsize;
00080             int i2=list(lrei,2);
00081             int j1=list(lrej,1)-xsize;
00082             int j2=list(lrej,2);
00083             Hess(i1,j1)+=locy(i2).u_bar[j2-1];
00084           }
00085         }
00086       }
00087       else
00088       {
00089         for (i=1;i<=num_local_re;i++)
00090         {
00091           int lrei=lre_index(i);
00092           for (j=1;j<=num_local_re;j++)
00093           {
00094             int lrej=lre_index(j);
00095             int i1=list(lrei,1)-xsize;
00096             int i2=list(lrei,2);
00097             int j1=list(lrej,1)-xsize;
00098             int j2=list(lrej,2);
00099 
00100             if (i1<=j1)
00101             {
00102               sparse_count++;
00103               (*sparse_triplet2)((*sparse_iterator)(sparse_count))
00104                 +=locy(i2).u_bar[j2-1];
00105             }
00106           }
00107         }
00108       }
00109     }
00110   }
00111 
00112   // Now the adjoint code
00113 
00114   if (num_local_re > 0)
00115   {
00116     int j;
00117     if (hesstype==3)
00118     {
00119       for (i=1;i<=num_local_re;i++)
00120       {
00121         int lrei=lre_index(i);
00122         for (j=1;j<=num_local_re;j++)
00123         {
00124           int lrej=lre_index(j);
00125           int i1=list(lrei,1)-xsize;
00126           int i2=list(lrei,2);
00127           int j1=list(lrej,1)-xsize;
00128           int j2=list(lrej,2);
00129           if (i1>=j1)
00130           {
00131             //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00132             locy(i2).get_u_bar_tilde()[j2-1]=(*bHessadjoint)(i1,j1);
00133           }
00134         }
00135       }
00136     }
00137     else if (hesstype==4)
00138     {
00139       if (pmin->lapprox->sparse_hessian_flag==0)
00140       {
00141         for (i=1;i<=num_local_re;i++)
00142         {
00143           int lrei=lre_index(i);
00144           for (j=1;j<=num_local_re;j++)
00145           {
00146             int lrej=lre_index(j);
00147             int i1=list(lrei,1)-xsize;
00148             int i2=list(lrei,2);
00149             int j1=list(lrej,1)-xsize;
00150             int j2=list(lrej,2);
00151             {
00152               //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00153               locy(i2).get_u_bar_tilde()[j2-1]=Hessadjoint(i1,j1);
00154             }
00155           }
00156         }
00157       }
00158       else
00159       {
00160         dcompressed_triplet * vsparse_triplet_adjoint =
00161           pmin->lapprox->vsparse_triplet_adjoint;
00162         for (i=1;i<=num_local_re;i++)
00163         {
00164           int lrei=lre_index(i);
00165           for (j=1;j<=num_local_re;j++)
00166           {
00167             int lrej=lre_index(j);
00168             int i1=list(lrei,1)-xsize;
00169             int i2=list(lrei,2);
00170             int j1=list(lrej,1)-xsize;
00171             int j2=list(lrej,2);
00172             {
00173               //(*bHess)(i1,j1)+=locy(i2).u_bar[j2-1];
00174               //locy(i2).get_u_bar_tilde()[j2-1]=Hessadjoint(i1,j1);
00175               if (i1<=j1)
00176               {
00177                 //(*sparse_triplet2)((*sparse_iterator)(sparse_count_adjoint))
00178                 //  +=locy(i2).u_bar[j2-1];
00179                 sparse_count_adjoint++;
00180                 locy(i2).get_u_bar_tilde()[j2-1]=
00181                   (*vsparse_triplet_adjoint)
00182                     ((*sparse_iterator)(sparse_count_adjoint));
00183               }
00184             }
00185           }
00186         }
00187       }
00188     }
00189   }
00190 
00191   df1b2variable::passnumber=2;
00192   df1b2_gradcalc1();
00193 
00194   df1b2variable::passnumber=3;
00195   df1b2_gradcalc1();
00196 
00197   f1b2gradlist->reset();
00198   f1b2gradlist->list.initialize();
00199   f1b2gradlist->list2.initialize();
00200   f1b2gradlist->list3.initialize();
00201   f1b2gradlist->nlist.initialize();
00202   f1b2gradlist->nlist2.initialize();
00203   f1b2gradlist->nlist3.initialize();
00204 
00205 
00206   for (i=1;i<=num_local_fe;i++)
00207   {
00208     int lfei=lfe_index(i);
00209     int i1=list(lfei,1);
00210     //if (lfei - list(lfei,2))
00211      // cout << lfei << " " <<  list(lfei,2) << endl;
00212     local_dtemp(i1)+=*locy(lfei).get_u_tilde();
00213 #if !defined(OPT_LIB)
00214     int mmin=xadjoint.indexmin();
00215     int mmax=xadjoint.indexmax();
00216     if (i1<mmin || i1> mmax)
00217     {
00218        cerr << "this can't happen" << endl;
00219       ad_exit(1);
00220     }
00221 #endif
00222   }
00223 
00224   for (i=1;i<=num_local_re;i++)
00225   {
00226     int i1=list(lre_index(i),1)-xsize;
00227     int i2=list(lre_index(i),2);
00228     uadjoint(i1)+=*locy(i2).get_u_tilde();
00229   }
00230 
00231   f1b2gradlist->reset();
00232   f1b2gradlist->list.initialize();
00233   f1b2gradlist->list2.initialize();
00234   f1b2gradlist->list3.initialize();
00235   f1b2gradlist->nlist.initialize();
00236   f1b2gradlist->nlist2.initialize();
00237   f1b2gradlist->nlist3.initialize();
00238   funnel_init_var::num_vars=0;
00239   funnel_init_var::num_active_parameters=0;
00240   funnel_init_var::num_inactive_vars=0;
00241 }
00242 #endif // if defined(USE_LAPLACE)