this article is intended to illustrate how to implement
a nontrivial iAD Model Builder vectorized function for 
GPU programming using opencl. Only first derivatives (dvariable type)
are considered.   The function to be vectorized is the
log of the negative binomial density.

The code for the log of the negative binomial density looks like

  double log_negbinomial_density(double x,double mu,double tau)
  {
    if (tau<=1.0)
    {
      cerr << "tau <=1 in log_negbinomial_density " << endl;
      ad_exit(1);
    }
    double r=mu/(1.e-120+(tau-1.0));
    return gammln(x+r)-gammln(r)-gammln(x+1)
      +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
  }

We see that most of the work is done by the gammln function
(the log of the gamma function.  So frist get some C code to
evaluate that.

double gammln(double xx)
{
  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;
  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) 
  {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}

This is simple moderately accurate version of the log of the gamma
function.  Keep in mind that it is more important that the derivative
of the function is accurate than it is that the function itself is
overly accurate i.e the derivative should accurate reflect the
function actually being calculated.

Now the gammln function has only one independent variable so that
the derviatives can be calculated easily by forward mode 
automatic differentiation (AD).
To make this simple we need to define a struct df11variable 
which will hold the value of variables and the derivative.  


struct df11variable
{
  double v[2]; // v[0] holds value, v[1] holds derivative
};

Here  is a first go at the version of gammln with forward
mode AD for the derivatives. Later we will need to interface
the derivatives with the negbin function. Note that so far 
nothing has been said about GPU stuff. The idea is that we
can produce correct working C code on the PC and then just
moldify it with a bit of boilerplate so that it gives us a working
vectorized  opencl function.

struct df11variable df11gammln(struct df11variable xx)
{
  struct  df11variable x,tmp,ser,z;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};

  xx.v[1]=1.0;   // initialize 
  int j;
  x.v[0]=xx.v[0]-1.0;
  x.v[1]=xx.v[1];

  tmp.v[1]=x.v[1];
  tmp.v[0]=x.v[0]+5.5;
  

  tmp.v[1] -= ( x.v[1]*log(tmp.v[0]) + (x.v[0]+0.5)/tmp.v[0]*tmp.v[1] );
  tmp.v[0] -= (x.v[0]+0.5)*log(tmp.v[0]);


  ser.v[1]=0.0;
  ser.v[0]=1.0;

  for (j=0;j<=5;j++) 
  {
    x.v[0] += 1.0;
    ser.v[1] -= cof[j]/(x.v[0]*x.v[0])*x.v[1];
    ser.v[0] += cof[j]/x.v[0];
  }
 

  z.v[0]= -tmp.v[0] + log(2.50662827465*ser.v[0]);
  z.v[1]= -tmp.v[1]+ 1.0/ser.v[0]*ser.v[1];

  return z;
}


If this code were written in C++ we would overload all the
functions to carry the derivative information along automatically,
but in C we need to do this by hand.  The idea is that if we 
have a line of code

     x.v[0] = f(z.v[0]);

we then include

     x.v[1] = f'(z.v[0])*z.v[1];

where f'(z.v[0]) is the derivative of f.
For a function of two variables such as

     x.v[0] = g(y.z[0],z.v[0]);

the corresponding derivative code is

     x.v[1] = g_1(y.z[0],z.v[0])*y.v[1] + g_2(y.z[0],z.v[0])*z.v[1];

where g_1 is the partial derivatvie of g wrt its first variable, y.

There is one gotcha here. for code where the same variable appears
on both sides you must do the derivative calculations first or the 
code wil be incorrect, so the line

  tmp -= (x+0.5)*log(tmp);

is replaced by

  tmp.v[1] -= ( x.v[1]*log(tmp.v[0]) + (x.v[0]+0.5)/tmp.v[0]*tmp.v[1] );
  tmp.v[0] -= (x.v[0]+0.5)*log(tmp.v[0]);






Testing the code

