ADMB Documentation  11.1.2246
 All Classes Files Functions Variables Typedefs Friends Defines
combv.cpp
Go to the documentation of this file.
00001 
00011 #include <fvar.hpp>
00012 
00013 dvariable gammln(const dvariable& xx);
00014 dvariable factln(const dvariable& n);
00015 double factln(double n);
00016 
00024 dvariable log_comb(const dvariable& n, double k)
00025 {
00026   return factln(n)-factln(k)-factln(n-k);
00027 }
00028 
00036 dvariable log_comb(const dvariable& n, const dvariable& k)
00037 {
00038   return factln(n)-factln(k)-factln(n-k);
00039 }
00040 
00041 
00049 dvariable log_comb(double n, const dvariable& k)
00050 {
00051   return factln(n)-factln(k)-factln(n-k);
00052 }
00053 
00059 dvariable factln(const dvariable& n)
00060 {
00061   return gammln(n+1.0);
00062 }
00063 
00070 dvar_vector log_comb(const dvar_vector& n, const dvar_vector& r)
00071 {
00072     RETURN_ARRAYS_INCREMENT();
00073     int mmin=n.indexmin();
00074     int mmax=n.indexmax();
00075     if (mmin != r.indexmin() || mmax != r.indexmax())
00076     {
00077       cerr << "Incompatible array bounds in function "
00078         "dvar_vector log_comb(const dvar_vector& n, const dvector& r)" << endl;
00079       ad_exit(1);
00080     }
00081     dvar_vector tmp(mmin,mmax);
00082     for (int i=mmin;i<=mmax;i++)
00083     {
00084       tmp(i)=log_comb(n(i),r(i));
00085     }
00086     RETURN_ARRAYS_DECREMENT();
00087     return tmp;
00088 }
00089 
00096 dvar_vector log_comb(const dvector& n, const dvar_vector& r)
00097 {
00098     RETURN_ARRAYS_INCREMENT();
00099     int mmin=n.indexmin();
00100     int mmax=n.indexmax();
00101     if (mmin != r.indexmin() || mmax != r.indexmax())
00102     {
00103       cerr << "Incompatible array bounds in function "
00104         "dvar_vector log_comb(const dvar_vector& n, const dvector& r)" << endl;
00105       ad_exit(1);
00106     }
00107     dvar_vector tmp(mmin,mmax);
00108     for (int i=mmin;i<=mmax;i++)
00109     {
00110       tmp(i)=log_comb(n(i),r(i));
00111     }
00112     RETURN_ARRAYS_DECREMENT();
00113     return tmp;
00114 }
00115 
00122 dvar_vector log_comb(const dvar_vector& n, const dvector& r)
00123 {
00124     RETURN_ARRAYS_INCREMENT();
00125     int mmin=n.indexmin();
00126     int mmax=n.indexmax();
00127     if (mmin != r.indexmin() || mmax != r.indexmax())
00128     {
00129       cerr << "Incompatible array bounds in function "
00130         "dvar_vector log_comb(const dvar_vector& n, const dvector& r)" << endl;
00131       ad_exit(1);
00132     }
00133     dvar_vector tmp(mmin,mmax);
00134     for (int i=mmin;i<=mmax;i++)
00135     {
00136       tmp(i)=log_comb(n(i),r(i));
00137     }
00138     RETURN_ARRAYS_DECREMENT();
00139     return tmp;
00140 }
00141 
00148 dvar_vector log_comb(const dvariable& n, const dvector& r)
00149 {
00150     RETURN_ARRAYS_INCREMENT();
00151     int mmin=r.indexmin();
00152     int mmax=r.indexmax();
00153     dvar_vector tmp(mmin,mmax);
00154     for (int i=mmin;i<=mmax;i++)
00155     {
00156       tmp(i)=log_comb(n,r(i));
00157     }
00158     RETURN_ARRAYS_DECREMENT();
00159     return tmp;
00160 }
00161 
00168 dvar_vector log_comb(const dvariable& n, const dvar_vector& r)
00169 {
00170     RETURN_ARRAYS_INCREMENT();
00171     int mmin=r.indexmin();
00172     int mmax=r.indexmax();
00173     dvar_vector tmp(mmin,mmax);
00174     for (int i=mmin;i<=mmax;i++)
00175     {
00176       tmp(i)=log_comb(n,r(i));
00177     }
00178     RETURN_ARRAYS_DECREMENT();
00179     return tmp;
00180 }
00181 
00187 dvar_vector factln(const dvar_vector& r)
00188 {
00189     RETURN_ARRAYS_INCREMENT();
00190     int mmin=r.indexmin();
00191     int mmax=r.indexmax();
00192     dvar_vector tmp(mmin,mmax);
00193     for (int i=mmin;i<=mmax;i++)
00194     {
00195       tmp(i)=gammln(r(i)+1.0);
00196     }
00197     RETURN_ARRAYS_DECREMENT();
00198     return tmp;
00199 }
00200 
00207 dvar_vector log_comb(double n, const dvar_vector& r)
00208 {
00209     RETURN_ARRAYS_INCREMENT();
00210     int mmin=r.indexmin();
00211     int mmax=r.indexmax();
00212     dvar_vector tmp(mmin,mmax);
00213     for (int i=mmin;i<=mmax;i++)
00214     {
00215       tmp(i)=log_comb(n,r(i));
00216     }
00217     RETURN_ARRAYS_DECREMENT();
00218     return tmp;
00219 }
00220 
00226 dvar_vector gammln(const dvar_vector& v)
00227 {
00228   RETURN_ARRAYS_INCREMENT();
00229   int mmin=v.indexmin();
00230   int mmax=v.indexmax();
00231   dvar_vector tmp(mmin,mmax);
00232   for (int i=mmin;i<=mmax;i++)
00233   {
00234     tmp(i)=gammln(v(i));
00235   }
00236   RETURN_ARRAYS_DECREMENT();
00237   return tmp;
00238 }
00239 
00240 static dvariable gammlnguts(const prevariable _z)
00241 {
00242   double  z = value(_z);
00243   //double zdot=1.0;
00244   //const double lpi =1.1447298858494001741434272;
00245   //const double pi =3.1415926535897932384626432;
00246   const double lpp =0.9189385332046727417803297;
00247   int n=7;
00248   const double c[9]={0.99999999999980993,
00249     676.5203681218851,
00250     -1259.1392167224028,
00251      771.32342877765313,
00252     -176.61502916214059,
00253     12.507343278686905,
00254      -0.13857109526572012,
00255     9.9843695780195716e-6,
00256     1.5056327351493116e-7};
00257   z-=1.0;
00258   double x=c[0];
00259   double xdot=0.0;
00260   for (int i=1;i<=n+1;i++)
00261   {
00262     double zinv=1.0/(z+i);
00263     x+=c[i]*zinv;
00264     //xdot-=c[i]/square(z+i)*zdot;  since zdot=1.0
00265     xdot-=c[i]*square(zinv);
00266   }
00267   double t=z+n+0.5;
00268   //return lpp + (z+0.5)*log(t) -t + log(x);
00269   double ans= lpp + (z+0.5)*log(t) -t + log(x);
00270   //double tdot=zdot;
00271   //double ansdot=zdot*log(t) + (z+0.5)/t*tdot -tdot +xdot/x;
00272   // since tdot=1.0
00273   // since zdot=1.0
00274   double ansdot=log(t) + (z+0.5)/t -1.0 +xdot/x;
00275   dvariable u;
00276   u.v->x=ans;
00277   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation,
00278     &(u.v->x), &(_z.v->x), ansdot );
00279   return(u);
00280 }
00281 
00282 dvariable gammln(const prevariable& z)
00283 {
00284   const double lpi =1.1447298858494001741434272;
00285   const double pi =3.1415926535897932384626432;
00286   if (z<0.5)
00287   {
00288     return lpi - log(sin(pi*z)) - gammlnguts(1.0-z);
00289   }
00290   else
00291   {
00292     return gammlnguts(z);
00293   }
00294 }