matrixjpl.c

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

#include <math.h>
#include <malloc.h>
#include <stddef.h>
#include <stdio.h>
#include "matrixjpl.h"

/*-------------------------------------------------------------------------*/

void vextract(int n, int col, double **mat, double *vec)
{ 
	int i;
	
        for(i=0; i<n; i++)
		vec[i] = mat[i][col];
}

/*-------------------------------------------------------------------------*/


void vinsert(int n, int col, double *vec, double **mat)
{ 
	int i;

	for(i=0; i<n; i++)
		mat[i][col] = vec[i];
}

/*-------------------------------------------------------------------------*/


void ivextract(int n, int col, int **mat, int *vec)
{ 
	int i;

	for(i=0; i<n; i++)
		vec[i] = mat[i][col];
}

/*-------------------------------------------------------------------------*/



void vcopy(int n, double *vec1, double *vec2)
{ 
	int i;

	for(i=0; i<n; i++)
		vec2[i] = vec1[i];
}

/*-------------------------------------------------------------------------*/



void vprint(int n, double *vec)
{
int i;

for(i=0; i<n; i++)
printf(" %8.3lf \n",vec[i]);

printf("\n");

}


/*-------------------------------------------------------------------------*/


void ivprint(int n, int *vec)
{
int i;

for(i=0; i<n; i++)
printf(" %8d \n",vec[i]);

printf("\n");

}


/*-------------------------------------------------------------------------*/


void mprint(int n, int m, double **vec)
{
int i, j;

for(i=0; i<n; i++){
	for(j=0; j<m; j++)
	printf(" %8.3lf ",vec[i][j]);
	printf(" \n");
}

printf("\n");


}


/*-------------------------------------------------------------------------*/

void imprint(int n, int m, int **vec)
{
int i, j;

for(i=0; i<n; i++){
	for(j=0; j<m; j++)
	printf(" %8d ",vec[i][j]);
	printf(" \n");
}

printf("\n");

}

/*-------------------------------------------------------------------------*/


/* In-place Permutations 

   permute:    OUT[i]       = IN[perm[i]]     i = 0 .. N-1
   invpermute: OUT[perm[i]] = IN[i]           i = 0 .. N-1

   PERM is an index map, i.e. a vector which contains a permutation of
   the integers 0 .. N-1.

   From Knuth "Sorting and Searching", Volume 3 (3rd ed), Section 5.2
   Exercise 10 (answers), p 617
*/

void permute(int n, int *p, double *data)
{
  int i, k, pk;
  double r1, t;
  
  for (i = 0; i < n; i++)
    {
      k = p[i];
      
      while (k > i) 
        k = p[k];
      
      if (k < i)
        continue ;
      
      /* Now have k == i, i.e the least in its cycle */
      
      pk = p[k];
      
      if (pk == i)
        continue ;
      
      /* shuffle the elements of the cycle */
      
      t = data[i];
      
      while (pk != i)
        {
          r1 = data[pk];
          data[k] = r1;
          k = pk;
          pk = p[k];
        };
      
      data[k] = t;
    }
}


/*-------------------------------------------------------------------------*/

void permute_inverse (int n, int *p, double *data)
{
  int i, k, pk;

  double r1, t;

  for (i = 0; i < n; i++)
    {
      k = p[i];
          
      while (k > i) 
        k = p[k];

      if (k < i)
        continue ;
      
      /* Now have k == i, i.e the least in its cycle */

      pk = p[k];

      if (pk == i)
        continue ;
      
      /* shuffle the elements of the cycle in the inverse direction */
      
      t = data[k];

      while (pk != i)
        {
          r1 = data[pk];
          data[pk] = t;
          k = pk;
          pk = p[k];
          t = r1;
        };
      
      data[pk] = t;
    }

}


/*-------------------------------------------------------------------------*/


void matmulc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{ int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]*y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] * y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] * x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] * y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] * x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("matmulc: nonconformability of matrices \n");
   }

}


/*-------------------------------------------------------------------------*/



void matdivc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{ int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]/y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] / y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] / x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] / y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] / x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("matdivc: nonconformability of matrices \n");
   }

}


/*-------------------------------------------------------------------------*/



void mataddc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{ int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i]+y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] + y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] + x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] + y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] + x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("mataddc: nonconformability of matrices \n");
   }

}

/*-------------------------------------------------------------------------*/


void matsubc(double z[], double x[], double y[], int size,
             int m1, int m2, int n1, int n2)
{ int i, j;
  int flag = 0;
  
  /* case of two equal matrices */
  if (m1 == m2  &&  n1 == n2){
  for (i=0; i <= size-1; i++){
  z[i] = x[i] - y[i];
  }
    flag = 1;
  }
  if (m1 == m2 && n2 == 1) {
  size = 0;
  for (j=0; j<= n1-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = x[size] - y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (m1 == m2 && n1 == 1) {
  size = 0;
  for (j=0; j<= n2-1; j++){
  for (i=0; i<= m1-1; i++){
  z[size] = y[size] - x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  if (n1 == n2 && m2 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m1-1; j++){
  z[size] = x[size] - y[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
  if (n1 == n2 && m1 == 1) {
  size = 0;
  for (i=0; i<= n1-1; i++){
  for (j=0; j<= m2-1; j++){
  z[size] = y[size] - x[i];
  size = size + 1;
  }
  }
      flag = 1;
  }
  
   if (flag == 0){
      printf("matsubc: nonconformability of matrices \n");
   }

}


/*-------------------------------------------------------------------------*/


void meanvar(double data[], int n, double *ave, double *svar)
{ 
int j;
double s;

*ave = (*svar) = 0.0;
for(j=0; j<n; j++) *ave += data[j];
*ave/= n;
for(j=0; j<n; j++){
   s = data[j]-(*ave);
   *svar += s*s;
   }
   *svar /= (n-1);
}

NUMERICS_EXPORT BOOL choldcmp(double **a, int n, double *p)
{
    int i, j, k;
    double sum;
    
    *p = sqrt(a[0][0]);

    for(i = 1; i < n; i++)
	{
    
		a[i][0] /= *p;
    
	}
    
	for(i = 1; i < n; i++)
	{
        
		sum = a[i][i];
        
		for(k = 0; k < i; ++k)
		
		{
            
			sum -= a[i][k] * a[i][k];
        }


        if(sum <= 0.0)
		
		{
            return FALSE;