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