Here is the source code for testing the derivatives and verifying that the
function value is correct when compared to the double version of gammln
You shiould be able to compile this just by grabbing the stuff between the

      ***************************************

and putting it in a file named say testgamm.c


*****************************************************************
*****************************************************************
*****************************************************************
#include <math.h>
#include <stdio.h>

double gammln(double xx)
{
  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;
  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) 
  {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}

struct df11variable
{
  double v[2]; // v[0] holds value, v[1] holds derivative
};

struct df11variable df11gammln(struct df11variable xx)
{
  struct  df11variable x,tmp,ser,z;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};

  xx.v[1]=1.0;   // initialize 
  int j;
  x.v[0]=xx.v[0]-1.0;
  x.v[1]=xx.v[1];

  tmp.v[1]=x.v[1];
  tmp.v[0]=x.v[0]+5.5;
  

  tmp.v[1] -= ( x.v[1]*log(tmp.v[0]) + (x.v[0]+0.5)/tmp.v[0]*tmp.v[1] );
  tmp.v[0] -= (x.v[0]+0.5)*log(tmp.v[0]);


  ser.v[1]=0.0;
  ser.v[0]=1.0;

  for (j=0;j<=5;j++) 
  {
    x.v[0] += 1.0;
    ser.v[1] -= cof[j]/(x.v[0]*x.v[0])*x.v[1];
    ser.v[0] += cof[j]/x.v[0];
  }
 

  z.v[0]= -tmp.v[0] + log(2.50662827465*ser.v[0]);
  z.v[1]= -tmp.v[1]+ 1.0/ser.v[0]*ser.v[1];

  return z;
}

main()
{
  struct df11variable x,z,z1,z2;
  double xsave,diff,lfact,lcfact;
  double delta=1.e-6;
  do
  {
    printf("enter value for x\n");
    scanf("%lf",&(x.v[0]));
    lcfact=gammln(x.v[0]);
    z=df11gammln(x);
    xsave=x.v[0];
    x.v[0]=xsave+delta;
    z1=df11gammln(x);
    x.v[0]=xsave-delta;
    z2=df11gammln(x);
    diff=(z1.v[0]-z2.v[0])/(2.0*delta);
    lfact=z.v[0];
    printf("finit diff %lf   forward AD %lf df fvalue %lf df cvalue %lf \n",
       diff,z.v[1],lfact,lcfact);
  }
  while(1);
}

*****************************************************************
*****************************************************************
*****************************************************************

when the program is run it asks for a value of x to input. 
say you put in 6.18. 

then the output is

   finit diff 1.738236   forward AD 1.738236 df fvalue 5.097499 df cvalue 5.097499

You are seeing the finite difference derivative followed by the forward mode AD
derivative and then the value of the function for the double and df11variable
versions.

Extending the code to the log of the negative binomail density.

The negative binomial density is a function of two variables, the mean ,mu, and the overdispersion, tau.  To do forward model AD we need a struct which
holds the two partial derivatives.


 struct df12variable
 {  
   double v[3]; // v[0] holds value, 
                //v[1] holds derivative wrt first variable
                //v[2] holds derivative wrt second variable
 };

