ADMB Documentation  11.1.2545
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f21.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f21.cpp 2535 2014-10-31 12:28:33Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00012 #include <df1b2fun.h>
00013 
00014 #ifndef OPT_LIB
00015   #include <cassert>
00016   #include <climits>
00017 #endif
00018 
00019 void ad_read_pass1_der_values(void);
00020 void read_pass1_1_dv(void);
00021 void read_pass1_2_dv(void);
00022 void read_pass1_3_dv(void);
00023 
00028 int df1b2_gradlist::write_pass1(const df1b2variable * _px,
00029   df1b2variable * pz, double df, double d2f, double d3f)
00030 {
00031   ADUNCONST(df1b2variable*,px)
00032   ncount++;
00033  /*
00034   if (28126<ncount && 28129> ncount)
00035     cout << "trap" << endl;
00036   */
00037 #if defined(CHECK_COUNT)
00038   if (ncount >= ncount_check)
00039     ncount_checker(ncount,ncount_check);
00040 #endif
00041   int _nvar=df1b2variable::nvar;
00042 #ifndef OPT_LIB
00043   assert(_nvar > 0);
00044 #endif
00045   size_t nvar= (size_t)_nvar;
00046 
00047   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)
00048    // +sizeof(char*)
00049     +4*sizeof(double)+nvar*sizeof(double);
00050 // string identifier debug stuff
00051 #if defined(SAFE_ALL)
00052   char ids[]="HU";
00053   int slen=strlen(ids);
00054   total_bytes+=slen;
00055 #endif
00056 
00057   list.check_buffer_size(total_bytes);
00058   void * tmpptr=list.bptr;
00059 #if defined(SAFE_ALL)
00060   memcpy(list,ids,slen);
00061 #endif
00062 // end of string identifier debug stuff
00063 
00064   memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00065   memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00066 
00067   // replace this with the derivative values
00068   //memcpy(list,&pf,sizeof(char*));
00069 
00070   memcpy(list,&df,sizeof(double));
00071   memcpy(list,&d2f,sizeof(double));
00072   memcpy(list,&d3f,sizeof(double));
00073 
00074 
00075 
00076   memcpy(list,px->get_u(),sizeof(double));
00077   //list.bptr+=sizeof(double);
00078   memcpy(list,px->get_u_dot(),nvar*(int)sizeof(double));
00079   //list.bptr+=nvar*sizeof(double);
00080   // ***** write  record size
00081   nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00082   nlist.bptr->pf=(ADrfptr)(&ad_read_pass1_der_values);
00083   ++nlist;
00084   return 0;
00085 }
00086 
00091 void ad_read_pass1_der_values(void)
00092 {
00093   switch(df1b2variable::passnumber)
00094   {
00095   case 1:
00096     read_pass1_1_dv();
00097     break;
00098   case 2:
00099     read_pass1_2_dv();
00100     break;
00101   case 3:
00102     read_pass1_3_dv();
00103     break;
00104   default:
00105     cerr << "illegal value for df1b2variable::pass = "
00106          << df1b2variable::passnumber << endl;
00107     exit(1);
00108   }
00109 }
00110 
00115 void read_pass1_1_dv(void)
00116 {
00117   // We are going backword for bptr and forward for bptr2
00118   // the current entry+2 in bptr is the size of the record i.e
00119   // points to the next record
00120   int _nvar=df1b2variable::nvar;
00121 #ifndef OPT_LIB
00122   assert(_nvar > 0);
00123 #endif
00124   size_t nvar= (size_t)_nvar;
00125   fixed_smartlist & nlist=f1b2gradlist->nlist;
00126   test_smartlist& list=f1b2gradlist->list;
00127    // nlist-=sizeof(int);
00128   // get record size
00129   int num_bytes=nlist.bptr->numbytes;
00130   // backup the size of the record
00131   list-=num_bytes;
00132   //list.bptr-=num_bytes;
00133   list.saveposition(); // save pointer to beginning of record;
00134   // save the pointer to the beginning of the record
00135   double xu;
00136   double * xdot;
00137   //df1b2_header x,z;
00138   //df1b2function1 * pf;
00139 
00140   // get info from tape1
00141   // get info from tape1
00142 #if defined(SAFE_ALL)
00143   checkidentiferstring("HU",list);
00144 #endif
00145   df1b2_header * px=(df1b2_header *) list.bptr;
00146   list.bptr+=sizeof(df1b2_header);
00147   df1b2_header * pz=(df1b2_header *) list.bptr;
00148   list.bptr+=sizeof(df1b2_header);
00149   double df=*(double*) list.bptr;
00150   list.bptr+=sizeof(double);
00151   double d2f=*(double*) list.bptr;
00152   list.bptr+=sizeof(double);
00153   //double d3f=*(double*) list.bptr;
00154   list.bptr+=sizeof(double);
00155   memcpy(&xu,list.bptr,sizeof(double));
00156   list.bptr+=sizeof(double);
00157   xdot=(double*)list.bptr;
00158   list.restoreposition(); // save pointer to beginning of record;
00159 
00160   // Do first reverse paSS calculations
00161   // ****************************************************************
00162   // turn this off if no third derivatives are calculated
00163   // if (!no_third_derivatives)
00164   // {
00165   // save for second reverse pass
00166   // save identifier 1
00167      fixed_smartlist2& nlist2=f1b2gradlist->nlist2;
00168      test_smartlist& list2=f1b2gradlist->list2;
00169 
00170   size_t total_bytes=2*nvar*sizeof(double);
00171 // string identifier debug stuff
00172 #if defined(SAFE_ALL)
00173   char ids[]="JE";
00174   int slen=strlen(ids);
00175   total_bytes+=slen;
00176 #endif
00177 
00178   list2.check_buffer_size(total_bytes);
00179   void * tmpptr2=list2.bptr;
00180 
00181 #if defined(SAFE_ALL)
00182   memcpy(list2,ids,slen);
00183 #endif
00184 
00185    memcpy(list2,pz->get_u_bar(),nvar*(int)sizeof(double));
00186    memcpy(list2,pz->get_u_dot_bar(),nvar*(int)sizeof(double));
00187    *nlist2.bptr=adptr_diff(list2.bptr,tmpptr2);
00188    ++nlist2;
00189   // }
00190   //
00191   // ****************************************************************
00192 #if defined(PRINT_DERS)
00193  print_derivatives(pz,"z");
00194  print_derivatives(px,"x");
00195 #endif
00196 
00197   //double df=(pf->df)(xu);
00198   //double d2f=(pf->d2f)(xu);
00199   //double d3f=(pf->d3f)(xu);
00200 
00201   for (size_t i=0;i<nvar;i++)
00202   {
00203     //px->u_bar[i]+=(pf->df)(xu)* pz->u_bar[i];
00204     px->u_bar[i]+=df * pz->u_bar[i];
00205   }
00206   for (size_t i=0;i<nvar;i++)
00207   {
00208     //px->u_bar[i]+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_bar[i];
00209     px->u_bar[i]+=d2f*xdot[i]*pz->u_dot_bar[i];
00210   }
00211   for (size_t i=0;i<nvar;i++)
00212   {
00213     //px->u_dot_bar[i]+=(pf->df)(xu)*pz->u_dot_bar[i];
00214     px->u_dot_bar[i]+=df*pz->u_dot_bar[i];
00215   }
00216 
00217   // !!!!!!!!!!!!!!!!!!!!!!
00218   for (size_t i=0;i<nvar;i++)
00219   {
00220     pz->u_bar[i]=0;
00221   }
00222   for (size_t i=0;i<nvar;i++)
00223   {
00224     pz->u_dot_bar[i]=0;
00225   }
00226 
00227 #if defined(PRINT_DERS)
00228  print_derivatives(pz,"z");
00229  print_derivatives(px,"x");
00230 #endif
00231 }
00232 
00237 void read_pass1_2_dv(void)
00238 {
00239   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00240   // We are going forward for bptr and backword for bptr2
00241   //
00242   // list 1
00243   //
00244   int _nvar=df1b2variable::nvar;
00245 #ifndef OPT_LIB
00246   assert(_nvar > 0);
00247 #endif
00248   size_t nvar= (size_t)_nvar;
00249   test_smartlist & list=f1b2gradlist->list;
00250 
00251   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)+sizeof(char*)
00252     +sizeof(double)+nvar*sizeof(double);
00253 // string identifier debug stuff
00254 #if defined(SAFE_ALL)
00255   char ids[]="HU";
00256   int slen=strlen(ids);
00257   total_bytes+=slen;
00258 #endif
00259 
00260   list.check_buffer_size(total_bytes);
00261 // end of string identifier debug stuff
00262 
00263   list.saveposition(); // save pointer to beginning of record;
00264   fixed_smartlist & nlist=f1b2gradlist->nlist;
00265    // nlist-=sizeof(int);
00266   // get record size
00267   int num_bytes=nlist.bptr->numbytes;
00268     // nlist+=nlist_record_size;
00269   //
00270   // list 2
00271   //
00272   test_smartlist & list2=f1b2gradlist->list2;
00273   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00274   // get record size
00275   int num_bytes2=*nlist2.bptr;
00276   --nlist2;
00277   // backup the size of the record
00278   list2-=num_bytes2;
00279   list2.saveposition();
00280   // save the pointer to the beginning of the record
00281   // bptr and bptr2 now both point to the beginning of their records
00282 
00283   double xu;
00284   double * xdot;
00285   //df1b2_header x,z;
00286   //df1b2function1 * pf;
00287 
00288   // get info from tape1
00289   // get info from tape1
00290 #if defined(SAFE_ALL)
00291   checkidentiferstring("HU",list);
00292   checkidentiferstring("JE",list2);
00293 #endif
00294   df1b2_header * px=(df1b2_header *) list.bptr;
00295   list.bptr+=sizeof(df1b2_header);
00296   df1b2_header * pz=(df1b2_header *) list.bptr;
00297   list.bptr+=sizeof(df1b2_header);
00298 
00299   //pf=*(df1b2function1 **) list.bptr;
00300 
00301   double df=*(double*) list.bptr;
00302   list.bptr+=sizeof(double);
00303   double d2f=*(double*) list.bptr;
00304   list.bptr+=sizeof(double);
00305   double d3f=*(double*) list.bptr;
00306   list.bptr+=sizeof(double);
00307   //list.bptr+=sizeof(char*);
00308   memcpy(&xu,list.bptr,sizeof(double));
00309   list.bptr+=sizeof(double);
00310   xdot=(double*)list.bptr;
00311   list.restoreposition(num_bytes); // save pointer to beginning of record;
00312 
00313   double* zbar = (double*)list2.bptr;
00314   double* zdotbar = (double*)(list2.bptr+nvar*sizeof(double));
00315 
00316   double* x_tilde=px->get_u_tilde();
00317   double* x_dot_tilde=px->get_u_dot_tilde();
00318   double* x_bar_tilde=px->get_u_bar_tilde();
00319   double* x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00320   double* z_bar_tilde=pz->get_u_bar_tilde();
00321   double* z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00322 #if defined(PRINT_DERS)
00323  print_derivatives(pz,"z");
00324  print_derivatives(px,"x");
00325 #endif
00326   // Do second "reverse-reverse" pass calculations
00327 
00328   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00329   for (size_t i=0;i<nvar;i++)
00330   {
00331     z_bar_tilde[i]=0;
00332     z_dot_bar_tilde[i]=0;
00333   }
00334 
00335   //double df=(pf->df)(xu);
00336   //double d2f=(pf->d2f)(xu);
00337   //double d3f=(pf->d3f)(xu);
00338   for (size_t i=0;i<nvar;i++)
00339   {
00340     //*x_tilde+=(pf->d2f)(xu)*zbar[i]*x_bar_tilde[i];
00341     *x_tilde+=d2f*zbar[i]*x_bar_tilde[i];
00342     //z_bar_tilde[i]+=(pf->df)(xu)*x_bar_tilde[i];
00343     z_bar_tilde[i]+=df*x_bar_tilde[i];
00344   }
00345 
00346   for (size_t i=0;i<nvar;i++)
00347   {
00348     //*x_tilde+=(pf->d2f)(xu)*zdotbar[i]*x_dot_bar_tilde[i];
00349     *x_tilde+=d2f*zdotbar[i]*x_dot_bar_tilde[i];
00350     //z_dot_bar_tilde[i]+=(pf->df)(xu)*x_dot_bar_tilde[i];
00351     z_dot_bar_tilde[i]+=df*x_dot_bar_tilde[i];
00352   }
00353 
00354   for (size_t i=0;i<nvar;i++)
00355   {
00356     //x_dot_tilde[i]+=(pf->d2f)(xu)*zdotbar[i]*x_bar_tilde[i];
00357     //z_dot_bar_tilde[i]+=(pf->d2f)(xu)*xdot[i]*x_bar_tilde[i];
00358     //*x_tilde+=(pf->d3f)(xu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00359     x_dot_tilde[i]+=d2f*zdotbar[i]*x_bar_tilde[i];
00360     z_dot_bar_tilde[i]+=d2f*xdot[i]*x_bar_tilde[i];
00361     *x_tilde+=d3f*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00362   }
00363   list2.restoreposition();
00364 #if defined(PRINT_DERS)
00365  print_derivatives(pz,"z");
00366  print_derivatives(px,"x");
00367 #endif
00368 }
00369 
00374 void read_pass1_3_dv(void)
00375 {
00376   // We are going backword for bptr and forward for bptr2
00377   // the current entry+2 in bptr is the size of the record i.e
00378   // points to the next record
00379   int _nvar=df1b2variable::nvar;
00380 #ifndef OPT_LIB
00381   assert(_nvar > 0);
00382 #endif
00383   size_t nvar= (size_t)_nvar;
00384   fixed_smartlist & nlist=f1b2gradlist->nlist;
00385   test_smartlist& list=f1b2gradlist->list;
00386    // nlist-=sizeof(int);
00387   // get record size
00388   int num_bytes=nlist.bptr->numbytes;
00389   // backup the size of the record
00390   list-=num_bytes;
00391   list.saveposition(); // save pointer to beginning of record;
00392   // save the pointer to the beginning of the record
00393   double xu;
00394   double * xdot;
00395   //df1b2_header x,z;
00396   //df1b2function1 * pf;
00397 
00398   // get info from tape1
00399 #if defined(SAFE_ALL)
00400   checkidentiferstring("HU",list);
00401 #endif
00402   df1b2_header * px=(df1b2_header *) list.bptr;
00403   list.bptr+=sizeof(df1b2_header);
00404   df1b2_header * pz=(df1b2_header *) list.bptr;
00405   list.bptr+=sizeof(df1b2_header);
00406   //pf=*(df1b2function1 **) list.bptr;
00407   //list.bptr+=sizeof(char*);
00408   double df=*(double*) list.bptr;
00409   list.bptr+=sizeof(double);
00410   double d2f=*(double*) list.bptr;
00411   list.bptr+=sizeof(double);
00412   //double d3f=*(double*) list.bptr;
00413   list.bptr+=sizeof(double);
00414   memcpy(&xu,list.bptr,sizeof(double));
00415   list.bptr+=sizeof(double);
00416   xdot=(double*)list.bptr;
00417   list.restoreposition(); // save pointer to beginning of record;
00418 
00419 #if defined(PRINT_DERS)
00420  print_derivatives(pz,"z");
00421  print_derivatives(px,"x");
00422 #endif
00423 
00424   //double df=(pf->df)(xu);
00425   //double d2f=(pf->d2f)(xu);
00426   //*(px->u_tilde)+=(pf->df)(xu)* *(pz->u_tilde);
00427   *(px->u_tilde)+=df * *(pz->u_tilde);
00428   for (size_t i=0;i<nvar;i++)
00429   {
00430     //*(px->u_tilde)+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_tilde[i];
00431     *(px->u_tilde)+=d2f*xdot[i]*pz->u_dot_tilde[i];
00432   }
00433   for (size_t i=0;i<nvar;i++)
00434   {
00435     //px->u_dot_tilde[i]+=(pf->df)(xu)*pz->u_dot_tilde[i];
00436     px->u_dot_tilde[i]+=df*pz->u_dot_tilde[i];
00437   }
00438   *(pz->u_tilde)=0;
00439   for (size_t i=0;i<nvar;i++)
00440   {
00441     pz->u_dot_tilde[i]=0;
00442   }
00443 #if defined(PRINT_DERS)
00444  print_derivatives(pz,"z");
00445  print_derivatives(px,"x");
00446 #endif
00447 }