ADMB Documentation  11.1.1958
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f22.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f22.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  */
00012 #include "df1b2test.h"
00013 
00018   df1b2function2c::df1b2function2c(double (*_f)(double,double),
00019     double (*_df)(double,double),
00020     double (*_d2f)(double,double),
00021     double (*_d3f)(double,double),
00022     const adstring & _s)
00023   {
00024     f=_f;
00025     df=_df;
00026     d2f=_d2f;
00027     d3f=_d3f;
00028     funname=_s;
00029   }
00030 
00035   df1b2variable df1b2function2c::operator () (const df1b2variable& _x,
00036     double y)
00037   {
00038     ADUNCONST(df1b2variable,x)
00039     df1b2variable z;
00040     double xu=*x.get_u();
00041     double * xd=x.get_u_dot();
00042     double * zd=z.get_u_dot();
00043     *z.get_u()=(*f)(xu,y);
00044     double dfx=(*df)(xu,y);
00045     for (int i=0;i<df1b2variable::nvar;i++)
00046     {
00047       *zd++ =dfx * *xd++ ;
00048     }
00049 
00050     // WRITE WHATEVER ON TAPE
00051     if (!df1b2_gradlist::no_derivatives)
00052       f1b2gradlist->write_pass1c(&x,y,&z,this);
00053     return z;
00054   }
00055 
00056 void ad_read_pass2c(void);
00057 
00062  int df1b2_gradlist::write_pass1c(const df1b2variable * _px,
00063    double py,df1b2variable * pz,df1b2function2c * pf)
00064  {
00065    ADUNCONST(df1b2variable*,px)
00066    ncount++;
00067 #if defined(CHECK_COUNT)
00068   if (ncount >= ncount_check)
00069     ncount_checker(ncount,ncount_check);
00070 #endif
00071    int nvar=df1b2variable::nvar;
00072 
00073    int total_bytes=2*sizeof(df1b2_header)+sizeof(char*)
00074      +(nvar+2)*sizeof(double);
00075 // string identifier debug stuff
00076 #if defined(SAFE_ALL)
00077   char ids[]="BY";
00078   int slen=strlen(ids);
00079   total_bytes+=slen;
00080 #endif
00081   list.check_buffer_size(total_bytes);
00082   void * tmpptr=list.bptr;
00083 #if defined(SAFE_ALL)
00084   memcpy(list,ids,slen);
00085 #endif
00086 // end of string identifier debug stuff
00087 
00088    memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00089    memcpy(list,&py,sizeof(double));
00090    memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00091    memcpy(list,&pf,sizeof(char *));
00092    //*(char**)(list.bptr)=(char*)pf;
00093    memcpy(list,px->get_u(),sizeof(double));
00094    //memcpy(list,py->get_u(),sizeof(double));
00095    memcpy(list,px->get_u_dot(),nvar*sizeof(double));
00096    //memcpy(list,py->get_u_dot(),nvar*sizeof(double));
00097    // ***** write  record size
00098    nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00099    nlist.bptr->pf=(ADrfptr)(&ad_read_pass2c);
00100    ++nlist;
00101    return 0;
00102  }
00103 
00104 
00105 void read_pass2_1c(void);
00106 void read_pass2_2c(void);
00107 void read_pass2_3c(void);
00108 
00113 void ad_read_pass2c(void)
00114 {
00115   switch(df1b2variable::passnumber)
00116   {
00117   case 1:
00118     read_pass2_1c();
00119     break;
00120   case 2:
00121     read_pass2_2c();
00122     break;
00123   case 3:
00124     read_pass2_3c();
00125     break;
00126   default:
00127     cerr << "illegal value for df1b2variable::pass = "
00128          << df1b2variable::passnumber << endl;
00129     exit(1);
00130   }
00131 }
00132 
00137 void read_pass2_1c(void)
00138 {
00139   // We are going backword for bptr and nbptr
00140   // and  forward for bptr2 and nbptr2
00141   // the current entry+2 in bptr is the size of the record i.e
00142   // points to the next record
00143   //char * bptr=f1b2gradlist->bptr;
00144   //char * bptr2=f1b2gradlist2->bptr;
00145   int nvar=df1b2variable::nvar;
00146   test_smartlist& list=f1b2gradlist->list;
00147   //f1b2gradlist->nlist-=sizeof(int);
00148   int num_bytes=f1b2gradlist->nlist.bptr->numbytes;
00149   list-=num_bytes;
00150   list.saveposition(); // save pointer to beginning of record;
00151   double xu, yu = 0;
00152   //ad_dstar xdot,ydot;
00153   df1b2function2c * pf;
00154 
00155   // get info from tape1
00156 #if defined(SAFE_ARRAYS)
00157   checkidentiferstring("BY",f1b2gradlist->list);
00158 #endif
00159   char * bptr=f1b2gradlist->list.bptr;
00160   df1b2_header * px=(df1b2_header *) bptr;
00161   bptr+=sizeof(df1b2_header);
00162   //double * py=(double *) bptr;
00163   bptr+=sizeof(double);
00164   df1b2_header * pz=(df1b2_header *) bptr;
00165   bptr+=sizeof(df1b2_header);
00166   pf=*(df1b2function2c **) bptr;
00167   bptr+=sizeof(char*);
00168   memcpy(&xu,bptr,sizeof(double));
00169   bptr+=sizeof(double);
00170   //memcpy(&yu,bptr,sizeof(double));
00171   //bptr+=sizeof(double);
00172   double * xdot=(double*)bptr;
00173   //bptr+=nvar*sizeof(double);
00174   //double * ydot=(double*)bptr;
00175 
00176   list.restoreposition(); // save pointer to beginning of record;
00177 
00178   // ****************************************************************
00179   // turn this off if no third derivatives are calculated
00180   // if (!no_third_derivatives)
00181   // {
00182   // save for second reverse pass
00183   // save identifier 1
00184      test_smartlist & list2 = f1b2gradlist->list2;
00185 
00186 
00187    int total_bytes=2*nvar*sizeof(double);
00188 // string identifier debug stuff
00189 #if defined(SAFE_ALL)
00190   char ids[]="FW";
00191   int slen=strlen(ids);
00192   total_bytes+=slen;
00193 #endif
00194   list2.check_buffer_size(total_bytes);
00195   void * tmpptr=list2.bptr;
00196 #if defined(SAFE_ALL)
00197   memcpy(list2,ids,slen);
00198 #endif
00199 
00200      fixed_smartlist2 & nlist2 = f1b2gradlist->nlist2;
00201      memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00202      memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00203      *nlist2.bptr=adptr_diff(list2.bptr,tmpptr);
00204      ++nlist2;
00205   // }
00206   //
00207   // ****************************************************************
00208 #if defined(PRINT_DERS)
00209  print_derivatives(pz,"z");
00210  print_derivatives(px,"x");
00211 #endif
00212 #if defined(__DERCHECK__)
00213   if (derchecker)
00214   if (derchecker->node_number)
00215   {
00216     if (derchecker->counter == derchecker->node_number)
00217     {
00218       switch (derchecker->pass_number) // increment the variable of interest
00219       {
00220       case 2:
00221         switch(derchecker->vartype)
00222         {
00223         case 1:
00224           if (!derchecker->dotflag)
00225             px->u_bar[derchecker->index-1]+=derchecker->delta;
00226           else
00227             px->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00228           break;
00229         case 2:
00230         case 3:
00231           if (!derchecker->dotflag)
00232             pz->u_bar[derchecker->index-1]+=derchecker->delta;
00233           else
00234             pz->u_dot_bar[derchecker->index-1]+=derchecker->delta;
00235           break;
00236         default:
00237           cerr << "Invalid index value for dercheck_index was "
00238                << derchecker->index << endl;
00239           break;
00240         }
00241         break;
00242       case 3:
00243         switch(derchecker->vartype)
00244         {
00245         case 1:
00246           if (!derchecker->dotflag)
00247             px->u_bar[derchecker->index-1]-=derchecker->delta;
00248           else
00249             px->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00250           break;
00251         case 2:
00252         case 3:
00253           if (!derchecker->dotflag)
00254             pz->u_bar[derchecker->index-1]-=derchecker->delta;
00255           else
00256             pz->u_dot_bar[derchecker->index-1]-=derchecker->delta;
00257           break;
00258         default:
00259           cerr << "Invalid index value for dercheck_index was "
00260                << derchecker->index << endl;
00261           break;
00262         }
00263         break;
00264       }
00265     }
00266   }
00267 #endif
00268 
00269   // Do first reverse pass calculations
00270   int i;
00271   for (i=0;i<nvar;i++)
00272   {
00273     //double df=(pf->df)(xu,yu);
00274     px->u_bar[i]+=(pf->df)(xu,yu)*pz->u_bar[i];
00275   }
00276   for (i=0;i<nvar;i++)
00277   {
00278     px->u_bar[i]+=(pf->d2f)(xu,yu)*xdot[i]*pz->u_dot_bar[i];
00279 #if defined(ADDEBUG_PRINT)
00280     cout << px->u_bar[i] << " " << pz->u_dot_bar[i] << " " << addebug_count
00281          << endl;
00282 #endif
00283   }
00284   for (i=0;i<nvar;i++)
00285   {
00286     //px->u_dot_bar[i]+=(pf->df1)(*(px->u),*(py->u))*pz->u_dot_bar[i];
00287     px->u_dot_bar[i]+=(pf->df)(xu,yu)*pz->u_dot_bar[i];
00288 #if defined(ADDEBUG_PRINT)
00289     cout << px->u_dot_bar[i] << " " << addebug_count << endl;
00290     cout << pz->u_dot_bar[i] << " " << addebug_count << endl;
00291 #endif
00292   }
00293 
00294 
00295   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
00296   for (i=0;i<nvar;i++)
00297   {
00298     pz->u_bar[i]=0;
00299   }
00300   for (i=0;i<nvar;i++)
00301   {
00302     pz->u_dot_bar[i]=0;
00303   }
00304 
00305 #if defined(PRINT_DERS)
00306  print_derivatives(pz,"z");
00307  print_derivatives(px,"x");
00308 #endif
00309 }
00310 
00315 void read_pass2_2c(void)
00316 {
00317   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00318   // We are going forward for bptr and backword for bptr2
00319   //
00320   // list 1
00321   //
00322   int nvar=df1b2variable::nvar;
00323   test_smartlist & list=f1b2gradlist->list;
00324 
00325  int total_bytes=2*sizeof(df1b2_header)+sizeof(char*)
00326    +(nvar+2)*sizeof(double);
00327 // string identifier debug stuff
00328 #if defined(SAFE_ALL)
00329   char ids[]="BY";
00330   int slen=strlen(ids);
00331   total_bytes+=slen;
00332 #endif
00333   list.check_buffer_size(total_bytes);
00334 // end of string identifier debug stuff
00335 
00336   list.saveposition(); // save pointer to beginning of record;
00337   fixed_smartlist & nlist=f1b2gradlist->nlist;
00338    // nlist-=sizeof(int);
00339   // get record size
00340   int num_bytes=nlist.bptr->numbytes;
00341     // nlist+=nlist_record_size;
00342   //
00343   // list 2
00344   //
00345   test_smartlist & list2=f1b2gradlist->list2;
00346   fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00347   // get record size
00348   int num_bytes2=*nlist2.bptr;
00349   --nlist2;
00350   // backup the size of the record
00351   list2-=num_bytes2;
00352   list2.saveposition(); // save pointer to beginning of record;
00353   // save the pointer to the beginning of the record
00354   // bptr and bptr2 now both point to the beginning of their records
00355 
00356   double xu,yu;
00357   double * xdot;
00358   //df1b2_header x,z;
00359   df1b2function2c * pf;
00360 
00361   // get info from tape1
00362   // get info from tape1
00363 #if defined(SAFE_ARRAYS)
00364   checkidentiferstring("BY",list);
00365   checkidentiferstring("FW",list2);
00366 #endif
00367   df1b2_header * px=(df1b2_header *) list.bptr;
00368   list.bptr+=sizeof(df1b2_header);
00369   //double * py=(double *) list.bptr;
00370   list.bptr+=sizeof(double);
00371   df1b2_header * pz=(df1b2_header *) list.bptr;
00372   list.bptr+=sizeof(df1b2_header);
00373   pf=*(df1b2function2c **) list.bptr;
00374   list.bptr+=sizeof(char*);
00375   memcpy(&xu,list.bptr,sizeof(double));
00376   list.bptr+=sizeof(double);
00377   memcpy(&yu,list.bptr,sizeof(double));
00378   list.bptr+=sizeof(double);
00379   xdot=(double*)list.bptr;
00380   list.restoreposition(num_bytes); // save pointer to beginning of record;
00381 
00382   double * zbar;
00383   double * zdotbar;
00384 
00385   zbar=(double*)list2.bptr;
00386   zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00387   list2.restoreposition(); // save pointer to beginning of record;
00388 
00389   double * x_tilde=px->get_u_tilde();
00390   double * x_dot_tilde=px->get_u_dot_tilde();
00391   double * x_bar_tilde=px->get_u_bar_tilde();
00392   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00393   double * z_bar_tilde=pz->get_u_bar_tilde();
00394   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00395   // Do second "reverse-reverse" pass calculations
00396 #if defined(PRINT_DERS)
00397  print_derivatives(pz,"z");
00398  print_derivatives(px,"x");
00399 #endif
00400 
00401   int i;
00402 
00403   for (i=0;i<nvar;i++)
00404   {
00405     z_bar_tilde[i]=0;
00406     z_dot_bar_tilde[i]=0;
00407   }
00408 
00409   // start with x and add y
00410   for (i=0;i<nvar;i++)
00411   {
00412     *x_tilde+=(pf->d2f)(xu,yu)*zbar[i]*x_bar_tilde[i];
00413     z_bar_tilde[i]+=(pf->df)(xu,yu)*x_bar_tilde[i];
00414   }
00415 
00416   for (i=0;i<nvar;i++)
00417   {
00418     *x_tilde+=(pf->d2f)(xu,yu)*zdotbar[i]*x_dot_bar_tilde[i];
00419     z_dot_bar_tilde[i]+=(pf->df)(xu,yu)*x_dot_bar_tilde[i];
00420   }
00421 
00422   for (i=0;i<nvar;i++)
00423   {
00424     x_dot_tilde[i]+=(pf->d2f)(xu,yu)*zdotbar[i]*x_bar_tilde[i];
00425     z_dot_bar_tilde[i]+=(pf->d2f)(xu,yu)*xdot[i]*x_bar_tilde[i];
00426     *x_tilde+=(pf->d3f)(xu,yu)*xdot[i]*zdotbar[i]*x_bar_tilde[i];
00427   }
00428 
00429 #if defined(__DERCHECK__)
00430   if (derchecker->node_number)
00431   {
00432     if (derchecker->counter == derchecker->node_number)
00433     {
00434       if (derchecker->pass_number==1) // increment the variable of interest
00435       {
00436         switch(derchecker->vartype)
00437         {
00438         case 1:
00439           if (!derchecker->dotflag)
00440             derchecker->der_value=
00441               px->u_bar_tilde[derchecker->index-1];
00442           else
00443             derchecker->der_value=
00444               px->u_dot_bar_tilde[derchecker->index-1];
00445           break;
00446         case 2:
00447           if (!derchecker->dotflag)
00448             derchecker->der_value=
00449               py->u_bar_tilde[derchecker->index-1];
00450           else
00451             derchecker->der_value=
00452               py->u_dot_bar_tilde[derchecker->index-1];
00453           break;
00454         case 3:
00455           if (!derchecker->dotflag)
00456             derchecker->der_value=
00457               pz->u_bar_tilde[derchecker->index-1];
00458           else
00459             derchecker->der_value=
00460               pz->u_dot_bar_tilde[derchecker->index-1];
00461           break;
00462         default:
00463           cerr << "Invalid index value for dercheck_index was "
00464                << derchecker->index << endl;
00465         }
00466       }
00467     }
00468   }
00469 #endif
00470 #if defined(PRINT_DERS)
00471  print_derivatives(pz,"z");
00472  print_derivatives(px,"x");
00473 #endif
00474 }
00475 
00480 void read_pass2_3c(void)
00481 {
00482   // We are going backword for bptr and forward for bptr2
00483   // the current entry+2 in bptr is the size of the record i.e
00484   // points to the next record
00485   int nvar=df1b2variable::nvar;
00486   fixed_smartlist & nlist=f1b2gradlist->nlist;
00487   test_smartlist& list=f1b2gradlist->list;
00488    // nlist-=sizeof(int);
00489   // get record size
00490   int num_bytes=nlist.bptr->numbytes;
00491   // backup the size of the record
00492   list-=num_bytes;
00493   list.saveposition(); // save pointer to beginning of record;
00494   // save the pointer to the beginning of the record
00495   double xu;
00496   double yu;
00497   double * xdot;
00498   //df1b2_header x,z;
00499   df1b2function2c * pf;
00500 
00501   // get info from tape1
00502   // get info from tape1
00503 #if defined(SAFE_ARRAYS)
00504   checkidentiferstring("BY",list);
00505 #endif
00506   df1b2_header * px=(df1b2_header *) list.bptr;
00507   list.bptr+=sizeof(df1b2_header);
00508   //double * py=(double *) list.bptr;
00509   list.bptr+=sizeof(double);
00510   df1b2_header * pz=(df1b2_header *) list.bptr;
00511   list.bptr+=sizeof(df1b2_header);
00512   pf=*(df1b2function2c **) list.bptr;
00513   list.bptr+=sizeof(char*);
00514   memcpy(&xu,list.bptr,sizeof(double));
00515   list.bptr+=sizeof(double);
00516   memcpy(&yu,list.bptr,sizeof(double));
00517   list.bptr+=sizeof(double);
00518   xdot=(double*)list.bptr;
00519   list.restoreposition();
00520   int i;
00521 #if defined(PRINT_DERS)
00522  print_derivatives(pz,"z");
00523  print_derivatives(px,"x");
00524 #endif
00525 
00526   *(px->u_tilde)+=(pf->df)(xu,yu)* *(pz->u_tilde);
00527   for (i=0;i<nvar;i++)
00528   {
00529     *(px->u_tilde)+=(pf->d2f)(xu,yu)*xdot[i]*pz->u_dot_tilde[i];
00530   }
00531   for (i=0;i<nvar;i++)
00532   {
00533     px->u_dot_tilde[i]+=(pf->df)(xu,yu)*pz->u_dot_tilde[i];
00534   }
00535   *(pz->u_tilde)=0;
00536   for (i=0;i<nvar;i++)
00537   {
00538     pz->u_dot_tilde[i]=0;
00539   }
00540 #if defined(PRINT_DERS)
00541  print_derivatives(pz,"z");
00542  print_derivatives(px,"x");
00543 #endif
00544 }
00545 
00546 double AD_pow_1(double x,double y);
00547 double AD_pow_11(double x,double y);
00548 double AD_pow_111(double x,double y);
00549 
00550 static df1b2function2c AD_pow_c2(::pow,AD_pow_1,AD_pow_11,AD_pow_111,"pow_vc");
00551 
00552 /*
00553 df1b2variable pow(const df1b2variable& x,double y)
00554 {
00555   return AD_pow_c2(x,y);
00556 }
00557 */