Here is the code. Note that while there are two independent variables,
most of the calculations are in calls to gammln whihc depends on only one
independent variable. This fact can be used to simplify the code in that
in the gammln part we only need to calculate one derivative.

 struct df12variablelog_negbinomial_density(double x,struct df12variable mu,
   struct df12variabledouble tau)
  {
    struct df12variable r;
    struct df12variable z;
    struct df11variable r1;
    struct df11variable r2;
    struct df11variable z1;
    struct df11variable z2;
    mu.v[1]=1.0;
    mu.v[2]=0.0;
    tau.v[1]=1.0;
    tau.v[2]=0.0;
    if (tau.v[0]<=1.0)
    {
      cerr << "tau <=1 in log_negbinomial_density " << endl;
      ad_exit(1);
    }
    r.v[0]=mu.v[0]/(1.e-120+(tau.v[0]-1.0));
    r.v[1]=mu.v[1]/(1.e-120+(tau.v[0]-1.0)); //der wrt mu
    r.v[2]=-mu.v[0]/((1.e-120+(tau.v[0]-1.0))*(1.e-120+(tau.v[0]-1.0)))*
              tau.v[2]; //der wrt tau
    // return gammln(x+r)-gammln(r)-gammln(x+1)
    //  +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
    r1.v[0]=r.v[0]+x;
    z1=df11gammln(r1); 
    r2.v[0]=r.v[0];
    z2=df11gammln(r2); 
    z.v[0]=0.0;
    z.v[1]=0.0;
    z.v[2]=0.0;
    // do it step by step
    z.v[0]+=z1.v[0]; // this is z+=gammln(x+r)
    z.v[1]+=z1.v[1]*r.v[1]; // chain rule for der wrt mu
    z.v[2]+=z1.v[1]*r.v[2]; // chain rule for der wrt tau

    z.v[0]-=z2.v[0]; // this is z-=gammln(r)
    z.v[1]-=z2.v[1]*r.v[1]; // chain rule for der wrt mu
    z.v[2]-=z2.v[1]*r.v[2]; // chain rule for der wrt tau

    z.v[0]-=gammln(x+1.0); // constant term so no ders

    z.v[0]+= r.v[0]*log(r.v[0]); // this is z+=r*log(r)
    z.v[1]+= (log(r.v[0])+1.0)*r.v[1]; // chain rule for der wrt mu
    z.v[2]+= (log(r.v[0])+1.0)*r.v[2]; // chain rule for der wrt tau

    z.v[0]+= x*log(mu.v[0]);
    z.v[1]+= x/mu.v[0]*mu.v[1];; // chain rule for der wrt mu
                               // note that mu.v[2]=0.0
    z.v[0]-=(r.v[0]+x)*log(r.v[0]+mu.v[0]);
    z.v[1]-=log(r.v[0]+mu.v[0])*r.v[1]
    z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[1];
    z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*mu.v[1];
    z.v[2]-=log(r.v[0]+mu.v[0])*r.v[2]
    z.v[2]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[2];
    return z;
  }

The complete code needed to test the derivatives for the
log_negatvie binomial is here.

*****************************************************************
*****************************************************************
*****************************************************************


#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double gammln(double xx)
{
  double x,tmp,ser;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};
  int j;
  x=xx-1.0;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.0;
  for (j=0;j<=5;j++) 
  {
    x += 1.0;
    ser += cof[j]/x;
  }
  return -tmp+log(2.50662827465*ser);
}
double log_negbinomial_density(double x,double mu,double tau)
{
  if (tau<=1.0)
  {
    printf("tau <=1 in log_negbinomial_density \n");
    exit(1);
  }
  double r=mu/(1.e-120+(tau-1.0));
  return gammln(x+r)-gammln(r)-gammln(x+1)
    +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
}



struct df11variable
{
  double v[2]; // v[0] holds value, v[1] holds derivative
};

struct df11variable df11gammln(struct df11variable xx)
{
  struct  df11variable x,tmp,ser,z;
  static double cof[6]={76.18009173,-86.50532033,24.01409822,
    -1.231739516,0.120858003e-2,-0.536382e-5};

  xx.v[1]=1.0;   // initialize 
  int j;
  x.v[0]=xx.v[0]-1.0;
  x.v[1]=xx.v[1];

  tmp.v[1]=x.v[1];
  tmp.v[0]=x.v[0]+5.5;
  

  tmp.v[1] -= ( x.v[1]*log(tmp.v[0]) + (x.v[0]+0.5)/tmp.v[0]*tmp.v[1] );
  tmp.v[0] -= (x.v[0]+0.5)*log(tmp.v[0]);


  ser.v[1]=0.0;
  ser.v[0]=1.0;

