ADMB Documentation  11.1.2192
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fun.cpp 2118 2014-06-18 22:44:43Z 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 /*
00140 static void stupid_xxx(int){;}
00141 static void stupid_xxx(void *){;}
00142 */
00143 
00144 typedef void (**ADprfptr)(void);
00145 typedef void (*ADrfptr)(void);
00146 
00151 void df1b2_gradcalc1(void)
00152 {
00153   //smartlist & list=f1b2gradlist->list;
00154   fixed_smartlist & nlist=f1b2gradlist->nlist;
00155   int ncount=f1b2gradlist->ncount;
00156   //ADrfptr rf2;
00157   int xcount=0;
00158   int tmpcount;
00159   int tcount=f1b2gradlist->ncount;
00160 
00161   //int check_pool_flag3=0;
00162   switch (df1b2variable::passnumber)
00163   {
00164   case 1:
00165     f1b2gradlist->list.save_end();
00166     f1b2gradlist->nlist.save_end();
00167   case 3:
00168 #if defined(__DERCHECK__)
00169     //  derchecker->counter=f1b2gradlist->ncount;
00170     mydercheckercounter=f1b2gradlist->ncount;
00171 #endif
00172     f1b2gradlist->list.set_reverse();
00173     f1b2gradlist->list.restore_end();
00174     f1b2gradlist->nlist.set_reverse();
00175     f1b2gradlist->nlist.restore_end();
00176 
00177     if (df1b2variable::passnumber==3)
00178     {
00179       f1b2gradlist->nlist3.save_end();
00180       f1b2gradlist->nlist3.set_reverse();
00181       f1b2gradlist->nlist3.restore_end();
00182       f1b2gradlist->list3.save_end();
00183       f1b2gradlist->list3.set_reverse();
00184       f1b2gradlist->list3.restore_end();
00185     }
00186     tmpcount=ncount;
00187     do
00188     {
00189       tmpcount--;
00190      /*
00191       if (!(tmpcount %100))
00192       {
00193         cout << "B " << tmpcount << endl;
00194       }
00195       if (tmpcount < 2)
00196       {
00197         cout << "B " << tmpcount << endl;
00198       }
00199       */
00200       //nlist-=sizeof(char*);
00201       --nlist;
00202       if (nlist.eof())
00203         break;
00204       ADrfptr rf=nlist.bptr->pf;
00205       (*rf)();
00206      /*
00207       if (check_pool_flag3)
00208       {
00209         if (xcount > 722)
00210         {
00211           cout << xcount << "  ";
00212           //check_pool_depths();
00213         }
00214       }
00215      */
00216 #if defined(__DERCHECK__)
00217         //derchecker->counter--;
00218       mydercheckercounter--;
00219        // stupid_xxx(derchecker->counter);
00220 #endif
00221         xcount++;
00222           tcount--;
00223        if (xcount > 99999999) cout << xcount << endl;
00224        //if (tcount == 6599 )
00225         //  cout << tcount << endl;
00226       /*
00227         if (initial_df1b2params::current_phase==2)
00228         {
00229           tcount--;
00230           if (rf2 != (nlist.buffer+122)->pf || xcount > 12488)
00231              cout << tcount << " " << xcount << " "
00232                   << (nlist.buffer+122)->pf << endl;
00233         }
00234        */
00235     }
00236     while(nlist.bptr>=nlist.buffer);
00237     break;
00238   case 2:
00239     {
00240       f1b2gradlist->list2.save_end();
00241       f1b2gradlist->list2.set_reverse();
00242       f1b2gradlist->list2.restore_end();
00243       f1b2gradlist->nlist2.save_end();
00244       f1b2gradlist->nlist2.set_reverse();
00245       f1b2gradlist->nlist2.restore_end();
00246 #if defined(__DERCHECK__)
00247       // derchecker->counter=0;
00248       mydercheckercounter=0;
00249 #endif
00250       f1b2gradlist->list.set_forward();
00251       f1b2gradlist->list.rewind();
00252       f1b2gradlist->nlist.set_forward();
00253       f1b2gradlist->nlist.rewind();
00254 
00255       //nlist.reset();  // set pointer to beginning of list
00256       --(f1b2gradlist->nlist2);  // backup one record
00257       int icount=0;
00258       do
00259       {
00260 #if defined(__DERCHECK__)
00261         // derchecker->counter++;
00262         mydercheckercounter++;
00263 #endif
00264         //ADrfptr rf=*ADprfptr(nlist.bptr);
00265         ADrfptr rf=nlist.bptr->pf;
00266         (*rf)();
00267         ++nlist;
00268 
00269        /*
00270         if (initial_df1b2params::current_phase==2)
00271         {
00272           if (icount >=4579)
00273           {
00274             cout << icount << endl;
00275           }
00276         }
00277 
00278         if (icount>28115)
00279         {
00280           cout << "icount = " << icount << endl;
00281         }
00282         */
00283       }
00284       while(++icount<ncount);
00285       // need this to get pointer right for step 3?
00286        // nlist-=sizeof(int); // set nlist to point to second record;
00287       //nlist++;
00288 
00289       break;
00290     }
00291 
00292   default:
00293      cerr << "illega value for df1b2variable::passnumber "
00294       " value = " << df1b2variable::passnumber << endl;
00295   }
00296 }
00297 
00298 
00299 
00300 double nsin(double x) {return -sin(x);}
00301 double ncos(double x) {return -cos(x);}
00302 
00303 double ADmult_add_fun(double x,double y);
00304 df1b2function1 ADf1b2_sin(::sin,::cos,::nsin,::ncos,"sin");
00305 df1b2function1 ADf1b2_cos(::cos,::nsin,::ncos,::sin);
00306 df1b2function1 ADf1b2_exp(::exp,::exp,::exp,::exp,"exp");
00307 
00312 double AD_df1_atan(double x)
00313 {
00314   return double(1.0)/(1+square(x));
00315 }
00316 
00321 double AD_df2_atan(double x)
00322 {
00323   return double(-2.0)*x/square(1+square(x));
00324 }
00325 
00330 double AD_df1_tan(double x)
00331 {
00332   return double(1.0)+square(tan(x));
00333 }
00334 
00339 double AD_df2_tan(double x)
00340 {
00341   double y=tan(x);
00342   return double(2.0)*y*(double(1.0)+square(y));
00343 }
00344 
00349 double AD_df3_atan(double x)
00350 {
00351   return double(-2.0)/square(double(1)+square(x))
00352          + double(12.0)*square(x)/cube(double(1)+square(x));
00353 }
00354 
00359 double AD_df3_tan(double x)
00360 {
00361   double y=square(tan(x));
00362 
00363   return double(2.0) * (double(1.0)+double(3.0)*y) * (double(1) + y);
00364 }
00365 
00366 df1b2function1 ADf1b2_tan(::tan,::AD_df1_tan,::AD_df2_tan ,::AD_df3_tan,"tan");
00367 
00368 df1b2function1 ADf1b2_atan(::atan,::AD_df1_atan,::AD_df2_atan ,::AD_df3_atan,
00369   "atan");
00370 
00375 double AD_arg_inv(double x)
00376 {
00377   return double(1.0)/x;
00378 }
00379 
00384 double AD_minus_arg_inv2(double x)
00385 {
00386   return double(-1.0)/(x*x);
00387 }
00388 
00393 double AD_2arg_inv3(double x)
00394 {
00395   return double(2.0)/(x*x*x);
00396 }
00397 
00402 double AD_minus6_arg_inv4(double x)
00403 {
00404   return double(-6.0)/(x*x*x*x);
00405 }
00406 
00407 df1b2function1 ADf1b2_inv(AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3,
00408   AD_minus6_arg_inv4);
00409 
00410 df1b2function1 ADf1b2_log(::log,AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3);
00411 
00412 /*
00413 df1b2variable sin(const df1b2variable& x)
00414 {
00415   return ADf1b2_sin(x);
00416 }
00417 */
00418 
00423 df1b2variable atan(const df1b2variable& x)
00424 {
00425   return ADf1b2_atan(x);
00426 }
00427 
00432 df1b2variable tan(const df1b2variable& x)
00433 {
00434   return ADf1b2_tan(x);
00435 }
00436 
00441 df1b2variable& df1b2variable::operator *= (const df1b2variable& v)
00442 {
00443   return *this=*this*v;
00444 }
00445 
00450 df1b2variable& df1b2variable::operator /= (const df1b2variable& v)
00451 {
00452   return *this=*this/v;
00453 }
00454 
00455 
00456 /*
00457 df1b2variable cos(const df1b2variable& x)
00458 {
00459   cout << "cos not implemented yet" << endl;
00460   ad_exit(1);
00461   return ADf1b2_sin(x);
00462 }
00463 */
00464 
00469 df1b2variable exp(const df1b2variable& x)
00470 {
00471   return ADf1b2_exp(x);
00472 }
00473 
00478 df1b2variable mfexp(const df1b2variable& x)
00479 {
00480   double b=60;
00481   if (value(x)<=b && value(x)>=-b)
00482   {
00483     return ADf1b2_exp(x);
00484   }
00485   else if (value(x)>b)
00486   {
00487     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00488   }
00489   else
00490   {
00491     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00492   }
00493 }
00494 
00499 df1b2variable mfexp(const df1b2variable& x,double b)
00500 {
00501   if (value(x)<=b && value(x)>=-b)
00502   {
00503     return ADf1b2_exp(x);
00504   }
00505   else if (value(x)>b)
00506   {
00507     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00508   }
00509   else
00510   {
00511     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00512   }
00513 }
00514 
00519 df1b2variable log(const df1b2variable& x)
00520 {
00521   return ADf1b2_log(x);
00522 }
00523 
00528 df1b2variable inv(const df1b2variable& x)
00529 {
00530   return ADf1b2_inv(x);
00531 }
00532 
00537 double ADproduct_fun(double x,double y)
00538 {
00539   return x*y;
00540 }
00541 
00546 double ADmult_add_fun(double x,double y)
00547 {
00548   return x*y+x;
00549 }
00550 
00555 double ADdiv_fun(double x,double y)
00556 {
00557   return x/y;
00558 }
00559 
00564 double ADsum_fun(double x,double y)
00565 {
00566   return x+y;
00567 }
00568 
00573 double ADdiff_fun(double x,double y)
00574 {
00575   return x-y;
00576 }
00577 
00582 double ADzero_fun(double x,double y)
00583 {
00584   return 0.0;
00585 }
00586 
00591 double ADzero_fun(double x)
00592 {
00593   return 0.0;
00594 }
00595 
00600 double AD1_fun(double x)
00601 {
00602   return 1.0;
00603 }
00604 
00609 double AD1_fun(double x,double y)
00610 {
00611   return 1.0;
00612 }
00613 
00618 double ADm1_fun(double x,double y)
00619 {
00620   return -1.0;
00621 }
00622 
00627 double AD_id(double x)
00628 {
00629   return x;
00630 }
00631 
00636 double ADfirst_arg(double x,double y)
00637 {
00638   return x;
00639 }
00640 
00645 double ADsecond_arg_plus1(double x,double y)
00646 {
00647   return y+1;
00648 }
00649 
00654 double ADsecond_arg(double x,double y)
00655 {
00656   return y;
00657 }
00658 
00663 double AD_div_1(double x,double y)
00664 {
00665   return 1.0/y;
00666 }
00667 
00672 double AD_div_2(double x,double y)
00673 {
00674   return -x/(y*y);
00675 }
00676 
00681 double AD_div_22(double x,double y)
00682 {
00683   return 2.0*x/(y*y*y);
00684 }
00685 
00690 double AD_div_122(double x,double y)
00691 {
00692   return 2.0/(y*y*y);
00693 }
00694 
00699 double AD_div_12(double x,double y)
00700 {
00701   return -1.0/(y*y);
00702 }
00703 
00708 double AD_div_11(double x,double y)
00709 {
00710   return 0.0;
00711 }
00712 
00717 double AD_div_111(double x,double y)
00718 {
00719   return 0.0;
00720 }
00721 
00726 double AD_div_112(double x,double y)
00727 {
00728   return 0.0;
00729 }
00730 
00735 double AD_div_222(double x,double y)
00736 {
00737   return -6.0*x/(y*y*y*y);
00738 }
00739 
00740 df1b2function2 ADf1b2_div(ADdiv_fun,
00741  AD_div_1,AD_div_2,
00742  AD_div_11,AD_div_12,AD_div_22,
00743  AD_div_111,AD_div_112,AD_div_122,
00744  AD_div_222);
00745 
00746 df1b2function2 ADf1b2_mult_add(ADmult_add_fun,
00747   ADsecond_arg_plus1,ADfirst_arg,
00748   ADzero_fun,AD1_fun,ADzero_fun,
00749   ADzero_fun, ADzero_fun, ADzero_fun,
00750   ADzero_fun,"mult_add");
00751 
00756 df1b2variable mult_add(const df1b2variable& x,const df1b2variable& y)
00757 {
00758   return ADf1b2_mult_add(x,y);
00759 }
00760 
00761 
00762 df1b2function2 ADf1b2_product(ADproduct_fun,
00763   ADsecond_arg,ADfirst_arg,
00764   ADzero_fun,AD1_fun,ADzero_fun,
00765   ADzero_fun, ADzero_fun, ADzero_fun,
00766   ADzero_fun,"product");
00767 
00768 df1b2function2 ADf1b2_diff(ADdiff_fun,
00769   AD1_fun,ADm1_fun,
00770   ADzero_fun,ADzero_fun,ADzero_fun,
00771   ADzero_fun,ADzero_fun,ADzero_fun,
00772   ADzero_fun);
00773 
00778 double ADsquare_fun(double x)
00779 {
00780   return x*x;
00781 }
00782 
00787 double ADthree_square_fun(double x)
00788 {
00789   return 3.0*x*x;
00790 }
00791 
00796 double ADcube_fun(double x)
00797 {
00798   return x*x*x;
00799 }
00800 
00805 double ADtwo_id_fun(double x)
00806 {
00807   return 2.0*x;
00808 }
00809 
00814 double ADsix_id_fun(double x)
00815 {
00816   return 6.0*x;
00817 }
00818 
00823 double ADsix_fun(double x)
00824 {
00825   return 6.0;
00826 }
00827 
00832 double ADtwo_fun(double x)
00833 {
00834   return 2.0;
00835 }
00836 
00837 
00838 df1b2function1 ADf1b2_square(ADsquare_fun,ADtwo_id_fun,ADtwo_fun,ADzero_fun,
00839   "square");
00840 
00845 df1b2variable square(const df1b2variable& x)
00846 {
00847   return ADf1b2_square(x);
00848 }
00849 
00850 df1b2function1 ADf1b2_cube(ADcube_fun,ADthree_square_fun,ADsix_id_fun,ADsix_fun,
00851   "cube");
00852 
00857 df1b2variable cube(const df1b2variable& x)
00858 {
00859   return ADf1b2_cube(x);
00860 }
00861 
00862 /*
00863 df1b2variable fourth(const df1b2variable& _z)
00864 {
00865   ADUNCONST(df1b2variable,z)
00866   df1b2variable u;
00867   double x=*z.get_u();
00868   tmp.get_u()=x3*x.get_u();
00869   double dfx=4.0*cube(x);
00870   double d2fx=12.0*square(x);
00871   double d3fx=24.0*x;
00872 
00873   double * xd=z.get_u_dot();
00874   double * zd=u.get_u_dot();
00875   *u.get_u()=f;
00876   for (int i=0;i<df1b2variable::nvar;i++)
00877   {
00878     *zd++ =dfx * *xd++;
00879   }
00880 
00881   if (!df1b2_gradlist::no_derivatives)
00882     f1b2gradlist->write_pass1(&z,&u,dfx,d2f,d3f);
00883   return(u);
00884 }
00885 */
00886 
00887 /*
00888 df1b2variable operator * (const df1b2variable& x,const df1b2variable& y)
00889 {
00890   return ADf1b2_product(x,y);
00891 }
00892 
00893 
00894 df1b2variable operator / (const df1b2variable& x,const df1b2variable& y)
00895 {
00896   return ADf1b2_div(x,y);
00897 }
00898 */
00899 
00900 
00901 /*
00902 df1b2variable operator - (const df1b2variable& x,const df1b2variable& y)
00903 {
00904   return ADf1b2_diff(x,y);
00905 }
00906 */
00907 
00908 /*
00909 df1b2variable operator * (double x,const df1b2variable& y)
00910 {
00911   return ADf1b2_product(x,y);
00912 }
00913 */
00914 /*
00915 df1b2variable operator * (const df1b2variable& x,double y)
00916 {
00917   return ADf1b2_product(y,x);
00918 }
00919 */
00920 
00921 // the boilerplate for defining the + operator
00922 
00923 df1b2function2 ADf1b2_sum(ADsum_fun,
00924   AD1_fun,AD1_fun,
00925   ADzero_fun,ADzero_fun,ADzero_fun,
00926   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun,"sum");
00927 /*
00928 mydf1b2function2 ADf1b2_mysum(ADsum_fun,
00929   AD1_fun,AD1_fun,
00930   ADzero_fun,ADzero_fun,ADzero_fun,
00931   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun);
00932 */
00933 
00934 df1b2function1 ADf1b2_assign(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00935 
00936 /*
00937 df1b2function1 ADf1b2_pluseq(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00938 
00939 df1b2variable& df1b2variable::operator += (const df1b2variable& y)
00940 {
00941   return ADf1b2_pluseq(*this,y,1); // 1 so that special function is used
00942 }
00943 
00944 
00945 df1b2variable operator + (const df1b2variable& x,const df1b2variable& y)
00946 {
00947   return ADf1b2_sum(x,y);
00948 }
00949 */
00950 
00951 
00952 
00953 /*
00954 df1b2variable mysum(const df1b2variable& x,const df1b2variable& y)
00955 {
00956   return ADf1b2_mysum(x,y);
00957 }
00958 */
00959 
00960 /*
00961 df1b2variable operator + (double x,const df1b2variable& y)
00962 {
00963   return ADf1b2_sum(x,y);
00964 }
00965 */
00966 /*
00967 df1b2variable operator + (const df1b2variable& x,double y)
00968 {
00969   return ADf1b2_sum(x,y);
00970 }
00971 */