dirirnd.c

Fit a multivariate gaussian mixture by a cross-entropy method. Cross-Entropy is a powerfull tool to

#include "mex.h"
#include "randomlib.h"





/* dirirnd mex file 

   usage in matlab  D = dirirnd(A);
 
   Input 

          A      = (vector (N x M1,...,Mn))
  
   Output
  
         D       = Dirichlet samples

   Author : S閎astien PARIS ([e-mail])
 
   Example 
   
   A = ceil(3*rand(3 , 3 , 2)) + 1;
   D = dirirnd(A);
   sum(D)

   In matlab command, type  mex -output dirirnd.dll dirirnd.c matrixjpl.c randomlib.c to compile
 
*/
        

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{

  double *D, *A; 
	  
  register double sumD , invsum;

  const double one = 1.0, zero = 0.0;

  int *dimsA;

  int i, j , t , v , N , M=1 , numdimsA;
  
  /* Check input */

  if(nrhs != 1)
	  
	{

     mexErrMsgTxt("D = dirichlet(A): need (N x M1 , ... x Mn) entry.");
	
	}


    A         = mxGetPr(prhs[0]);

	numdimsA  = mxGetNumberOfDimensions(prhs[0]);

	dimsA     = mxGetDimensions(prhs[0]);

	N         = dimsA[0];

	for (i=1 ; i<numdimsA ; i++)
	{
        M *= dimsA[i];
	}
      

 
/* Output */


	plhs[0] = mxCreateNumericArray(numdimsA, dimsA , mxDOUBLE_CLASS, mxREAL);
	
    D       = mxGetPr(plhs[0]);


    for(j = 0 ; j < M ; j++)
	
	{

     t       = j*N;

	 sumD    = zero;

     for(i = 0; i<N; i++)
		 
	 {
         
		 v         = i + t;
      
		 D[v]      = sgamma(A[v]);
	  
		 sumD     += D[v];

	 }

	 invsum = one/sumD;

	 for(i = 0 ; i < N ; i++) 

	 {
		 D[i + t] *=  invsum;
	 }
	
	}



}