/************************************************************************* * * * 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 "YapInterface.h" #include "config.h" #include #if defined(__MINGW32__) || _MSC_VER #include #endif #if HAVE_STRING_H #include #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) memmove((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) memmove((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; // function KahanSum(input) double c = 0.0; // A running compensation for lost low-order bits. for (i = 0; i < mat[MAT_SIZE]; i++) { double y = data[i] - c; // So far, so good: c is zero. double t = sum + y; // Alas, sum is big, y small, so low-order digits of y are lost. c = (t - sum) - y; // (t - sum) cancels the high-order part of y; // subtracting y recovers negative (low part of y) sum = t; // Algebraically, c should always be zero. Beware // overly-aggressive optimizing compilers! } 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; } static YAP_Bool get_float_from_address(void) { YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1); YAP_Int off = YAP_IntOfTerm(YAP_ARG2); return YAP_Unify(YAP_ARG3, YAP_MkFloatTerm(fp[off])); } static YAP_Bool set_float_from_address(void) { YAP_Float *fp = (YAP_Float *)YAP_IntOfTerm(YAP_ARG1); YAP_Int off = YAP_IntOfTerm(YAP_ARG2); YAP_Float f = YAP_FloatOfTerm(YAP_ARG3); fp[off] = f; 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); YAP_UserCPredicate("get_float_from_address", get_float_from_address, 3); YAP_UserCPredicate("set_float_from_address", set_float_from_address, 3); } #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