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