#include<math.h>

void llclust (double *Y, double *lambda, double *theta,  double *a, 
               int *n, int *m, double *sm) 
{
   int i, j, n0, n1;
   double lam;
   for(j = 0 ; j < *m ; j++ ) {

      n0=0; n1=0;  
      for(i=0 ; i< *n ; i++ )     { 
         if      (Y[ i*(*m) + j ]==1 ) {n1++; }
         else if (Y[ i*(*m) + j ]==0 ) {n0++; }
                                   }
      lam= lgamma(a[0]+a[1]) - lgamma(a[0]) -lgamma(a[1]) +
           lgamma(a[0]+n1) + lgamma(a[1]+n0) - lgamma(a[0]+a[1]+n1+n0) - 
           n1*log(theta[j]) - n0*log((1.0-theta[j])) ;

      *sm= *sm + lam +log( exp(-lam)+exp(lambda[j] ))- log(1.0+exp( lambda[j])) ;
                              }
}