  for (j=0;j<=5;j++) 
  {
    x.v[0] += 1.0;
    ser.v[1] -= cof[j]/(x.v[0]*x.v[0])*x.v[1];
    ser.v[0] += cof[j]/x.v[0];
  }
 

  z.v[0]= -tmp.v[0] + log(2.50662827465*ser.v[0]);
  z.v[1]= -tmp.v[1]+ 1.0/ser.v[0]*ser.v[1];

  return z;
}
struct df12variable
{
  double v[3]; // v[0] holds value,
              //v[1] holds derivative wrt first variable
              //v[2] holds derivative wrt second variable
};



struct df12variable df12log_negbinomial_density(double x,struct df12variable mu,
  struct df12variable tau)
{
  struct df12variable r;
  struct df12variable z;
  struct df11variable r1;
  struct df11variable r2;
  struct df11variable z1;
  struct df11variable z2;
  mu.v[1]=1.0;
  mu.v[2]=0.0;
  tau.v[1]=0.0;
  tau.v[2]=1.0;
  if (tau.v[0]<=1.0)
  {
    printf("tau <=1 in log_negbinomial_density \n");
    exit(1);
  }
  r.v[0]=mu.v[0]/(1.e-120+(tau.v[0]-1.0));
  r.v[1]=mu.v[1]/(1.e-120+(tau.v[0]-1.0)); //der wrt mu
  r.v[2]=-mu.v[0]/((1.e-120+(tau.v[0]-1.0))*(1.e-120+(tau.v[0]-1.0)))*
            tau.v[2]; //der wrt tau
  // return gammln(x+r)-gammln(r)-gammln(x+1)
  //  +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
  r1.v[0]=r.v[0]+x;
  z1=df11gammln(r1); 
  r2.v[0]=r.v[0];
  z2=df11gammln(r2); 
  z.v[0]=0.0;
  z.v[1]=0.0;
  z.v[2]=0.0;
  // do it step by step
  z.v[0]+=z1.v[0]; // this is z+=gammln(x+r)
  z.v[1]+=z1.v[1]*r.v[1]; // chain rule for der wrt mu
  z.v[2]+=z1.v[1]*r.v[2]; // chain rule for der wrt tau

  z.v[0]-=z2.v[0]; // this is z-=gammln(r)
  z.v[1]-=z2.v[1]*r.v[1]; // chain rule for der wrt mu
  z.v[2]-=z2.v[1]*r.v[2]; // chain rule for der wrt tau

  z.v[0]-=gammln(x+1.0); // constant term so no ders

  z.v[0]+= r.v[0]*log(r.v[0]); // this is z+=r*log(r)
  z.v[1]+= (log(r.v[0])+1.0)*r.v[1]; // chain rule for der wrt mu
  z.v[2]+= (log(r.v[0])+1.0)*r.v[2]; // chain rule for der wrt tau

  z.v[0]+= x*log(mu.v[0]);
  z.v[1]+= x/mu.v[0]*mu.v[1];; // chain rule for der wrt mu
                             // note that mu.v[2]=0.0
  z.v[0]-=(r.v[0]+x)*log(r.v[0]+mu.v[0]);
  z.v[1]-=log(r.v[0]+mu.v[0])*r.v[1];
  z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[1];
  z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*mu.v[1];
  z.v[2]-=log(r.v[0]+mu.v[0])*r.v[2];
  z.v[2]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[2];
  return z;
}

