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