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