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