ADMB Documentation: vbetacf.cpp (* $Id: vbetacf.cpp 789 2010-10-05 01:01:09Z johnoel $) Source File

Go to the documentation of this file.
00001 
00007 #include <fvar.hpp>
00008 #include <math.h>
00009 #define EPS 1.0e-9
00010 #define FPMIN 1.0e-30
00011 
00024 dvariable betacf(const dvariable& _a, const dvariable& _b, const dvariable& _x,
00025   int MAXIT)
00026 {
00027   double qab,qam,qap;
00028   double a=value(_a);
00029   double b=value(_b);
00030   double x=value(_x);
00031 
00032   qab=a+b;
00033   qap=a+1.0;
00034   qam=a-1.0;
00035   dvector c1(0,MAXIT);
00036   dvector c(1,MAXIT);
00037   dvector d1(0,MAXIT);
00038   dvector d(1,MAXIT);
00039   dvector del(1,MAXIT);
00040   dvector h1(0,MAXIT);
00041   dvector h(1,MAXIT);
00042   dvector aa(1,MAXIT);
00043   dvector aa1(1,MAXIT);
00044   c1(0)=1.0;
00045   d1(0)=1.0/(1.0-qab*x/qap);
00046   h1(0)=d1(0);
00047 
00048   int m = 1;
00049   for (; m <= MAXIT; m++)
00050   {
00051     int i=m;
00052     int m2=2*m;
00053     aa(i)=m*(b-m)*x/((qam+m2)*(a+m2));
00054     d(i)=1.0/(1.0+aa(i)*d1(i-1));
00055     c(i)=1.0+aa(i)/c1(i-1);
00056     h(i) = h1(i-1)*d(i)*c(i);
00057     aa1(i) = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00058     d1(i)=1.0/(1.0+aa1(i)*d(i));
00059     c1(i)=1.0+aa1(i)/c(i);
00060     del(i)=d1(i)*c1(i);
00061     h1(i) = h(i)*del(i);
00062     if (fabs(del(i)-1.0) < EPS) break;
00063   }
00064   if (m > MAXIT)
00065   {
00066     cerr << "a or b too big, or MAXIT too small in cumulative beta function"
00067       " routine" << endl;
00068     m=MAXIT;
00069   }
00070   int mmax=m;
00071   dvariable hh;
00072   value(hh)=h1(mmax);
00073 
00074 
00075   dvector dfc1(0,MAXIT);
00076   dvector dfc(1,MAXIT);
00077   dvector dfd1(0,MAXIT);
00078   dvector dfd(1,MAXIT);
00079   dvector dfh1(0,MAXIT);
00080   dvector dfh(1,MAXIT);
00081   dvector dfaa(1,MAXIT);
00082   dvector dfaa1(1,MAXIT);
00083   dvector dfdel(1,MAXIT);
00084 
00085   dfc1.initialize();
00086   dfc.initialize();
00087   dfaa1.initialize();
00088   dfaa.initialize();
00089   dfd1.initialize();
00090   dfd.initialize();
00091   dfh1.initialize();
00092   dfh.initialize();
00093   dfdel.initialize();
00094   dfh1(mmax)=1.0;
00095   double dfqab=0.0;
00096   double dfqam=0.0;
00097   double dfqap=0.0;
00098   double dfa=0.0;
00099   double dfb=0.0;
00100   double dfx=0.0;
00101 
00102   for (m=mmax;m>=1;m--)
00103   {
00104    /*
00105     int i=m;
00106     m2=2*m;
00107     aa(i)=m*(b-m)*x/((qam+m2)*(a+m2));
00108     d(i)=1.0/(1.0+aa(i)*d1(i-1));
00109     c(i)=1.0+aa(i)/c1(i-1);
00110     h(i) = h1(i-1)*d(i)*c(i);
00111     aa1(i) = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00112     d1(i)=1.0/(1.0+aa1(i)*d(i));
00113     c1(i)=1.0+aa1(i)/c(i);
00114     del(i)=d1(i)*c1(i);
00115     h1(i) = h(i)*del(i);
00116    */
00117 
00118     int i=m;
00119     int m2=2*m;
00120 
00121     //h1(i) = h(i)*del(i);
00122 
00123     dfh(i)+=dfh1(i)*del(i);
00124     dfdel(i)+=dfh1(i)*h(i);
00125     dfh1(i)=0.0;
00126 
00127     //del(i)=d1(i)*c1(i);
00128 
00129     dfd1(i)+=dfdel(i)*c1(i);
00130     dfc1(i)+=dfdel(i)*d1(i);
00131     dfdel(i)=0.0;
00132 
00133     //c1(i)=1.0+aa1(i)/c(i);
00134 
00135     dfaa1(i)+=dfc1(i)/c(i);
00136     dfc(i)-=dfc1(i)*aa1(i)/(c(i)*c(i));
00137     dfc1(i)=0.0;
00138 
00139     //d1(i)=1.0/(1.0+aa1(i)*d(i));
00140     double sq=square(d1(i));
00141     dfaa1(i)-=dfd1(i)*sq*d(i);
00142     dfd(i)-=dfd1(i)*sq*aa1(i);
00143     dfd1(i)=0.0;
00144 
00145     //aa1(i) = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
00146     dfx -= dfaa1(i) *
00147      (a+m)*(qab+m)/((a+m2)*(qap+m2));
00148 
00149     dfa += dfaa1(i) * aa1(i)* (1.0/(a+m) - 1.0/(a+m2));
00150     dfqab += dfaa1(i) * aa1(i)/(qab+m);
00151     dfqap += dfaa1(i) * aa1(i)* (-1.0/(qap+m2));
00152     dfaa1(i)=0.0;
00153 
00154 
00155     //h(i) = h1(i-1)*d(i)*c(i);
00156     dfh1(i-1)+=dfh(i)*d(i)*c(i);
00157     dfd(i)+=dfh(i)*h1(i-1)*c(i);
00158     dfc(i)+=dfh(i)*h1(i-1)*d(i);
00159     dfh(i)=0.0;
00160 
00161     //c(i)=1.0+aa(i)/c1(i-1);
00162     dfaa(i)+=dfc(i)/c1(i-1);
00163     dfc1(i-1)-=dfc(i)*aa(i)/square(c1(i-1));
00164     dfc(i)=0.0;
00165 
00166 
00167     //d(i)=1.0/(1.0+aa(i)*d1(i-1));
00168     dfaa(i)-=dfd(i)*square(d(i))*d1(i-1);
00169     dfd1(i-1)-=dfd(i)*square(d(i))*aa(i);
00170     dfd(i)=0.0;
00171 
00172     //aa(i)=m*(b-m)*x/((qam+m2)*(a+m2));
00173     dfx+=dfaa(i)*
00174       m*(b-m)/((qam+m2)*(a+m2));
00175 
00176     dfb+=dfaa(i)*
00177       m*x/((qam+m2)*(a+m2));
00178 
00179 
00180     dfa-=dfaa(i)*aa(i)/(a+m2);
00181     dfqam-=dfaa(i)*aa(i)/(qam+m2);
00182     dfaa(i)=0.0;
00183   }
00184   /*
00185   c1(0)=1.0;
00186   d1(0)=1.0/(1.0-qab*x/qap);
00187   h1(0)=d1(0);
00188  */
00189 
00190   //h1(0)=d1(0);
00191   dfd1(0)+=dfh1(0);
00192   dfh1(0)=0.0;
00193 
00194   //d1(0)=1.0/(1.0-qab*x/qap);
00195   double sq1=square(d1(0))/qap;
00196   dfx+=dfd1(0)*sq1*qab;
00197   dfqab+=dfd1(0)*sq1*x;
00198   dfqap-=dfd1(0)*sq1*qab*x/qap;
00199   dfd1(0)=0.0;
00200 
00201   /*
00202   qab=a+b;
00203   qap=a+1.0;
00204   qam=a-1.0;
00205  */
00206 
00207   //qam=a-1.0;
00208   dfa+=dfqam;
00209 
00210   //qap=a+1.0;
00211   dfa+=dfqap;
00212 
00213   //qab=a+b;
00214   dfa+=dfqab;
00215   dfb+=dfqab;
00216 
00217 
00218   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation3ind,
00219     &(value(hh)) ,&(value(_a)),dfa ,&(value(_b)),dfb ,&(value(_x)),dfx);
00220 
00221 
00222   return hh;
00223 }
00224 
00225 #undef MAXIT
00226 #undef EPS
00227 #undef FPMIN