ADMB Documentation  11.1.2382
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fn3.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fn3.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 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   int i;
00074   ostream * derout;
00075   derout=&cout;
00076   //*derout << derprintcount << " " << endl;
00077   if (bflag)
00078   {
00079     *derout << endl << "           --------------------------------- " << endl;
00080   }
00081   *derout << "pass " << df1b2variable::passnumber;
00082   *derout << "  variable " << s << "  address "
00083           << int(px->u) << endl;
00084   *derout << "u\t\t = " << *px->u << " ";
00085   *derout << endl;
00086 
00087   *derout << "udot\t\t = ";
00088   for (i=0;i<df1b2variable::nvar;i++)
00089     *derout <<  px->u_dot[i] << " ";
00090   *derout << endl;
00091 
00092   *derout << "u_bar\t\t = ";
00093   for (i=0;i<df1b2variable::nvar;i++)
00094     *derout <<  px->u_bar[i] << " ";
00095   *derout << endl;
00096 
00097   *derout << "u_dot_bar\t = ";
00098   for (i=0;i<df1b2variable::nvar;i++)
00099     *derout <<  px->u_dot_bar[i] << " ";
00100   *derout << endl;
00101 
00102   if (df1b2variable::passnumber>1)
00103   {
00104     *derout << "u_tilde\t\t = " << *px->u_tilde << " ";
00105     *derout << endl;
00106 
00107     *derout << "u_dot_tilde\t = ";
00108     for (i=0;i<df1b2variable::nvar;i++)
00109       *derout <<  px->u_dot_tilde[i] << " ";
00110     *derout << endl;
00111 
00112     *derout << "u_bar_tilde\t = ";
00113     for (i=0;i<df1b2variable::nvar;i++)
00114       *derout <<  px->u_bar_tilde[i] << " ";
00115     *derout << endl;
00116 
00117     *derout << "u_dot_bar_tilde\t = ";
00118     for (i=0;i<df1b2variable::nvar;i++)
00119       *derout <<  px->u_dot_bar_tilde[i] << " ";
00120     *derout << endl;
00121   }
00122   *derout << endl;
00123 }
00124 #endif
00125 
00130 df1b2variable& df1b2variable::operator += (const df1b2variable& _x)
00131 {
00132   ADUNCONST(df1b2variable,x)
00133   double * xd=x.get_u_dot();
00134   double * zd=get_u_dot();
00135   *get_u()+=*x.get_u();
00136   for (int i=0;i<df1b2variable::nvar;i++)
00137   {
00138     *zd++ += *xd++;
00139   }
00140 
00141   // WRITE WHATEVER ON TAPE
00142   //df1b2tape->set_tapeinfo_header(&x,&z,this,xd);
00143   // save stuff for first reverse pass
00144   // need &x, &z, this,
00145   if (!df1b2_gradlist::no_derivatives)
00146     f1b2gradlist->write_pass1_pluseq(&x,this);
00147   return *this;
00148 }
00149 void ad_read_pass1_plus_eq(void);
00150 
00155 int df1b2_gradlist::write_pass1_pluseq(const df1b2variable * _px,
00156   df1b2variable * pz)
00157 {
00158   ncount++;
00159 #if defined(CHECK_COUNT)
00160   if (ncount >= ncount_check)
00161     ncount_checker(ncount,ncount_check);
00162 #endif
00163   //int nvar=df1b2variable::nvar;
00164   ADUNCONST(df1b2variable*,px)
00165   fixed_smartlist & nlist=f1b2gradlist->nlist;
00166   test_smartlist& list=f1b2gradlist->list;
00167 
00168   int total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header);
00169 #if defined(SAFE_ALL)
00170   char ids[]="JK";
00171   int slen=strlen(ids);
00172   total_bytes+=slen;
00173 #endif
00174   list.check_buffer_size(total_bytes);
00175   void * tmpptr=list.bptr;
00176 #if defined(SAFE_ALL)
00177   memcpy(list,ids,slen);
00178 #endif
00179 
00180   memcpy(list,(df1b2_header*)(px),sizeof(df1b2_header));
00181   memcpy(list,(df1b2_header*)(pz),sizeof(df1b2_header));
00182 
00183   // ***** write  record size
00184   nlist.bptr->numbytes=adptr_diff(list.bptr,tmpptr);
00185   nlist.bptr->pf=(ADrfptr)(&ad_read_pass1_plus_eq);
00186   ++nlist;
00187   return 0;
00188 }
00189 
00194 void ad_read_pass1_plus_eq(void)
00195 {
00196   switch(df1b2variable::passnumber)
00197   {
00198   case 1:
00199     read_pass1_plus_eq_1();
00200     break;
00201   case 2:
00202     read_pass1_plus_eq_2();
00203     break;
00204   case 3:
00205     read_pass1_plus_eq_3();
00206     break;
00207   default:
00208     cerr << "illegal value for df1b2variable::pass = "
00209          << df1b2variable::passnumber << endl;
00210     exit(1);
00211   }
00212 }
00213 
00218 void read_pass1_plus_eq_1(void)
00219 {
00220   // We are going backword for bptr and forward for bptr2
00221   // the current entry+2 in bptr is the size of the record i.e
00222   // points to the next record
00223   int nvar=df1b2variable::nvar;
00224   fixed_smartlist & nlist=f1b2gradlist->nlist;
00225   test_smartlist& list=f1b2gradlist->list;
00226    // nlist-=sizeof(int);
00227   // get record size
00228   int num_bytes=nlist.bptr->numbytes;
00229   // backup the size of the record
00230   list-=num_bytes;
00231   list.saveposition(); // save pointer to beginning of record;
00232   // save the pointer to the beginning of the record
00233 #if defined(SAFE_ALL)
00234   checkidentiferstring("JK",list);
00235 #endif
00236 
00237   // get info from tape1
00238   df1b2_header * px=(df1b2_header *) list.bptr;
00239   list.bptr+=sizeof(df1b2_header);
00240   df1b2_header * pz=(df1b2_header *) list.bptr;
00241 
00242   list.restoreposition(); // save pointer to beginning of record;
00243   int i;
00244 
00245   // Do first reverse paSS calculations
00246   // ****************************************************************
00247   // turn this off if no third derivatives are calculated
00248   // if (!no_third_derivatives)
00249   // {
00250   // save for second reverse pass
00251   // save identifier 1
00252   //   fixed_smartlist2& nlist2=f1b2gradlist->nlist2;
00253   //   test_smartlist& list2=f1b2gradlist->list2;
00254   //int total_bytes=2*nvar*sizeof(double);
00255 // string identifier debug stuff
00256 #if defined(SAFE_ALL)
00257   //char ids[]="IL";
00258   //int slen=strlen(ids);
00259   //total_bytes+=slen;
00260 #endif
00261   //list2.check_buffer_size(total_bytes);
00262   //void * tmpptr2=list2.bptr;
00263 #if defined(SAFE_ALL)
00264   //memcpy(list2,ids,slen);
00265 #endif
00266      //memcpy(list2,pz->get_u_bar(),nvar*sizeof(double));
00267      //memcpy(list2,pz->get_u_dot_bar(),nvar*sizeof(double));
00268      //*nlist2.bptr=adptr_diff(list2.bptr,tmpptr2);
00269      //nlist2++;
00270   // }
00271   //
00272   // ****************************************************************
00273 #if defined(PRINT_DERS)
00274  print_derivatives(px,"x");
00275  print_derivatives(pz,"z");
00276 #endif
00277 
00278   for (i=0;i<nvar;i++)
00279   {
00280     px->u_bar[i]+=pz->u_bar[i];
00281   }
00282   for (i=0;i<nvar;i++)
00283   {
00284     px->u_dot_bar[i]+=pz->u_dot_bar[i];
00285   }
00286 #if defined(PRINT_DERS)
00287  print_derivatives(px,"x");
00288  print_derivatives(pz,"z");
00289 #endif
00290 }
00291 
00296 void read_pass1_plus_eq_2(void)
00297 {
00298   //const int nlist_record_size=sizeof(int)+sizeof(char*);
00299   // We are going forward for bptr and backword for bptr2
00300   //
00301   // list 1
00302   //
00303   int nvar=df1b2variable::nvar;
00304   test_smartlist & list=f1b2gradlist->list;
00305 
00306   int total_bytes=sizeof(df1b2_header)+sizeof(df1b2_header);
00307 #if defined(SAFE_ALL)
00308   char ids[]="JK";
00309   int slen=strlen(ids);
00310   total_bytes+=slen;
00311 #endif
00312   list.check_buffer_size(total_bytes);
00313 
00314   list.saveposition(); // save pointer to beginning of record;
00315   fixed_smartlist & nlist=f1b2gradlist->nlist;
00316    // nlist-=sizeof(int);
00317   // get record size
00318   int num_bytes=nlist.bptr->numbytes;
00319     // nlist+=nlist_record_size;
00320   //
00321   // list 2
00322   //
00323   //test_smartlist & list2=f1b2gradlist->list2;
00324   //fixed_smartlist2 & nlist2=f1b2gradlist->nlist2;
00325   // get record size
00326   //int num_bytes2=*nlist2.bptr;
00327   //nlist2--;
00328   // backup the size of the record
00329   //list2-=num_bytes2;
00330   //list2.saveposition(); // save pointer to beginning of record;
00331   // save the pointer to the beginning of the record
00332   // bptr and bptr2 now both point to the beginning of their records
00333 #if defined(SAFE_ALL)
00334   checkidentiferstring("JK",list);
00335   //checkidentiferstring("IL",list2);
00336 #endif
00337 
00338 
00339   // get info from tape1
00340   df1b2_header * px=(df1b2_header *) list.bptr;
00341   list.bptr+=sizeof(df1b2_header);
00342   df1b2_header * pz=(df1b2_header *) list.bptr;
00343 
00344   list.restoreposition(num_bytes); // save pointer to beginning of record;
00345 
00346 
00347   //double * zbar=(double*)list2.bptr;
00348   //double * zdotbar=(double*)(list2.bptr+nvar*sizeof(double));
00349 
00350   double * x_bar_tilde=px->get_u_bar_tilde();
00351   double * x_dot_bar_tilde=px->get_u_dot_bar_tilde();
00352   double * z_bar_tilde=pz->get_u_bar_tilde();
00353   double * z_dot_bar_tilde=pz->get_u_dot_bar_tilde();
00354   // Do second "reverse-reverse" pass calculations
00355   int i;
00356   for (i=0;i<nvar;i++)
00357   {
00358     z_bar_tilde[i]+=x_bar_tilde[i];
00359   }
00360 
00361   for (i=0;i<nvar;i++)
00362   {
00363     z_dot_bar_tilde[i]+=x_dot_bar_tilde[i];
00364   }
00365   //list2.restoreposition(); // save pointer to beginning of record;
00366 #if defined(PRINT_DERS)
00367  print_derivatives(px,"x");
00368  print_derivatives(pz,"z");
00369 #endif
00370 }
00371 
00376 void read_pass1_plus_eq_3(void)
00377 {
00378   // We are going backword for bptr and forward for bptr2
00379   // the current entry+2 in bptr is the size of the record i.e
00380   // points to the next record
00381   int nvar=df1b2variable::nvar;
00382   fixed_smartlist & nlist=f1b2gradlist->nlist;
00383   test_smartlist& list=f1b2gradlist->list;
00384    // nlist-=sizeof(int);
00385   // get record size
00386   int num_bytes=nlist.bptr->numbytes;
00387   // backup the size of the record
00388   list-=num_bytes;
00389   list.saveposition(); // save pointer to beginning of record;
00390   // save the pointer to the beginning of the record
00391 
00392 #if defined(SAFE_ALL)
00393   checkidentiferstring("JK",list);
00394 #endif
00395   // get info from tape1
00396   df1b2_header * px=(df1b2_header *) list.bptr;
00397   list.bptr+=sizeof(df1b2_header);
00398   df1b2_header * pz=(df1b2_header *) list.bptr;
00399 
00400   list.restoreposition(); // save pointer to beginning of record;
00401 
00402   int i;
00403   *(px->u_tilde)+=*pz->u_tilde;
00404   for (i=0;i<nvar;i++)
00405   {
00406     px->u_dot_tilde[i]+=pz->u_dot_tilde[i];
00407   }
00408 #if defined(PRINT_DERS)
00409  print_derivatives(px,"x");
00410  print_derivatives(pz,"z");
00411 #endif
00412 }
00413 
00418 df1b2variable fabs(const df1b2variable& x)
00419 {
00420   if (value(x)>=0.0)
00421     return x;
00422   else
00423     return -x;
00424 }
00425 
00430 df1b2vector fabs(const df1b2vector& t1)
00431 {
00432    df1b2vector tmp(t1.indexmin(),t1.indexmax());
00433    for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
00434    {
00435      tmp(i)=fabs(t1(i));
00436    }
00437 
00438    return(tmp);
00439 }
00440 
00445 df1b2variable max(const df1b2vector& t1)
00446 {
00447    df1b2variable tmp;
00448    int mmin=t1.indexmin();
00449    int mmax=t1.indexmax();
00450    tmp=t1(mmin);
00451    for (int i=mmin+1; i<=mmax; i++)
00452    {
00453      if (value(tmp)<value(t1(i))) tmp=t1(i);
00454    }
00455    return(tmp);
00456 }