ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fun.cpp 2613 2014-11-10 21:23:04Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00013 extern "C"
00014 {
00015   int kbhit(void) { return 0;}
00016 }
00017 
00018 class df1b2_gradlist;
00019 //df1b2_gradlist * f1b2gradlist = NULL;
00020 
00021 //int df1b2variable::noallocate=0;
00022 
00023 //initial_df1b2params ** initial_df1b2params::varsptr
00024 //  =new P_INITIAL_DF1B2PARAMS[1000];
00025 //int initial_df1b2params::num_initial_df1b2params=0;         // array
00026 
00027 //int initial_df1b2params::num_initial_df1b2params_sav=0;         // array
00028 //initial_df1b2params ** initial_df1b2params::varsptr_sav=0;
00029 
00030 //int initial_df1b2params::current_phase=0;
00031 
00032 #if defined(__DERCHECK__)
00033 int mydercheckercounter=0;
00034 #endif
00035 
00040 void df1b2variable::save_adpool_pointer(void)
00041 {
00042   if (adpool_stack_pointer> adpool_stack_size-1)
00043   {
00044     cerr << "overflow in save_adpool_pointer" << endl;
00045     ad_exit(1);
00046   }
00047   adpool_stack[adpool_stack_pointer]=pool;
00048   adpool_nvar_stack[adpool_stack_pointer++]=nvar;
00049 }
00050 
00055 void df1b2variable::restore_adpool_pointer(void)
00056 {
00057   if (adpool_stack_pointer<=0)
00058   {
00059     cerr << "underflow in save_adpool_pointer" << endl;
00060     ad_exit(1);
00061   }
00062   pool=adpool_stack[--adpool_stack_pointer];
00063   nvar=adpool_nvar_stack[adpool_stack_pointer];
00064 }
00065 
00070 void initial_df1b2params::save_varsptr(void)
00071 {
00072   initial_df1b2params::varsptr_sav=initial_df1b2params::varsptr;
00073 
00074   varsptr=new P_INITIAL_DF1B2PARAMS[1000];
00075   if (varsptr == 0)
00076   {
00077     cerr << "error allocating memory for "
00078       "initial_df1b2params::varsptr  " << endl;
00079     ad_exit(1);
00080   }
00081 
00082   num_initial_df1b2params_sav=num_initial_df1b2params;
00083   num_initial_df1b2params =0;
00084 }
00085 
00090 void initial_df1b2params::restore_varsptr(void)
00091 {
00092   if (num_initial_df1b2params == 0)
00093   {
00094     if (varsptr)
00095     {
00096       delete [] varsptr;
00097       varsptr=0;
00098     }
00099     varsptr=initial_df1b2params::varsptr_sav;
00100     varsptr_sav=0;
00101 
00102     num_initial_df1b2params= num_initial_df1b2params_sav;
00103     num_initial_df1b2params_sav=0;
00104   }
00105   else
00106   {
00107     if (num_initial_df1b2params+num_initial_df1b2params_sav
00108      > 1000)
00109     {
00110       cerr << "maximum numver of iitial_df1b2params is 1000 "
00111            << endl;
00112       ad_exit(1);
00113     }
00114 
00115     for (int i=0;i<num_initial_df1b2params_sav;i++)
00116     {
00117       varsptr[i+num_initial_df1b2params]=varsptr_sav[i];
00118     }
00119     num_initial_df1b2params+=num_initial_df1b2params_sav;
00120     delete varsptr_sav;
00121     varsptr_sav=0;
00122     num_initial_df1b2params_sav=0;
00123   }
00124 }
00125 
00130 initial_df1b2params::initial_df1b2params(void) : ind_index(0)
00131 {
00132   scalefactor=0.0;
00133   phase_start=1;
00134   phase_save=1;
00135   add_to_list();
00136 }
00137 
00138 typedef void (**ADprfptr)(void);
00139 typedef void (*ADrfptr)(void);
00140 
00145 void df1b2_gradcalc1(void)
00146 {
00147   //smartlist & list=f1b2gradlist->list;
00148   fixed_smartlist & nlist=f1b2gradlist->nlist;
00149   int ncount=f1b2gradlist->ncount;
00150   //ADrfptr rf2;
00151   int xcount=0;
00152   int tmpcount;
00153   int tcount=f1b2gradlist->ncount;
00154 
00155   //int check_pool_flag3=0;
00156   switch (df1b2variable::passnumber)
00157   {
00158   case 1:
00159     f1b2gradlist->list.save_end();
00160     f1b2gradlist->nlist.save_end();
00161   case 3:
00162 #if defined(__DERCHECK__)
00163     //  derchecker->counter=f1b2gradlist->ncount;
00164     mydercheckercounter=f1b2gradlist->ncount;
00165 #endif
00166     f1b2gradlist->list.set_reverse();
00167     f1b2gradlist->list.restore_end();
00168     f1b2gradlist->nlist.set_reverse();
00169     f1b2gradlist->nlist.restore_end();
00170 
00171     if (df1b2variable::passnumber==3)
00172     {
00173       f1b2gradlist->nlist3.save_end();
00174       f1b2gradlist->nlist3.set_reverse();
00175       f1b2gradlist->nlist3.restore_end();
00176       f1b2gradlist->list3.save_end();
00177       f1b2gradlist->list3.set_reverse();
00178       f1b2gradlist->list3.restore_end();
00179     }
00180     tmpcount=ncount;
00181     do
00182     {
00183       tmpcount--;
00184      /*
00185       if (!(tmpcount %100))
00186       {
00187         cout << "B " << tmpcount << endl;
00188       }
00189       if (tmpcount < 2)
00190       {
00191         cout << "B " << tmpcount << endl;
00192       }
00193       */
00194       //nlist-=sizeof(char*);
00195       --nlist;
00196       if (nlist.eof())
00197         break;
00198       ADrfptr rf=nlist.bptr->pf;
00199       (*rf)();
00200      /*
00201       if (check_pool_flag3)
00202       {
00203         if (xcount > 722)
00204         {
00205           cout << xcount << "  ";
00206           //check_pool_depths();
00207         }
00208       }
00209      */
00210 #if defined(__DERCHECK__)
00211         //derchecker->counter--;
00212       mydercheckercounter--;
00213 #endif
00214         xcount++;
00215           tcount--;
00216        if (xcount > 99999999) cout << xcount << endl;
00217        //if (tcount == 6599 )
00218         //  cout << tcount << endl;
00219       /*
00220         if (initial_df1b2params::current_phase==2)
00221         {
00222           tcount--;
00223           if (rf2 != (nlist.buffer+122)->pf || xcount > 12488)
00224              cout << tcount << " " << xcount << " "
00225                   << (nlist.buffer+122)->pf << endl;
00226         }
00227        */
00228     }
00229     while(nlist.bptr>=nlist.buffer);
00230     break;
00231   case 2:
00232     {
00233       f1b2gradlist->list2.save_end();
00234       f1b2gradlist->list2.set_reverse();
00235       f1b2gradlist->list2.restore_end();
00236       f1b2gradlist->nlist2.save_end();
00237       f1b2gradlist->nlist2.set_reverse();
00238       f1b2gradlist->nlist2.restore_end();
00239 #if defined(__DERCHECK__)
00240       // derchecker->counter=0;
00241       mydercheckercounter=0;
00242 #endif
00243       f1b2gradlist->list.set_forward();
00244       f1b2gradlist->list.rewind();
00245       f1b2gradlist->nlist.set_forward();
00246       f1b2gradlist->nlist.rewind();
00247 
00248       //nlist.reset();  // set pointer to beginning of list
00249       --(f1b2gradlist->nlist2);  // backup one record
00250       int icount=0;
00251       do
00252       {
00253 #if defined(__DERCHECK__)
00254         // derchecker->counter++;
00255         mydercheckercounter++;
00256 #endif
00257         //ADrfptr rf=*ADprfptr(nlist.bptr);
00258         ADrfptr rf=nlist.bptr->pf;
00259         (*rf)();
00260         ++nlist;
00261 
00262        /*
00263         if (initial_df1b2params::current_phase==2)
00264         {
00265           if (icount >=4579)
00266           {
00267             cout << icount << endl;
00268           }
00269         }
00270 
00271         if (icount>28115)
00272         {
00273           cout << "icount = " << icount << endl;
00274         }
00275         */
00276       }
00277       while(++icount<ncount);
00278       // need this to get pointer right for step 3?
00279        // nlist-=sizeof(int); // set nlist to point to second record;
00280       //nlist++;
00281 
00282       break;
00283     }
00284 
00285   default:
00286      cerr << "illega value for df1b2variable::passnumber "
00287       " value = " << df1b2variable::passnumber << endl;
00288   }
00289 }
00290 
00291 
00292 
00293 double nsin(double x) {return -sin(x);}
00294 double ncos(double x) {return -cos(x);}
00295 
00296 double ADmult_add_fun(double x,double y);
00297 df1b2function1 ADf1b2_sin(::sin,::cos,::nsin,::ncos,"sin");
00298 df1b2function1 ADf1b2_cos(::cos,::nsin,::ncos,::sin);
00299 df1b2function1 ADf1b2_exp(::exp,::exp,::exp,::exp,"exp");
00300 
00305 double AD_df1_atan(double x)
00306 {
00307   return double(1.0)/(1+square(x));
00308 }
00309 
00314 double AD_df2_atan(double x)
00315 {
00316   return double(-2.0)*x/square(1+square(x));
00317 }
00318 
00323 double AD_df1_tan(double x)
00324 {
00325   return double(1.0)+square(tan(x));
00326 }
00327 
00332 double AD_df2_tan(double x)
00333 {
00334   double y=tan(x);
00335   return double(2.0)*y*(double(1.0)+square(y));
00336 }
00337 
00342 double AD_df3_atan(double x)
00343 {
00344   return double(-2.0)/square(double(1)+square(x))
00345          + double(12.0)*square(x)/cube(double(1)+square(x));
00346 }
00347 
00352 double AD_df3_tan(double x)
00353 {
00354   double y=square(tan(x));
00355 
00356   return double(2.0) * (double(1.0)+double(3.0)*y) * (double(1) + y);
00357 }
00358 
00359 df1b2function1 ADf1b2_tan(::tan,::AD_df1_tan,::AD_df2_tan ,::AD_df3_tan,"tan");
00360 
00361 df1b2function1 ADf1b2_atan(::atan,::AD_df1_atan,::AD_df2_atan ,::AD_df3_atan,
00362   "atan");
00363 
00368 double AD_arg_inv(double x)
00369 {
00370   return double(1.0)/x;
00371 }
00372 
00377 double AD_minus_arg_inv2(double x)
00378 {
00379   return double(-1.0)/(x*x);
00380 }
00381 
00386 double AD_2arg_inv3(double x)
00387 {
00388   return double(2.0)/(x*x*x);
00389 }
00390 
00395 double AD_minus6_arg_inv4(double x)
00396 {
00397   return double(-6.0)/(x*x*x*x);
00398 }
00399 
00400 df1b2function1 ADf1b2_inv(AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3,
00401   AD_minus6_arg_inv4);
00402 
00403 df1b2function1 ADf1b2_log(::log,AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3);
00404 
00405 /*
00406 df1b2variable sin(const df1b2variable& x)
00407 {
00408   return ADf1b2_sin(x);
00409 }
00410 */
00411 
00416 df1b2variable atan(const df1b2variable& x)
00417 {
00418   return ADf1b2_atan(x);
00419 }
00420 
00425 df1b2variable tan(const df1b2variable& x)
00426 {
00427   return ADf1b2_tan(x);
00428 }
00429 
00434 df1b2variable& df1b2variable::operator *= (const df1b2variable& v)
00435 {
00436   return *this=*this*v;
00437 }
00438 
00443 df1b2variable& df1b2variable::operator /= (const df1b2variable& v)
00444 {
00445   return *this=*this/v;
00446 }
00447 
00448 
00449 /*
00450 df1b2variable cos(const df1b2variable& x)
00451 {
00452   cout << "cos not implemented yet" << endl;
00453   ad_exit(1);
00454   return ADf1b2_sin(x);
00455 }
00456 */
00457 
00462 df1b2variable exp(const df1b2variable& x)
00463 {
00464   return ADf1b2_exp(x);
00465 }
00466 
00471 df1b2variable mfexp(const df1b2variable& x)
00472 {
00473   double b=60;
00474   if (value(x)<=b && value(x)>=-b)
00475   {
00476     return ADf1b2_exp(x);
00477   }
00478   else if (value(x)>b)
00479   {
00480     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00481   }
00482   else
00483   {
00484     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00485   }
00486 }
00487 
00492 df1b2variable mfexp(const df1b2variable& x,double b)
00493 {
00494   if (value(x)<=b && value(x)>=-b)
00495   {
00496     return ADf1b2_exp(x);
00497   }
00498   else if (value(x)>b)
00499   {
00500     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00501   }
00502   else
00503   {
00504     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00505   }
00506 }
00507 
00512 df1b2variable log(const df1b2variable& x)
00513 {
00514   return ADf1b2_log(x);
00515 }
00516 
00521 df1b2variable inv(const df1b2variable& x)
00522 {
00523   return ADf1b2_inv(x);
00524 }
00525 
00530 double ADproduct_fun(double x,double y)
00531 {
00532   return x*y;
00533 }
00534 
00539 double ADmult_add_fun(double x,double y)
00540 {
00541   return x*y+x;
00542 }
00543 
00548 double ADdiv_fun(double x,double y)
00549 {
00550   return x/y;
00551 }
00552 
00557 double ADsum_fun(double x,double y)
00558 {
00559   return x+y;
00560 }
00561 
00566 double ADdiff_fun(double x,double y)
00567 {
00568   return x-y;
00569 }
00570 
00575 double ADzero_fun(double x,double y)
00576 {
00577   return 0.0;
00578 }
00579 
00584 double ADzero_fun(double x)
00585 {
00586   return 0.0;
00587 }
00588 
00593 double AD1_fun(double x)
00594 {
00595   return 1.0;
00596 }
00597 
00602 double AD1_fun(double x,double y)
00603 {
00604   return 1.0;
00605 }
00606 
00611 double ADm1_fun(double x,double y)
00612 {
00613   return -1.0;
00614 }
00615 
00620 double AD_id(double x)
00621 {
00622   return x;
00623 }
00624 
00629 double ADfirst_arg(double x,double y)
00630 {
00631   return x;
00632 }
00633 
00638 double ADsecond_arg_plus1(double x,double y)
00639 {
00640   return y+1;
00641 }
00642 
00647 double ADsecond_arg(double x,double y)
00648 {
00649   return y;
00650 }
00651 
00656 double AD_div_1(double x,double y)
00657 {
00658   return 1.0/y;
00659 }
00660 
00665 double AD_div_2(double x,double y)
00666 {
00667   return -x/(y*y);
00668 }
00669 
00674 double AD_div_22(double x,double y)
00675 {
00676   return 2.0*x/(y*y*y);
00677 }
00678 
00683 double AD_div_122(double x,double y)
00684 {
00685   return 2.0/(y*y*y);
00686 }
00687 
00692 double AD_div_12(double x,double y)
00693 {
00694   return -1.0/(y*y);
00695 }
00696 
00701 double AD_div_11(double x,double y)
00702 {
00703   return 0.0;
00704 }
00705 
00710 double AD_div_111(double x,double y)
00711 {
00712   return 0.0;
00713 }
00714 
00719 double AD_div_112(double x,double y)
00720 {
00721   return 0.0;
00722 }
00723 
00728 double AD_div_222(double x,double y)
00729 {
00730   return -6.0*x/(y*y*y*y);
00731 }
00732 
00733 df1b2function2 ADf1b2_div(ADdiv_fun,
00734  AD_div_1,AD_div_2,
00735  AD_div_11,AD_div_12,AD_div_22,
00736  AD_div_111,AD_div_112,AD_div_122,
00737  AD_div_222);
00738 
00739 df1b2function2 ADf1b2_mult_add(ADmult_add_fun,
00740   ADsecond_arg_plus1,ADfirst_arg,
00741   ADzero_fun,AD1_fun,ADzero_fun,
00742   ADzero_fun, ADzero_fun, ADzero_fun,
00743   ADzero_fun,"mult_add");
00744 
00749 df1b2variable mult_add(const df1b2variable& x,const df1b2variable& y)
00750 {
00751   return ADf1b2_mult_add(x,y);
00752 }
00753 
00754 
00755 df1b2function2 ADf1b2_product(ADproduct_fun,
00756   ADsecond_arg,ADfirst_arg,
00757   ADzero_fun,AD1_fun,ADzero_fun,
00758   ADzero_fun, ADzero_fun, ADzero_fun,
00759   ADzero_fun,"product");
00760 
00761 df1b2function2 ADf1b2_diff(ADdiff_fun,
00762   AD1_fun,ADm1_fun,
00763   ADzero_fun,ADzero_fun,ADzero_fun,
00764   ADzero_fun,ADzero_fun,ADzero_fun,
00765   ADzero_fun);
00766 
00771 double ADsquare_fun(double x)
00772 {
00773   return x*x;
00774 }
00775 
00780 double ADthree_square_fun(double x)
00781 {
00782   return 3.0*x*x;
00783 }
00784 
00789 double ADcube_fun(double x)
00790 {
00791   return x*x*x;
00792 }
00793 
00798 double ADtwo_id_fun(double x)
00799 {
00800   return 2.0*x;
00801 }
00802 
00807 double ADsix_id_fun(double x)
00808 {
00809   return 6.0*x;
00810 }
00811 
00816 double ADsix_fun(double x)
00817 {
00818   return 6.0;
00819 }
00820 
00825 double ADtwo_fun(double x)
00826 {
00827   return 2.0;
00828 }
00829 
00830 
00831 df1b2function1 ADf1b2_square(ADsquare_fun,ADtwo_id_fun,ADtwo_fun,ADzero_fun,
00832   "square");
00833 
00838 df1b2variable square(const df1b2variable& x)
00839 {
00840   return ADf1b2_square(x);
00841 }
00842 
00843 df1b2function1 ADf1b2_cube(ADcube_fun,ADthree_square_fun,ADsix_id_fun,ADsix_fun,
00844   "cube");
00845 
00850 df1b2variable cube(const df1b2variable& x)
00851 {
00852   return ADf1b2_cube(x);
00853 }
00854 
00855 /*
00856 df1b2variable fourth(const df1b2variable& _z)
00857 {
00858   ADUNCONST(df1b2variable,z)
00859   df1b2variable u;
00860   double x=*z.get_u();
00861   tmp.get_u()=x3*x.get_u();
00862   double dfx=4.0*cube(x);
00863   double d2fx=12.0*square(x);
00864   double d3fx=24.0*x;
00865 
00866   double * xd=z.get_u_dot();
00867   double * zd=u.get_u_dot();
00868   *u.get_u()=f;
00869   for (int i=0;i<df1b2variable::nvar;i++)
00870   {
00871     *zd++ =dfx * *xd++;
00872   }
00873 
00874   if (!df1b2_gradlist::no_derivatives)
00875     f1b2gradlist->write_pass1(&z,&u,dfx,d2f,d3f);
00876   return(u);
00877 }
00878 */
00879 
00880 /*
00881 df1b2variable operator * (const df1b2variable& x,const df1b2variable& y)
00882 {
00883   return ADf1b2_product(x,y);
00884 }
00885 
00886 
00887 df1b2variable operator / (const df1b2variable& x,const df1b2variable& y)
00888 {
00889   return ADf1b2_div(x,y);
00890 }
00891 */
00892 
00893 
00894 /*
00895 df1b2variable operator - (const df1b2variable& x,const df1b2variable& y)
00896 {
00897   return ADf1b2_diff(x,y);
00898 }
00899 */
00900 
00901 /*
00902 df1b2variable operator * (double x,const df1b2variable& y)
00903 {
00904   return ADf1b2_product(x,y);
00905 }
00906 */
00907 /*
00908 df1b2variable operator * (const df1b2variable& x,double y)
00909 {
00910   return ADf1b2_product(y,x);
00911 }
00912 */
00913 
00914 // the boilerplate for defining the + operator
00915 
00916 df1b2function2 ADf1b2_sum(ADsum_fun,
00917   AD1_fun,AD1_fun,
00918   ADzero_fun,ADzero_fun,ADzero_fun,
00919   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun,"sum");
00920 /*
00921 mydf1b2function2 ADf1b2_mysum(ADsum_fun,
00922   AD1_fun,AD1_fun,
00923   ADzero_fun,ADzero_fun,ADzero_fun,
00924   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun);
00925 */
00926 
00927 df1b2function1 ADf1b2_assign(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00928 
00929 /*
00930 df1b2function1 ADf1b2_pluseq(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00931 
00932 df1b2variable& df1b2variable::operator += (const df1b2variable& y)
00933 {
00934   return ADf1b2_pluseq(*this,y,1); // 1 so that special function is used
00935 }
00936 
00937 
00938 df1b2variable operator + (const df1b2variable& x,const df1b2variable& y)
00939 {
00940   return ADf1b2_sum(x,y);
00941 }
00942 */
00943 
00944 
00945 
00946 /*
00947 df1b2variable mysum(const df1b2variable& x,const df1b2variable& y)
00948 {
00949   return ADf1b2_mysum(x,y);
00950 }
00951 */
00952 
00953 /*
00954 df1b2variable operator + (double x,const df1b2variable& y)
00955 {
00956   return ADf1b2_sum(x,y);
00957 }
00958 */
00959 /*
00960 df1b2variable operator + (const df1b2variable& x,double y)
00961 {
00962   return ADf1b2_sum(x,y);
00963 }
00964 */