main()
{
  struct df12variable x,y,z,z1,z2,zz1,zz2;
  double xsave,ysave,diff,diffy,lfact,lcfact;
  double delta=1.e-6;
  //do
  {
    x.v[0]=4.2;
    x.v[1]=0.0;
    x.v[2]=1.0;
    y.v[0]=9.0;
    y.v[1]=1.0;
    y.v[2]=0.0;
    //printf("enter value for x\n");
    //scanf("%lf",&(x.v[0]));
    lcfact=log_negbinomial_density(12.0,x.v[0],y.v[0]);
    z=df12log_negbinomial_density(12.0,x,y);
    xsave=x.v[0];
    x.v[0]=xsave+delta;
    z1=df12log_negbinomial_density(12.0,x,y);
    x.v[0]=xsave-delta;
    z2=df12log_negbinomial_density(12.0,x,y);
    diff=(z1.v[0]-z2.v[0])/(2.0*delta);
    lfact=z.v[0];
    x.v[0]=xsave;

    ysave=y.v[0];
    y.v[0]=ysave+delta;
    zz1=df12log_negbinomial_density(12.0,x,y);
    y.v[0]=ysave-delta;
    zz2=df12log_negbinomial_density(12.0,x,y);
    diffy=(zz1.v[0]-zz2.v[0])/(2.0*delta);
    lfact=z.v[0];
    y.v[0]=ysave;


    printf("finit diff %lf   forward AD %lf df fvalue %lf df cvalue %lf \n",
       diff,z.v[1],lfact,lcfact);
    printf("finit diff %lf   forward AD %lf df fvalue %lf df cvalue %lf \n",
       diffy,z.v[2],lfact,lcfact);
  }
  //while(1);
}

*****************************************************************
*****************************************************************
*****************************************************************

The next step is to vecotrize the function and put it in a form
suitable for running on the GPU. we do this in a TPL file
putting the new code in the GLOBALS_SECTION. This will allow
simple derivative checking.
Here is the TPL file
 
*****************************************************************
*****************************************************************
*****************************************************************

DATA_SECTION
  int n
 !! n=5;
  vector obs(1,n)
 !! obs.fill_seqadd(1,1);

PARAMETER_SECTION
  init_vector mu(1,n)
  init_vector tau(1,n)
  objective_function_value f
PROCEDURE_SECTION
  dvar_vector z=my_log_negbinomial_density(obs,mu,tau);
  f=sum(z);
