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