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