ADMB Documentation  11.1.2490
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f21.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f21.cpp 2329 2014-09-12 00:30:49Z 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 
00043   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)
00044    // +sizeof(char*)
00045     +4*sizeof(double)+nvar*sizeof(double);
00046 // string identifier debug stuff
00047 #if defined(SAFE_ALL)
00048   char ids[]="HU";
00049   int slen=strlen(ids);
00050   total_bytes+=slen;
00051 #endif
00052 
00053 #ifndef OPT_LIB
00054   assert(total_bytes <= (size_t)INT_MAX);
00055 #endif
00056 
00057   list.check_buffer_size((int)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   fixed_smartlist & nlist=f1b2gradlist->nlist;
00122   test_smartlist& list=f1b2gradlist->list;
00123    // nlist-=sizeof(int);
00124   // get record size
00125   int num_bytes=nlist.bptr->numbytes;
00126   // backup the size of the record
00127   list-=num_bytes;
00128   //list.bptr-=num_bytes;
00129   list.saveposition(); // save pointer to beginning of record;
00130   // save the pointer to the beginning of the record
00131   double xu;
00132   double * xdot;
00133   //df1b2_header x,z;
00134   //df1b2function1 * pf;
00135 
00136   // get info from tape1
00137   // get info from tape1
00138 #if defined(SAFE_ALL)
00139   checkidentiferstring("HU",list);
00140 #endif
00141   df1b2_header * px=(df1b2_header *) list.bptr;
00142   list.bptr+=sizeof(df1b2_header);
00143   df1b2_header * pz=(df1b2_header *) list.bptr;
00144   list.bptr+=sizeof(df1b2_header);
00145   double df=*(double*) list.bptr;
00146   list.bptr+=sizeof(double);
00147   double d2f=*(double*) list.bptr;
00148   list.bptr+=sizeof(double);
00149   //double d3f=*(double*) list.bptr;
00150   list.bptr+=sizeof(double);
00151   memcpy(&xu,list.bptr,sizeof(double));
00152   list.bptr+=sizeof(double);
00153   xdot=(double*)list.bptr;
00154   list.restoreposition(); // save pointer to beginning of record;
00155   int i;
00156 
00157   // Do first reverse paSS calculations
00158   // ****************************************************************
00159   // turn this off if no third derivatives are calculated
00160   // if (!no_third_derivatives)
00161   // {
00162   // save for second reverse pass
00163   // save identifier 1
00164      fixed_smartlist2& nlist2=f1b2gradlist->nlist2;
00165      test_smartlist& list2=f1b2gradlist->list2;
00166 
00167   size_t total_bytes=2*nvar*sizeof(double);
00168 // string identifier debug stuff
00169 #if defined(SAFE_ALL)
00170   char ids[]="JE";
00171   int slen=strlen(ids);
00172   total_bytes+=slen;
00173 #endif
00174 
00175 #ifndef OPT_LIB
00176   assert(total_bytes <= (size_t)INT_MAX);
00177 #endif
00178 
00179   list2.check_buffer_size((int)total_bytes);
00180   void * tmpptr2=list2.bptr;
00181 
00182 #if defined(SAFE_ALL)
00183   memcpy(list2,ids,slen);
00184 #endif
00185 
00186    memcpy(list2,pz->get_u_bar(),nvar*(int)sizeof(double));
00187    memcpy(list2,pz->get_u_dot_bar(),nvar*(int)sizeof(double));
00188    *nlist2.bptr=adptr_diff(list2.bptr,tmpptr2);
00189    ++nlist2;
00190   // }
00191   //
00192   // ****************************************************************
00193 #if defined(PRINT_DERS)
00194  print_derivatives(pz,"z");
00195  print_derivatives(px,"x");
00196 #endif
00197 
00198   //double df=(pf->df)(xu);
00199   //double d2f=(pf->d2f)(xu);
00200   //double d3f=(pf->d3f)(xu);
00201 
00202   for (i=0;i<nvar;i++)
00203   {
00204     //px->u_bar[i]+=(pf->df)(xu)* pz->u_bar[i];
00205     px->u_bar[i]+=df * pz->u_bar[i];
00206   }
00207   for (i=0;i<nvar;i++)
00208   {
00209     //px->u_bar[i]+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_bar[i];
00210     px->u_bar[i]+=d2f*xdot[i]*pz->u_dot_bar[i];
00211   }
00212   for (i=0;i<nvar;i++)
00213   {
00214     //px->u_dot_bar[i]+=(pf->df)(xu)*pz->u_dot_bar[i];
00215     px->u_dot_bar[i]+=df*pz->u_dot_bar[i];
00216   }
00217 
00218   // !!!!!!!!!!!!!!!!!!!!!!
00219   for (i=0;i<nvar;i++)
00220   {
00221     pz->u_bar[i]=0;
00222   }
00223   for (i=0;i<nvar;i++)
00224   {
00225     pz->u_dot_bar[i]=0;
00226   }
00227 
00228 #if defined(PRINT_DERS)
00229  print_derivatives(pz,"z");
00230  print_derivatives(px,"x");
00231 #endif
00232 }
00233 
00238 void read_pass1_2_dv(void)
00239 {
00240   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00241   // We are going forward for bptr and backword for bptr2
00242   //
00243   // list 1
00244   //
00245   int nvar=df1b2variable::nvar;
00246   test_smartlist & list=f1b2gradlist->list;
00247 
00248   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header)+sizeof(char*)
00249     +sizeof(double)+nvar*sizeof(double);
00250 // string identifier debug stuff
00251 #if defined(SAFE_ALL)
00252   char ids[]="HU";
00253   int slen=strlen(ids);
00254   total_bytes+=slen;
00255 #endif
00256 
00257 #ifndef OPT_LIB
00258   assert(total_bytes <= (size_t)INT_MAX);
00259 #endif
00260 
00261   list.check_buffer_size((int)total_bytes);
00262 // end of string identifier debug stuff
00263 
00264   list.saveposition(); // save pointer to beginning of record;
00265   fixed_smartlist & nlist=f1b2gradlist->nlist;
00266    // nlist-=sizeof(int);
00267   // get record size
00268   int num_bytes=nlist.bptr->numbytes;
00269     // nlist+=nlist_record_size;
00270   //
00271   // list 2
00272   //
00273   test_smartlist & list2=f1b2gradlist->list2;
00274   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00275   // get record size
00276   int num_bytes2=*nlist2.bptr;
00277   --nlist2;
00278   // backup the size of the record
00279   list2-=num_bytes2;
00280   list2.saveposition();
00281   // save the pointer to the beginning of the record
00282   // bptr and bptr2 now both point to the beginning of their records
00283 
00284   double xu;
00285   double * xdot;
00286   //df1b2_header x,z;
00287   //df1b2function1 * pf;
00288 
00289   // get info from tape1
00290   // get info from tape1
00291 #if defined(SAFE_ALL)
00292   checkidentiferstring("HU",list);
00293   checkidentiferstring("JE",list2);
00294 #endif
00295   df1b2_header * px=(df1b2_header *) list.bptr;
00296   list.bptr+=sizeof(df1b2_header);
00297   df1b2_header * pz=(df1b2_header *) list.bptr;
00298   list.bptr+=sizeof(df1b2_header);
00299 
00300   //pf=*(df1b2function1 **) list.bptr;
00301 
00302   double df=*(double*) list.bptr;
00303   list.bptr+=sizeof(double);
00304   double d2f=*(double*) list.bptr;
00305   list.bptr+=sizeof(double);
00306   double d3f=*(double*) list.bptr;
00307   list.bptr+=sizeof(double);
00308   //list.bptr+=sizeof(char*);
00309   memcpy(&xu,list.bptr,sizeof(double));
00310   list.bptr+=sizeof(double);
00311   xdot=(double*)list.bptr;
00312   list.restoreposition(num_bytes); // save pointer to beginning of record;
00313 
00314   double* zbar = (double*)list2.bptr;
00315   double* zdotbar = (double*)(list2.bptr+nvar*sizeof(double));
00316 
00317   double* x_tilde=px->get_u_tilde();
00318   double* x_dot_tilde=px->get_u_dot_tilde();
00319   double* x_bar_tilde=px->get_u_bar_tilde();
00320   double* x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00321   double* z_bar_tilde=pz->get_u_bar_tilde();
00322   double* z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00323 #if defined(PRINT_DERS)
00324  print_derivatives(pz,"z");
00325  print_derivatives(px,"x");
00326 #endif
00327   // Do second "reverse-reverse" pass calculations
00328   int i;
00329 
00330   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00331   for (i=0;i<nvar;i++)
00332   {
00333     z_bar_tilde[i]=0;
00334     z_dot_bar_tilde[i]=0;
00335   }
00336 
00337   //double df=(pf->df)(xu);
00338   //double d2f=(pf->d2f)(xu);
00339   //double d3f=(pf->d3f)(xu);
00340   for (i=0;i<nvar;i++)
00341   {
00342     //*x_tilde+=(pf->d2f)(xu)*zbar[i]*x_bar_tilde[i];
00343     *x_tilde+=d2f*zbar[i]*x_bar_tilde[i];
00344     //z_bar_tilde[i]+=(pf->df)(xu)*x_bar_tilde[i];
00345     z_bar_tilde[i]+=df*x_bar_tilde[i];
00346   }
00347 
00348   for (i=0;i<nvar;i++)
00349   {
00350     //*x_tilde+=(pf->d2f)(xu)*zdotbar[i]*x_dot_bar_tilde[i];
00351     *x_tilde+=d2f*zdotbar[i]*x_dot_bar_tilde[i];
00352     //z_dot_bar_tilde[i]+=(pf->df)(xu)*x_dot_bar_tilde[i];
00353     z_dot_bar_tilde[i]+=df*x_dot_bar_tilde[i];
00354   }
00355 
00356   for (i=0;i<nvar;i++)
00357   {
00358     //x_dot_tilde[i]+=(pf->d2f)(xu)*zdotbar[i]*x_bar_tilde[i];
00359     //z_dot_bar_tilde[i]+=(pf->d2f)(xu)*xdot[i]*x_bar_tilde[i];
00360     //*x_tilde+=(pf->d3f)(xu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00361     x_dot_tilde[i]+=d2f*zdotbar[i]*x_bar_tilde[i];
00362     z_dot_bar_tilde[i]+=d2f*xdot[i]*x_bar_tilde[i];
00363     *x_tilde+=d3f*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00364   }
00365   list2.restoreposition();
00366 #if defined(PRINT_DERS)
00367  print_derivatives(pz,"z");
00368  print_derivatives(px,"x");
00369 #endif
00370 }
00371 
00376 void read_pass1_3_dv(void)
00377 {
00378   // We are going backword for bptr and forward for bptr2
00379   // the current entry+2 in bptr is the size of the record i.e
00380   // points to the next record
00381   int nvar=df1b2variable::nvar;
00382   fixed_smartlist & nlist=f1b2gradlist->nlist;
00383   test_smartlist& list=f1b2gradlist->list;
00384    // nlist-=sizeof(int);
00385   // get record size
00386   int num_bytes=nlist.bptr->numbytes;
00387   // backup the size of the record
00388   list-=num_bytes;
00389   list.saveposition(); // save pointer to beginning of record;
00390   // save the pointer to the beginning of the record
00391   double xu;
00392   double * xdot;
00393   //df1b2_header x,z;
00394   //df1b2function1 * pf;
00395 
00396   // get info from tape1
00397 #if defined(SAFE_ALL)
00398   checkidentiferstring("HU",list);
00399 #endif
00400   df1b2_header * px=(df1b2_header *) list.bptr;
00401   list.bptr+=sizeof(df1b2_header);
00402   df1b2_header * pz=(df1b2_header *) list.bptr;
00403   list.bptr+=sizeof(df1b2_header);
00404   //pf=*(df1b2function1 **) list.bptr;
00405   //list.bptr+=sizeof(char*);
00406   double df=*(double*) list.bptr;
00407   list.bptr+=sizeof(double);
00408   double d2f=*(double*) list.bptr;
00409   list.bptr+=sizeof(double);
00410   //double d3f=*(double*) list.bptr;
00411   list.bptr+=sizeof(double);
00412   memcpy(&xu,list.bptr,sizeof(double));
00413   list.bptr+=sizeof(double);
00414   xdot=(double*)list.bptr;
00415   list.restoreposition(); // save pointer to beginning of record;
00416   int i;
00417 
00418 #if defined(PRINT_DERS)
00419  print_derivatives(pz,"z");
00420  print_derivatives(px,"x");
00421 #endif
00422 
00423   //double df=(pf->df)(xu);
00424   //double d2f=(pf->d2f)(xu);
00425   //*(px->u_tilde)+=(pf->df)(xu)* *(pz->u_tilde);
00426   *(px->u_tilde)+=df * *(pz->u_tilde);
00427   for (i=0;i<nvar;i++)
00428   {
00429     //*(px->u_tilde)+=(pf->d2f)(xu)*xdot[i]*pz->u_dot_tilde[i];
00430     *(px->u_tilde)+=d2f*xdot[i]*pz->u_dot_tilde[i];
00431   }
00432   for (i=0;i<nvar;i++)
00433   {
00434     //px->u_dot_tilde[i]+=(pf->df)(xu)*pz->u_dot_tilde[i];
00435     px->u_dot_tilde[i]+=df*pz->u_dot_tilde[i];
00436   }
00437   *(pz->u_tilde)=0;
00438   for (i=0;i<nvar;i++)
00439   {
00440     pz->u_dot_tilde[i]=0;
00441   }
00442 #if defined(PRINT_DERS)
00443  print_derivatives(pz,"z");
00444  print_derivatives(px,"x");
00445 #endif
00446 }