GLOBALS_SECTION

  #include <admodel.h>
 extern "C"
 {
 
 double mygammln(double xx)
 {
   double x,tmp,ser;
   static double cof[6]={76.18009173,-86.50532033,24.01409822,
     -1.231739516,0.120858003e-2,-0.536382e-5};
   int j;
   x=xx-1.0;
   tmp=x+5.5;
   tmp -= (x+0.5)*log(tmp);
   ser=1.0;
   for (j=0;j<=5;j++) 
   {
     x += 1.0;
     ser += cof[j]/x;
   }
   return -tmp+log(2.50662827465*ser);
 }
 double mylog_negbinomial_density(double x,double mu,double tau)
 {
   if (tau<=1.0)
   {
     printf("tau <=1 in log_negbinomial_density \n");
     exit(1);
   }
   double r=mu/(1.e-120+(tau-1.0));
   return mygammln(x+r)-gammln(r)-gammln(x+1)
     +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
 }
 

  struct df11variable
  {
    double v[2]; // v[0] holds value, v[1] holds derivative
  };
  
  struct df11variable df11gammln(struct df11variable xx)
  {
    struct  df11variable x,tmp,ser,z;
    static double cof[6]={76.18009173,-86.50532033,24.01409822,
      -1.231739516,0.120858003e-2,-0.536382e-5};
  
    xx.v[1]=1.0;   // initialize 
    int j;
    x.v[0]=xx.v[0]-1.0;
    x.v[1]=xx.v[1];
  
    tmp.v[1]=x.v[1];
    tmp.v[0]=x.v[0]+5.5;
    
  
    tmp.v[1] -= ( x.v[1]*log(tmp.v[0]) + (x.v[0]+0.5)/tmp.v[0]*tmp.v[1] );
    tmp.v[0] -= (x.v[0]+0.5)*log(tmp.v[0]);
  
  
    ser.v[1]=0.0;
    ser.v[0]=1.0;
  
    for (j=0;j<=5;j++) 
    {
      x.v[0] += 1.0;
      ser.v[1] -= cof[j]/(x.v[0]*x.v[0])*x.v[1];
      ser.v[0] += cof[j]/x.v[0];
    }
   
  
    z.v[0]= -tmp.v[0] + log(2.50662827465*ser.v[0]);
    z.v[1]= -tmp.v[1]+ 1.0/ser.v[0]*ser.v[1];
  
    return z;
  }
  struct df12variable
  {
    double v[3]; // v[0] holds value,
                //v[1] holds derivative wrt first variable
                //v[2] holds derivative wrt second variable
  };
  
  
  
  struct df12variable df12log_negbinomial_density(double x,struct df12variable mu,
    struct df12variable tau)
  {
    struct df12variable r;
    struct df12variable z;
    struct df11variable r1;
    struct df11variable r2;
    struct df11variable z1;
    struct df11variable z2;
    mu.v[1]=1.0;
    mu.v[2]=0.0;
    tau.v[1]=0.0;
    tau.v[2]=1.0;
    if (tau.v[0]<=1.0)
    {
      printf("tau <=1 in log_negbinomial_density \n");
      exit(1);
    }
    r.v[0]=mu.v[0]/(1.e-120+(tau.v[0]-1.0));
    r.v[1]=mu.v[1]/(1.e-120+(tau.v[0]-1.0)); //der wrt mu
    r.v[2]=-mu.v[0]/((1.e-120+(tau.v[0]-1.0))*(1.e-120+(tau.v[0]-1.0)))*
              tau.v[2]; //der wrt tau
    // return mygammln(x+r)-gammln(r)-gammln(x+1)
    //  +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
    r1.v[0]=r.v[0]+x;
    z1=df11gammln(r1); 
    r2.v[0]=r.v[0];
    z2=df11gammln(r2); 
    z.v[0]=0.0;
    z.v[1]=0.0;
    z.v[2]=0.0;
    // do it step by step
    z.v[0]+=z1.v[0]; // this is z+=gammln(x+r)
    z.v[1]+=z1.v[1]*r.v[1]; // chain rule for der wrt mu
    z.v[2]+=z1.v[1]*r.v[2]; // chain rule for der wrt tau
  
    z.v[0]-=z2.v[0]; // this is z-=gammln(r)
    z.v[1]-=z2.v[1]*r.v[1]; // chain rule for der wrt mu
    z.v[2]-=z2.v[1]*r.v[2]; // chain rule for der wrt tau
  
    z.v[0]-=mygammln(x+1.0); // constant term so no ders
  
    z.v[0]+= r.v[0]*log(r.v[0]); // this is z+=r*log(r)
    z.v[1]+= (log(r.v[0])+1.0)*r.v[1]; // chain rule for der wrt mu
    z.v[2]+= (log(r.v[0])+1.0)*r.v[2]; // chain rule for der wrt tau
  
    z.v[0]+= x*log(mu.v[0]);
    z.v[1]+= x/mu.v[0]*mu.v[1];; // chain rule for der wrt mu
                               // note that mu.v[2]=0.0
    z.v[0]-=(r.v[0]+x)*log(r.v[0]+mu.v[0]);
    z.v[1]-=log(r.v[0]+mu.v[0])*r.v[1];
    z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[1];
    z.v[1]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*mu.v[1];
    z.v[2]-=log(r.v[0]+mu.v[0])*r.v[2];
    z.v[2]-=(r.v[0]+x)/(r.v[0]+mu.v[0])*r.v[2];
    return z;
  }
  
  } // extern "C"


  void DF_dvdvnegbin(void);
  
    dvar_vector my_log_negbinomial_density(dvector& x,dvar_vector& mu,dvar_vector& tau){
      int mmin=x.indexmin();
      int mmax=x.indexmax();
      if (mmin !=mu.indexmin() || mmin != tau.indexmin() ||
          mmax !=tau.indexmax() || mmax != tau.indexmax())
      {
        cerr << "shape error in log_negbinomial_density" << endl;
        ad_exit(1);
      }
      dvar_vector tmp(mmin,mmax);
      dvector dfx(mmin,mmax);
      dvector dfy(mmin,mmax);
      if (USE_GPU==0)
      {
        for (int i=mmin;i<=mmax;i++)
        {
          struct df12variable z,xmu,xtau;
          xmu.v[0]=value(mu(i));
          xtau.v[0]=value(tau(i));
          z=df12log_negbinomial_density(x(i),xmu,xtau);
          value(tmp(i))=z.v[0];
          dfx(i)=z.v[1];
          dfy(i)=z.v[2];
        }
      }
      else
      {
        dvector cmu=value(mu);
        dvector ctau=value(tau);
        ocl_vec_negbinomial_density(&(x(mmin)),&(cmu(mmin)),&(ctau(mmin)),
          &(dfx(mmin)),&(dfy(mmin)));
      }
  
      save_identifier_string("tg");
      tmp.save_dvar_vector_position();
      mu.save_dvar_vector_position();
      tau.save_dvar_vector_position();
      save_identifier_string("ez");
      dfx.save_dvector_value();
      dfx.save_dvector_position();
      dfy.save_dvector_value();
      dfy.save_dvector_position();
      save_identifier_string("ey");
      gradient_structure::GRAD_STACK1->
                set_gradient_stack(DF_dvdvnegbin);
      return tmp; 
    }
    dvar_vector ocl_vec_negbinomial_density(dvector& x,dvar_vector& mu,dvar_vector& tau){
      int mmin=x.indexmin();
      int mmax=x.indexmax();
      if (mmin !=mu.indexmin() || mmin != tau.indexmin() ||
          mmax !=tau.indexmax() || mmax != tau.indexmax())
      {
        cerr << "shape error in log_negbinomial_density" << endl;
        ad_exit(1);
      }
      dvar_vector tmp(mmin,mmax);
      dvector dfx(mmin,mmax);
      dvector dfy(mmin,mmax);
  
      dvector cmu=value(mu);
      dvector ctau=value(tau);
      ocl_vec_negbinomial_density_kernel(&(x(mmin)),&(cmu(mmin)),&(ctau(mmin)),
        &(dfx(mmin)),&(dfy(mmin)));

    save_identifier_string("tg");
    tmp.save_dvar_vector_position();
    mu.save_dvar_vector_position();
    tau.save_dvar_vector_position();
    save_identifier_string("ez");
    dfx.save_dvector_value();
    dfx.save_dvector_position();
    dfy.save_dvector_value();
    dfy.save_dvector_position();
    save_identifier_string("ey");
    gradient_structure::GRAD_STACK1->
              set_gradient_stack(DF_dvdvnegbin);
    return tmp; 
  }
  
  void DF_dvdvnegbin(void)
  {
    // int ierr=fsetpos(gradient_structure::get_fp(),&filepos);
    verify_identifier_string("ey");
    dvector_position dfy_pos=restore_dvector_position();
    dvector dfy=restore_dvector_value(dfy_pos);
    dvector_position dfx_pos=restore_dvector_position();
    dvector dfx=restore_dvector_value(dfx_pos);
    verify_identifier_string("ez");
  
    dvar_vector_position tau_pos=restore_dvar_vector_position();
    dvar_vector_position mu_pos=restore_dvar_vector_position();
    dvar_vector_position tmp_pos=restore_dvar_vector_position();
  
    verify_identifier_string("tg");
    dvector dfvtmp=restore_dvar_vector_derivatives(tmp_pos);
  
    dvector dfmu(dfvtmp.indexmin(),dfvtmp.indexmax());
    dvector dftau(dfvtmp.indexmin(),dfvtmp.indexmax());
  
    for (int i=dfvtmp.indexmin();i<=dfvtmp.indexmax();i++)
    {
      dfmu(i)=dfx(i)*dfvtmp(i);
      dftau(i)=dfy(i)*dfvtmp(i);
    }
    dfmu.save_dvector_derivatives(mu_pos);
    dftau.save_dvector_derivatives(tau_pos);
  }