ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fn3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fn3.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  */
00011 #include <df1b2fun.h>
00012 void read_pass1_plus_eq_1(void);
00013 void read_pass1_plus_eq_2(void);
00014 void read_pass1_plus_eq_3(void);
00015 //#define ADDEBUG_PRINT
00016 #if defined(ADDEBUG_PRINT)
00017   extern int addebug_count;
00018 #endif
00019 //#define PRINT_DERS
00020 
00021 
00022 #if defined(PRINT_DERS)
00023 
00027 void print_derivatives(const adstring& s, double f, double df,
00028   double d2f,double d3f,int bflag)
00029 {
00030   ostream * derout;
00031   derout=&cout;
00032   if (bflag)
00033   {
00034     *derout << "           ---------------------------------------- " << endl;
00035   }
00036   *derout << "Function: " << s << " " << "f = " << f
00037           << " df = " << df << " d2f = " << d2f << " d3f = " << d3f << endl;
00038 }
00039 
00044 void print_derivatives(const adstring& s, double f, double df1,
00045   double df2,double df11,double df12, double df22,
00046   double df111, double df112, double df122, double df222,int bflag)
00047 {
00048   ostream * derout;
00049   derout=&cout;
00050   if (bflag)
00051   {
00052     *derout << endl << "           --------------------------------- " << endl;
00053   }
00054   *derout << "Function: " << s << " " << "f = " << f
00055           << " df1 = " << df1
00056           << " df2 = " << df2
00057           << " df11 = " << df11
00058           << " df12 = " << df12
00059           << " df22 = " << df22
00060           << " df111 = " << df111
00061           << " df112 = " << df112
00062           << " df122 = " << df122
00063           << " df222 = " << df222 << endl;
00064 }
00065 
00070 void print_derivatives(df1b2_header * px,const char * s,
00071   int bflag)
00072 {
00073   ostream * derout;
00074   derout=&cout;
00075   //*derout << derprintcount << " " << endl;
00076   if (bflag)
00077   {
00078     *derout << endl << "           --------------------------------- " << endl;
00079   }
00080   *derout << "pass " << df1b2variable::passnumber;
00081   *derout << "  variable " << s << "  address "
00082           << int(px->u) << endl;
00083   *derout << "u\t\t = " << *px->u << " ";
00084   *derout << endl;
00085 
00086   *derout << "udot\t\t = ";
00087   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00088     *derout <<  px->u_dot[i] << " ";
00089   *derout << endl;
00090 
00091   *derout << "u_bar\t\t = ";
00092   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00093     *derout <<  px->u_bar[i] << " ";
00094   *derout << endl;
00095 
00096   *derout << "u_dot_bar\t = ";
00097   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00098     *derout <<  px->u_dot_bar[i] << " ";
00099   *derout << endl;
00100 
00101   if (df1b2variable::passnumber>1)
00102   {
00103     *derout << "u_tilde\t\t = " << *px->u_tilde << " ";
00104     *derout << endl;
00105 
00106     *derout << "u_dot_tilde\t = ";
00107     for (unsigned int i=0;i<df1b2variable::nvar;i++)
00108       *derout <<  px->u_dot_tilde[i] << " ";
00109     *derout << endl;
00110 
00111     *derout << "u_bar_tilde\t = ";
00112     for (unsigned int i=0;i<df1b2variable::nvar;i++)
00113       *derout <<  px->u_bar_tilde[i] << " ";
00114     *derout << endl;
00115 
00116     *derout << "u_dot_bar_tilde\t = ";
00117     for (unsigned int i=0;i<df1b2variable::nvar;i++)
00118       *derout <<  px->u_dot_bar_tilde[i] << " ";
00119     *derout << endl;
00120   }
00121   *derout << endl;
00122 }
00123 #endif
00124 
00129 df1b2variable& df1b2variable::operator += (const df1b2variable& _x)
00130 {
00131   ADUNCONST(df1b2variable,x)
00132   double * xd=x.get_u_dot();
00133   double * zd=get_u_dot();
00134   *get_u()+=*x.get_u();
00135   for (unsigned int i=0;i<df1b2variable::nvar;i++)
00136   {
00137     *zd++ += *xd++;
00138   }
00139 
00140   // WRITE WHATEVER ON TAPE
00141   //df1b2tape->set_tapeinfo_header(&x,&z,this,xd);
00142   // save stuff for first reverse pass
00143   // need &x, &z, this,
00144   if (!df1b2_gradlist::no_derivatives)
00145     f1b2gradlist->write_pass1_pluseq(&x,this);
00146   return *this;
00147 }
00148 void ad_read_pass1_plus_eq(void);
00149 
00154 int df1b2_gradlist::write_pass1_pluseq(const df1b2variable * _px,
00155   df1b2variable * pz)
00156 {
00157   ncount++;
00158 #if defined(CHECK_COUNT)
00159   if (ncount >= ncount_check)
00160     ncount_checker(ncount,ncount_check);
00161 #endif
00162   //int nvar=df1b2variable::nvar;
00163   ADUNCONST(df1b2variable*,px)
00164   fixed_smartlist & nlist=f1b2gradlist->nlist;
00165   test_smartlist& list=f1b2gradlist->list;
00166 
00167   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header);
00168 #if defined(SAFE_ALL)
00169   char ids[]="JK";
00170   int slen=strlen(ids);
00171   total_bytes+=slen;
00172 #endif
00173   list.check_buffer_size(total_bytes);
00174   void * tmpptr=list.bptr;
00175 #if defined(SAFE_ALL)
00176   memcpy(list,ids,slen);
00177 #endif
00178 
00179   memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00180   memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00181 
00182   // ***** write  record size
00183   nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00184   nlist.bptr->pf=(ADrfptr)(&ad_read_pass1_plus_eq);
00185   ++nlist;
00186   return 0;
00187 }
00188 
00193 void ad_read_pass1_plus_eq(void)
00194 {
00195   switch(df1b2variable::passnumber)
00196   {
00197   case 1:
00198     read_pass1_plus_eq_1();
00199     break;
00200   case 2:
00201     read_pass1_plus_eq_2();
00202     break;
00203   case 3:
00204     read_pass1_plus_eq_3();
00205     break;
00206   default:
00207     cerr << "illegal value for df1b2variable::pass = "
00208          << df1b2variable::passnumber << endl;
00209     exit(1);
00210   }
00211 }
00212 
00217 void read_pass1_plus_eq_1(void)
00218 {
00219   // We are going backword for bptr and forward for bptr2
00220   // the current entry+2 in bptr is the size of the record i.e
00221   // points to the next record
00222   unsigned int nvar=df1b2variable::nvar;
00223   fixed_smartlist & nlist=f1b2gradlist->nlist;
00224   test_smartlist& list=f1b2gradlist->list;
00225    // nlist-=sizeof(int);
00226   // get record size
00227   int num_bytes=nlist.bptr->numbytes;
00228   // backup the size of the record
00229   list-=num_bytes;
00230   list.saveposition(); // save pointer to beginning of record;
00231   // save the pointer to the beginning of the record
00232 #if defined(SAFE_ALL)
00233   checkidentiferstring("JK",list);
00234 #endif
00235 
00236   // get info from tape1
00237   df1b2_header * px=(df1b2_header *) list.bptr;
00238   list.bptr+=sizeof(df1b2_header);
00239   df1b2_header * pz=(df1b2_header *) list.bptr;
00240 
00241   list.restoreposition(); // save pointer to beginning of record;
00242 
00243   // Do first reverse paSS calculations
00244   // ****************************************************************
00245   // turn this off if no third derivatives are calculated
00246   // if (!no_third_derivatives)
00247   // {
00248   // save for second reverse pass
00249   // save identifier 1
00250   //   fixed_smartlist2& nlist2=f1b2gradlist->nlist2;
00251   //   test_smartlist& list2=f1b2gradlist->list2;
00252   //int total_bytes=2*nvar*sizeof(double);
00253 // string identifier debug stuff
00254 #if defined(SAFE_ALL)
00255   //char ids[]="IL";
00256   //int slen=strlen(ids);
00257   //total_bytes+=slen;
00258 #endif
00259   //list2.check_buffer_size(total_bytes);
00260   //void * tmpptr2=list2.bptr;
00261 #if defined(SAFE_ALL)
00262   //memcpy(list2,ids,slen);
00263 #endif
00264      //memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00265      //memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00266      //*nlist2.bptr=adptr_diff(list2.bptr,tmpptr2);
00267      //nlist2++;
00268   // }
00269   //
00270   // ****************************************************************
00271 #if defined(PRINT_DERS)
00272  print_derivatives(px,"x");
00273  print_derivatives(pz,"z");
00274 #endif
00275 
00276   for (unsigned int i=0;i<nvar;i++)
00277   {
00278     px->u_bar[i]+=pz->u_bar[i];
00279   }
00280   for (unsigned int i=0;i<nvar;i++)
00281   {
00282     px->u_dot_bar[i]+=pz->u_dot_bar[i];
00283   }
00284 #if defined(PRINT_DERS)
00285  print_derivatives(px,"x");
00286  print_derivatives(pz,"z");
00287 #endif
00288 }
00289 
00294 void read_pass1_plus_eq_2(void)
00295 {
00296   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00297   // We are going forward for bptr and backword for bptr2
00298   //
00299   // list 1
00300   //
00301   unsigned int nvar=df1b2variable::nvar;
00302   test_smartlist & list=f1b2gradlist->list;
00303 
00304   size_t total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header);
00305 #if defined(SAFE_ALL)
00306   char ids[]="JK";
00307   int slen=strlen(ids);
00308   total_bytes+=slen;
00309 #endif
00310   list.check_buffer_size(total_bytes);
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 #if defined(SAFE_ALL)
00332   checkidentiferstring("JK",list);
00333   //checkidentiferstring("IL",list2);
00334 #endif
00335 
00336   // get info from tape1
00337   df1b2_header * px=(df1b2_header *) list.bptr;
00338   list.bptr+=sizeof(df1b2_header);
00339   df1b2_header * pz=(df1b2_header *) list.bptr;
00340 
00341   list.restoreposition(num_bytes); // save pointer to beginning of record;
00342 
00343 
00344   //double * zbar=(double*)list2.bptr;
00345   //double * zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00346 
00347   double * x_bar_tilde=px->get_u_bar_tilde();
00348   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00349   double * z_bar_tilde=pz->get_u_bar_tilde();
00350   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00351   // Do second "reverse-reverse" pass calculations
00352   for (unsigned int i=0;i<nvar;i++)
00353   {
00354     z_bar_tilde[i]+=x_bar_tilde[i];
00355   }
00356 
00357   for (unsigned int i=0;i<nvar;i++)
00358   {
00359     z_dot_bar_tilde[i]+=x_dot_bar_tilde[i];
00360   }
00361   //list2.restoreposition(); // save pointer to beginning of record;
00362 #if defined(PRINT_DERS)
00363  print_derivatives(px,"x");
00364  print_derivatives(pz,"z");
00365 #endif
00366 }
00367 
00372 void read_pass1_plus_eq_3(void)
00373 {
00374   // We are going backword for bptr and forward for bptr2
00375   // the current entry+2 in bptr is the size of the record i.e
00376   // points to the next record
00377   unsigned int nvar=df1b2variable::nvar;
00378   fixed_smartlist & nlist=f1b2gradlist->nlist;
00379   test_smartlist& list=f1b2gradlist->list;
00380    // nlist-=sizeof(int);
00381   // get record size
00382   int num_bytes=nlist.bptr->numbytes;
00383   // backup the size of the record
00384   list-=num_bytes;
00385   list.saveposition(); // save pointer to beginning of record;
00386   // save the pointer to the beginning of the record
00387 
00388 #if defined(SAFE_ALL)
00389   checkidentiferstring("JK",list);
00390 #endif
00391   // get info from tape1
00392   df1b2_header * px=(df1b2_header *) list.bptr;
00393   list.bptr+=sizeof(df1b2_header);
00394   df1b2_header * pz=(df1b2_header *) list.bptr;
00395 
00396   list.restoreposition(); // save pointer to beginning of record;
00397 
00398   *(px->u_tilde)+=*pz->u_tilde;
00399   for (unsigned int i=0;i<nvar;i++)
00400   {
00401     px->u_dot_tilde[i]+=pz->u_dot_tilde[i];
00402   }
00403 #if defined(PRINT_DERS)
00404  print_derivatives(px,"x");
00405  print_derivatives(pz,"z");
00406 #endif
00407 }
00408 
00413 df1b2variable fabs(const df1b2variable& x)
00414 {
00415   if (value(x)>=0.0)
00416     return x;
00417   else
00418     return -x;
00419 }
00420 
00425 df1b2vector fabs(const df1b2vector& t1)
00426 {
00427    df1b2vector tmp(t1.indexmin(),t1.indexmax());
00428    for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00429    {
00430      tmp(i)=fabs(t1(i));
00431    }
00432 
00433    return(tmp);
00434 }
00435 
00440 df1b2variable max(const df1b2vector& t1)
00441 {
00442    df1b2variable tmp;
00443    int mmin=t1.indexmin();
00444    int mmax=t1.indexmax();
00445    tmp=t1(mmin);
00446    for (int i=mmin+1; i<=mmax; i++)
00447    {
00448      if (value(tmp)<value(t1(i))) tmp=t1(i);
00449    }
00450    return(tmp);
00451 }