ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f27.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f27.cpp 2601 2014-11-08 22:34:58Z 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 #endif
00016 
00017 void ncount_checker(int ncount,int ncount_check);
00018 
00019 void ad_read_pass1x(void);
00020 void ad_read_pass2x(void);
00021 void read_pass1_3x(void);
00022 // should inline this
00023 void read_pass1_1x(void);
00024 void read_pass1_2x(void);
00025 
00030 int df1b2_gradlist::write_pass1x(const df1b2variable * _px,
00031   df1b2variable * pz, df1b2function1 * pf)
00032 {
00033   ADUNCONST(df1b2variable*,px)
00034   ncount++;
00035 #if defined(CHECK_COUNT)
00036   if (ncount >= ncount_check)
00037     ncount_checker(ncount,ncount_check);
00038 #endif
00039   unsigned int nvar=df1b2variable::nvar;
00040 
00041   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)+sizeof(char*)
00042     +sizeof(double)+nvar*sizeof(double);
00043 // string identifier debug stuff
00044 #if defined(SAFE_ALL)
00045   char ids[]="CX";
00046   int slen=strlen(ids);
00047   total_bytes+=slen;
00048 #endif
00049   list.check_buffer_size(total_bytes);
00050   void * tmpptr=list.bptr;
00051 #if defined(SAFE_ALL)
00052   memcpy(list,ids,slen);
00053 #endif
00054 // end of string identifier debug stuff
00055 
00056   memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00057   memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00058   //list.bptr+=sizeof(df1b2_header);
00059   memcpy(list,&pf,sizeof(char*));
00060   //*(char**)(list.bptr)=(char*)pf;
00061   //list.bptr+=sizeof(char*);
00062   memcpy(list,px->get_u(),sizeof(double));
00063   //list.bptr+=sizeof(double);
00064   memcpy(list,px->get_u_dot(),nvar*sizeof(double));
00065   //list.bptr+=nvar*sizeof(double);
00066   // ***** write  record size
00067   nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00068   nlist.bptr->pf=(ADrfptr)(&ad_read_pass1x);
00069   ++nlist;
00070   return 0;
00071 }
00072 
00077 void ad_read_pass1x(void)
00078 {
00079   switch(df1b2variable::passnumber)
00080   {
00081   case 1:
00082     read_pass1_1x();
00083     break;
00084   case 2:
00085     read_pass1_2x();
00086     break;
00087   case 3:
00088     read_pass1_3x();
00089     break;
00090   default:
00091     cerr << "illegal value for df1b2variable::pass = "
00092          << df1b2variable::passnumber << endl;
00093     exit(1);
00094   }
00095 }
00096 
00101 void read_pass1_1x(void)
00102 {
00103   // We are going backword for bptr and forward for bptr2
00104   // the current entry+2 in bptr is the size of the record i.e
00105   // points to the next record
00106   unsigned int nvar=df1b2variable::nvar;
00107   fixed_smartlist & nlist=f1b2gradlist->nlist;
00108   test_smartlist& list=f1b2gradlist->list;
00109    // nlist-=sizeof(int);
00110   // get record size
00111   int num_bytes=nlist.bptr->numbytes;
00112   // backup the size of the record
00113   list-=num_bytes;
00114   //list.bptr-=num_bytes;
00115   list.saveposition(); // save pointer to beginning of record;
00116   // save the pointer to the beginning of the record
00117   double xu;
00118   double * xdot;
00119   //df1b2_header x,z;
00120   df1b2function1 * pf;
00121 
00122   // get info from tape1
00123   // get info from tape1
00124 #if defined(SAFE_ARRAYS)
00125   checkidentiferstring("CX",list);
00126 #endif
00127   df1b2_header * px=(df1b2_header *) list.bptr;
00128   list.bptr+=sizeof(df1b2_header);
00129   df1b2_header * pz=(df1b2_header *) list.bptr;
00130   list.bptr+=sizeof(df1b2_header);
00131   pf=*(df1b2function1 **) list.bptr;
00132   list.bptr+=sizeof(char*);
00133   memcpy(&xu,list.bptr,sizeof(double));
00134   list.bptr+=sizeof(double);
00135   xdot=(double*)list.bptr;
00136   list.restoreposition(); // save pointer to beginning of record;
00137 
00138   // Do first reverse paSS calculations
00139   // ****************************************************************
00140   // turn this off if no third derivatives are calculated
00141   // if (!no_third_derivatives)
00142   // {
00143   // save for second reverse pass
00144   // save identifier 1
00145      fixed_smartlist2& nlist2=f1b2gradlist->nlist2;
00146      test_smartlist& list2=f1b2gradlist->list2;
00147 
00148 
00149   size_t total_bytes=2*nvar*sizeof(double);
00150 // string identifier debug stuff
00151 #if defined(SAFE_ALL)
00152   char ids[]="DU";
00153   int slen=strlen(ids);
00154   total_bytes+=slen;
00155 #endif
00156 
00157   list2.check_buffer_size(total_bytes);
00158   void * tmpptr2=list2.bptr;
00159 
00160 #if defined(SAFE_ALL)
00161   memcpy(list2,ids,slen);
00162 #endif
00163 
00164    memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00165    memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00166    *nlist2.bptr=adptr_diff(list2.bptr,tmpptr2);
00167    ++nlist2;
00168   // }
00169   //
00170   // ****************************************************************
00171 #if defined(PRINT_DERS)
00172  print_derivatives(pf->funname,(pf->df)(xu),(pf->df)(xu),(pf->d2f)(xu),
00173   (pf->d3f)(xu),1);
00174  print_derivatives(pz,"z");
00175  print_derivatives(px,"x");
00176 #endif
00177 
00178   double df=(pf->df)(xu);
00179   double d2f=(pf->d2f)(xu);
00180   //double d3f=(pf->d3f)(xu);
00181 
00182   for (size_t i=0;i<nvar;i++)
00183   {
00184     //px->u_bar[i]+=(pf->df)(xu)* pz->u_bar[i];
00185     px->u_bar[i]+=df * pz->u_bar[i];
00186   }
00187   for (size_t i=0;i<nvar;i++)
00188   {
00189     //px->u_bar[i]+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_bar[i];
00190     px->u_bar[i]+=d2f*xdot[i]*pz->u_dot_bar[i];
00191   }
00192   for (size_t i=0;i<nvar;i++)
00193   {
00194     //px->u_dot_bar[i]+=(pf->df)(xu)*pz->u_dot_bar[i];
00195     px->u_dot_bar[i]+=df*pz->u_dot_bar[i];
00196   }
00197 
00198   // !!!!!!!!!!!!!!!!!!!!!!
00199   for (size_t i=0;i<nvar;i++)
00200   {
00201     pz->u_bar[i]=0;
00202   }
00203   for (size_t i=0;i<nvar;i++)
00204   {
00205     pz->u_dot_bar[i]=0;
00206   }
00207 
00208 #if defined(PRINT_DERS)
00209  print_derivatives(pz,"z");
00210  print_derivatives(px,"x");
00211 #endif
00212 }
00213 
00218 void read_pass1_2x(void)
00219 {
00220   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00221   // We are going forward for bptr and backword for bptr2
00222   //
00223   // list 1
00224   //
00225   unsigned int nvar=df1b2variable::nvar;
00226   test_smartlist & list=f1b2gradlist->list;
00227 
00228   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)+sizeof(char*)
00229     +sizeof(double)+nvar*sizeof(double);
00230 // string identifier debug stuff
00231 #if defined(SAFE_ALL)
00232   char ids[]="CX";
00233   int slen=strlen(ids);
00234   total_bytes+=slen;
00235 #endif
00236   list.check_buffer_size(total_bytes);
00237 // end of string identifier debug stuff
00238 
00239   list.saveposition(); // save pointer to beginning of record;
00240   fixed_smartlist & nlist=f1b2gradlist->nlist;
00241    // nlist-=sizeof(int);
00242   // get record size
00243   int num_bytes=nlist.bptr->numbytes;
00244     // nlist+=nlist_record_size;
00245   //
00246   // list 2
00247   //
00248   test_smartlist & list2=f1b2gradlist->list2;
00249   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00250   // get record size
00251   int num_bytes2=*nlist2.bptr;
00252   --nlist2;
00253   // backup the size of the record
00254   list2-=num_bytes2;
00255   list2.saveposition();
00256   // save the pointer to the beginning of the record
00257   // bptr and bptr2 now both point to the beginning of their records
00258 
00259   double xu;
00260   double * xdot;
00261   //df1b2_header x,z;
00262   df1b2function1 * pf;
00263 
00264   // get info from tape1
00265   // get info from tape1
00266 #if defined(SAFE_ARRAYS)
00267   checkidentiferstring("CX",list);
00268   checkidentiferstring("DU",list2);
00269 #endif
00270   df1b2_header * px=(df1b2_header *) list.bptr;
00271   list.bptr+=sizeof(df1b2_header);
00272   df1b2_header * pz=(df1b2_header *) list.bptr;
00273   list.bptr+=sizeof(df1b2_header);
00274   pf=*(df1b2function1 **) list.bptr;
00275   list.bptr+=sizeof(char*);
00276   memcpy(&xu,list.bptr,sizeof(double));
00277   list.bptr+=sizeof(double);
00278   xdot=(double*)list.bptr;
00279   list.restoreposition(num_bytes); // save pointer to beginning of record;
00280 
00281   double * zbar;
00282   double * zdotbar;
00283 
00284   zbar=(double*)list2.bptr;
00285   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00286 
00287   double * x_tilde=px->get_u_tilde();
00288   double * x_dot_tilde=px->get_u_dot_tilde();
00289   double * x_bar_tilde=px->get_u_bar_tilde();
00290   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00291   double * z_bar_tilde=pz->get_u_bar_tilde();
00292   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00293 #if defined(PRINT_DERS)
00294  print_derivatives(pf->funname,(pf->df)(xu),(pf->df)(xu),(pf->d2f)(xu),
00295   (pf->d3f)(xu),1);
00296  print_derivatives(pz,"z");
00297  print_derivatives(px,"x");
00298 #endif
00299   // Do second "reverse-reverse" pass calculations
00300 
00301   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00302   for (size_t i=0;i<nvar;i++)
00303   {
00304     z_bar_tilde[i]=0;
00305     z_dot_bar_tilde[i]=0;
00306   }
00307 
00308   double df=(pf->df)(xu);
00309   double d2f=(pf->d2f)(xu);
00310   double d3f=(pf->d3f)(xu);
00311   for (size_t i=0;i<nvar;i++)
00312   {
00313     //*x_tilde+=(pf->d2f)(xu)*zbar[i]*x_bar_tilde[i];
00314     *x_tilde+=d2f*zbar[i]*x_bar_tilde[i];
00315     //z_bar_tilde[i]+=(pf->df)(xu)*x_bar_tilde[i];
00316     z_bar_tilde[i]+=df*x_bar_tilde[i];
00317   }
00318 
00319   for (size_t i=0;i<nvar;i++)
00320   {
00321     //*x_tilde+=(pf->d2f)(xu)*zdotbar[i]*x_dot_bar_tilde[i];
00322     *x_tilde+=d2f*zdotbar[i]*x_dot_bar_tilde[i];
00323     //z_dot_bar_tilde[i]+=(pf->df)(xu)*x_dot_bar_tilde[i];
00324     z_dot_bar_tilde[i]+=df*x_dot_bar_tilde[i];
00325   }
00326 
00327   for (size_t i=0;i<nvar;i++)
00328   {
00329     //x_dot_tilde[i]+=(pf->d2f)(xu)*zdotbar[i]*x_bar_tilde[i];
00330     //z_dot_bar_tilde[i]+=(pf->d2f)(xu)*xdot[i]*x_bar_tilde[i];
00331     //*x_tilde+=(pf->d3f)(xu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00332     x_dot_tilde[i]+=d2f*zdotbar[i]*x_bar_tilde[i];
00333     z_dot_bar_tilde[i]+=d2f*xdot[i]*x_bar_tilde[i];
00334     *x_tilde+=d3f*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00335   }
00336   list2.restoreposition();
00337 #if defined(PRINT_DERS)
00338  print_derivatives(pz,"z");
00339  print_derivatives(px,"x");
00340 #endif
00341 }
00342 
00347 void read_pass1_3x(void)
00348 {
00349   // We are going backword for bptr and forward for bptr2
00350   // the current entry+2 in bptr is the size of the record i.e
00351   // points to the next record
00352   unsigned int nvar=df1b2variable::nvar;
00353   fixed_smartlist & nlist=f1b2gradlist->nlist;
00354   test_smartlist& list=f1b2gradlist->list;
00355    // nlist-=sizeof(int);
00356   // get record size
00357   int num_bytes=nlist.bptr->numbytes;
00358   // backup the size of the record
00359   list-=num_bytes;
00360   list.saveposition(); // save pointer to beginning of record;
00361   // save the pointer to the beginning of the record
00362   double xu;
00363   double * xdot;
00364   //df1b2_header x,z;
00365   df1b2function1 * pf;
00366 
00367   // get info from tape1
00368 #if defined(SAFE_ARRAYS)
00369   checkidentiferstring("CX",list);
00370 #endif
00371   df1b2_header * px=(df1b2_header *) list.bptr;
00372   list.bptr+=sizeof(df1b2_header);
00373   df1b2_header * pz=(df1b2_header *) list.bptr;
00374   list.bptr+=sizeof(df1b2_header);
00375   pf=*(df1b2function1 **) list.bptr;
00376   list.bptr+=sizeof(char*);
00377   memcpy(&xu,list.bptr,sizeof(double));
00378   list.bptr+=sizeof(double);
00379   xdot=(double*)list.bptr;
00380   list.restoreposition(); // save pointer to beginning of record;
00381 
00382 #if defined(PRINT_DERS)
00383  print_derivatives(pf->funname,(pf->df)(xu),(pf->df)(xu),(pf->d2f)(xu),
00384   (pf->d3f)(xu),1);
00385  print_derivatives(pz,"z");
00386  print_derivatives(px,"x");
00387 #endif
00388 
00389   double df=(pf->df)(xu);
00390   double d2f=(pf->d2f)(xu);
00391   //*(px->u_tilde)+=(pf->df)(xu)* *(pz->u_tilde);
00392   *(px->u_tilde)+=df * *(pz->u_tilde);
00393   for (size_t i=0;i<nvar;i++)
00394   {
00395     //*(px->u_tilde)+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_tilde[i];
00396     *(px->u_tilde)+=d2f*xdot[i]*pz->u_dot_tilde[i];
00397   }
00398   for (size_t i=0;i<nvar;i++)
00399   {
00400     //px->u_dot_tilde[i]+=(pf->df)(xu)*pz->u_dot_tilde[i];
00401     px->u_dot_tilde[i]+=df*pz->u_dot_tilde[i];
00402   }
00403   *(pz->u_tilde)=0;
00404   for (size_t i=0;i<nvar;i++)
00405   {
00406     pz->u_dot_tilde[i]=0;
00407   }
00408 #if defined(PRINT_DERS)
00409  print_derivatives(pz,"z");
00410  print_derivatives(px,"x");
00411 #endif
00412 }