3053 lines
		
	
	
		
			69 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
			
		
		
	
	
			3053 lines
		
	
	
		
			69 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
/*************************************************************************
 | 
						|
*									 *
 | 
						|
*	 YAP Prolog 							 *
 | 
						|
*									 *
 | 
						|
*	Yap Prolog was developed at NCCUP - Universidade do Porto	 *
 | 
						|
*									 *
 | 
						|
* Copyright L.Damas, V.S.Costa and Universidade do Porto 1985-1997	 *
 | 
						|
*									 *
 | 
						|
**************************************************************************
 | 
						|
*									 *
 | 
						|
* File:		matrix.c						 *
 | 
						|
* Last rev:								 *
 | 
						|
* mods:									 *
 | 
						|
* comments:	numerical arrays		                         *
 | 
						|
*									 *
 | 
						|
*************************************************************************/
 | 
						|
 | 
						|
#include "config.h"
 | 
						|
#include "YapInterface.h"
 | 
						|
#include <math.h>
 | 
						|
#if defined(__MINGW32__) || _MSC_VER
 | 
						|
#include <windows.h>
 | 
						|
#endif
 | 
						|
#if HAVE_STRING_H
 | 
						|
#include <string.h>
 | 
						|
#endif
 | 
						|
 | 
						|
/*
 | 
						|
  A matrix is something of the form
 | 
						|
 | 
						|
  TYPE = {int,double}
 | 
						|
  #DIMS = an int
 | 
						|
  DIM1
 | 
						|
  ...
 | 
						|
  DIMn
 | 
						|
  DATA in C format.
 | 
						|
 | 
						|
  floating point matrixes may need to be aligned, so we always have an
 | 
						|
  extra element at the end.
 | 
						|
*/
 | 
						|
 | 
						|
/* maximal number of dimensions, 1024 should be enough */
 | 
						|
#define MAX_DIMS 1024
 | 
						|
 | 
						|
typedef enum {
 | 
						|
  INT_MATRIX,
 | 
						|
  FLOAT_MATRIX
 | 
						|
} mat_data_type;
 | 
						|
 | 
						|
typedef enum {
 | 
						|
  MAT_TYPE=0,
 | 
						|
  MAT_NDIMS=1,
 | 
						|
  MAT_SIZE=2,
 | 
						|
  MAT_ALIGN=3,
 | 
						|
  MAT_DIMS=4,
 | 
						|
} mat_type;
 | 
						|
 | 
						|
typedef enum {
 | 
						|
  MAT_PLUS=0,
 | 
						|
  MAT_SUB=1,
 | 
						|
  MAT_TIMES=2,
 | 
						|
  MAT_DIV=3,
 | 
						|
  MAT_IDIV=4,
 | 
						|
  MAT_ZDIV=5,
 | 
						|
  MAT_LOG=6,
 | 
						|
  MAT_EXP=7
 | 
						|
} op_type;
 | 
						|
 | 
						|
static long int *
 | 
						|
matrix_long_data(int *mat, int ndims)
 | 
						|
{
 | 
						|
  return (long int *)(mat+(MAT_DIMS+ndims));
 | 
						|
}
 | 
						|
 | 
						|
static double *
 | 
						|
matrix_double_data(int *mat, int ndims)
 | 
						|
{
 | 
						|
  return (double *)(mat+(MAT_DIMS+ndims));
 | 
						|
}
 | 
						|
 | 
						|
static unsigned int
 | 
						|
