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