/************************************************************************* * * * 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} BASE = integer #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_BASE=1, MAT_NDIMS=2, MAT_SIZE=3, MAT_ALIGN=4, MAT_DIMS=5, } 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; YAP_Functor FunctorM; YAP_Atom AtomC; 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++) { int v; pos /= mat[MAT_DIMS+i]; v = indx[i]-mat[MAT_BASE]; if (v >= mat[MAT_DIMS+i]) { return off; } off += pos*v; } 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_BASE] = 0; 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_BASE] = 0; 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 YAP_Bool 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 YAP_Bool scan_dims_args(int ndims, YAP_Term tl, int dims[MAX_DIMS]) { int i; for (i = 0; i < ndims; i++) { YAP_Term th; int d; th = YAP_ArgOfTerm(2+i, tl); if (!YAP_IsIntTerm(th)) { /* ERROR */ return FALSE; } d = YAP_IntOfTerm(th); if (d < 0) { /* ERROR */ return FALSE; } dims[i] = d; } return TRUE; } static YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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_int_list2(int nelems, int base, 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]+base),tf); if (tf == tn) { /* error */ return tn; } } return tf; } static YAP_Term mk_rep_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),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 YAP_Bool 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 YAP_Bool matrix_set2(void) { int dims[MAX_DIMS], *mat; YAP_Term tf, t = YAP_ARG1; mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1,t)); if (!mat) { /* Error */ return FALSE; } if (!scan_dims_args(mat[MAT_NDIMS], t, dims)) { /* Error */ return FALSE; } tf = YAP_ARG2; 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool do_matrix_access2(void) { int dims[MAX_DIMS], *mat; YAP_Term tf, t = YAP_ARG1; mat = (int *)YAP_BlobOfTerm(YAP_ArgOfTerm(1, t)); if (!mat) { /* Error */ return FALSE; } if (!scan_dims_args(mat[MAT_NDIMS], t, dims)) { /* Error */ return FALSE; } tf = matrix_access(mat, dims); return YAP_Unify(tf, YAP_ARG2); } static YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool matrix_set_base(void) { int *mat; mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { /* Error */ return FALSE; } mat[MAT_BASE] = YAP_IntOfTerm(YAP_ARG2); return TRUE; } static YAP_Bool 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 YAP_Bool matrix_dims3(void) { int *mat; YAP_Term tf, tof; mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { /* Error */ return FALSE; } tf = mk_int_list(mat[MAT_NDIMS],mat+MAT_DIMS); tof = mk_rep_int_list(mat[MAT_NDIMS],mat[MAT_BASE]); return YAP_Unify(YAP_ARG2, tf) && YAP_Unify(YAP_ARG3, tof); } static YAP_Bool 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 YAP_Bool 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 YAP_Bool matrix_type(void) { int *mat; YAP_Term tf; mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { /* not an error, it may be called on a term matrix */ return FALSE; } if (mat[MAT_TYPE] == INT_MATRIX) { tf = YAP_MkIntTerm(0); } else { tf = YAP_MkIntTerm(1); } return YAP_Unify(YAP_ARG2, tf); } static YAP_Bool 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 YAP_Bool 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_list2(mat[MAT_NDIMS], mat[MAT_BASE], 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool matrix_log_all2(void) { int *mat; mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { /* Error */ return FALSE; } if (mat[MAT_TYPE] == INT_MATRIX) { YAP_Term out; long int *data = matrix_long_data(mat, mat[MAT_NDIMS]); double *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((double)data[i]); } if (YAP_IsVarTerm(YAP_ARG2)) { return YAP_Unify(YAP_ARG2, out); } } 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool matrix_exp_all2(void) { int *mat; mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { /* Error */ return FALSE; } if (mat[MAT_TYPE] == INT_MATRIX) { YAP_Term out; long int *data = matrix_long_data(mat, mat[MAT_NDIMS]); double *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((double)data[i]); } if (YAP_IsVarTerm(YAP_ARG2)) { return YAP_Unify(YAP_ARG2, out); } } 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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_SUB) { int i; for (i = 0; i < mat[MAT_SIZE]; i++) { ndata[i] = num - data[i]; } } 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_SUB: { int i; for (i = 0; i < mat[MAT_SIZE]; i++) { ndata[i] = num-data[i]; } } 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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]; memset(ndata, 0, 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 YAP_Bool 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 YAP_Bool 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 YAP_Bool 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); } static YAP_Bool matrix_m(void) { int ndims, i, size; YAP_Term tm, *tp; int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1); if (!mat) { return YAP_Unify(YAP_ARG1, YAP_ARG2); } ndims = mat[MAT_NDIMS]; size = mat[MAT_SIZE]; tm = YAP_MkNewApplTerm(FunctorM, 5); tp = YAP_ArgsOfTerm(tm); tp[0] = mk_int_list(ndims,mat+MAT_DIMS); tp[1] = YAP_MkIntTerm(ndims); tp[2] = YAP_MkIntTerm(size); tp[3] = mk_rep_int_list(ndims,mat[MAT_BASE]); tp[4] = YAP_MkNewApplTerm(YAP_MkFunctor(AtomC, size), size); tp = YAP_ArgsOfTerm(tp[3]); if (mat[MAT_TYPE] == INT_MATRIX) { long int *data; /* in case the matrix moved */ mat = (int *)YAP_BlobOfTerm(YAP_ARG1); data = matrix_long_data(mat,ndims); for (i=0; i< mat[MAT_SIZE]; i++) { tp[i] = YAP_MkIntTerm(data[i]); } } else { double *data; /* in case the matrix moved */ mat = (int *)YAP_BlobOfTerm(YAP_ARG1); data = matrix_double_data(mat,ndims); for (i=0; i< mat[MAT_SIZE]; i++) { tp[i] = YAP_MkFloatTerm(data[i]); } } return YAP_Unify(YAP_ARG2, tm); } static YAP_Bool is_matrix(void) { YAP_Term t = YAP_ARG1; int *mat = (int *)YAP_BlobOfTerm(t); if (!mat) { if (!YAP_IsApplTerm(t)) return FALSE; return YAP_FunctorOfTerm(t) == FunctorM; } return TRUE; } X_API void init_matrix( void ); X_API void init_matrix(void) { AtomC = YAP_LookupAtom("c"); FunctorM = YAP_MkFunctor(YAP_LookupAtom("$matrix"), 5); 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("matrixn_set", matrix_set, 3); YAP_UserCPredicate("matrix_set", matrix_set2, 2); YAP_UserCPredicate("matrix_set_all", matrix_set_all, 2); YAP_UserCPredicate("matrix_add", matrix_add, 3); YAP_UserCPredicate("matrixn_get", do_matrix_access, 3); YAP_UserCPredicate("matrixn_get", do_matrix_access2, 2); 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("matrixn_to_list", matrix_to_list, 2); YAP_UserCPredicate("matrixn_set_base", matrix_set_base, 2); YAP_UserCPredicate("matrixn_dims", matrix_dims, 2); YAP_UserCPredicate("matrixn_dims", matrix_dims3, 3); YAP_UserCPredicate("matrixn_ndims", matrix_ndims, 2); YAP_UserCPredicate("matrixn_size", matrix_size, 2); YAP_UserCPredicate("matrix_type_as_number", matrix_type, 2); YAP_UserCPredicate("matrixn_arg_to_offset", matrix_arg_to_offset, 3); YAP_UserCPredicate("matrixn_offset_to_arg", matrix_offset_to_arg, 3); YAP_UserCPredicate("matrixn_max", matrix_max, 2); YAP_UserCPredicate("matrixn_maxarg", matrix_maxarg, 2); YAP_UserCPredicate("matrixn_min", matrix_min, 2); YAP_UserCPredicate("matrixn_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("matrixn_to_logs", matrix_log_all2,2); YAP_UserCPredicate("matrixn_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); YAP_UserCPredicate("matrix_m", matrix_m, 2); YAP_UserCPredicate("matrix", is_matrix, 1); } #ifdef _WIN32 int WINAPI 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