matrix_get_offset(int *mat, int* indx)
 | 
						|
{
 | 
						|
  unsigned int i, pos = mat[MAT_SIZE], off = 0;
 | 
						|
 | 
						|
  /* find where we are */
 | 
						|
  for (i = 0; i < mat[MAT_NDIMS]; i++) {
 | 
						|
    pos /= mat[MAT_DIMS+i];
 | 
						|
    if (indx[i] >= mat[MAT_DIMS+i]) {
 | 
						|
      return off;
 | 
						|
    }
 | 
						|
    off += pos*indx[i];
 | 
						|
  }
 | 
						|
  return off;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_get_index(int *mat, unsigned int offset, int* indx)
 | 
						|
{
 | 
						|
  unsigned int i, pos = mat[MAT_SIZE];
 | 
						|
 | 
						|
  /* find where we are */
 | 
						|
  for (i = 0; i < mat[MAT_NDIMS]; i++) {
 | 
						|
    pos /= mat[MAT_DIMS+i];
 | 
						|
    indx[i] = offset / pos;
 | 
						|
    offset = offset % pos;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_next_index(int *dims, int ndims, int* indx)
 | 
						|
{
 | 
						|
  unsigned int i;
 | 
						|
 | 
						|
  /* find where we are */
 | 
						|
  for (i = ndims; i >0; ) {
 | 
						|
    i--;
 | 
						|
    indx[i]++;
 | 
						|
    if (indx[i]!=dims[i]) return;
 | 
						|
    indx[i] = 0;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
new_int_matrix(int ndims, int dims[], long int data[])
 | 
						|
{
 | 
						|
  unsigned int sz;
 | 
						|
  unsigned int i, nelems=1;
 | 
						|
  YAP_Term blob;
 | 
						|
  int *mat;
 | 
						|
  long int *bdata;
 | 
						|
  int idims[MAX_DIMS];
 | 
						|
 | 
						|
  /* in case we don't have enough room and need to shift the stack, we can't
 | 
						|
     really afford to keep a pointer to the global */
 | 
						|
  for (i=0;i< ndims;i++) {
 | 
						|
    idims[i] = dims[i];
 | 
						|
    nelems *= dims[i];
 | 
						|
  }
 | 
						|
  sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+nelems*sizeof(long int))/sizeof(YAP_CELL);
 | 
						|
  blob = YAP_MkBlobTerm(sz);
 | 
						|
  if (blob == YAP_TermNil())  {
 | 
						|
    return blob;
 | 
						|
  }
 | 
						|
  mat = (int *)YAP_BlobOfTerm(blob);
 | 
						|
  mat[MAT_TYPE] = INT_MATRIX;
 | 
						|
  mat[MAT_NDIMS] = ndims;
 | 
						|
  mat[MAT_SIZE] = nelems;
 | 
						|
  for (i=0;i< ndims;i++) {
 | 
						|
    mat[MAT_DIMS+i] = idims[i];
 | 
						|
  }
 | 
						|
  bdata = matrix_long_data(mat,ndims);
 | 
						|
  if (data)
 | 
						|
    memcpy((void *)bdata,(void *)data,sizeof(double)*nelems);
 | 
						|
  return blob;
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
new_float_matrix(int ndims, int dims[], double data[])
 | 
						|
{
 | 
						|
   unsigned int sz;
 | 
						|
  unsigned int i, nelems=1;
 | 
						|
  YAP_Term blob;
 | 
						|
  int *mat;
 | 
						|
  double *bdata;
 | 
						|
  int idims[MAX_DIMS];
 | 
						|
 | 
						|
  /* in case we don't have enough room and need to shift the stack, we can't
 | 
						|
     really afford to keep a pointer to the global */
 | 
						|
  for (i=0;i< ndims;i++) {
 | 
						|
    idims[i] = dims[i];
 | 
						|
    nelems *= dims[i];
 | 
						|
  }
 | 
						|
  sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+(nelems+1)*sizeof(double)+(sizeof(YAP_CELL)-1))/sizeof(YAP_CELL);
 | 
						|
  blob = YAP_MkBlobTerm(sz);
 | 
						|
  if (blob == YAP_TermNil())
 | 
						|
    return blob;
 | 
						|
  mat = YAP_BlobOfTerm(blob);
 | 
						|
  mat[MAT_TYPE] = FLOAT_MATRIX;
 | 
						|
  mat[MAT_NDIMS] = ndims;
 | 
						|
  mat[MAT_SIZE] = nelems;
 | 
						|
  for (i=0;i< ndims;i++) {
 | 
						|
    mat[MAT_DIMS+i] = idims[i];
 | 
						|
  }
 | 
						|
  bdata = matrix_double_data(mat,ndims);
 | 
						|
  if (data)
 | 
						|
    memcpy((void *)bdata,(void *)data,sizeof(double)*nelems);
 | 
						|
  return blob;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i = 0; i < ndims; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
    int d;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tl)) {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    th = YAP_HeadOfTerm(tl);
 | 
						|
    if (!YAP_IsIntTerm(th)) {
 | 
						|
      /* ERROR */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    d = YAP_IntOfTerm(th);
 | 
						|
    if (d < 0) {
 | 
						|
      /* ERROR */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    dims[i] = d;
 | 
						|
    tl = YAP_TailOfTerm(tl);
 | 
						|
  }
 | 
						|
  if (tl != YAP_TermNil()) {
 | 
						|
    /* ERROR */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
cp_int_matrix(YAP_Term tl,YAP_Term matrix)
 | 
						|
{
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(matrix);
 | 
						|
  int i, nelems = mat[MAT_SIZE];
 | 
						|
  long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  for (i = 0; i < nelems; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
    int d;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tl)) {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    th = YAP_HeadOfTerm(tl);
 | 
						|
    if (!YAP_IsIntTerm(th)) {
 | 
						|
      /* ERROR */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    d = YAP_IntOfTerm(th);
 | 
						|
    j[i] = d;
 | 
						|
    tl = YAP_TailOfTerm(tl);
 | 
						|
  }
 | 
						|
  if (tl != YAP_TermNil()) {
 | 
						|
    /* ERROR */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
cp_float_matrix(YAP_Term tl,YAP_Term matrix)
 | 
						|
{
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(matrix);
 | 
						|
  int i, nelems = mat[MAT_SIZE];
 | 
						|
  double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  for (i = 0; i < nelems; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
    double d;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tl)) {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
    th = YAP_HeadOfTerm(tl);
 | 
						|
    if (YAP_IsIntTerm(th)) {
 | 
						|
      d = YAP_IntOfTerm(th);
 | 
						|
    } else if (!YAP_IsFloatTerm(th)) {
 | 
						|
      /* ERROR */
 | 
						|
      return FALSE;
 | 
						|
    } else {
 | 
						|
      d = YAP_FloatOfTerm(th);
 | 
						|
    }
 | 
						|
    j[i] = d;
 | 
						|
    tl = YAP_TailOfTerm(tl);
 | 
						|
  }
 | 
						|
  if (tl != YAP_TermNil()) {
 | 
						|
    /* ERROR */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
 | 
						|
static int
 | 
						|
set_int_matrix(YAP_Term matrix,long int set)
 | 
						|
{
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(matrix);
 | 
						|
  int i, nelems = mat[MAT_SIZE];
 | 
						|
  long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  for (i = 0; i < nelems; i++) {
 | 
						|
    j[i] = set;
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
set_float_matrix(YAP_Term matrix,double set)
 | 
						|
{
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(matrix);
 | 
						|
  int i, nelems = mat[MAT_SIZE];
 | 
						|
  double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  for (i = 0; i < nelems; i++) {
 | 
						|
    j[i] = set;
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
new_ints_matrix(void)
 | 
						|
{
 | 
						|
  int ndims = YAP_IntOfTerm(YAP_ARG1);
 | 
						|
  YAP_Term tl = YAP_ARG2, out;
 | 
						|
  int dims[MAX_DIMS];
 | 
						|
  
 | 
						|
  if (!scan_dims(ndims, tl, dims))
 | 
						|
    return FALSE;
 | 
						|
  out = new_int_matrix(ndims, dims, NULL);
 | 
						|
  if (out == YAP_TermNil())
 | 
						|
    return FALSE;
 | 
						|
  if (!cp_int_matrix(YAP_ARG3,out))
 | 
						|
    return FALSE;
 | 
						|
  return YAP_Unify(YAP_ARG4, out);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
new_ints_matrix_set(void)
 | 
						|
{
 | 
						|
  int ndims = YAP_IntOfTerm(YAP_ARG1);
 | 
						|
  YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
 | 
						|
  int dims[MAX_DIMS];
 | 
						|
  long int set;
 | 
						|
  
 | 
						|
  if (!YAP_IsIntTerm(tset)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  set = YAP_IntOfTerm(tset);
 | 
						|
  if (!scan_dims(ndims, tl, dims))
 | 
						|
    return FALSE;
 | 
						|
  out = new_int_matrix(ndims, dims, NULL);
 | 
						|
  if (!set_int_matrix(out,set))
 | 
						|
    return FALSE;
 | 
						|
  return YAP_Unify(YAP_ARG4, out);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
new_floats_matrix(void)
 | 
						|
{
 | 
						|
  int ndims = YAP_IntOfTerm(YAP_ARG1);
 | 
						|
  YAP_Term tl = YAP_ARG2, out;
 | 
						|
  int dims[MAX_DIMS];
 | 
						|
  if (!scan_dims(ndims, tl, dims))
 | 
						|
    return FALSE;
 | 
						|
  out = new_float_matrix(ndims, dims, NULL);
 | 
						|
  if (out == YAP_TermNil())
 | 
						|
    return FALSE;
 | 
						|
  if (!cp_float_matrix(YAP_ARG3,out))
 | 
						|
    return FALSE;
 | 
						|
  return YAP_Unify(YAP_ARG4, out);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
new_floats_matrix_set(void)
 | 
						|
{
 | 
						|
  int ndims = YAP_IntOfTerm(YAP_ARG1);
 | 
						|
  YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
 | 
						|
  int dims[MAX_DIMS];
 | 
						|
  double set;
 | 
						|
  
 | 
						|
  if (!YAP_IsFloatTerm(tset)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  set = YAP_FloatOfTerm(tset);
 | 
						|
  if (!scan_dims(ndims, tl, dims))
 | 
						|
    return FALSE;
 | 
						|
  out = new_float_matrix(ndims, dims, NULL);
 | 
						|
  if (!set_float_matrix(out,set))
 | 
						|
    return FALSE;
 | 
						|
  return YAP_Unify(YAP_ARG4, out);
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
float_matrix_to_list(int *mat) {
 | 
						|
  double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  /* prepare for worst case with double taking two cells */
 | 
						|
  if (YAP_RequiresExtraStack(6*mat[MAT_SIZE])) {
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
  }
 | 
						|
  return YAP_FloatsToList(data, mat[MAT_SIZE]);
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
mk_int_list(int nelems, int *data)
 | 
						|
{
 | 
						|
  YAP_Term tn = YAP_TermNil();
 | 
						|
  YAP_Term tf = tn;
 | 
						|
  int i = 0;
 | 
						|
 | 
						|
  for (i = nelems-1; i>= 0; i--) {
 | 
						|
    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
 | 
						|
    if (tf == tn) {
 | 
						|
      /* error */
 | 
						|
      return tn;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return tf;
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
mk_long_list(int nelems, long int *data)
 | 
						|
{
 | 
						|
  YAP_Term tn = YAP_TermNil();
 | 
						|
  YAP_Term tf = tn;
 | 
						|
  int i = 0;
 | 
						|
 | 
						|
  for (i = nelems-1; i>= 0; i--) {
 | 
						|
    tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
 | 
						|
    if (tf == tn) {
 | 
						|
      /* error */
 | 
						|
      return tn;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return tf;
 | 
						|
}
 | 
						|
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
long_matrix_to_list(int *mat) {
 | 
						|
  long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  /* prepare for worst case with longs evrywhere (3cells + 1) */
 | 
						|
  if (YAP_RequiresExtraStack(5*mat[MAT_SIZE])) {
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
  }
 | 
						|
  return mk_long_list(mat[MAT_SIZE], data);
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
matrix_access(int *mat, int *indx)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  if (mat[MAT_TYPE]==FLOAT_MATRIX)
 | 
						|
    return YAP_MkFloatTerm((matrix_double_data(mat,mat[MAT_NDIMS]))[off]);
 | 
						|
  else
 | 
						|
    return YAP_MkIntTerm((matrix_long_data(mat,mat[MAT_NDIMS]))[off]);
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_float_set(int *mat, int *indx, double nval)
 | 
						|
{
 | 
						|
  unsigned int off = 0;
 | 
						|
 | 
						|
  off = matrix_get_offset(mat, indx);
 | 
						|
  (matrix_double_data(mat,mat[MAT_NDIMS]))[off] = nval;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_set(int *mat, int *indx, long int nval)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  (matrix_long_data(mat,mat[MAT_NDIMS]))[off] = nval;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_float_set_all(int *mat, double nval)
 | 
						|
{
 | 
						|
  int i;
 | 
						|
  double *data = matrix_double_data(mat,mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  for (i = 0; i< mat[MAT_SIZE]; i++)
 | 
						|
    data[i] = nval;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_set_all(int *mat, long int nval)
 | 
						|
{
 | 
						|
  int i;
 | 
						|
  long int *data = matrix_long_data(mat,mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  if (nval == 0) {
 | 
						|
    memset((void *)data,0,sizeof(long int)*mat[MAT_SIZE]);
 | 
						|
  } else {
 | 
						|
    for (i = 0; i< mat[MAT_SIZE]; i++)
 | 
						|
      data[i] = nval;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_float_add(int *mat, int *indx, double nval)
 | 
						|
{
 | 
						|
  unsigned int off;
 | 
						|
  double *dat = matrix_double_data(mat,mat[MAT_NDIMS]);
 | 
						|
 | 
						|
  off = matrix_get_offset(mat, indx);
 | 
						|
  dat[off] += nval;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_add(int *mat, int *indx, long int nval)
 | 
						|
{
 | 
						|
  long int *dat = matrix_long_data(mat,mat[MAT_NDIMS]);
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  dat[off] += nval;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_inc(int *mat, int *indx)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  if (mat[MAT_TYPE]==FLOAT_MATRIX)
 | 
						|
    (matrix_double_data(mat,mat[MAT_NDIMS])[off])++;
 | 
						|
  else
 | 
						|
    ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])++;
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_dec(int *mat, int *indx)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  if (mat[MAT_TYPE]==FLOAT_MATRIX)
 | 
						|
    (matrix_double_data(mat,mat[MAT_NDIMS])[off])--;
 | 
						|
  else
 | 
						|
    ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])--;
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
matrix_inc2(int *mat, int *indx)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  if (mat[MAT_TYPE]==FLOAT_MATRIX) {
 | 
						|
    double *data  = matrix_double_data(mat,mat[MAT_NDIMS]);
 | 
						|
    double d = data[off];
 | 
						|
    d++;
 | 
						|
    data[off] = d;
 | 
						|
    return YAP_MkFloatTerm(d);
 | 
						|
  } else {
 | 
						|
    long int *data  = matrix_long_data(mat,mat[MAT_NDIMS]);
 | 
						|
    long int d = data[off];
 | 
						|
    d++;
 | 
						|
    data[off] = d;
 | 
						|
    return YAP_MkIntTerm(d);
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static YAP_Term
 | 
						|
matrix_dec2(int *mat, int *indx)
 | 
						|
{
 | 
						|
  unsigned int off = matrix_get_offset(mat, indx);
 | 
						|
  if (mat[MAT_TYPE]==FLOAT_MATRIX) {
 | 
						|
    double *data  = matrix_double_data(mat,mat[MAT_NDIMS]);
 | 
						|
    double d = data[off];
 | 
						|
    d--;
 | 
						|
    data[off] = d;
 | 
						|
    return YAP_MkFloatTerm(d);
 | 
						|
  } else {
 | 
						|
    long int *data  = matrix_long_data(mat,mat[MAT_NDIMS]);
 | 
						|
    long int d = data[off];
 | 
						|
    d--;
 | 
						|
    data[off] = d;
 | 
						|
    return YAP_MkIntTerm(d);
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
 | 
						|
static int
 | 
						|
matrix_set(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  tf = YAP_ARG3;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_long_set(mat, dims, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_long_set(mat, dims, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_float_set(mat, dims, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_float_set(mat, dims, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_set_all(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  tf = YAP_ARG2;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_long_set_all(mat, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_long_set_all(mat, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_float_set_all(mat, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_float_set_all(mat, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_add(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  tf = YAP_ARG3;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_long_add(mat, dims, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_long_add(mat, dims, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    if (YAP_IsIntTerm(tf)) {
 | 
						|
      matrix_float_add(mat, dims, YAP_IntOfTerm(tf));
 | 
						|
    } else if (YAP_IsFloatTerm(tf)) {
 | 
						|
      matrix_float_add(mat, dims, YAP_FloatOfTerm(tf));
 | 
						|
    } else {
 | 
						|
      /* Error */
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
do_matrix_access(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  tf = matrix_access(mat, dims);
 | 
						|
  return YAP_Unify(tf, YAP_ARG3);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
do_matrix_inc(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  matrix_inc(mat, dims);
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
do_matrix_dec(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  matrix_dec(mat, dims);
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
do_matrix_inc2(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return
 | 
						|
    YAP_Unify(matrix_inc2(mat, dims), YAP_ARG3);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
do_matrix_dec2(void)
 | 
						|
{
 | 
						|
  int dims[MAX_DIMS], *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return
 | 
						|
    YAP_Unify(matrix_dec2(mat, dims), YAP_ARG3);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_to_list(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX)
 | 
						|
    tf = long_matrix_to_list(mat);
 | 
						|
  else
 | 
						|
    tf = float_matrix_to_list(mat);
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_dims(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  tf = mk_int_list(mat[MAT_NDIMS],mat+MAT_DIMS);
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_size(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_SIZE]));
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_ndims(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_NDIMS]));
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_type(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    tf = YAP_MkIntTerm(0);
 | 
						|
  } else {
 | 
						|
    tf = YAP_MkIntTerm(1);
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_arg_to_offset(void)
 | 
						|
{
 | 
						|
  int indx[MAX_DIMS], *mat;
 | 
						|
  unsigned int off;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, indx)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  off = matrix_get_offset(mat, indx);
 | 
						|
 | 
						|
  return YAP_Unify(YAP_ARG3, YAP_MkIntTerm(off));
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_offset_to_arg(void)
 | 
						|
{
 | 
						|
  int indx[MAX_DIMS], *mat;
 | 
						|
  unsigned int off;
 | 
						|
  YAP_Term ti, tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (!YAP_IsIntTerm(ti = YAP_ARG2)) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  off = YAP_IntOfTerm(ti);
 | 
						|
  matrix_get_index(mat, off, indx);
 | 
						|
  tf = mk_int_list(mat[MAT_NDIMS], indx);
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
static unsigned int
 | 
						|
scan_max_long(int sz, long int *data)
 | 
						|
{
 | 
						|
  int i, off=0;
 | 
						|
  long int max= data[0];
 | 
						|
  for (i=1; i<sz; i++) {
 | 
						|
    if (data[i]>max) {
 | 
						|
      off=i;
 | 
						|
      max = data[i];
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return off;
 | 
						|
}
 | 
						|
 | 
						|
static unsigned int
 | 
						|
scan_max_float(int sz, double *data)
 | 
						|
{
 | 
						|
  int i, off=0;
 | 
						|
  double max= data[0];
 | 
						|
  for (i=1; i<sz; i++) {
 | 
						|
    if (data[i]>max) {
 | 
						|
      max = data[i];
 | 
						|
      off=i;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return off;
 | 
						|
}
 | 
						|
 | 
						|
static unsigned int
 | 
						|
scan_min_long(int sz, long int *data)
 | 
						|
{
 | 
						|
  int i, off=0;
 | 
						|
  long int max= data[0];
 | 
						|
  for (i=1; i<sz; i++) {
 | 
						|
    if (data[i]<max) {
 | 
						|
      max = data[i];
 | 
						|
      off=i;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return off;
 | 
						|
}
 | 
						|
 | 
						|
static unsigned int
 | 
						|
scan_min_float(int sz, double *data)
 | 
						|
{
 | 
						|
  int i, off=0;
 | 
						|
  double max= data[0];
 | 
						|
  for (i=1; i<sz; i++) {
 | 
						|
    if (data[i]<max) {
 | 
						|
      max = data[i];
 | 
						|
      off=i;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return off;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_max(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  unsigned int off;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_max_long(mat[MAT_SIZE], data);
 | 
						|
    tf = YAP_MkIntTerm(data[off]);
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_max_float(mat[MAT_SIZE], data);
 | 
						|
    tf = YAP_MkFloatTerm(data[off]);
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_maxarg(void)
 | 
						|
{
 | 
						|
  int indx[MAX_DIMS], *mat;
 | 
						|
  unsigned int off;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_max_long(mat[MAT_SIZE], data);
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_max_float(mat[MAT_SIZE], data);
 | 
						|
  }
 | 
						|
  matrix_get_index(mat, off, indx);
 | 
						|
  tf = mk_int_list(mat[MAT_NDIMS], indx);
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_min(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  unsigned int off;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_min_long(mat[MAT_SIZE], data);
 | 
						|
    tf = YAP_MkIntTerm(data[off]);
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_min_float(mat[MAT_SIZE], data);
 | 
						|
    tf = YAP_MkFloatTerm(data[off]);
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_log_all(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    return FALSE;
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    int i;
 | 
						|
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      data[i] = log(data[i]);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_log_all2(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    return FALSE;
 | 
						|
  } else {
 | 
						|
    YAP_Term out;
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
 | 
						|
    int i;
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    if (!YAP_IsVarTerm(YAP_ARG2)) {
 | 
						|
      out = YAP_ARG2;
 | 
						|
    } else {
 | 
						|
      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
 | 
						|
      if (out == YAP_TermNil())
 | 
						|
	return FALSE;
 | 
						|
    }
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(out);
 | 
						|
    ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      ndata[i] = log(data[i]);
 | 
						|
    }
 | 
						|
    if (YAP_IsVarTerm(YAP_ARG2)) {
 | 
						|
      return YAP_Unify(YAP_ARG2, out);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_exp_all(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    return FALSE;
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    int i;
 | 
						|
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      data[i] = exp(data[i]);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_exp2_all(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    return FALSE;
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    int i;
 | 
						|
    double max = data[0];
 | 
						|
 | 
						|
    for (i=1; i< mat[MAT_SIZE]; i++) {
 | 
						|
      if (data[i] > max) max = data[i];
 | 
						|
    }
 | 
						|
    
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      data[i] = exp(data[i]-max);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_exp_all2(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    return FALSE;
 | 
						|
  } else {
 | 
						|
    YAP_Term out;
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
 | 
						|
    int i;
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    if (!YAP_IsVarTerm(YAP_ARG2)) {
 | 
						|
      out = YAP_ARG2;
 | 
						|
    } else {
 | 
						|
      out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
 | 
						|
      if (out == YAP_TermNil())
 | 
						|
	return FALSE;
 | 
						|
    }
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(out);
 | 
						|
    ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      ndata[i] = exp(data[i]);
 | 
						|
    }
 | 
						|
    if (YAP_IsVarTerm(YAP_ARG2)) {
 | 
						|
      return YAP_Unify(YAP_ARG2, out);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return TRUE;
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_minarg(void)
 | 
						|
{
 | 
						|
  int indx[MAX_DIMS], *mat;
 | 
						|
  unsigned int off;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_min_long(mat[MAT_SIZE], data);
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    off = scan_min_float(mat[MAT_SIZE], data);
 | 
						|
  }
 | 
						|
  matrix_get_index(mat, off, indx);
 | 
						|
  tf = mk_int_list(mat[MAT_NDIMS], indx);
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_sum(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
 | 
						|
    int i;
 | 
						|
    long int sum = 0;
 | 
						|
 | 
						|
    for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
      sum += data[i];
 | 
						|
    }
 | 
						|
    tf = YAP_MkIntTerm(sum);
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
 | 
						|
    int i;
 | 
						|
    double sum = 0.0;
 | 
						|
 | 
						|
    for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
      sum += data[i];
 | 
						|
    }
 | 
						|
    tf = YAP_MkFloatTerm(sum);
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG2, tf);
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_int_lines(int total,int nlines,long int *mat0,long int *matf)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<ncols;i++) {
 | 
						|
    long int sum = 0;
 | 
						|
    int j;
 | 
						|
 | 
						|
    for (j=i;j<total;j+=ncols) {
 | 
						|
      sum += mat0[j];
 | 
						|
    }
 | 
						|
    matf[i] = sum;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_double_lines(int total,int nlines,double *mat0,double *matf)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<ncols;i++) {
 | 
						|
    double sum = 0;
 | 
						|
    int j;
 | 
						|
 | 
						|
    for (j=i;j<total;j+=ncols) {
 | 
						|
      sum += mat0[j];
 | 
						|
    }
 | 
						|
    matf[i] = sum;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_agg_lines(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  YAP_Term top = YAP_ARG2;
 | 
						|
  op_type op;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* create a new array without first dimension */
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    tf = new_int_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
    data = matrix_long_data(mat, dims);
 | 
						|
    ndata = matrix_long_data(nmat, dims-1);
 | 
						|
    if (op == MAT_PLUS) {
 | 
						|
      add_int_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
 | 
						|
    } else
 | 
						|
      return FALSE;
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    tf = new_float_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    data = matrix_double_data(mat, dims);
 | 
						|
    ndata = matrix_double_data(nmat, dims-1);
 | 
						|
    if (op == MAT_PLUS) {
 | 
						|
      add_double_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
 | 
						|
    } else
 | 
						|
      return FALSE;
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3,tf);
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_int_cols(int total,int nlines,long int *mat0,long int *matf)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i, j = 0;
 | 
						|
  for (i=0;i<nlines;i++) {
 | 
						|
    long int sum = 0;
 | 
						|
    int max = (i+1)*ncols;
 | 
						|
 | 
						|
    for (;j<max;j++) {
 | 
						|
      sum += mat0[j];
 | 
						|
    }
 | 
						|
    matf[i] = sum;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_double_cols(int total,int nlines,double *mat0,double *matf)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i, j = 0;
 | 
						|
  for (i=0;i<nlines;i++) {
 | 
						|
    double sum = 0;
 | 
						|
    int max = (i+1)*ncols;
 | 
						|
 | 
						|
    for (;j<max;j++) {
 | 
						|
      sum += mat0[j];
 | 
						|
    }
 | 
						|
    matf[i] = sum;
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_agg_cols(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf;
 | 
						|
  YAP_Term top = YAP_ARG2;
 | 
						|
  op_type op;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* create a new array without first dimension */
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    tf = new_int_matrix(1,mat+MAT_DIMS,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
    data = matrix_long_data(mat, dims);
 | 
						|
    ndata = matrix_long_data(nmat, 1);
 | 
						|
    if (op == MAT_PLUS) {
 | 
						|
      add_int_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
 | 
						|
    } else
 | 
						|
      return FALSE;
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    tf = new_float_matrix(1,mat+MAT_DIMS,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
    data = matrix_double_data(mat, dims);
 | 
						|
    ndata = matrix_double_data(nmat, 1);
 | 
						|
    if (op == MAT_PLUS) {
 | 
						|
      add_double_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
 | 
						|
    } else
 | 
						|
      return FALSE;
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3,tf);
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
div_int_by_lines(int total,int nlines,long int *mat1,long int *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = ((double)mat1[i])/mat2[i%ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
div_int_by_dlines(int total,int nlines,long int *mat1,double *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i]/mat2[i%ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
div_float_long_by_lines(int total,int nlines,double *mat1,long int *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i]/mat2[i%ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
div_float_by_lines(int total,int nlines,double *mat1,double *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int ncols = total/nlines, i;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i]/mat2[i%ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_op_to_lines(void)
 | 
						|
{
 | 
						|
  int *mat1, *mat2;
 | 
						|
  YAP_Term top = YAP_ARG3;
 | 
						|
  op_type op;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat1) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
  if (!mat2) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* create a new array without first dimension */
 | 
						|
  if (mat1[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data1;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;   
 | 
						|
    data1 = matrix_long_data(mat1, dims);
 | 
						|
 | 
						|
    if (mat2[MAT_TYPE] == INT_MATRIX) {
 | 
						|
      long int *data2 = matrix_long_data(mat2, dims-1);
 | 
						|
      if (op == MAT_DIV) {
 | 
						|
	double *ndata;
 | 
						|
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_double_data(nmat, dims);
 | 
						|
	div_int_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
 | 
						|
      double *data2 = matrix_double_data(mat2, dims-1);
 | 
						|
      if (op == MAT_DIV) {
 | 
						|
	double *ndata;
 | 
						|
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_double_data(nmat, dims);
 | 
						|
	div_int_by_dlines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data1, *ndata;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    data1 = matrix_double_data(mat1, dims);
 | 
						|
    tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL); 
 | 
						|
    nmat = YAP_BlobOfTerm(tf);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    ndata = matrix_double_data(nmat, dims);
 | 
						|
    if (mat2[MAT_TYPE] == INT_MATRIX) {
 | 
						|
      long int *data2 = matrix_long_data(mat2, dims-1);
 | 
						|
      if (op == MAT_DIV) {
 | 
						|
	div_float_long_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
 | 
						|
      double *data2 = matrix_double_data(mat2, dims-1);
 | 
						|
      if (op == MAT_DIV) {
 | 
						|
	div_float_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG4,tf);
 | 
						|
}
 | 
						|
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_add_data(long int *nmat, int siz, long int mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]+mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_add_data(double *nmat, int siz, long int mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]+mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_double_add_data(double *nmat, int siz, double mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]+mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_sub_data(long int *nmat, int siz, long int mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]-mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_sub_data(double *nmat, int siz, long int mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]-mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_rsub_data(double *nmat, int siz, double mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat2[i]-mat1[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_double_sub_data(double *nmat, int siz, double mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]-mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_mult_data(long int *nmat, int siz, long int mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]*mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_mult_data(double *nmat, int siz, long int mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]*mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_double_mult_data(double *nmat, int siz, double mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]*mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_div_data(long int *nmat, int siz, long int mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_div_data(double *nmat, int siz, long int mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_div2_data(double *nmat, int siz, double mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_double_div_data(double *nmat, int siz, double mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_zdiv_data(long int *nmat, int siz, long int mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    if (mat1[i] == 0)
 | 
						|
      nmat[i] = 0;
 | 
						|
    else
 | 
						|
      nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_zdiv_data(double *nmat, int siz, long int mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    if (mat1[i] == 0)
 | 
						|
      nmat[i] = 0;
 | 
						|
    else
 | 
						|
      nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_long_double_zdiv2_data(double *nmat, int siz, double mat1[], long int mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    if (mat1[i] == 0.0)
 | 
						|
      nmat[i] = 0;
 | 
						|
    else
 | 
						|
      nmat[i] = mat1[i]/mat2[i];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
matrix_double_zdiv_data(double *nmat, int siz, double mat1[], double mat2[])
 | 
						|
{
 | 
						|
  int i;
 | 
						|
 | 
						|
  for (i=0; i< siz; i++) {
 | 
						|
    if (mat1[i] == 0.0) {
 | 
						|
      nmat[i] = 0.0;
 | 
						|
    } else {
 | 
						|
      nmat[i] = mat1[i]/mat2[i];
 | 
						|
    }
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_op(void)
 | 
						|
{
 | 
						|
  int *mat1, *mat2;
 | 
						|
  YAP_Term top = YAP_ARG3;
 | 
						|
  op_type op;
 | 
						|
  YAP_Term tf = YAP_ARG4;
 | 
						|
  int create = TRUE;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat1) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
  if (!mat2) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (tf == YAP_ARG1 || tf == YAP_ARG2) {
 | 
						|
    create = FALSE;
 | 
						|
  }
 | 
						|
  if (mat1[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data1;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;   
 | 
						|
    data1 = matrix_long_data(mat1, dims);
 | 
						|
 | 
						|
    if (mat2[MAT_TYPE] == INT_MATRIX) {
 | 
						|
      long int *data2;
 | 
						|
      long int *ndata;
 | 
						|
 | 
						|
      if (create)
 | 
						|
	tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
      if (tf == YAP_TermNil()) {
 | 
						|
	return FALSE;
 | 
						|
      } else {
 | 
						|
	/* there may have been an overflow */
 | 
						|
	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
	data1 = matrix_long_data(mat1, dims);
 | 
						|
	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
	data2 = matrix_long_data(mat2, dims);
 | 
						|
      }
 | 
						|
      nmat = YAP_BlobOfTerm(tf);
 | 
						|
      ndata = matrix_long_data(nmat, dims);
 | 
						|
      switch (op) {
 | 
						|
      case MAT_PLUS:
 | 
						|
	matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_SUB:
 | 
						|
	matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_TIMES:
 | 
						|
	matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_DIV:
 | 
						|
	matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_ZDIV:
 | 
						|
	matrix_long_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      default:
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
 | 
						|
      double *data2;
 | 
						|
      double *ndata;
 | 
						|
 | 
						|
      if (create)
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
      if (tf == YAP_TermNil()) {
 | 
						|
	return FALSE;
 | 
						|
      } else {
 | 
						|
	/* there may have been an overflow */
 | 
						|
	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
	data1 = matrix_long_data(mat1, dims);
 | 
						|
	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
	data2  = matrix_double_data(mat2, dims);
 | 
						|
      }
 | 
						|
      nmat = YAP_BlobOfTerm(tf);
 | 
						|
      ndata = matrix_double_data(nmat, dims);
 | 
						|
      switch (op) {
 | 
						|
      case MAT_PLUS:
 | 
						|
	matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_SUB:
 | 
						|
	matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_TIMES:
 | 
						|
	matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_DIV:
 | 
						|
	matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_ZDIV:
 | 
						|
	matrix_long_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      default:
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data1;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;   
 | 
						|
    data1 = matrix_double_data(mat1, dims);
 | 
						|
 | 
						|
    if (mat2[MAT_TYPE] == INT_MATRIX) {
 | 
						|
      long int *data2;
 | 
						|
      double *ndata;
 | 
						|
 | 
						|
      if (create)
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
      if (tf == YAP_TermNil()) {
 | 
						|
	return FALSE;
 | 
						|
      } else {
 | 
						|
	/* there may have been an overflow */
 | 
						|
	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
	data1 = matrix_double_data(mat1, dims);
 | 
						|
	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
	data2  = matrix_long_data(mat2, dims);
 | 
						|
      }
 | 
						|
      nmat = YAP_BlobOfTerm(tf);
 | 
						|
      ndata = matrix_double_data(nmat, dims);
 | 
						|
      switch (op) {
 | 
						|
      case MAT_PLUS:
 | 
						|
	matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
 | 
						|
	break;
 | 
						|
      case MAT_SUB:
 | 
						|
	matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_TIMES:
 | 
						|
	matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
 | 
						|
	break;
 | 
						|
      case MAT_DIV:
 | 
						|
	matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_ZDIV:
 | 
						|
	matrix_long_double_zdiv2_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      default:
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
 | 
						|
      double *data2;
 | 
						|
      double *ndata;
 | 
						|
 | 
						|
      if (create)
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
      if (tf == YAP_TermNil()) {
 | 
						|
	return FALSE;
 | 
						|
      } else {
 | 
						|
	/* there may have been an overflow */
 | 
						|
	mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
	data1 = matrix_double_data(mat1, dims);
 | 
						|
	mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
	data2  = matrix_double_data(mat2, dims);
 | 
						|
      }
 | 
						|
      nmat = YAP_BlobOfTerm(tf);
 | 
						|
      ndata = matrix_double_data(nmat, dims);
 | 
						|
      switch (op) {
 | 
						|
      case MAT_PLUS:
 | 
						|
	matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_SUB:
 | 
						|
	matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_TIMES:
 | 
						|
	matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_DIV:
 | 
						|
	matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      case MAT_ZDIV:
 | 
						|
	matrix_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
 | 
						|
	break;
 | 
						|
      default:
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG4,tf);
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_int_by_cols(int total,int nlines,long int *mat1,long int *mat2,long int *ndata)
 | 
						|
{
 | 
						|
  int i, ncols = total/nlines;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i] + mat2[i/ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_int_by_dcols(int total,int nlines,long int *mat1,double *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int i, ncols = total/nlines;
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i] + mat2[i/ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static void
 | 
						|
add_double_by_cols(int total,int nlines,double *mat1,double *mat2,double *ndata)
 | 
						|
{
 | 
						|
  int i;
 | 
						|
  int ncols = total/nlines;
 | 
						|
 | 
						|
  for (i=0;i<total;i++) {
 | 
						|
    ndata[i] = mat1[i] + mat2[i/ncols];
 | 
						|
  }
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_op_to_cols(void)
 | 
						|
{
 | 
						|
  int *mat1, *mat2;
 | 
						|
  YAP_Term top = YAP_ARG3;
 | 
						|
  op_type op;
 | 
						|
  YAP_Term tf;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat1) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
 | 
						|
  if (!mat2) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (mat1[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data1;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;   
 | 
						|
    data1 = matrix_long_data(mat1, dims);
 | 
						|
 | 
						|
    if (mat2[MAT_TYPE] == INT_MATRIX) {
 | 
						|
      long int *data2 = matrix_long_data(mat2, 1);
 | 
						|
      if (op == MAT_PLUS) {
 | 
						|
	long int *ndata;
 | 
						|
 | 
						|
	tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_long_data(nmat, dims);
 | 
						|
	add_int_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
 | 
						|
      double *data2 = matrix_double_data(mat2, 1);
 | 
						|
      if (op == MAT_PLUS) {
 | 
						|
	double *ndata;
 | 
						|
 | 
						|
	tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_double_data(nmat, dims);
 | 
						|
	add_int_by_dcols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data1, *data2, *ndata;
 | 
						|
    int dims = mat1[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
 | 
						|
    if (mat2[MAT_TYPE] != FLOAT_MATRIX)
 | 
						|
      return FALSE;
 | 
						|
    tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    nmat = YAP_BlobOfTerm(tf);
 | 
						|
    data1 = matrix_double_data(mat1, dims);
 | 
						|
    data2 = matrix_double_data(mat2, 1);
 | 
						|
    ndata = matrix_double_data(nmat, dims);
 | 
						|
    if (op == MAT_PLUS) {
 | 
						|
      add_double_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
 | 
						|
    } else
 | 
						|
      return FALSE;
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG4,tf);
 | 
						|
}
 | 
						|
 | 
						|
static int
 | 
						|
matrix_op_to_all(void)
 | 
						|
{
 | 
						|
  int *mat;
 | 
						|
  YAP_Term tf = 0;
 | 
						|
  YAP_Term top = YAP_ARG2;
 | 
						|
  op_type op;
 | 
						|
  int create = FALSE;
 | 
						|
 | 
						|
  if (!YAP_IsIntTerm(top)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  op = YAP_IntOfTerm(top);
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  if (YAP_IsVarTerm(YAP_ARG4)) {
 | 
						|
    create = TRUE;
 | 
						|
  } 
 | 
						|
  /* create a new array with same dimensions */
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
    YAP_Term tnum = YAP_ARG3;
 | 
						|
 | 
						|
    if (YAP_IsIntTerm(tnum)) {
 | 
						|
      long int num;
 | 
						|
      long int *ndata;
 | 
						|
      
 | 
						|
      num = YAP_IntOfTerm(tnum);
 | 
						|
      data = matrix_long_data(mat, dims);
 | 
						|
      if (create) {
 | 
						|
	tf = new_int_matrix(dims,mat+(MAT_DIMS),NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_long_data(nmat, dims);
 | 
						|
      } else {
 | 
						|
	nmat = mat;
 | 
						|
	ndata = data; 
 | 
						|
      }
 | 
						|
      if (op == MAT_PLUS) {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] + num;
 | 
						|
	} 
 | 
						|
      } else if (op == MAT_TIMES) {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] * num;
 | 
						|
	}
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
    } else if (YAP_IsFloatTerm(tnum)) {
 | 
						|
      double num;
 | 
						|
      double *ndata;
 | 
						|
      
 | 
						|
      
 | 
						|
      num = YAP_FloatOfTerm(tnum);
 | 
						|
      if (create) {
 | 
						|
	tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
 | 
						|
	if (tf == YAP_TermNil())
 | 
						|
	  return FALSE;
 | 
						|
	nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
	ndata = matrix_double_data(nmat, dims);
 | 
						|
      } else {
 | 
						|
	return FALSE;
 | 
						|
      }
 | 
						|
      data = matrix_long_data(mat, dims);
 | 
						|
      if (op == MAT_PLUS) {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] + num;
 | 
						|
	}
 | 
						|
      } else if (op == MAT_TIMES) {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] * num;
 | 
						|
	}
 | 
						|
      } else if (op == MAT_DIV) {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] / num;
 | 
						|
	}
 | 
						|
      }
 | 
						|
    } else {
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
    int dims = mat[MAT_NDIMS];
 | 
						|
    int *nmat;
 | 
						|
    YAP_Term tnum = YAP_ARG3;
 | 
						|
    double num;
 | 
						|
 | 
						|
    if (YAP_IsFloatTerm(tnum)) {
 | 
						|
      num = YAP_FloatOfTerm(tnum);
 | 
						|
    } else if (!YAP_IntOfTerm(tnum)) {
 | 
						|
      return FALSE;
 | 
						|
    } else {
 | 
						|
      if (!create)
 | 
						|
	return FALSE;
 | 
						|
      num = (double)YAP_IntOfTerm(tnum);
 | 
						|
    }
 | 
						|
    data = matrix_double_data(mat, dims);
 | 
						|
    if (create) {
 | 
						|
      tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
 | 
						|
      if (tf == YAP_TermNil())
 | 
						|
	return FALSE;
 | 
						|
      nmat = (int *)YAP_BlobOfTerm(tf);
 | 
						|
      ndata = matrix_double_data(nmat, dims);
 | 
						|
    } else {
 | 
						|
      nmat = mat;
 | 
						|
      ndata = data;
 | 
						|
    }
 | 
						|
    switch(op) {
 | 
						|
    case MAT_PLUS:
 | 
						|
      {
 | 
						|
	int i;
 | 
						|
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] + num;
 | 
						|
	}
 | 
						|
      }
 | 
						|
      break;
 | 
						|
    case MAT_TIMES:
 | 
						|
      {
 | 
						|
	int i;
 | 
						|
	
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] * num;
 | 
						|
	}
 | 
						|
      }
 | 
						|
      break;
 | 
						|
    case MAT_DIV:
 | 
						|
      {
 | 
						|
	int i;
 | 
						|
      
 | 
						|
	for (i = 0; i < mat[MAT_SIZE]; i++) {
 | 
						|
	  ndata[i] = data[i] / num;
 | 
						|
	}
 | 
						|
      } 
 | 
						|
      break;
 | 
						|
    default:
 | 
						|
      return FALSE;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  if (create)
 | 
						|
    return YAP_Unify(YAP_ARG4,tf);
 | 
						|
  return YAP_Unify(YAP_ARG4,YAP_ARG1);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, build a new reordered matrix to follow
 | 
						|
   the new order
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_transpose(void)
 | 
						|
{
 | 
						|
  int ndims, i, *dims, *dimsn;
 | 
						|
  int conv[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
 | 
						|
  YAP_Term tconv, tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(ndims,mat+MAT_DIMS,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
  } else {
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(ndims,mat+MAT_DIMS,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
  }
 | 
						|
  /* just in case there was an overflow */
 | 
						|
  mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
  nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  dimsn = nmat+MAT_DIMS;
 | 
						|
  /* we now have our target matrix, let us grab our conversion matrix */
 | 
						|
  tconv = YAP_ARG2;
 | 
						|
  for (i=0; i < ndims; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
    long int j;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tconv))
 | 
						|
      return FALSE;
 | 
						|
    th = YAP_HeadOfTerm(tconv);
 | 
						|
    if (!YAP_IsIntTerm(th))
 | 
						|
      return FALSE;
 | 
						|
    conv[i] = j = YAP_IntOfTerm(th);
 | 
						|
    dimsn[i] = dims[j];
 | 
						|
    tconv = YAP_TailOfTerm(tconv);
 | 
						|
  }
 | 
						|
  /*
 | 
						|
    we now got all the dimensions set up, so what we need to do
 | 
						|
    next is to copy the elements to the new matrix.
 | 
						|
  */
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data = matrix_long_data(mat,ndims);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      long int x = data[i];
 | 
						|
      int j;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0; j < ndims; j++)  {
 | 
						|
	nindx[j] = indx[conv[j]];
 | 
						|
      }
 | 
						|
      matrix_long_set(nmat, nindx, x);
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data = matrix_double_data(mat,ndims);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      double x = data[i];
 | 
						|
      long j;
 | 
						|
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0; j < ndims; j++) 
 | 
						|
	nindx[j] = indx[conv[j]];
 | 
						|
      matrix_float_set(nmat, nindx, x);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, fold one of the dimensions of the
 | 
						|
   matrix on one of the elements
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_select(void)
 | 
						|
{
 | 
						|
  int ndims, i, j, newdims, prdim, leftarg, *dims, indx[MAX_DIMS];
 | 
						|
  int nindx[MAX_DIMS];
 | 
						|
  YAP_Term tpdim, tdimarg, tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tpdim = YAP_ARG2;
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  if (!YAP_IsIntTerm(tpdim)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  prdim = YAP_IntOfTerm(tpdim);
 | 
						|
  tdimarg = YAP_ARG3;
 | 
						|
  if (!YAP_IsIntTerm(tdimarg)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  leftarg = YAP_IntOfTerm(tdimarg);
 | 
						|
  for (i=0, j=0; i< ndims; i++) {
 | 
						|
    if (i != prdim) {
 | 
						|
      nindx[j]= dims[i];
 | 
						|
      j++;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  newdims = ndims-1;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      int j,k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(nmat, i, indx);
 | 
						|
      for (j = 0, k=0; j < newdims; j++,k++)  {
 | 
						|
	if (j == prdim) {
 | 
						|
	  nindx[k] = leftarg;
 | 
						|
	  k++;
 | 
						|
	}
 | 
						|
	nindx[k]= indx[j];
 | 
						|
      }
 | 
						|
      if (k == prdim) {
 | 
						|
	  nindx[k] = leftarg;
 | 
						|
      }
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, nindx)];
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      int j,k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(nmat, i, indx);
 | 
						|
      for (j = 0, k=0; j < newdims; j++,k++)  {
 | 
						|
	if (j == prdim) {
 | 
						|
	  nindx[k] = leftarg;
 | 
						|
	  k++;
 | 
						|
	}
 | 
						|
	nindx[k]= indx[j];
 | 
						|
      }
 | 
						|
      if (k == prdim) {
 | 
						|
	  nindx[k] = leftarg;
 | 
						|
      }
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, nindx)];
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG4, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of N-1 dims, get the first dimension
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_column(void)
 | 
						|
{
 | 
						|
  int size, i, ndims, newdims[1];
 | 
						|
  int indx[MAX_DIMS];
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  YAP_Term tconv, tf;
 | 
						|
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tconv = YAP_ARG2;
 | 
						|
  for (i=1; i < ndims; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tconv))
 | 
						|
      return FALSE;
 | 
						|
    th = YAP_HeadOfTerm(tconv);
 | 
						|
    if (!YAP_IsIntTerm(th))
 | 
						|
      return FALSE;
 | 
						|
    indx[i] = YAP_IntOfTerm(th);
 | 
						|
    tconv = YAP_TailOfTerm(tconv);
 | 
						|
  }
 | 
						|
  if (tconv != YAP_TermNil())
 | 
						|
    return FALSE;
 | 
						|
  newdims[0] = size = mat[MAT_DIMS];
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(1, newdims, NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,1);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0; i< size; i++) {
 | 
						|
      indx[0]=i;
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, indx)];
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(1,newdims,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,1);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0; i< size; i++) {
 | 
						|
      indx[0]=i;
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, indx)];
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, sum out one of the dimensions
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_sum_out(void)
 | 
						|
{
 | 
						|
  int ndims, i, j, newdims, prdim;
 | 
						|
  int indx[MAX_DIMS], nindx[MAX_DIMS];
 | 
						|
  YAP_Term tpdim, tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tpdim = YAP_ARG2;
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  if (!YAP_IsIntTerm(tpdim)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  prdim = YAP_IntOfTerm(tpdim);
 | 
						|
  newdims = ndims-1;
 | 
						|
  for (i=0, j=0; i< ndims; i++) {
 | 
						|
    if (i != prdim) {
 | 
						|
      nindx[j]= (mat+MAT_DIMS)[i];
 | 
						|
      j++;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (j != prdim) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] += data[i];
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0.0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (j != prdim) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] += data[i];
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, sum out one of the dimensions
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_sum_out_several(void)
 | 
						|
{
 | 
						|
  int ndims, i, *dims, newdims;
 | 
						|
  int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
 | 
						|
  YAP_Term tf, tconv;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tconv = YAP_ARG2;
 | 
						|
  for (i=0, newdims=0; i < ndims; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tconv))
 | 
						|
      return FALSE;
 | 
						|
    th = YAP_HeadOfTerm(tconv);
 | 
						|
    if (!YAP_IsIntTerm(th))
 | 
						|
      return FALSE;
 | 
						|
    conv[i] = YAP_IntOfTerm(th);
 | 
						|
    if (!conv[i]) {
 | 
						|
      nindx[newdims++] = dims[i];
 | 
						|
    }
 | 
						|
    tconv = YAP_TailOfTerm(tconv);
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (!conv[j]) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0.0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (!conv[j]) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, sum out one of the dimensions
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_sum_out_logs(void)
 | 
						|
{
 | 
						|
  int ndims, i, j, *dims, newdims, prdim;
 | 
						|
  int nindx[MAX_DIMS];
 | 
						|
  YAP_Term tpdim, tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tpdim = YAP_ARG2;
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  if (!YAP_IsIntTerm(tpdim)) {
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  prdim = YAP_IntOfTerm(tpdim);
 | 
						|
  newdims = ndims-1;
 | 
						|
  for (i=0, j=0; i< ndims; i++) {
 | 
						|
    if (i != prdim) {
 | 
						|
      nindx[j]= dims[i];
 | 
						|
      j++;
 | 
						|
    }
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
    int d = 1, j = 0, dd = 1;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    while (j < prdim) {
 | 
						|
      d = d*dims[j];
 | 
						|
      j++;
 | 
						|
    }
 | 
						|
    dd = d*dims[prdim];
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) {
 | 
						|
      int j = i % d + (i/dd)*d;
 | 
						|
      ndata[j] = exp(data[i]);
 | 
						|
    }
 | 
						|
    for (; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j = i % d + (i/dd)*d;
 | 
						|
      ndata[j] += exp(data[i]);
 | 
						|
    }
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      ndata[i] = log(ndata[i]);
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
    int d = 1, j = 0, dd = 1;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    
 | 
						|
    j = ndims-1;
 | 
						|
    while (j > prdim) {
 | 
						|
      d = d*dims[j];
 | 
						|
      j--;
 | 
						|
    }
 | 
						|
    dd = d*dims[prdim];
 | 
						|
    bzero(ndata, sizeof(double)*nmat[MAT_SIZE]);
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      YAP_Int k = i % d + (i/dd)*d;
 | 
						|
      ndata[k] += exp(data[i]);
 | 
						|
    }
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      ndata[i] = log(ndata[i]);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, sum out one of the dimensions
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_sum_out_logs_several(void)
 | 
						|
{
 | 
						|
  int ndims, i, *dims, newdims;
 | 
						|
  int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
 | 
						|
  YAP_Term tf, tconv;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  /* we now have our target matrix, let us grab our conversion arguments */
 | 
						|
  tconv = YAP_ARG2;
 | 
						|
  for (i=0, newdims=0; i < ndims; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tconv))
 | 
						|
      return FALSE;
 | 
						|
    th = YAP_HeadOfTerm(tconv);
 | 
						|
    if (!YAP_IsIntTerm(th))
 | 
						|
      return FALSE;
 | 
						|
    conv[i] = YAP_IntOfTerm(th);
 | 
						|
    if (!conv[i]) {
 | 
						|
      nindx[newdims++] = dims[i];
 | 
						|
    }
 | 
						|
    tconv = YAP_TailOfTerm(tconv);
 | 
						|
  }
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (!conv[j]) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    /* create a new matrix with smaller size */
 | 
						|
    for (i=0;i<nmat[MAT_SIZE];i++) 
 | 
						|
      ndata[i] = 0.0;
 | 
						|
    for (i=0; i< mat[MAT_SIZE]; i++) {
 | 
						|
      int j, k;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      for (j = 0, k=0; j < ndims; j++)  {
 | 
						|
	if (!conv[j]) {
 | 
						|
	  nindx[k++]= indx[j];
 | 
						|
	}
 | 
						|
      }
 | 
						|
      ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
 | 
						|
    }
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      ndata[i] = log(ndata[i]);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, build a matrix to follow
 | 
						|
   the new order
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_expand(void)
 | 
						|
{
 | 
						|
  int ndims, i, *dims, newdims=0, olddims = 0;
 | 
						|
  int new[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
 | 
						|
  YAP_Term tconv, tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  /* we now have our target matrix, let us grab our conversion matrix */
 | 
						|
  tconv = YAP_ARG2;
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  for (i=0; i < MAX_DIMS; i++) {
 | 
						|
    YAP_Term th;
 | 
						|
    long int j;
 | 
						|
 | 
						|
    if (!YAP_IsPairTerm(tconv)) {
 | 
						|
      if (tconv != YAP_TermNil())
 | 
						|
	return FALSE;
 | 
						|
      break;
 | 
						|
    }
 | 
						|
    th = YAP_HeadOfTerm(tconv);
 | 
						|
    if (!YAP_IsIntTerm(th))
 | 
						|
      return FALSE;
 | 
						|
    newdims++;
 | 
						|
    j = YAP_IntOfTerm(th);
 | 
						|
    if (j==0) {
 | 
						|
      new[i] = 0;
 | 
						|
      nindx[i] = dims[olddims];
 | 
						|
      olddims++;
 | 
						|
    } else {
 | 
						|
      new[i] = 1;
 | 
						|
      nindx[i] = j;
 | 
						|
    }
 | 
						|
    tconv = YAP_TailOfTerm(tconv);
 | 
						|
  }
 | 
						|
  if (olddims != ndims)
 | 
						|
    return FALSE;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,newdims);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      int j,k=0;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(nmat, i, indx);
 | 
						|
      for (j = 0; j < newdims; j++)  {
 | 
						|
	if (!new[j])
 | 
						|
	  nindx[k++] = indx[j];
 | 
						|
      }
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, nindx)];
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(newdims,nindx,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,newdims);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i < newdims; i++)
 | 
						|
      indx[i] = 0;
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      int j,k=0;
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      for (j = 0; j < newdims; j++)  {
 | 
						|
	if (!new[j])
 | 
						|
	  nindx[k++] = indx[j];
 | 
						|
      }
 | 
						|
      ndata[i] = data[matrix_get_offset(mat, nindx)];
 | 
						|
      matrix_next_index(nmat+MAT_DIMS, newdims, indx);
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG3, tf);
 | 
						|
}
 | 
						|
 | 
						|
/* given a matrix M and a set of dims, build contract a matrix to follow
 | 
						|
   the new order
 | 
						|
*/ 
 | 
						|
static int
 | 
						|
matrix_set_all_that_disagree(void)
 | 
						|
{
 | 
						|
  int ndims, i, *dims;
 | 
						|
  int indx[MAX_DIMS];
 | 
						|
  YAP_Term tf;
 | 
						|
  int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
 | 
						|
  int dim = YAP_IntOfTerm(YAP_ARG2);
 | 
						|
  int pos = YAP_IntOfTerm(YAP_ARG3);
 | 
						|
 | 
						|
  if (!mat) {
 | 
						|
    /* Error */
 | 
						|
    return FALSE;
 | 
						|
  }
 | 
						|
  ndims = mat[MAT_NDIMS];
 | 
						|
  dims = mat+MAT_DIMS;
 | 
						|
  if (mat[MAT_TYPE] == INT_MATRIX) {
 | 
						|
    long int *data, *ndata, val;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_int_matrix(ndims,dims,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_long_data(mat,ndims);
 | 
						|
    ndata = matrix_long_data(nmat,ndims);
 | 
						|
    if (!YAP_IsIntTerm(YAP_ARG4))
 | 
						|
      return FALSE;
 | 
						|
    val = YAP_IntOfTerm(YAP_ARG4);
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      if (indx[dim] != pos)
 | 
						|
	ndata[i] = val;
 | 
						|
      else
 | 
						|
	ndata[i] = data[i];
 | 
						|
    }
 | 
						|
  } else {
 | 
						|
    double *data, *ndata, val;
 | 
						|
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    tf = new_float_matrix(ndims,dims,NULL);
 | 
						|
    if (tf == YAP_TermNil())
 | 
						|
      return FALSE;
 | 
						|
    /* in case the matrix moved */
 | 
						|
    mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
 | 
						|
    nmat = (int *)YAP_BlobOfTerm(tf);  
 | 
						|
    data = matrix_double_data(mat,ndims);
 | 
						|
    ndata = matrix_double_data(nmat,ndims);
 | 
						|
    if (YAP_IsFloatTerm(YAP_ARG4))
 | 
						|
      val = YAP_FloatOfTerm(YAP_ARG4);
 | 
						|
    else if (YAP_IsIntTerm(YAP_ARG4))
 | 
						|
      val = YAP_IntOfTerm(YAP_ARG4);
 | 
						|
    else
 | 
						|
      return FALSE;
 | 
						|
    /* create a new matrix with the same size */
 | 
						|
    for (i=0; i< nmat[MAT_SIZE]; i++) {
 | 
						|
      
 | 
						|
      /*
 | 
						|
	not very efficient, we could try to take advantage of the fact
 | 
						|
	that we usually only change an index at a time
 | 
						|
      */
 | 
						|
      matrix_get_index(mat, i, indx);
 | 
						|
      if (indx[dim] != pos)
 | 
						|
	ndata[i] = val;
 | 
						|
      else
 | 
						|
	ndata[i] = data[i];
 | 
						|
    }
 | 
						|
  }
 | 
						|
  return YAP_Unify(YAP_ARG5, tf);
 | 
						|
}
 | 
						|
 | 
						|
void PROTO(init_matrix, (void));
 | 
						|
 | 
						|
void
 | 
						|
init_matrix(void)
 | 
						|
{
 | 
						|
  YAP_UserCPredicate("new_ints_matrix", new_ints_matrix, 4);
 | 
						|
  YAP_UserCPredicate("new_ints_matrix_set", new_ints_matrix_set, 4);
 | 
						|
  YAP_UserCPredicate("new_floats_matrix", new_floats_matrix, 4);
 | 
						|
  YAP_UserCPredicate("new_floats_matrix_set", new_floats_matrix_set, 4);
 | 
						|
  YAP_UserCPredicate("matrix_set", matrix_set, 3);
 | 
						|
  YAP_UserCPredicate("matrix_set_all", matrix_set_all, 2);
 | 
						|
  YAP_UserCPredicate("matrix_add", matrix_add, 3);
 | 
						|
  YAP_UserCPredicate("matrix_get", do_matrix_access, 3);
 | 
						|
  YAP_UserCPredicate("matrix_inc", do_matrix_inc, 2);
 | 
						|
  YAP_UserCPredicate("matrix_dec", do_matrix_dec, 2);
 | 
						|
  YAP_UserCPredicate("matrix_inc", do_matrix_inc2, 3);
 | 
						|
  YAP_UserCPredicate("matrix_dec", do_matrix_dec2, 3);
 | 
						|
  YAP_UserCPredicate("matrix_to_list", matrix_to_list, 2);
 | 
						|
  YAP_UserCPredicate("matrix_dims", matrix_dims, 2);
 | 
						|
  YAP_UserCPredicate("matrix_ndims", matrix_ndims, 2);
 | 
						|
  YAP_UserCPredicate("matrix_size", matrix_size, 2);
 | 
						|
  YAP_UserCPredicate("matrix_type_as_number", matrix_type, 2);
 | 
						|
  YAP_UserCPredicate("matrix_arg_to_offset", matrix_arg_to_offset, 3);
 | 
						|
  YAP_UserCPredicate("matrix_offset_to_arg", matrix_offset_to_arg, 3);
 | 
						|
  YAP_UserCPredicate("matrix_max", matrix_max, 2);
 | 
						|
  YAP_UserCPredicate("matrix_maxarg", matrix_maxarg, 2);
 | 
						|
  YAP_UserCPredicate("matrix_min", matrix_min, 2);
 | 
						|
  YAP_UserCPredicate("matrix_minarg", matrix_minarg, 2);
 | 
						|
  YAP_UserCPredicate("matrix_sum", matrix_sum, 2);
 | 
						|
  YAP_UserCPredicate("matrix_shuffle", matrix_transpose, 3);
 | 
						|
  YAP_UserCPredicate("matrix_expand", matrix_expand, 3);
 | 
						|
  YAP_UserCPredicate("matrix_select", matrix_select, 4);
 | 
						|
  YAP_UserCPredicate("matrix_column", matrix_column, 3);
 | 
						|
  YAP_UserCPredicate("matrix_to_logs", matrix_log_all,1);
 | 
						|
  YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1);
 | 
						|
  YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1);
 | 
						|
  YAP_UserCPredicate("matrix_to_logs", matrix_log_all2,2);
 | 
						|
  YAP_UserCPredicate("matrix_to_exps", matrix_exp_all2, 2);
 | 
						|
  YAP_UserCPredicate("matrix_sum_out", matrix_sum_out, 3);
 | 
						|
  YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3);
 | 
						|
  YAP_UserCPredicate("matrix_sum_logs_out", matrix_sum_out_logs, 3);
 | 
						|
  YAP_UserCPredicate("matrix_sum_logs_out_several", matrix_sum_out_logs_several, 3);
 | 
						|
  YAP_UserCPredicate("matrix_set_all_that_disagree", matrix_set_all_that_disagree, 5);
 | 
						|
  YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
 | 
						|
  YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
 | 
						|
  YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);
 | 
						|
  YAP_UserCPredicate("do_matrix_op_to_all", matrix_op_to_all, 4);
 | 
						|
  YAP_UserCPredicate("do_matrix_op_to_lines", matrix_op_to_lines, 4);
 | 
						|
  YAP_UserCPredicate("do_matrix_op_to_cols", matrix_op_to_cols, 4);
 | 
						|
}
 | 
						|
 | 
						|
#ifdef _WIN32
 | 
						|
 | 
						|
int WINAPI PROTO(win_matrixs, (HANDLE, DWORD, LPVOID));
 | 
						|
 | 
						|
int WINAPI win_matrixs(HANDLE hinst, DWORD reason, LPVOID reserved)
 | 
						|
{
 | 
						|
  switch (reason) 
 | 
						|
    {
 | 
						|
    case DLL_PROCESS_ATTACH:
 | 
						|
      break;
 | 
						|
    case DLL_PROCESS_DETACH:
 | 
						|
      break;
 | 
						|
    case DLL_THREAD_ATTACH:
 | 
						|
      break;
 | 
						|
    case DLL_THREAD_DETACH:
 | 
						|
      break;
 | 
						|
    }
 | 
						|
  return 1;
 | 
						|
}
 | 
						|
#endif
 |