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