ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2fun.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 
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 int mydercheckercounter=0;
00034 
00039 void df1b2variable::save_adpool_pointer(void)
00040 {
00041   if (adpool_stack_pointer> adpool_stack_size-1)
00042   {
00043     cerr << "overflow in save_adpool_pointer" << endl;
00044     ad_exit(1);
00045   }
00046   adpool_stack[adpool_stack_pointer]=pool;
00047   adpool_nvar_stack[adpool_stack_pointer++]=nvar;
00048 }
00049 
00054 void df1b2variable::restore_adpool_pointer(void)
00055 {
00056   if (adpool_stack_pointer<=0)
00057   {
00058     cerr << "underflow in save_adpool_pointer" << endl;
00059     ad_exit(1);
00060   }
00061   pool=adpool_stack[--adpool_stack_pointer];
00062   nvar=adpool_nvar_stack[adpool_stack_pointer];
00063 }
00064 
00069 void initial_df1b2params::save_varsptr(void)
00070 {
00071   initial_df1b2params::varsptr_sav=initial_df1b2params::varsptr;
00072 
00073   varsptr=new P_INITIAL_DF1B2PARAMS[1000];
00074   if (varsptr == 0)
00075   {
00076     cerr << "error allocating memory for "
00077       "initial_df1b2params::varsptr  " << endl;
00078     ad_exit(1);
00079   }
00080 
00081   num_initial_df1b2params_sav=num_initial_df1b2params;
00082   num_initial_df1b2params =0;
00083 }
00084 
00089 void initial_df1b2params::restore_varsptr(void)
00090 {
00091   if (num_initial_df1b2params == 0)
00092   {
00093     if (varsptr)
00094     {
00095       delete [] varsptr;
00096       varsptr=0;
00097     }
00098     varsptr=initial_df1b2params::varsptr_sav;
00099     varsptr_sav=0;
00100 
00101     num_initial_df1b2params= num_initial_df1b2params_sav;
00102     num_initial_df1b2params_sav=0;
00103   }
00104   else
00105   {
00106     if (num_initial_df1b2params+num_initial_df1b2params_sav
00107      > 1000)
00108     {
00109       cerr << "maximum numver of iitial_df1b2params is 1000 "
00110            << endl;
00111       ad_exit(1);
00112     }
00113 
00114     for (int i=0;i<num_initial_df1b2params_sav;i++)
00115     {
00116       varsptr[i+num_initial_df1b2params]=varsptr_sav[i];
00117     }
00118     num_initial_df1b2params+=num_initial_df1b2params_sav;
00119     delete varsptr_sav;
00120     varsptr_sav=0;
00121     num_initial_df1b2params_sav=0;
00122   }
00123 }
00124 
00129 initial_df1b2params::initial_df1b2params(void) : ind_index(0)
00130 {
00131   scalefactor=0.0;
00132   phase_start=1;
00133   phase_save=1;
00134   add_to_list();
00135 }
00136 
00137 /*
00138 static void stupid_xxx(int){;}
00139 static void stupid_xxx(void *){;}
00140 */
00141 
00142 typedef void (**ADprfptr)(void);
00143 typedef void (*ADrfptr)(void);
00144 
00149 void df1b2_gradcalc1(void)
00150 {
00151   //smartlist & list=f1b2gradlist->list;
00152   fixed_smartlist & nlist=f1b2gradlist->nlist;
00153   int ncount=f1b2gradlist->ncount;
00154   //ADrfptr rf2;
00155   int xcount=0;
00156   int tmpcount;
00157   int tcount=f1b2gradlist->ncount;
00158 
00159   //int check_pool_flag3=0;
00160   switch (df1b2variable::passnumber)
00161   {
00162   case 1:
00163     f1b2gradlist->list.save_end();
00164     f1b2gradlist->nlist.save_end();
00165   case 3:
00166 #if defined(__DERCHECK__)
00167     //  derchecker->counter=f1b2gradlist->ncount;
00168       mydercheckercounter=f1b2gradlist->ncount;
00169 #   endif
00170     f1b2gradlist->list.set_reverse();
00171     f1b2gradlist->list.restore_end();
00172     f1b2gradlist->nlist.set_reverse();
00173     f1b2gradlist->nlist.restore_end();
00174 
00175     if (df1b2variable::passnumber==3)
00176     {
00177       f1b2gradlist->nlist3.save_end();
00178       f1b2gradlist->nlist3.set_reverse();
00179       f1b2gradlist->nlist3.restore_end();
00180       f1b2gradlist->list3.save_end();
00181       f1b2gradlist->list3.set_reverse();
00182       f1b2gradlist->list3.restore_end();
00183     }
00184     tmpcount=ncount;
00185     do
00186     {
00187       tmpcount--;
00188      /*
00189       if (!(tmpcount %100))
00190       {
00191         cout << "B " << tmpcount << endl;
00192       }
00193       if (tmpcount < 2)
00194       {
00195         cout << "B " << tmpcount << endl;
00196       }
00197       */
00198       //nlist-=sizeof(char*);
00199       --nlist;
00200       if (nlist.eof())
00201         break;
00202       ADrfptr rf=nlist.bptr->pf;
00203       (*rf)();
00204      /*
00205       if (check_pool_flag3)
00206       {
00207         if (xcount > 722)
00208         {
00209           cout << xcount << "  ";
00210           //check_pool_depths();
00211         }
00212       }
00213      */
00214 #if   defined(__DERCHECK__)
00215         //derchecker->counter--;
00216       mydercheckercounter--;
00217        // stupid_xxx(derchecker->counter);
00218 #     endif
00219         xcount++;
00220           tcount--;
00221        if (xcount > 99999999) cout << xcount << endl;
00222        //if (tcount == 6599 )
00223         //  cout << tcount << endl;
00224       /*
00225         if (initial_df1b2params::current_phase==2)
00226         {
00227           tcount--;
00228           if (rf2 != (nlist.buffer+122)->pf || xcount > 12488)
00229              cout << tcount << " " << xcount << " "
00230                   << (nlist.buffer+122)->pf << endl;
00231         }
00232        */
00233     }
00234     while(nlist.bptr>=nlist.buffer);
00235     break;
00236   case 2:
00237     {
00238       f1b2gradlist->list2.save_end();
00239       f1b2gradlist->list2.set_reverse();
00240       f1b2gradlist->list2.restore_end();
00241       f1b2gradlist->nlist2.save_end();
00242       f1b2gradlist->nlist2.set_reverse();
00243       f1b2gradlist->nlist2.restore_end();
00244 #if   defined(__DERCHECK__)
00245       //  derchecker->counter=0;
00246       mydercheckercounter=0;
00247 #     endif
00248       f1b2gradlist->list.set_forward();
00249       f1b2gradlist->list.rewind();
00250       f1b2gradlist->nlist.set_forward();
00251       f1b2gradlist->nlist.rewind();
00252 
00253       //nlist.reset();  // set pointer to beginning of list
00254       --(f1b2gradlist->nlist2);  // backup one record
00255       int icount=0;
00256       do
00257       {
00258 #if     defined(__DERCHECK__)
00259          // derchecker->counter++;
00260           mydercheckercounter++;
00261 #       endif
00262         //ADrfptr rf=*ADprfptr(nlist.bptr);
00263         ADrfptr rf=nlist.bptr->pf;
00264         (*rf)();
00265         ++nlist;
00266 
00267        /*
00268         if (initial_df1b2params::current_phase==2)
00269         {
00270           if (icount >=4579)
00271           {
00272             cout << icount << endl;
00273           }
00274         }
00275 
00276         if (icount>28115)
00277         {
00278           cout << "icount = " << icount << endl;
00279         }
00280         */
00281       }
00282       while(++icount<ncount);
00283       // need this to get pointer right for step 3?
00284        // nlist-=sizeof(int); // set nlist to point to second record;
00285       //nlist++;
00286 
00287       break;
00288     }
00289 
00290   default:
00291      cerr << "illega value for df1b2variable::passnumber "
00292       " value = " << df1b2variable::passnumber << endl;
00293   }
00294 }
00295 
00296 
00297 
00298 double nsin(double x) {return -sin(x);}
00299 double ncos(double x) {return -cos(x);}
00300 
00301 double ADmult_add_fun(double x,double y);
00302 df1b2function1 ADf1b2_sin(::sin,::cos,::nsin,::ncos,"sin");
00303 df1b2function1 ADf1b2_cos(::cos,::nsin,::ncos,::sin);
00304 df1b2function1 ADf1b2_exp(::exp,::exp,::exp,::exp,"exp");
00305 
00310 double AD_df1_atan(double x)
00311 {
00312   return double(1.0)/(1+square(x));
00313 }
00314 
00319 double AD_df2_atan(double x)
00320 {
00321   return double(-2.0)*x/square(1+square(x));
00322 }
00323 
00328 double AD_df1_tan(double x)
00329 {
00330   return double(1.0)+square(tan(x));
00331 }
00332 
00337 double AD_df2_tan(double x)
00338 {
00339   double y=tan(x);
00340   return double(2.0)*y*(double(1.0)+square(y));
00341 }
00342 
00347 double AD_df3_atan(double x)
00348 {
00349   return double(-2.0)/square(double(1)+square(x))
00350          + double(12.0)*square(x)/cube(double(1)+square(x));
00351 }
00352 
00357 double AD_df3_tan(double x)
00358 {
00359   double y=square(tan(x));
00360 
00361   return double(2.0) * (double(1.0)+double(3.0)*y) * (double(1) + y);
00362 }
00363 
00364 df1b2function1 ADf1b2_tan(::tan,::AD_df1_tan,::AD_df2_tan ,::AD_df3_tan,"tan");
00365 
00366 df1b2function1 ADf1b2_atan(::atan,::AD_df1_atan,::AD_df2_atan ,::AD_df3_atan,
00367   "atan");
00368 
00373 double AD_arg_inv(double x)
00374 {
00375   return double(1.0)/x;
00376 }
00377 
00382 double AD_minus_arg_inv2(double x)
00383 {
00384   return double(-1.0)/(x*x);
00385 }
00386 
00391 double AD_2arg_inv3(double x)
00392 {
00393   return double(2.0)/(x*x*x);
00394 }
00395 
00400 double AD_minus6_arg_inv4(double x)
00401 {
00402   return double(-6.0)/(x*x*x*x);
00403 }
00404 
00405 df1b2function1 ADf1b2_inv(AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3,
00406   AD_minus6_arg_inv4);
00407 
00408 df1b2function1 ADf1b2_log(::log,AD_arg_inv,AD_minus_arg_inv2,AD_2arg_inv3);
00409 
00410 /*
00411 df1b2variable sin(const df1b2variable& x)
00412 {
00413   return ADf1b2_sin(x);
00414 }
00415 */
00416 
00421 df1b2variable atan(const df1b2variable& x)
00422 {
00423   return ADf1b2_atan(x);
00424 }
00425 
00430 df1b2variable tan(const df1b2variable& x)
00431 {
00432   return ADf1b2_tan(x);
00433 }
00434 
00439 df1b2variable& df1b2variable::operator *= (const df1b2variable& v)
00440 {
00441   return *this=*this*v;
00442 }
00443 
00448 df1b2variable& df1b2variable::operator /= (const df1b2variable& v)
00449 {
00450   return *this=*this/v;
00451 }
00452 
00453 
00454 /*
00455 df1b2variable cos(const df1b2variable& x)
00456 {
00457   cout << "cos not implemented yet" << endl;
00458   ad_exit(1);
00459   return ADf1b2_sin(x);
00460 }
00461 */
00462 
00467 df1b2variable exp(const df1b2variable& x)
00468 {
00469   return ADf1b2_exp(x);
00470 }
00471 
00476 df1b2variable mfexp(const df1b2variable& x)
00477 {
00478   double b=60;
00479   if (value(x)<=b && value(x)>=-b)
00480   {
00481     return ADf1b2_exp(x);
00482   }
00483   else if (value(x)>b)
00484   {
00485     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00486   }
00487   else
00488   {
00489     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00490   }
00491 }
00492 
00497 df1b2variable mfexp(const df1b2variable& x,double b)
00498 {
00499   if (value(x)<=b && value(x)>=-b)
00500   {
00501     return ADf1b2_exp(x);
00502   }
00503   else if (value(x)>b)
00504   {
00505     return exp(b)*(double(1.)+double(2.)*(x-b))/(double(1.)+x-b);
00506   }
00507   else
00508   {
00509     return exp(-b)*(double(1.)-x-b)/(double(1.)+double(2.)*(-x-b));
00510   }
00511 }
00512 
00517 df1b2variable log(const df1b2variable& x)
00518 {
00519   return ADf1b2_log(x);
00520 }
00521 
00526 df1b2variable inv(const df1b2variable& x)
00527 {
00528   return ADf1b2_inv(x);
00529 }
00530 
00535 double ADproduct_fun(double x,double y)
00536 {
00537   return x*y;
00538 }
00539 
00544 double ADmult_add_fun(double x,double y)
00545 {
00546   return x*y+x;
00547 }
00548 
00553 double ADdiv_fun(double x,double y)
00554 {
00555   return x/y;
00556 }
00557 
00562 double ADsum_fun(double x,double y)
00563 {
00564   return x+y;
00565 }
00566 
00571 double ADdiff_fun(double x,double y)
00572 {
00573   return x-y;
00574 }
00575 
00580 double ADzero_fun(double x,double y)
00581 {
00582   return 0.0;
00583 }
00584 
00589 double ADzero_fun(double x)
00590 {
00591   return 0.0;
00592 }
00593 
00598 double AD1_fun(double x)
00599 {
00600   return 1.0;
00601 }
00602 
00607 double AD1_fun(double x,double y)
00608 {
00609   return 1.0;
00610 }
00611 
00616 double ADm1_fun(double x,double y)
00617 {
00618   return -1.0;
00619 }
00620 
00625 double AD_id(double x)
00626 {
00627   return x;
00628 }
00629 
00634 double ADfirst_arg(double x,double y)
00635 {
00636   return x;
00637 }
00638 
00643 double ADsecond_arg_plus1(double x,double y)
00644 {
00645   return y+1;
00646 }
00647 
00652 double ADsecond_arg(double x,double y)
00653 {
00654   return y;
00655 }
00656 
00661 double AD_div_1(double x,double y)
00662 {
00663   return 1.0/y;
00664 }
00665 
00670 double AD_div_2(double x,double y)
00671 {
00672   return -x/(y*y);
00673 }
00674 
00679 double AD_div_22(double x,double y)
00680 {
00681   return 2.0*x/(y*y*y);
00682 }
00683 
00688 double AD_div_122(double x,double y)
00689 {
00690   return 2.0/(y*y*y);
00691 }
00692 
00697 double AD_div_12(double x,double y)
00698 {
00699   return -1.0/(y*y);
00700 }
00701 
00706 double AD_div_11(double x,double y)
00707 {
00708   return 0.0;
00709 }
00710 
00715 double AD_div_111(double x,double y)
00716 {
00717   return 0.0;
00718 }
00719 
00724 double AD_div_112(double x,double y)
00725 {
00726   return 0.0;
00727 }
00728 
00733 double AD_div_222(double x,double y)
00734 {
00735   return -6.0*x/(y*y*y*y);
00736 }
00737 
00738 df1b2function2 ADf1b2_div(ADdiv_fun,
00739  AD_div_1,AD_div_2,
00740  AD_div_11,AD_div_12,AD_div_22,
00741  AD_div_111,AD_div_112,AD_div_122,
00742  AD_div_222);
00743 
00744 df1b2function2 ADf1b2_mult_add(ADmult_add_fun,
00745   ADsecond_arg_plus1,ADfirst_arg,
00746   ADzero_fun,AD1_fun,ADzero_fun,
00747   ADzero_fun, ADzero_fun, ADzero_fun,
00748   ADzero_fun,"mult_add");
00749 
00754 df1b2variable mult_add(const df1b2variable& x,const df1b2variable& y)
00755 {
00756   return ADf1b2_mult_add(x,y);
00757 }
00758 
00759 
00760 df1b2function2 ADf1b2_product(ADproduct_fun,
00761   ADsecond_arg,ADfirst_arg,
00762   ADzero_fun,AD1_fun,ADzero_fun,
00763   ADzero_fun, ADzero_fun, ADzero_fun,
00764   ADzero_fun,"product");
00765 
00766 df1b2function2 ADf1b2_diff(ADdiff_fun,
00767   AD1_fun,ADm1_fun,
00768   ADzero_fun,ADzero_fun,ADzero_fun,
00769   ADzero_fun,ADzero_fun,ADzero_fun,
00770   ADzero_fun);
00771 
00776 double ADsquare_fun(double x)
00777 {
00778   return x*x;
00779 }
00780 
00785 double ADthree_square_fun(double x)
00786 {
00787   return 3.0*x*x;
00788 }
00789 
00794 double ADcube_fun(double x)
00795 {
00796   return x*x*x;
00797 }
00798 
00803 double ADtwo_id_fun(double x)
00804 {
00805   return 2.0*x;
00806 }
00807 
00812 double ADsix_id_fun(double x)
00813 {
00814   return 6.0*x;
00815 }
00816 
00821 double ADsix_fun(double x)
00822 {
00823   return 6.0;
00824 }
00825 
00830 double ADtwo_fun(double x)
00831 {
00832   return 2.0;
00833 }
00834 
00835 
00836 df1b2function1 ADf1b2_square(ADsquare_fun,ADtwo_id_fun,ADtwo_fun,ADzero_fun,
00837   "square");
00838 
00843 df1b2variable square(const df1b2variable& x)
00844 {
00845   return ADf1b2_square(x);
00846 }
00847 
00848 df1b2function1 ADf1b2_cube(ADcube_fun,ADthree_square_fun,ADsix_id_fun,ADsix_fun,
00849   "cube");
00850 
00855 df1b2variable cube(const df1b2variable& x)
00856 {
00857   return ADf1b2_cube(x);
00858 }
00859 
00860 /*
00861 df1b2variable fourth(const df1b2variable& _z)
00862 {
00863   ADUNCONST(df1b2variable,z)
00864   df1b2variable u;
00865   double x=*z.get_u();
00866   tmp.get_u()=x3*x.get_u();
00867   double dfx=4.0*cube(x);
00868   double d2fx=12.0*square(x);
00869   double d3fx=24.0*x;
00870 
00871   double * xd=z.get_u_dot();
00872   double * zd=u.get_u_dot();
00873   *u.get_u()=f;
00874   for (int i=0;i<df1b2variable::nvar;i++)
00875   {
00876     *zd++ =dfx * *xd++;
00877   }
00878 
00879   if (!df1b2_gradlist::no_derivatives)
00880     f1b2gradlist->write_pass1(&z,&u,dfx,d2f,d3f);
00881   return(u);
00882 }
00883 */
00884 
00885 /*
00886 df1b2variable operator * (const df1b2variable& x,const df1b2variable& y)
00887 {
00888   return ADf1b2_product(x,y);
00889 }
00890 
00891 
00892 df1b2variable operator / (const df1b2variable& x,const df1b2variable& y)
00893 {
00894   return ADf1b2_div(x,y);
00895 }
00896 */
00897 
00898 
00899 /*
00900 df1b2variable operator - (const df1b2variable& x,const df1b2variable& y)
00901 {
00902   return ADf1b2_diff(x,y);
00903 }
00904 */
00905 
00906 /*
00907 df1b2variable operator * (double x,const df1b2variable& y)
00908 {
00909   return ADf1b2_product(x,y);
00910 }
00911 */
00912 /*
00913 df1b2variable operator * (const df1b2variable& x,double y)
00914 {
00915   return ADf1b2_product(y,x);
00916 }
00917 */
00918 
00919 // the boilerplate for defining the + operator
00920 
00921 df1b2function2 ADf1b2_sum(ADsum_fun,
00922   AD1_fun,AD1_fun,
00923   ADzero_fun,ADzero_fun,ADzero_fun,
00924   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun,"sum");
00925 /*
00926 mydf1b2function2 ADf1b2_mysum(ADsum_fun,
00927   AD1_fun,AD1_fun,
00928   ADzero_fun,ADzero_fun,ADzero_fun,
00929   ADzero_fun, ADzero_fun, ADzero_fun, ADzero_fun);
00930 */
00931 
00932 df1b2function1 ADf1b2_assign(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00933 
00934 /*
00935 df1b2function1 ADf1b2_pluseq(AD_id,AD1_fun,ADzero_fun,ADzero_fun);
00936 
00937 df1b2variable& df1b2variable::operator += (const df1b2variable& y)
00938 {
00939   return ADf1b2_pluseq(*this,y,1); // 1 so that special function is used
00940 }
00941 
00942 
00943 df1b2variable operator + (const df1b2variable& x,const df1b2variable& y)
00944 {
00945   return ADf1b2_sum(x,y);
00946 }
00947 */
00948 
00949 
00950 
00951 /*
00952 df1b2variable mysum(const df1b2variable& x,const df1b2variable& y)
00953 {
00954   return ADf1b2_mysum(x,y);
00955 }
00956 */
00957 
00958 /*
00959 df1b2variable operator + (double x,const df1b2variable& y)
00960 {
00961   return ADf1b2_sum(x,y);
00962 }
00963 */
00964 /*
00965 df1b2variable operator + (const df1b2variable& x,double y)
00966 {
00967   return ADf1b2_sum(x,y);
00968 }
00969 */