1 /*************************************************************************
2 * *
3 * YAP Prolog *
4 * *
5 * Yap Prolog was developed at NCCUP - Universidade do Porto *
6 * *
7 * Copyright L.Damas, V.S.Costa and Universidade do Porto 1985-1997 *
8 * *
9 **************************************************************************
10 * *
11 * File: matrix.c *
12 * Last rev: *
13 * mods: *
14 * comments: numerical arrays *
15 * *
16 *************************************************************************/
17
18 #include "config.h"
19 #include "YapInterface.h"
20 #include <math.h>
21 #if defined(__MINGW32__) || _MSC_VER
22 #include <windows.h>
23 #endif
24 #if HAVE_STRING_H
25 #include <string.h>
26 #endif
27
28 /*
29 A matrix is something of the form
30
31 TYPE = {int,double}
32 #DIMS = an int
33 DIM1
34 ...
35 DIMn
36 DATA in C format.
37
38 floating point matrixes may need to be aligned, so we always have an
39 extra element at the end.
40 */
41
42 /* maximal number of dimensions, 1024 should be enough */
43 #define MAX_DIMS 1024
44
45 typedef enum {
46 INT_MATRIX,
47 FLOAT_MATRIX
48 } mat_data_type;
49
50 typedef enum {
51 MAT_TYPE=0,
52 MAT_NDIMS=1,
53 MAT_SIZE=2,
54 MAT_ALIGN=3,
55 MAT_DIMS=4,
56 } mat_type;
57
58 typedef enum {
59 MAT_PLUS=0,
60 MAT_SUB=1,
61 MAT_TIMES=2,
62 MAT_DIV=3,
63 MAT_IDIV=4,
64 MAT_ZDIV=5,
65 MAT_LOG=6,
66 MAT_EXP=7
67 } op_type;
68
69 static long int *
matrix_long_data(int * mat,int ndims)70 matrix_long_data(int *mat, int ndims)
71 {
72 return (long int *)(mat+(MAT_DIMS+ndims));
73 }
74
75 static double *
matrix_double_data(int * mat,int ndims)76 matrix_double_data(int *mat, int ndims)
77 {
78 return (double *)(mat+(MAT_DIMS+ndims));
79 }
80
81 static unsigned int
matrix_get_offset(int * mat,int * indx)82 matrix_get_offset(int *mat, int* indx)
83 {
84 unsigned int i, pos = mat[MAT_SIZE], off = 0;
85
86 /* find where we are */
87 for (i = 0; i < mat[MAT_NDIMS]; i++) {
88 pos /= mat[MAT_DIMS+i];
89 if (indx[i] >= mat[MAT_DIMS+i]) {
90 return off;
91 }
92 off += pos*indx[i];
93 }
94 return off;
95 }
96
97 static void
matrix_get_index(int * mat,unsigned int offset,int * indx)98 matrix_get_index(int *mat, unsigned int offset, int* indx)
99 {
100 unsigned int i, pos = mat[MAT_SIZE];
101
102 /* find where we are */
103 for (i = 0; i < mat[MAT_NDIMS]; i++) {
104 pos /= mat[MAT_DIMS+i];
105 indx[i] = offset / pos;
106 offset = offset % pos;
107 }
108 }
109
110 static void
matrix_next_index(int * dims,int ndims,int * indx)111 matrix_next_index(int *dims, int ndims, int* indx)
112 {
113 unsigned int i;
114
115 /* find where we are */
116 for (i = ndims; i >0; ) {
117 i--;
118 indx[i]++;
119 if (indx[i]!=dims[i]) return;
120 indx[i] = 0;
121 }
122 }
123
124 static YAP_Term
new_int_matrix(int ndims,int dims[],long int data[])125 new_int_matrix(int ndims, int dims[], long int data[])
126 {
127 unsigned int sz;
128 unsigned int i, nelems=1;
129 YAP_Term blob;
130 int *mat;
131 long int *bdata;
132 int idims[MAX_DIMS];
133
134 /* in case we don't have enough room and need to shift the stack, we can't
135 really afford to keep a pointer to the global */
136 for (i=0;i< ndims;i++) {
137 idims[i] = dims[i];
138 nelems *= dims[i];
139 }
140 sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+nelems*sizeof(long int))/sizeof(YAP_CELL);
141 blob = YAP_MkBlobTerm(sz);
142 if (blob == YAP_TermNil()) {
143 return blob;
144 }
145 mat = (int *)YAP_BlobOfTerm(blob);
146 mat[MAT_TYPE] = INT_MATRIX;
147 mat[MAT_NDIMS] = ndims;
148 mat[MAT_SIZE] = nelems;
149 for (i=0;i< ndims;i++) {
150 mat[MAT_DIMS+i] = idims[i];
151 }
152 bdata = matrix_long_data(mat,ndims);
153 if (data)
154 memcpy((void *)bdata,(void *)data,sizeof(double)*nelems);
155 return blob;
156 }
157
158 static YAP_Term
new_float_matrix(int ndims,int dims[],double data[])159 new_float_matrix(int ndims, int dims[], double data[])
160 {
161 unsigned int sz;
162 unsigned int i, nelems=1;
163 YAP_Term blob;
164 int *mat;
165 double *bdata;
166 int idims[MAX_DIMS];
167
168 /* in case we don't have enough room and need to shift the stack, we can't
169 really afford to keep a pointer to the global */
170 for (i=0;i< ndims;i++) {
171 idims[i] = dims[i];
172 nelems *= dims[i];
173 }
174 sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+(nelems+1)*sizeof(double)+(sizeof(YAP_CELL)-1))/sizeof(YAP_CELL);
175 blob = YAP_MkBlobTerm(sz);
176 if (blob == YAP_TermNil())
177 return blob;
178 mat = YAP_BlobOfTerm(blob);
179 mat[MAT_TYPE] = FLOAT_MATRIX;
180 mat[MAT_NDIMS] = ndims;
181 mat[MAT_SIZE] = nelems;
182 for (i=0;i< ndims;i++) {
183 mat[MAT_DIMS+i] = idims[i];
184 }
185 bdata = matrix_double_data(mat,ndims);
186 if (data)
187 memcpy((void *)bdata,(void *)data,sizeof(double)*nelems);
188 return blob;
189 }
190
191 static int
scan_dims(int ndims,YAP_Term tl,int dims[MAX_DIMS])192 scan_dims(int ndims, YAP_Term tl, int dims[MAX_DIMS])
193 {
194 int i;
195
196 for (i = 0; i < ndims; i++) {
197 YAP_Term th;
198 int d;
199
200 if (!YAP_IsPairTerm(tl)) {
201 return FALSE;
202 }
203 th = YAP_HeadOfTerm(tl);
204 if (!YAP_IsIntTerm(th)) {
205 /* ERROR */
206 return FALSE;
207 }
208 d = YAP_IntOfTerm(th);
209 if (d < 0) {
210 /* ERROR */
211 return FALSE;
212 }
213 dims[i] = d;
214 tl = YAP_TailOfTerm(tl);
215 }
216 if (tl != YAP_TermNil()) {
217 /* ERROR */
218 return FALSE;
219 }
220 return TRUE;
221 }
222
223 static int
cp_int_matrix(YAP_Term tl,YAP_Term matrix)224 cp_int_matrix(YAP_Term tl,YAP_Term matrix)
225 {
226 int *mat = (int *)YAP_BlobOfTerm(matrix);
227 int i, nelems = mat[MAT_SIZE];
228 long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
229
230 for (i = 0; i < nelems; i++) {
231 YAP_Term th;
232 int d;
233
234 if (!YAP_IsPairTerm(tl)) {
235 return FALSE;
236 }
237 th = YAP_HeadOfTerm(tl);
238 if (!YAP_IsIntTerm(th)) {
239 /* ERROR */
240 return FALSE;
241 }
242 d = YAP_IntOfTerm(th);
243 j[i] = d;
244 tl = YAP_TailOfTerm(tl);
245 }
246 if (tl != YAP_TermNil()) {
247 /* ERROR */
248 return FALSE;
249 }
250 return TRUE;
251 }
252
253 static int
cp_float_matrix(YAP_Term tl,YAP_Term matrix)254 cp_float_matrix(YAP_Term tl,YAP_Term matrix)
255 {
256 int *mat = (int *)YAP_BlobOfTerm(matrix);
257 int i, nelems = mat[MAT_SIZE];
258 double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
259
260 for (i = 0; i < nelems; i++) {
261 YAP_Term th;
262 double d;
263
264 if (!YAP_IsPairTerm(tl)) {
265 return FALSE;
266 }
267 th = YAP_HeadOfTerm(tl);
268 if (YAP_IsIntTerm(th)) {
269 d = YAP_IntOfTerm(th);
270 } else if (!YAP_IsFloatTerm(th)) {
271 /* ERROR */
272 return FALSE;
273 } else {
274 d = YAP_FloatOfTerm(th);
275 }
276 j[i] = d;
277 tl = YAP_TailOfTerm(tl);
278 }
279 if (tl != YAP_TermNil()) {
280 /* ERROR */
281 return FALSE;
282 }
283 return TRUE;
284 }
285
286
287 static int
set_int_matrix(YAP_Term matrix,long int set)288 set_int_matrix(YAP_Term matrix,long int set)
289 {
290 int *mat = (int *)YAP_BlobOfTerm(matrix);
291 int i, nelems = mat[MAT_SIZE];
292 long int *j = matrix_long_data(mat, mat[MAT_NDIMS]);
293
294 for (i = 0; i < nelems; i++) {
295 j[i] = set;
296 }
297 return TRUE;
298 }
299
300 static int
set_float_matrix(YAP_Term matrix,double set)301 set_float_matrix(YAP_Term matrix,double set)
302 {
303 int *mat = (int *)YAP_BlobOfTerm(matrix);
304 int i, nelems = mat[MAT_SIZE];
305 double *j = matrix_double_data(mat, mat[MAT_NDIMS]);
306
307 for (i = 0; i < nelems; i++) {
308 j[i] = set;
309 }
310 return TRUE;
311 }
312
313 static int
new_ints_matrix(void)314 new_ints_matrix(void)
315 {
316 int ndims = YAP_IntOfTerm(YAP_ARG1);
317 YAP_Term tl = YAP_ARG2, out;
318 int dims[MAX_DIMS];
319
320 if (!scan_dims(ndims, tl, dims))
321 return FALSE;
322 out = new_int_matrix(ndims, dims, NULL);
323 if (out == YAP_TermNil())
324 return FALSE;
325 if (!cp_int_matrix(YAP_ARG3,out))
326 return FALSE;
327 return YAP_Unify(YAP_ARG4, out);
328 }
329
330 static int
new_ints_matrix_set(void)331 new_ints_matrix_set(void)
332 {
333 int ndims = YAP_IntOfTerm(YAP_ARG1);
334 YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
335 int dims[MAX_DIMS];
336 long int set;
337
338 if (!YAP_IsIntTerm(tset)) {
339 return FALSE;
340 }
341 set = YAP_IntOfTerm(tset);
342 if (!scan_dims(ndims, tl, dims))
343 return FALSE;
344 out = new_int_matrix(ndims, dims, NULL);
345 if (!set_int_matrix(out,set))
346 return FALSE;
347 return YAP_Unify(YAP_ARG4, out);
348 }
349
350 static int
new_floats_matrix(void)351 new_floats_matrix(void)
352 {
353 int ndims = YAP_IntOfTerm(YAP_ARG1);
354 YAP_Term tl = YAP_ARG2, out;
355 int dims[MAX_DIMS];
356 if (!scan_dims(ndims, tl, dims))
357 return FALSE;
358 out = new_float_matrix(ndims, dims, NULL);
359 if (out == YAP_TermNil())
360 return FALSE;
361 if (!cp_float_matrix(YAP_ARG3,out))
362 return FALSE;
363 return YAP_Unify(YAP_ARG4, out);
364 }
365
366 static int
new_floats_matrix_set(void)367 new_floats_matrix_set(void)
368 {
369 int ndims = YAP_IntOfTerm(YAP_ARG1);
370 YAP_Term tl = YAP_ARG2, out, tset = YAP_ARG3;
371 int dims[MAX_DIMS];
372 double set;
373
374 if (!YAP_IsFloatTerm(tset)) {
375 return FALSE;
376 }
377 set = YAP_FloatOfTerm(tset);
378 if (!scan_dims(ndims, tl, dims))
379 return FALSE;
380 out = new_float_matrix(ndims, dims, NULL);
381 if (!set_float_matrix(out,set))
382 return FALSE;
383 return YAP_Unify(YAP_ARG4, out);
384 }
385
386 static YAP_Term
float_matrix_to_list(int * mat)387 float_matrix_to_list(int *mat) {
388 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
389 int i = 0;
390 YAP_Term tf = YAP_TermNil(), tnil = tf;
391
392 for (i = mat[MAT_SIZE]-1; i>= 0; i--) {
393 tf = YAP_MkPairTerm(YAP_MkFloatTerm(data[i]),tf);
394 if (tf == tnil) {
395 /* error */
396 return tnil;
397 }
398 }
399 return tf;
400 }
401
402 static YAP_Term
mk_int_list(int nelems,int * data)403 mk_int_list(int nelems, int *data)
404 {
405 YAP_Term tn = YAP_TermNil();
406 YAP_Term tf = tn;
407 int i = 0;
408
409 for (i = nelems-1; i>= 0; i--) {
410 tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
411 if (tf == tn) {
412 /* error */
413 return tn;
414 }
415 }
416 return tf;
417 }
418
419 static YAP_Term
mk_long_list(int nelems,long int * data)420 mk_long_list(int nelems, long int *data)
421 {
422 YAP_Term tn = YAP_TermNil();
423 YAP_Term tf = tn;
424 int i = 0;
425
426 for (i = nelems-1; i>= 0; i--) {
427 tf = YAP_MkPairTerm(YAP_MkIntTerm(data[i]),tf);
428 if (tf == tn) {
429 /* error */
430 return tn;
431 }
432 }
433 return tf;
434 }
435
436
437 static YAP_Term
long_matrix_to_list(int * mat)438 long_matrix_to_list(int *mat) {
439 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
440
441 return mk_long_list(mat[MAT_SIZE], data);
442 }
443
444 static YAP_Term
matrix_access(int * mat,int * indx)445 matrix_access(int *mat, int *indx)
446 {
447 unsigned int off = matrix_get_offset(mat, indx);
448 if (mat[MAT_TYPE]==FLOAT_MATRIX)
449 return YAP_MkFloatTerm((matrix_double_data(mat,mat[MAT_NDIMS]))[off]);
450 else
451 return YAP_MkIntTerm((matrix_long_data(mat,mat[MAT_NDIMS]))[off]);
452 }
453
454 static void
matrix_float_set(int * mat,int * indx,double nval)455 matrix_float_set(int *mat, int *indx, double nval)
456 {
457 unsigned int off = 0;
458
459 off = matrix_get_offset(mat, indx);
460 (matrix_double_data(mat,mat[MAT_NDIMS]))[off] = nval;
461 }
462
463 static void
matrix_long_set(int * mat,int * indx,long int nval)464 matrix_long_set(int *mat, int *indx, long int nval)
465 {
466 unsigned int off = matrix_get_offset(mat, indx);
467 (matrix_long_data(mat,mat[MAT_NDIMS]))[off] = nval;
468 }
469
470 static void
matrix_float_set_all(int * mat,double nval)471 matrix_float_set_all(int *mat, double nval)
472 {
473 int i;
474 double *data = matrix_double_data(mat,mat[MAT_NDIMS]);
475
476 for (i = 0; i< mat[MAT_SIZE]; i++)
477 data[i] = nval;
478 }
479
480 static void
matrix_long_set_all(int * mat,long int nval)481 matrix_long_set_all(int *mat, long int nval)
482 {
483 int i;
484 long int *data = matrix_long_data(mat,mat[MAT_NDIMS]);
485
486 if (nval == 0) {
487 memset((void *)data,0,sizeof(long int)*mat[MAT_SIZE]);
488 } else {
489 for (i = 0; i< mat[MAT_SIZE]; i++)
490 data[i] = nval;
491 }
492 }
493
494 static void
matrix_float_add(int * mat,int * indx,double nval)495 matrix_float_add(int *mat, int *indx, double nval)
496 {
497 unsigned int off;
498 double *dat = matrix_double_data(mat,mat[MAT_NDIMS]);
499
500 off = matrix_get_offset(mat, indx);
501 dat[off] += nval;
502 }
503
504 static void
matrix_long_add(int * mat,int * indx,long int nval)505 matrix_long_add(int *mat, int *indx, long int nval)
506 {
507 long int *dat = matrix_long_data(mat,mat[MAT_NDIMS]);
508 unsigned int off = matrix_get_offset(mat, indx);
509 dat[off] += nval;
510 }
511
512 static void
matrix_inc(int * mat,int * indx)513 matrix_inc(int *mat, int *indx)
514 {
515 unsigned int off = matrix_get_offset(mat, indx);
516 if (mat[MAT_TYPE]==FLOAT_MATRIX)
517 (matrix_double_data(mat,mat[MAT_NDIMS])[off])++;
518 else
519 ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])++;
520 }
521
522 static void
matrix_dec(int * mat,int * indx)523 matrix_dec(int *mat, int *indx)
524 {
525 unsigned int off = matrix_get_offset(mat, indx);
526 if (mat[MAT_TYPE]==FLOAT_MATRIX)
527 (matrix_double_data(mat,mat[MAT_NDIMS])[off])--;
528 else
529 ((matrix_long_data(mat,mat[MAT_NDIMS]))[off])--;
530 }
531
532 static YAP_Term
matrix_inc2(int * mat,int * indx)533 matrix_inc2(int *mat, int *indx)
534 {
535 unsigned int off = matrix_get_offset(mat, indx);
536 if (mat[MAT_TYPE]==FLOAT_MATRIX) {
537 double *data = matrix_double_data(mat,mat[MAT_NDIMS]);
538 double d = data[off];
539 d++;
540 data[off] = d;
541 return YAP_MkFloatTerm(d);
542 } else {
543 long int *data = matrix_long_data(mat,mat[MAT_NDIMS]);
544 long int d = data[off];
545 d++;
546 data[off] = d;
547 return YAP_MkIntTerm(d);
548 }
549 }
550
551 static YAP_Term
matrix_dec2(int * mat,int * indx)552 matrix_dec2(int *mat, int *indx)
553 {
554 unsigned int off = matrix_get_offset(mat, indx);
555 if (mat[MAT_TYPE]==FLOAT_MATRIX) {
556 double *data = matrix_double_data(mat,mat[MAT_NDIMS]);
557 double d = data[off];
558 d--;
559 data[off] = d;
560 return YAP_MkFloatTerm(d);
561 } else {
562 long int *data = matrix_long_data(mat,mat[MAT_NDIMS]);
563 long int d = data[off];
564 d--;
565 data[off] = d;
566 return YAP_MkIntTerm(d);
567 }
568 }
569
570
571 static int
matrix_set(void)572 matrix_set(void)
573 {
574 int dims[MAX_DIMS], *mat;
575 YAP_Term tf;
576
577 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
578 if (!mat) {
579 /* Error */
580 return FALSE;
581 }
582 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
583 /* Error */
584 return FALSE;
585 }
586 tf = YAP_ARG3;
587 if (mat[MAT_TYPE] == INT_MATRIX) {
588 if (YAP_IsIntTerm(tf)) {
589 matrix_long_set(mat, dims, YAP_IntOfTerm(tf));
590 } else if (YAP_IsFloatTerm(tf)) {
591 matrix_long_set(mat, dims, YAP_FloatOfTerm(tf));
592 } else {
593 /* Error */
594 return FALSE;
595 }
596 } else {
597 if (YAP_IsIntTerm(tf)) {
598 matrix_float_set(mat, dims, YAP_IntOfTerm(tf));
599 } else if (YAP_IsFloatTerm(tf)) {
600 matrix_float_set(mat, dims, YAP_FloatOfTerm(tf));
601 } else {
602 /* Error */
603 return FALSE;
604 }
605 }
606 return TRUE;
607 }
608
609 static int
matrix_set_all(void)610 matrix_set_all(void)
611 {
612 int *mat;
613 YAP_Term tf;
614
615 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
616 if (!mat) {
617 /* Error */
618 return FALSE;
619 }
620 tf = YAP_ARG2;
621 if (mat[MAT_TYPE] == INT_MATRIX) {
622 if (YAP_IsIntTerm(tf)) {
623 matrix_long_set_all(mat, YAP_IntOfTerm(tf));
624 } else if (YAP_IsFloatTerm(tf)) {
625 matrix_long_set_all(mat, YAP_FloatOfTerm(tf));
626 } else {
627 /* Error */
628 return FALSE;
629 }
630 } else {
631 if (YAP_IsIntTerm(tf)) {
632 matrix_float_set_all(mat, YAP_IntOfTerm(tf));
633 } else if (YAP_IsFloatTerm(tf)) {
634 matrix_float_set_all(mat, YAP_FloatOfTerm(tf));
635 } else {
636 /* Error */
637 return FALSE;
638 }
639 }
640 return TRUE;
641 }
642
643 static int
matrix_add(void)644 matrix_add(void)
645 {
646 int dims[MAX_DIMS], *mat;
647 YAP_Term tf;
648
649 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
650 if (!mat) {
651 /* Error */
652 return FALSE;
653 }
654 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
655 /* Error */
656 return FALSE;
657 }
658 tf = YAP_ARG3;
659 if (mat[MAT_TYPE] == INT_MATRIX) {
660 if (YAP_IsIntTerm(tf)) {
661 matrix_long_add(mat, dims, YAP_IntOfTerm(tf));
662 } else if (YAP_IsFloatTerm(tf)) {
663 matrix_long_add(mat, dims, YAP_FloatOfTerm(tf));
664 } else {
665 /* Error */
666 return FALSE;
667 }
668 } else {
669 if (YAP_IsIntTerm(tf)) {
670 matrix_float_add(mat, dims, YAP_IntOfTerm(tf));
671 } else if (YAP_IsFloatTerm(tf)) {
672 matrix_float_add(mat, dims, YAP_FloatOfTerm(tf));
673 } else {
674 /* Error */
675 return FALSE;
676 }
677 }
678 return TRUE;
679 }
680
681 static int
do_matrix_access(void)682 do_matrix_access(void)
683 {
684 int dims[MAX_DIMS], *mat;
685 YAP_Term tf;
686
687 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
688 if (!mat) {
689 /* Error */
690 return FALSE;
691 }
692 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
693 /* Error */
694 return FALSE;
695 }
696 tf = matrix_access(mat, dims);
697 return YAP_Unify(tf, YAP_ARG3);
698 }
699
700 static int
do_matrix_inc(void)701 do_matrix_inc(void)
702 {
703 int dims[MAX_DIMS], *mat;
704
705 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
706 if (!mat) {
707 /* Error */
708 return FALSE;
709 }
710 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
711 /* Error */
712 return FALSE;
713 }
714 matrix_inc(mat, dims);
715 return TRUE;
716 }
717
718 static int
do_matrix_dec(void)719 do_matrix_dec(void)
720 {
721 int dims[MAX_DIMS], *mat;
722
723 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
724 if (!mat) {
725 /* Error */
726 return FALSE;
727 }
728 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
729 /* Error */
730 return FALSE;
731 }
732 matrix_dec(mat, dims);
733 return TRUE;
734 }
735
736 static int
do_matrix_inc2(void)737 do_matrix_inc2(void)
738 {
739 int dims[MAX_DIMS], *mat;
740
741 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
742 if (!mat) {
743 /* Error */
744 return FALSE;
745 }
746 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
747 /* Error */
748 return FALSE;
749 }
750 return
751 YAP_Unify(matrix_inc2(mat, dims), YAP_ARG3);
752 }
753
754 static int
do_matrix_dec2(void)755 do_matrix_dec2(void)
756 {
757 int dims[MAX_DIMS], *mat;
758
759 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
760 if (!mat) {
761 /* Error */
762 return FALSE;
763 }
764 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, dims)) {
765 /* Error */
766 return FALSE;
767 }
768 return
769 YAP_Unify(matrix_dec2(mat, dims), YAP_ARG3);
770 }
771
772 static int
matrix_to_list(void)773 matrix_to_list(void)
774 {
775 int *mat;
776 YAP_Term tf;
777
778 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
779 if (!mat) {
780 /* Error */
781 return FALSE;
782 }
783 if (mat[MAT_TYPE] == INT_MATRIX)
784 tf = long_matrix_to_list(mat);
785 else
786 tf = float_matrix_to_list(mat);
787 return YAP_Unify(YAP_ARG2, tf);
788 }
789
790 static int
matrix_dims(void)791 matrix_dims(void)
792 {
793 int *mat;
794 YAP_Term tf;
795
796 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
797 if (!mat) {
798 /* Error */
799 return FALSE;
800 }
801 tf = mk_int_list(mat[MAT_NDIMS],mat+MAT_DIMS);
802 return YAP_Unify(YAP_ARG2, tf);
803 }
804
805 static int
matrix_size(void)806 matrix_size(void)
807 {
808 int *mat;
809
810 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
811 if (!mat) {
812 /* Error */
813 return FALSE;
814 }
815 return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_SIZE]));
816 }
817
818 static int
matrix_ndims(void)819 matrix_ndims(void)
820 {
821 int *mat;
822
823 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
824 if (!mat) {
825 /* Error */
826 return FALSE;
827 }
828 return YAP_Unify(YAP_ARG2, YAP_MkIntTerm(mat[MAT_NDIMS]));
829 }
830
831 static int
matrix_type(void)832 matrix_type(void)
833 {
834 int *mat;
835 YAP_Term tf;
836
837 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
838 if (!mat) {
839 /* Error */
840 return FALSE;
841 }
842 if (mat[MAT_TYPE] == INT_MATRIX) {
843 tf = YAP_MkIntTerm(0);
844 } else {
845 tf = YAP_MkIntTerm(1);
846 }
847 return YAP_Unify(YAP_ARG2, tf);
848 }
849
850 static int
matrix_arg_to_offset(void)851 matrix_arg_to_offset(void)
852 {
853 int indx[MAX_DIMS], *mat;
854 unsigned int off;
855
856 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
857 if (!mat) {
858 /* Error */
859 return FALSE;
860 }
861 if (!scan_dims(mat[MAT_NDIMS], YAP_ARG2, indx)) {
862 /* Error */
863 return FALSE;
864 }
865 off = matrix_get_offset(mat, indx);
866
867 return YAP_Unify(YAP_ARG3, YAP_MkIntTerm(off));
868 }
869
870 static int
matrix_offset_to_arg(void)871 matrix_offset_to_arg(void)
872 {
873 int indx[MAX_DIMS], *mat;
874 unsigned int off;
875 YAP_Term ti, tf;
876
877 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
878 if (!mat) {
879 /* Error */
880 return FALSE;
881 }
882 if (!YAP_IsIntTerm(ti = YAP_ARG2)) {
883 /* Error */
884 return FALSE;
885 }
886 off = YAP_IntOfTerm(ti);
887 matrix_get_index(mat, off, indx);
888 tf = mk_int_list(mat[MAT_NDIMS], indx);
889 return YAP_Unify(YAP_ARG3, tf);
890 }
891
892 static unsigned int
scan_max_long(int sz,long int * data)893 scan_max_long(int sz, long int *data)
894 {
895 int i, off=0;
896 long int max= data[0];
897 for (i=1; i<sz; i++) {
898 if (data[i]>max) {
899 off=i;
900 max = data[i];
901 }
902 }
903 return off;
904 }
905
906 static unsigned int
scan_max_float(int sz,double * data)907 scan_max_float(int sz, double *data)
908 {
909 int i, off=0;
910 double max= data[0];
911 for (i=1; i<sz; i++) {
912 if (data[i]>max) {
913 max = data[i];
914 off=i;
915 }
916 }
917 return off;
918 }
919
920 static unsigned int
scan_min_long(int sz,long int * data)921 scan_min_long(int sz, long int *data)
922 {
923 int i, off=0;
924 long int max= data[0];
925 for (i=1; i<sz; i++) {
926 if (data[i]<max) {
927 max = data[i];
928 off=i;
929 }
930 }
931 return off;
932 }
933
934 static unsigned int
scan_min_float(int sz,double * data)935 scan_min_float(int sz, double *data)
936 {
937 int i, off=0;
938 double max= data[0];
939 for (i=1; i<sz; i++) {
940 if (data[i]<max) {
941 max = data[i];
942 off=i;
943 }
944 }
945 return off;
946 }
947
948 static int
matrix_max(void)949 matrix_max(void)
950 {
951 int *mat;
952 unsigned int off;
953 YAP_Term tf;
954
955 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
956 if (!mat) {
957 /* Error */
958 return FALSE;
959 }
960 if (mat[MAT_TYPE] == INT_MATRIX) {
961 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
962 off = scan_max_long(mat[MAT_SIZE], data);
963 tf = YAP_MkIntTerm(data[off]);
964 } else {
965 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
966 off = scan_max_float(mat[MAT_SIZE], data);
967 tf = YAP_MkFloatTerm(data[off]);
968 }
969 return YAP_Unify(YAP_ARG2, tf);
970 }
971
972 static int
matrix_maxarg(void)973 matrix_maxarg(void)
974 {
975 int indx[MAX_DIMS], *mat;
976 unsigned int off;
977 YAP_Term tf;
978
979 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
980 if (!mat) {
981 /* Error */
982 return FALSE;
983 }
984 if (mat[MAT_TYPE] == INT_MATRIX) {
985 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
986 off = scan_max_long(mat[MAT_SIZE], data);
987 } else {
988 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
989 off = scan_max_float(mat[MAT_SIZE], data);
990 }
991 matrix_get_index(mat, off, indx);
992 tf = mk_int_list(mat[MAT_NDIMS], indx);
993 return YAP_Unify(YAP_ARG2, tf);
994 }
995
996 static int
matrix_min(void)997 matrix_min(void)
998 {
999 int *mat;
1000 unsigned int off;
1001 YAP_Term tf;
1002
1003 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1004 if (!mat) {
1005 /* Error */
1006 return FALSE;
1007 }
1008 if (mat[MAT_TYPE] == INT_MATRIX) {
1009 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1010 off = scan_min_long(mat[MAT_SIZE], data);
1011 tf = YAP_MkIntTerm(data[off]);
1012 } else {
1013 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1014 off = scan_min_float(mat[MAT_SIZE], data);
1015 tf = YAP_MkFloatTerm(data[off]);
1016 }
1017 return YAP_Unify(YAP_ARG2, tf);
1018 }
1019
1020 static int
matrix_log_all(void)1021 matrix_log_all(void)
1022 {
1023 int *mat;
1024
1025 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1026 if (!mat) {
1027 /* Error */
1028 return FALSE;
1029 }
1030 if (mat[MAT_TYPE] == INT_MATRIX) {
1031 return FALSE;
1032 } else {
1033 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1034 int i;
1035
1036 for (i=0; i< mat[MAT_SIZE]; i++) {
1037 data[i] = log(data[i]);
1038 }
1039 }
1040 return TRUE;
1041 }
1042
1043 static int
matrix_log_all2(void)1044 matrix_log_all2(void)
1045 {
1046 int *mat;
1047
1048 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1049 if (!mat) {
1050 /* Error */
1051 return FALSE;
1052 }
1053 if (mat[MAT_TYPE] == INT_MATRIX) {
1054 return FALSE;
1055 } else {
1056 YAP_Term out;
1057 double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
1058 int i;
1059 int *nmat;
1060
1061 if (!YAP_IsVarTerm(YAP_ARG2)) {
1062 out = YAP_ARG2;
1063 } else {
1064 out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
1065 if (out == YAP_TermNil())
1066 return FALSE;
1067 }
1068 nmat = (int *)YAP_BlobOfTerm(out);
1069 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1070 for (i=0; i< mat[MAT_SIZE]; i++) {
1071 ndata[i] = log(data[i]);
1072 }
1073 if (YAP_IsVarTerm(YAP_ARG2)) {
1074 return YAP_Unify(YAP_ARG2, out);
1075 }
1076 }
1077 return TRUE;
1078 }
1079
1080 static int
matrix_exp_all(void)1081 matrix_exp_all(void)
1082 {
1083 int *mat;
1084
1085 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1086 if (!mat) {
1087 /* Error */
1088 return FALSE;
1089 }
1090 if (mat[MAT_TYPE] == INT_MATRIX) {
1091 return FALSE;
1092 } else {
1093 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1094 int i;
1095
1096 for (i=0; i< mat[MAT_SIZE]; i++) {
1097 data[i] = exp(data[i]);
1098 }
1099 }
1100 return TRUE;
1101 }
1102
1103 static int
matrix_exp2_all(void)1104 matrix_exp2_all(void)
1105 {
1106 int *mat;
1107
1108 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1109 if (!mat) {
1110 /* Error */
1111 return FALSE;
1112 }
1113 if (mat[MAT_TYPE] == INT_MATRIX) {
1114 return FALSE;
1115 } else {
1116 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1117 int i;
1118 double max = data[0];
1119
1120 for (i=1; i< mat[MAT_SIZE]; i++) {
1121 if (data[i] > max) max = data[i];
1122 }
1123
1124 for (i=0; i< mat[MAT_SIZE]; i++) {
1125 data[i] = exp(data[i]-max);
1126 }
1127 }
1128 return TRUE;
1129 }
1130
1131 static int
matrix_exp_all2(void)1132 matrix_exp_all2(void)
1133 {
1134 int *mat;
1135
1136 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1137 if (!mat) {
1138 /* Error */
1139 return FALSE;
1140 }
1141 if (mat[MAT_TYPE] == INT_MATRIX) {
1142 return FALSE;
1143 } else {
1144 YAP_Term out;
1145 double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
1146 int i;
1147 int *nmat;
1148
1149 if (!YAP_IsVarTerm(YAP_ARG2)) {
1150 out = YAP_ARG2;
1151 } else {
1152 out = new_float_matrix(mat[MAT_NDIMS], mat+MAT_DIMS, NULL);
1153 if (out == YAP_TermNil())
1154 return FALSE;
1155 }
1156 nmat = (int *)YAP_BlobOfTerm(out);
1157 ndata = matrix_double_data(nmat, mat[MAT_NDIMS]);
1158 for (i=0; i< mat[MAT_SIZE]; i++) {
1159 ndata[i] = exp(data[i]);
1160 }
1161 if (YAP_IsVarTerm(YAP_ARG2)) {
1162 return YAP_Unify(YAP_ARG2, out);
1163 }
1164 }
1165 return TRUE;
1166 }
1167
1168 static int
matrix_minarg(void)1169 matrix_minarg(void)
1170 {
1171 int indx[MAX_DIMS], *mat;
1172 unsigned int off;
1173 YAP_Term tf;
1174
1175 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1176 if (!mat) {
1177 /* Error */
1178 return FALSE;
1179 }
1180 if (mat[MAT_TYPE] == INT_MATRIX) {
1181 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1182 off = scan_min_long(mat[MAT_SIZE], data);
1183 } else {
1184 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1185 off = scan_min_float(mat[MAT_SIZE], data);
1186 }
1187 matrix_get_index(mat, off, indx);
1188 tf = mk_int_list(mat[MAT_NDIMS], indx);
1189 return YAP_Unify(YAP_ARG2, tf);
1190 }
1191
1192 static int
matrix_sum(void)1193 matrix_sum(void)
1194 {
1195 int *mat;
1196 YAP_Term tf;
1197
1198 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1199 if (!mat) {
1200 /* Error */
1201 return FALSE;
1202 }
1203 if (mat[MAT_TYPE] == INT_MATRIX) {
1204 long int *data = matrix_long_data(mat, mat[MAT_NDIMS]);
1205 int i;
1206 long int sum = 0;
1207
1208 for (i = 0; i < mat[MAT_SIZE]; i++) {
1209 sum += data[i];
1210 }
1211 tf = YAP_MkIntTerm(sum);
1212 } else {
1213 double *data = matrix_double_data(mat, mat[MAT_NDIMS]);
1214 int i;
1215 double sum = 0.0;
1216
1217 for (i = 0; i < mat[MAT_SIZE]; i++) {
1218 sum += data[i];
1219 }
1220 tf = YAP_MkFloatTerm(sum);
1221 }
1222 return YAP_Unify(YAP_ARG2, tf);
1223 }
1224
1225 static void
add_int_lines(int total,int nlines,long int * mat0,long int * matf)1226 add_int_lines(int total,int nlines,long int *mat0,long int *matf)
1227 {
1228 int ncols = total/nlines, i;
1229 for (i=0;i<ncols;i++) {
1230 long int sum = 0;
1231 int j;
1232
1233 for (j=i;j<total;j+=ncols) {
1234 sum += mat0[j];
1235 }
1236 matf[i] = sum;
1237 }
1238 }
1239
1240 static void
add_double_lines(int total,int nlines,double * mat0,double * matf)1241 add_double_lines(int total,int nlines,double *mat0,double *matf)
1242 {
1243 int ncols = total/nlines, i;
1244 for (i=0;i<ncols;i++) {
1245 double sum = 0;
1246 int j;
1247
1248 for (j=i;j<total;j+=ncols) {
1249 sum += mat0[j];
1250 }
1251 matf[i] = sum;
1252 }
1253 }
1254
1255 static int
matrix_agg_lines(void)1256 matrix_agg_lines(void)
1257 {
1258 int *mat;
1259 YAP_Term tf;
1260 YAP_Term top = YAP_ARG2;
1261 op_type op;
1262
1263 if (!YAP_IsIntTerm(top)) {
1264 return FALSE;
1265 }
1266 op = YAP_IntOfTerm(top);
1267 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1268 if (!mat) {
1269 /* Error */
1270 return FALSE;
1271 }
1272 /* create a new array without first dimension */
1273 if (mat[MAT_TYPE] == INT_MATRIX) {
1274 long int *data, *ndata;
1275 int dims = mat[MAT_NDIMS];
1276 int *nmat;
1277
1278 tf = new_int_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
1279 if (tf == YAP_TermNil())
1280 return FALSE;
1281 nmat = (int *)YAP_BlobOfTerm(tf);
1282 data = matrix_long_data(mat, dims);
1283 ndata = matrix_long_data(nmat, dims-1);
1284 if (op == MAT_PLUS) {
1285 add_int_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
1286 } else
1287 return FALSE;
1288 } else {
1289 double *data, *ndata;
1290 int dims = mat[MAT_NDIMS];
1291 int *nmat;
1292
1293 tf = new_float_matrix(dims-1,mat+(MAT_DIMS+1),NULL);
1294 nmat = (int *)YAP_BlobOfTerm(tf);
1295 if (tf == YAP_TermNil())
1296 return FALSE;
1297 data = matrix_double_data(mat, dims);
1298 ndata = matrix_double_data(nmat, dims-1);
1299 if (op == MAT_PLUS) {
1300 add_double_lines(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
1301 } else
1302 return FALSE;
1303 }
1304 return YAP_Unify(YAP_ARG3,tf);
1305 }
1306
1307 static void
add_int_cols(int total,int nlines,long int * mat0,long int * matf)1308 add_int_cols(int total,int nlines,long int *mat0,long int *matf)
1309 {
1310 int ncols = total/nlines, i, j = 0;
1311 for (i=0;i<nlines;i++) {
1312 long int sum = 0;
1313 int max = (i+1)*ncols;
1314
1315 for (;j<max;j++) {
1316 sum += mat0[j];
1317 }
1318 matf[i] = sum;
1319 }
1320 }
1321
1322 static void
add_double_cols(int total,int nlines,double * mat0,double * matf)1323 add_double_cols(int total,int nlines,double *mat0,double *matf)
1324 {
1325 int ncols = total/nlines, i, j = 0;
1326 for (i=0;i<nlines;i++) {
1327 double sum = 0;
1328 int max = (i+1)*ncols;
1329
1330 for (;j<max;j++) {
1331 sum += mat0[j];
1332 }
1333 matf[i] = sum;
1334 }
1335 }
1336
1337 static int
matrix_agg_cols(void)1338 matrix_agg_cols(void)
1339 {
1340 int *mat;
1341 YAP_Term tf;
1342 YAP_Term top = YAP_ARG2;
1343 op_type op;
1344
1345 if (!YAP_IsIntTerm(top)) {
1346 return FALSE;
1347 }
1348 op = YAP_IntOfTerm(top);
1349 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
1350 if (!mat) {
1351 /* Error */
1352 return FALSE;
1353 }
1354 /* create a new array without first dimension */
1355 if (mat[MAT_TYPE] == INT_MATRIX) {
1356 long int *data, *ndata;
1357 int dims = mat[MAT_NDIMS];
1358 int *nmat;
1359
1360 tf = new_int_matrix(1,mat+MAT_DIMS,NULL);
1361 if (tf == YAP_TermNil())
1362 return FALSE;
1363 nmat = (int *)YAP_BlobOfTerm(tf);
1364 data = matrix_long_data(mat, dims);
1365 ndata = matrix_long_data(nmat, 1);
1366 if (op == MAT_PLUS) {
1367 add_int_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
1368 } else
1369 return FALSE;
1370 } else {
1371 double *data, *ndata;
1372 int dims = mat[MAT_NDIMS];
1373 int *nmat;
1374
1375 tf = new_float_matrix(1,mat+MAT_DIMS,NULL);
1376 if (tf == YAP_TermNil())
1377 return FALSE;
1378 nmat = (int *)YAP_BlobOfTerm(tf);
1379 data = matrix_double_data(mat, dims);
1380 ndata = matrix_double_data(nmat, 1);
1381 if (op == MAT_PLUS) {
1382 add_double_cols(mat[MAT_SIZE],mat[MAT_DIMS],data,ndata);
1383 } else
1384 return FALSE;
1385 }
1386 return YAP_Unify(YAP_ARG3,tf);
1387 }
1388
1389 static void
div_int_by_lines(int total,int nlines,long int * mat1,long int * mat2,double * ndata)1390 div_int_by_lines(int total,int nlines,long int *mat1,long int *mat2,double *ndata)
1391 {
1392 int ncols = total/nlines, i;
1393 for (i=0;i<total;i++) {
1394 ndata[i] = ((double)mat1[i])/mat2[i%ncols];
1395 }
1396 }
1397
1398 static void
div_int_by_dlines(int total,int nlines,long int * mat1,double * mat2,double * ndata)1399 div_int_by_dlines(int total,int nlines,long int *mat1,double *mat2,double *ndata)
1400 {
1401 int ncols = total/nlines, i;
1402 for (i=0;i<total;i++) {
1403 ndata[i] = mat1[i]/mat2[i%ncols];
1404 }
1405 }
1406
1407 static void
div_float_long_by_lines(int total,int nlines,double * mat1,long int * mat2,double * ndata)1408 div_float_long_by_lines(int total,int nlines,double *mat1,long int *mat2,double *ndata)
1409 {
1410 int ncols = total/nlines, i;
1411 for (i=0;i<total;i++) {
1412 ndata[i] = mat1[i]/mat2[i%ncols];
1413 }
1414 }
1415
1416 static void
div_float_by_lines(int total,int nlines,double * mat1,double * mat2,double * ndata)1417 div_float_by_lines(int total,int nlines,double *mat1,double *mat2,double *ndata)
1418 {
1419 int ncols = total/nlines, i;
1420 for (i=0;i<total;i++) {
1421 ndata[i] = mat1[i]/mat2[i%ncols];
1422 }
1423 }
1424
1425 static int
matrix_op_to_lines(void)1426 matrix_op_to_lines(void)
1427 {
1428 int *mat1, *mat2;
1429 YAP_Term top = YAP_ARG3;
1430 op_type op;
1431 YAP_Term tf;
1432
1433 if (!YAP_IsIntTerm(top)) {
1434 return FALSE;
1435 }
1436 op = YAP_IntOfTerm(top);
1437 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1438 if (!mat1) {
1439 /* Error */
1440 return FALSE;
1441 }
1442 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1443 if (!mat2) {
1444 /* Error */
1445 return FALSE;
1446 }
1447 /* create a new array without first dimension */
1448 if (mat1[MAT_TYPE] == INT_MATRIX) {
1449 long int *data1;
1450 int dims = mat1[MAT_NDIMS];
1451 int *nmat;
1452 data1 = matrix_long_data(mat1, dims);
1453
1454 if (mat2[MAT_TYPE] == INT_MATRIX) {
1455 long int *data2 = matrix_long_data(mat2, dims-1);
1456 if (op == MAT_DIV) {
1457 double *ndata;
1458
1459 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1460 if (tf == YAP_TermNil())
1461 return FALSE;
1462 nmat = YAP_BlobOfTerm(tf);
1463 ndata = matrix_double_data(nmat, dims);
1464 div_int_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1465 } else {
1466 return FALSE;
1467 }
1468 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1469 double *data2 = matrix_double_data(mat2, dims-1);
1470 if (op == MAT_DIV) {
1471 double *ndata;
1472
1473 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1474 if (tf == YAP_TermNil())
1475 return FALSE;
1476 nmat = YAP_BlobOfTerm(tf);
1477 ndata = matrix_double_data(nmat, dims);
1478 div_int_by_dlines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1479 } else {
1480 return FALSE;
1481 }
1482 } else {
1483 return FALSE;
1484 }
1485 } else {
1486 double *data1, *ndata;
1487 int dims = mat1[MAT_NDIMS];
1488 int *nmat;
1489
1490 data1 = matrix_double_data(mat1, dims);
1491 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1492 nmat = YAP_BlobOfTerm(tf);
1493 if (tf == YAP_TermNil())
1494 return FALSE;
1495 ndata = matrix_double_data(nmat, dims);
1496 if (mat2[MAT_TYPE] == INT_MATRIX) {
1497 long int *data2 = matrix_long_data(mat2, dims-1);
1498 if (op == MAT_DIV) {
1499 div_float_long_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1500 } else {
1501 return FALSE;
1502 }
1503 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1504 double *data2 = matrix_double_data(mat2, dims-1);
1505 if (op == MAT_DIV) {
1506 div_float_by_lines(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1507 } else {
1508 return FALSE;
1509 }
1510 } else {
1511 return FALSE;
1512 }
1513 }
1514 return YAP_Unify(YAP_ARG4,tf);
1515 }
1516
1517
1518 static void
matrix_long_add_data(long int * nmat,int siz,long int mat1[],long int mat2[])1519 matrix_long_add_data(long int *nmat, int siz, long int mat1[], long int mat2[])
1520 {
1521 int i;
1522
1523 for (i=0; i< siz; i++) {
1524 nmat[i] = mat1[i]+mat2[i];
1525 }
1526 }
1527
1528 static void
matrix_long_double_add_data(double * nmat,int siz,long int mat1[],double mat2[])1529 matrix_long_double_add_data(double *nmat, int siz, long int mat1[], double mat2[])
1530 {
1531 int i;
1532
1533 for (i=0; i< siz; i++) {
1534 nmat[i] = mat1[i]+mat2[i];
1535 }
1536 }
1537
1538 static void
matrix_double_add_data(double * nmat,int siz,double mat1[],double mat2[])1539 matrix_double_add_data(double *nmat, int siz, double mat1[], double mat2[])
1540 {
1541 int i;
1542
1543 for (i=0; i< siz; i++) {
1544 nmat[i] = mat1[i]+mat2[i];
1545 }
1546 }
1547
1548 static void
matrix_long_sub_data(long int * nmat,int siz,long int mat1[],long int mat2[])1549 matrix_long_sub_data(long int *nmat, int siz, long int mat1[], long int mat2[])
1550 {
1551 int i;
1552
1553 for (i=0; i< siz; i++) {
1554 nmat[i] = mat1[i]-mat2[i];
1555 }
1556 }
1557
1558 static void
matrix_long_double_sub_data(double * nmat,int siz,long int mat1[],double mat2[])1559 matrix_long_double_sub_data(double *nmat, int siz, long int mat1[], double mat2[])
1560 {
1561 int i;
1562
1563 for (i=0; i< siz; i++) {
1564 nmat[i] = mat1[i]-mat2[i];
1565 }
1566 }
1567
1568 static void
matrix_long_double_rsub_data(double * nmat,int siz,double mat1[],long int mat2[])1569 matrix_long_double_rsub_data(double *nmat, int siz, double mat1[], long int mat2[])
1570 {
1571 int i;
1572
1573 for (i=0; i< siz; i++) {
1574 nmat[i] = mat2[i]-mat1[i];
1575 }
1576 }
1577
1578 static void
matrix_double_sub_data(double * nmat,int siz,double mat1[],double mat2[])1579 matrix_double_sub_data(double *nmat, int siz, double mat1[], double mat2[])
1580 {
1581 int i;
1582
1583 for (i=0; i< siz; i++) {
1584 nmat[i] = mat1[i]-mat2[i];
1585 }
1586 }
1587
1588 static void
matrix_long_mult_data(long int * nmat,int siz,long int mat1[],long int mat2[])1589 matrix_long_mult_data(long int *nmat, int siz, long int mat1[], long int mat2[])
1590 {
1591 int i;
1592
1593 for (i=0; i< siz; i++) {
1594 nmat[i] = mat1[i]*mat2[i];
1595 }
1596 }
1597
1598 static void
matrix_long_double_mult_data(double * nmat,int siz,long int mat1[],double mat2[])1599 matrix_long_double_mult_data(double *nmat, int siz, long int mat1[], double mat2[])
1600 {
1601 int i;
1602
1603 for (i=0; i< siz; i++) {
1604 nmat[i] = mat1[i]*mat2[i];
1605 }
1606 }
1607
1608 static void
matrix_double_mult_data(double * nmat,int siz,double mat1[],double mat2[])1609 matrix_double_mult_data(double *nmat, int siz, double mat1[], double mat2[])
1610 {
1611 int i;
1612
1613 for (i=0; i< siz; i++) {
1614 nmat[i] = mat1[i]*mat2[i];
1615 }
1616 }
1617
1618 static void
matrix_long_div_data(long int * nmat,int siz,long int mat1[],long int mat2[])1619 matrix_long_div_data(long int *nmat, int siz, long int mat1[], long int mat2[])
1620 {
1621 int i;
1622
1623 for (i=0; i< siz; i++) {
1624 nmat[i] = mat1[i]/mat2[i];
1625 }
1626 }
1627
1628 static void
matrix_long_double_div_data(double * nmat,int siz,long int mat1[],double mat2[])1629 matrix_long_double_div_data(double *nmat, int siz, long int mat1[], double mat2[])
1630 {
1631 int i;
1632
1633 for (i=0; i< siz; i++) {
1634 nmat[i] = mat1[i]/mat2[i];
1635 }
1636 }
1637
1638 static void
matrix_long_double_div2_data(double * nmat,int siz,double mat1[],long int mat2[])1639 matrix_long_double_div2_data(double *nmat, int siz, double mat1[], long int mat2[])
1640 {
1641 int i;
1642
1643 for (i=0; i< siz; i++) {
1644 nmat[i] = mat1[i]/mat2[i];
1645 }
1646 }
1647
1648 static void
matrix_double_div_data(double * nmat,int siz,double mat1[],double mat2[])1649 matrix_double_div_data(double *nmat, int siz, double mat1[], double mat2[])
1650 {
1651 int i;
1652
1653 for (i=0; i< siz; i++) {
1654 nmat[i] = mat1[i]/mat2[i];
1655 }
1656 }
1657
1658 static void
matrix_long_zdiv_data(long int * nmat,int siz,long int mat1[],long int mat2[])1659 matrix_long_zdiv_data(long int *nmat, int siz, long int mat1[], long int mat2[])
1660 {
1661 int i;
1662
1663 for (i=0; i< siz; i++) {
1664 if (mat1[i] == 0)
1665 nmat[i] = 0;
1666 else
1667 nmat[i] = mat1[i]/mat2[i];
1668 }
1669 }
1670
1671 static void
matrix_long_double_zdiv_data(double * nmat,int siz,long int mat1[],double mat2[])1672 matrix_long_double_zdiv_data(double *nmat, int siz, long int mat1[], double mat2[])
1673 {
1674 int i;
1675
1676 for (i=0; i< siz; i++) {
1677 if (mat1[i] == 0)
1678 nmat[i] = 0;
1679 else
1680 nmat[i] = mat1[i]/mat2[i];
1681 }
1682 }
1683
1684 static void
matrix_long_double_zdiv2_data(double * nmat,int siz,double mat1[],long int mat2[])1685 matrix_long_double_zdiv2_data(double *nmat, int siz, double mat1[], long int mat2[])
1686 {
1687 int i;
1688
1689 for (i=0; i< siz; i++) {
1690 if (mat1[i] == 0.0)
1691 nmat[i] = 0;
1692 else
1693 nmat[i] = mat1[i]/mat2[i];
1694 }
1695 }
1696
1697 static void
matrix_double_zdiv_data(double * nmat,int siz,double mat1[],double mat2[])1698 matrix_double_zdiv_data(double *nmat, int siz, double mat1[], double mat2[])
1699 {
1700 int i;
1701
1702 for (i=0; i< siz; i++) {
1703 if (mat1[i] == 0.0) {
1704 nmat[i] = 0.0;
1705 } else {
1706 nmat[i] = mat1[i]/mat2[i];
1707 }
1708 }
1709 }
1710
1711 static int
matrix_op(void)1712 matrix_op(void)
1713 {
1714 int *mat1, *mat2;
1715 YAP_Term top = YAP_ARG3;
1716 op_type op;
1717 YAP_Term tf = YAP_ARG4;
1718 int create = TRUE;
1719
1720 if (!YAP_IsIntTerm(top)) {
1721 return FALSE;
1722 }
1723 op = YAP_IntOfTerm(top);
1724 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1725 if (!mat1) {
1726 /* Error */
1727 return FALSE;
1728 }
1729 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1730 if (!mat2) {
1731 /* Error */
1732 return FALSE;
1733 }
1734 if (tf == YAP_ARG1 || tf == YAP_ARG2) {
1735 create = FALSE;
1736 }
1737 if (mat1[MAT_TYPE] == INT_MATRIX) {
1738 long int *data1;
1739 int dims = mat1[MAT_NDIMS];
1740 int *nmat;
1741 data1 = matrix_long_data(mat1, dims);
1742
1743 if (mat2[MAT_TYPE] == INT_MATRIX) {
1744 long int *data2;
1745 long int *ndata;
1746
1747 if (create)
1748 tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
1749 if (tf == YAP_TermNil()) {
1750 return FALSE;
1751 } else {
1752 /* there may have been an overflow */
1753 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1754 data1 = matrix_long_data(mat1, dims);
1755 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1756 data2 = matrix_long_data(mat2, dims);
1757 }
1758 nmat = YAP_BlobOfTerm(tf);
1759 ndata = matrix_long_data(nmat, dims);
1760 switch (op) {
1761 case MAT_PLUS:
1762 matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1763 break;
1764 case MAT_SUB:
1765 matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1766 break;
1767 case MAT_TIMES:
1768 matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1769 break;
1770 case MAT_DIV:
1771 matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1772 break;
1773 case MAT_ZDIV:
1774 matrix_long_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1775 break;
1776 default:
1777 return FALSE;
1778 }
1779 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1780 double *data2;
1781 double *ndata;
1782
1783 if (create)
1784 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1785 if (tf == YAP_TermNil()) {
1786 return FALSE;
1787 } else {
1788 /* there may have been an overflow */
1789 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1790 data1 = matrix_long_data(mat1, dims);
1791 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1792 data2 = matrix_double_data(mat2, dims);
1793 }
1794 nmat = YAP_BlobOfTerm(tf);
1795 ndata = matrix_double_data(nmat, dims);
1796 switch (op) {
1797 case MAT_PLUS:
1798 matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1799 break;
1800 case MAT_SUB:
1801 matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1802 break;
1803 case MAT_TIMES:
1804 matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1805 break;
1806 case MAT_DIV:
1807 matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1808 break;
1809 case MAT_ZDIV:
1810 matrix_long_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1811 break;
1812 default:
1813 return FALSE;
1814 }
1815 } else {
1816 return FALSE;
1817 }
1818 } else {
1819 double *data1;
1820 int dims = mat1[MAT_NDIMS];
1821 int *nmat;
1822 data1 = matrix_double_data(mat1, dims);
1823
1824 if (mat2[MAT_TYPE] == INT_MATRIX) {
1825 long int *data2;
1826 double *ndata;
1827
1828 if (create)
1829 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1830 if (tf == YAP_TermNil()) {
1831 return FALSE;
1832 } else {
1833 /* there may have been an overflow */
1834 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1835 data1 = matrix_double_data(mat1, dims);
1836 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1837 data2 = matrix_long_data(mat2, dims);
1838 }
1839 nmat = YAP_BlobOfTerm(tf);
1840 ndata = matrix_double_data(nmat, dims);
1841 switch (op) {
1842 case MAT_PLUS:
1843 matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
1844 break;
1845 case MAT_SUB:
1846 matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
1847 break;
1848 case MAT_TIMES:
1849 matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
1850 break;
1851 case MAT_DIV:
1852 matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
1853 break;
1854 case MAT_ZDIV:
1855 matrix_long_double_zdiv2_data(ndata, mat1[MAT_SIZE], data1, data2);
1856 break;
1857 default:
1858 return FALSE;
1859 }
1860 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1861 double *data2;
1862 double *ndata;
1863
1864 if (create)
1865 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1866 if (tf == YAP_TermNil()) {
1867 return FALSE;
1868 } else {
1869 /* there may have been an overflow */
1870 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1871 data1 = matrix_double_data(mat1, dims);
1872 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1873 data2 = matrix_double_data(mat2, dims);
1874 }
1875 nmat = YAP_BlobOfTerm(tf);
1876 ndata = matrix_double_data(nmat, dims);
1877 switch (op) {
1878 case MAT_PLUS:
1879 matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
1880 break;
1881 case MAT_SUB:
1882 matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
1883 break;
1884 case MAT_TIMES:
1885 matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
1886 break;
1887 case MAT_DIV:
1888 matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
1889 break;
1890 case MAT_ZDIV:
1891 matrix_double_zdiv_data(ndata, mat1[MAT_SIZE], data1, data2);
1892 break;
1893 default:
1894 return FALSE;
1895 }
1896 } else {
1897 return FALSE;
1898 }
1899 }
1900 return YAP_Unify(YAP_ARG4,tf);
1901 }
1902
1903 static void
add_int_by_cols(int total,int nlines,long int * mat1,long int * mat2,long int * ndata)1904 add_int_by_cols(int total,int nlines,long int *mat1,long int *mat2,long int *ndata)
1905 {
1906 int i, ncols = total/nlines;
1907 for (i=0;i<total;i++) {
1908 ndata[i] = mat1[i] + mat2[i/ncols];
1909 }
1910 }
1911
1912 static void
add_int_by_dcols(int total,int nlines,long int * mat1,double * mat2,double * ndata)1913 add_int_by_dcols(int total,int nlines,long int *mat1,double *mat2,double *ndata)
1914 {
1915 int i, ncols = total/nlines;
1916 for (i=0;i<total;i++) {
1917 ndata[i] = mat1[i] + mat2[i/ncols];
1918 }
1919 }
1920
1921 static void
add_double_by_cols(int total,int nlines,double * mat1,double * mat2,double * ndata)1922 add_double_by_cols(int total,int nlines,double *mat1,double *mat2,double *ndata)
1923 {
1924 int i;
1925 int ncols = total/nlines;
1926
1927 for (i=0;i<total;i++) {
1928 ndata[i] = mat1[i] + mat2[i/ncols];
1929 }
1930 }
1931
1932 static int
matrix_op_to_cols(void)1933 matrix_op_to_cols(void)
1934 {
1935 int *mat1, *mat2;
1936 YAP_Term top = YAP_ARG3;
1937 op_type op;
1938 YAP_Term tf;
1939
1940 if (!YAP_IsIntTerm(top)) {
1941 return FALSE;
1942 }
1943 op = YAP_IntOfTerm(top);
1944 mat1 = (int *)YAP_BlobOfTerm(YAP_ARG1);
1945 if (!mat1) {
1946 /* Error */
1947 return FALSE;
1948 }
1949 mat2 = (int *)YAP_BlobOfTerm(YAP_ARG2);
1950 if (!mat2) {
1951 /* Error */
1952 return FALSE;
1953 }
1954 if (mat1[MAT_TYPE] == INT_MATRIX) {
1955 long int *data1;
1956 int dims = mat1[MAT_NDIMS];
1957 int *nmat;
1958 data1 = matrix_long_data(mat1, dims);
1959
1960 if (mat2[MAT_TYPE] == INT_MATRIX) {
1961 long int *data2 = matrix_long_data(mat2, 1);
1962 if (op == MAT_PLUS) {
1963 long int *ndata;
1964
1965 tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
1966 if (tf == YAP_TermNil())
1967 return FALSE;
1968 nmat = YAP_BlobOfTerm(tf);
1969 ndata = matrix_long_data(nmat, dims);
1970 add_int_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1971 } else {
1972 return FALSE;
1973 }
1974 } else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
1975 double *data2 = matrix_double_data(mat2, 1);
1976 if (op == MAT_PLUS) {
1977 double *ndata;
1978
1979 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1980 if (tf == YAP_TermNil())
1981 return FALSE;
1982 nmat = YAP_BlobOfTerm(tf);
1983 ndata = matrix_double_data(nmat, dims);
1984 add_int_by_dcols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
1985 } else {
1986 return FALSE;
1987 }
1988 } else {
1989 return FALSE;
1990 }
1991 } else {
1992 double *data1, *data2, *ndata;
1993 int dims = mat1[MAT_NDIMS];
1994 int *nmat;
1995
1996 if (mat2[MAT_TYPE] != FLOAT_MATRIX)
1997 return FALSE;
1998 tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
1999 if (tf == YAP_TermNil())
2000 return FALSE;
2001 nmat = YAP_BlobOfTerm(tf);
2002 data1 = matrix_double_data(mat1, dims);
2003 data2 = matrix_double_data(mat2, 1);
2004 ndata = matrix_double_data(nmat, dims);
2005 if (op == MAT_PLUS) {
2006 add_double_by_cols(mat1[MAT_SIZE],mat1[MAT_DIMS],data1,data2,ndata);
2007 } else
2008 return FALSE;
2009 }
2010 return YAP_Unify(YAP_ARG4,tf);
2011 }
2012
2013 static int
matrix_op_to_all(void)2014 matrix_op_to_all(void)
2015 {
2016 int *mat;
2017 YAP_Term tf = 0;
2018 YAP_Term top = YAP_ARG2;
2019 op_type op;
2020 int create = FALSE;
2021
2022 if (!YAP_IsIntTerm(top)) {
2023 return FALSE;
2024 }
2025 op = YAP_IntOfTerm(top);
2026 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2027 if (!mat) {
2028 /* Error */
2029 return FALSE;
2030 }
2031 if (YAP_IsVarTerm(YAP_ARG4)) {
2032 create = TRUE;
2033 }
2034 /* create a new array with same dimensions */
2035 if (mat[MAT_TYPE] == INT_MATRIX) {
2036 long int *data;
2037 int dims = mat[MAT_NDIMS];
2038 int *nmat;
2039 YAP_Term tnum = YAP_ARG3;
2040
2041 if (YAP_IsIntTerm(tnum)) {
2042 long int num;
2043 long int *ndata;
2044
2045 num = YAP_IntOfTerm(tnum);
2046 data = matrix_long_data(mat, dims);
2047 if (create) {
2048 tf = new_int_matrix(dims,mat+(MAT_DIMS),NULL);
2049 if (tf == YAP_TermNil())
2050 return FALSE;
2051 nmat = (int *)YAP_BlobOfTerm(tf);
2052 ndata = matrix_long_data(nmat, dims);
2053 } else {
2054 nmat = mat;
2055 ndata = data;
2056 }
2057 if (op == MAT_PLUS) {
2058 int i;
2059
2060 for (i = 0; i < mat[MAT_SIZE]; i++) {
2061 ndata[i] = data[i] + num;
2062 }
2063 } else if (op == MAT_TIMES) {
2064 int i;
2065
2066 for (i = 0; i < mat[MAT_SIZE]; i++) {
2067 ndata[i] = data[i] * num;
2068 }
2069 } else {
2070 return FALSE;
2071 }
2072 } else if (YAP_IsFloatTerm(tnum)) {
2073 double num;
2074 double *ndata;
2075
2076
2077 num = YAP_FloatOfTerm(tnum);
2078 if (create) {
2079 tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
2080 if (tf == YAP_TermNil())
2081 return FALSE;
2082 nmat = (int *)YAP_BlobOfTerm(tf);
2083 ndata = matrix_double_data(nmat, dims);
2084 } else {
2085 return FALSE;
2086 }
2087 data = matrix_long_data(mat, dims);
2088 if (op == MAT_PLUS) {
2089 int i;
2090
2091 for (i = 0; i < mat[MAT_SIZE]; i++) {
2092 ndata[i] = data[i] + num;
2093 }
2094 } else if (op == MAT_TIMES) {
2095 int i;
2096
2097 for (i = 0; i < mat[MAT_SIZE]; i++) {
2098 ndata[i] = data[i] * num;
2099 }
2100 } else if (op == MAT_DIV) {
2101 int i;
2102
2103 for (i = 0; i < mat[MAT_SIZE]; i++) {
2104 ndata[i] = data[i] / num;
2105 }
2106 }
2107 } else {
2108 return FALSE;
2109 }
2110 } else {
2111 double *data, *ndata;
2112 int dims = mat[MAT_NDIMS];
2113 int *nmat;
2114 YAP_Term tnum = YAP_ARG3;
2115 double num;
2116
2117 if (YAP_IsFloatTerm(tnum)) {
2118 num = YAP_FloatOfTerm(tnum);
2119 } else if (!YAP_IntOfTerm(tnum)) {
2120 return FALSE;
2121 } else {
2122 if (!create)
2123 return FALSE;
2124 num = (double)YAP_IntOfTerm(tnum);
2125 }
2126 data = matrix_double_data(mat, dims);
2127 if (create) {
2128 tf = new_float_matrix(dims,mat+(MAT_DIMS),NULL);
2129 if (tf == YAP_TermNil())
2130 return FALSE;
2131 nmat = (int *)YAP_BlobOfTerm(tf);
2132 ndata = matrix_double_data(nmat, dims);
2133 } else {
2134 nmat = mat;
2135 ndata = data;
2136 }
2137 switch(op) {
2138 case MAT_PLUS:
2139 {
2140 int i;
2141
2142 for (i = 0; i < mat[MAT_SIZE]; i++) {
2143 ndata[i] = data[i] + num;
2144 }
2145 }
2146 break;
2147 case MAT_TIMES:
2148 {
2149 int i;
2150
2151 for (i = 0; i < mat[MAT_SIZE]; i++) {
2152 ndata[i] = data[i] * num;
2153 }
2154 }
2155 break;
2156 case MAT_DIV:
2157 {
2158 int i;
2159
2160 for (i = 0; i < mat[MAT_SIZE]; i++) {
2161 ndata[i] = data[i] / num;
2162 }
2163 }
2164 break;
2165 default:
2166 return FALSE;
2167 }
2168 }
2169 if (create)
2170 return YAP_Unify(YAP_ARG4,tf);
2171 return YAP_Unify(YAP_ARG4,YAP_ARG1);
2172 }
2173
2174 /* given a matrix M and a set of dims, build a new reordered matrix to follow
2175 the new order
2176 */
2177 static int
matrix_transpose(void)2178 matrix_transpose(void)
2179 {
2180 int ndims, i, *dims, *dimsn;
2181 int conv[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
2182 YAP_Term tconv, tf;
2183 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2184 if (!mat) {
2185 /* Error */
2186 return FALSE;
2187 }
2188 ndims = mat[MAT_NDIMS];
2189 if (mat[MAT_TYPE] == INT_MATRIX) {
2190 /* create a new matrix with the same size */
2191 tf = new_int_matrix(ndims,mat+MAT_DIMS,NULL);
2192 if (tf == YAP_TermNil())
2193 return FALSE;
2194 } else {
2195 /* create a new matrix with the same size */
2196 tf = new_float_matrix(ndims,mat+MAT_DIMS,NULL);
2197 if (tf == YAP_TermNil())
2198 return FALSE;
2199 }
2200 /* just in case there was an overflow */
2201 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2202 nmat = (int *)YAP_BlobOfTerm(tf);
2203 dims = mat+MAT_DIMS;
2204 dimsn = nmat+MAT_DIMS;
2205 /* we now have our target matrix, let us grab our conversion matrix */
2206 tconv = YAP_ARG2;
2207 for (i=0; i < ndims; i++) {
2208 YAP_Term th;
2209 long int j;
2210
2211 if (!YAP_IsPairTerm(tconv))
2212 return FALSE;
2213 th = YAP_HeadOfTerm(tconv);
2214 if (!YAP_IsIntTerm(th))
2215 return FALSE;
2216 conv[i] = j = YAP_IntOfTerm(th);
2217 dimsn[i] = dims[j];
2218 tconv = YAP_TailOfTerm(tconv);
2219 }
2220 /*
2221 we now got all the dimensions set up, so what we need to do
2222 next is to copy the elements to the new matrix.
2223 */
2224 if (mat[MAT_TYPE] == INT_MATRIX) {
2225 long int *data = matrix_long_data(mat,ndims);
2226 /* create a new matrix with the same size */
2227 for (i=0; i< mat[MAT_SIZE]; i++) {
2228 long int x = data[i];
2229 int j;
2230 /*
2231 not very efficient, we could try to take advantage of the fact
2232 that we usually only change an index at a time
2233 */
2234 matrix_get_index(mat, i, indx);
2235 for (j = 0; j < ndims; j++) {
2236 nindx[j] = indx[conv[j]];
2237 }
2238 matrix_long_set(nmat, nindx, x);
2239 }
2240 } else {
2241 double *data = matrix_double_data(mat,ndims);
2242 /* create a new matrix with the same size */
2243 for (i=0; i< mat[MAT_SIZE]; i++) {
2244 double x = data[i];
2245 long j;
2246
2247 matrix_get_index(mat, i, indx);
2248 for (j = 0; j < ndims; j++)
2249 nindx[j] = indx[conv[j]];
2250 matrix_float_set(nmat, nindx, x);
2251 }
2252 }
2253 return YAP_Unify(YAP_ARG3, tf);
2254 }
2255
2256 /* given a matrix M and a set of dims, fold one of the dimensions of the
2257 matrix on one of the elements
2258 */
2259 static int
matrix_select(void)2260 matrix_select(void)
2261 {
2262 int ndims, i, j, *dims, newdims, prdim, leftarg;
2263 int indx[MAX_DIMS], nindx[MAX_DIMS];
2264 YAP_Term tpdim, tdimarg, tf;
2265 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2266 if (!mat) {
2267 /* Error */
2268 return FALSE;
2269 }
2270 /* we now have our target matrix, let us grab our conversion arguments */
2271 tpdim = YAP_ARG2;
2272 ndims = mat[MAT_NDIMS];
2273 dims = mat+MAT_DIMS;
2274 if (!YAP_IsIntTerm(tpdim)) {
2275 return FALSE;
2276 }
2277 prdim = YAP_IntOfTerm(tpdim);
2278 tdimarg = YAP_ARG3;
2279 if (!YAP_IsIntTerm(tdimarg)) {
2280 return FALSE;
2281 }
2282 leftarg = YAP_IntOfTerm(tdimarg);
2283 for (i=0, j=0; i< ndims; i++) {
2284 if (i != prdim) {
2285 nindx[j]= (mat+MAT_DIMS)[i];
2286 j++;
2287 }
2288 }
2289 newdims = ndims-1;
2290 if (mat[MAT_TYPE] == INT_MATRIX) {
2291 long int *data, *ndata;
2292
2293 /* create a new matrix with the same size */
2294 tf = new_int_matrix(newdims,nindx,NULL);
2295 if (tf == YAP_TermNil())
2296 return FALSE;
2297 /* in case the matrix moved */
2298 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2299 nmat = (int *)YAP_BlobOfTerm(tf);
2300 data = matrix_long_data(mat,ndims);
2301 ndata = matrix_long_data(nmat,newdims);
2302 /* create a new matrix with smaller size */
2303 for (i=0; i< nmat[MAT_SIZE]; i++) {
2304 int j,k;
2305 /*
2306 not very efficient, we could try to take advantage of the fact
2307 that we usually only change an index at a time
2308 */
2309 matrix_get_index(nmat, i, indx);
2310 for (j = 0, k=0; j < newdims; j++,k++) {
2311 if (j == prdim) {
2312 nindx[k] = leftarg;
2313 k++;
2314 }
2315 nindx[k]= indx[j];
2316 }
2317 if (k == prdim) {
2318 nindx[k] = leftarg;
2319 }
2320 ndata[i] = data[matrix_get_offset(mat, nindx)];
2321 }
2322 } else {
2323 double *data, *ndata;
2324
2325 /* create a new matrix with the same size */
2326 tf = new_float_matrix(newdims,nindx,NULL);
2327 if (tf == YAP_TermNil())
2328 return FALSE;
2329 /* in case the matrix moved */
2330 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2331 nmat = (int *)YAP_BlobOfTerm(tf);
2332 data = matrix_double_data(mat,ndims);
2333 ndata = matrix_double_data(nmat,newdims);
2334 /* create a new matrix with the same size */
2335 for (i=0; i< nmat[MAT_SIZE]; i++) {
2336 int j,k;
2337 /*
2338 not very efficient, we could try to take advantage of the fact
2339 that we usually only change an index at a time
2340 */
2341 matrix_get_index(nmat, i, indx);
2342 for (j = 0, k=0; j < newdims; j++,k++) {
2343 if (j == prdim) {
2344 nindx[k] = leftarg;
2345 k++;
2346 }
2347 nindx[k]= indx[j];
2348 }
2349 if (k == prdim) {
2350 nindx[k] = leftarg;
2351 }
2352 ndata[i] = data[matrix_get_offset(mat, nindx)];
2353 }
2354 }
2355 return YAP_Unify(YAP_ARG4, tf);
2356 }
2357
2358 /* given a matrix M and a set of N-1 dims, get the first dimension
2359 */
2360 static int
matrix_column(void)2361 matrix_column(void)
2362 {
2363 int size, i, ndims, newdims[1];
2364 int indx[MAX_DIMS];
2365 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2366 YAP_Term tconv, tf;
2367
2368 if (!mat) {
2369 /* Error */
2370 return FALSE;
2371 }
2372 ndims = mat[MAT_NDIMS];
2373 /* we now have our target matrix, let us grab our conversion arguments */
2374 tconv = YAP_ARG2;
2375 for (i=1; i < ndims; i++) {
2376 YAP_Term th;
2377
2378 if (!YAP_IsPairTerm(tconv))
2379 return FALSE;
2380 th = YAP_HeadOfTerm(tconv);
2381 if (!YAP_IsIntTerm(th))
2382 return FALSE;
2383 indx[i] = YAP_IntOfTerm(th);
2384 tconv = YAP_TailOfTerm(tconv);
2385 }
2386 if (tconv != YAP_TermNil())
2387 return FALSE;
2388 newdims[0] = size = mat[MAT_DIMS];
2389 if (mat[MAT_TYPE] == INT_MATRIX) {
2390 long int *data, *ndata;
2391
2392 /* create a new matrix with the same size */
2393 tf = new_int_matrix(1, newdims, NULL);
2394 if (tf == YAP_TermNil())
2395 return FALSE;
2396 /* in case the matrix moved */
2397 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2398 nmat = (int *)YAP_BlobOfTerm(tf);
2399 data = matrix_long_data(mat,ndims);
2400 ndata = matrix_long_data(nmat,1);
2401 /* create a new matrix with smaller size */
2402 for (i=0; i< size; i++) {
2403 indx[0]=i;
2404 ndata[i] = data[matrix_get_offset(mat, indx)];
2405 }
2406 } else {
2407 double *data, *ndata;
2408
2409 /* create a new matrix with the same size */
2410 tf = new_float_matrix(1,newdims,NULL);
2411 if (tf == YAP_TermNil())
2412 return FALSE;
2413 /* in case the matrix moved */
2414 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2415 nmat = (int *)YAP_BlobOfTerm(tf);
2416 data = matrix_double_data(mat,ndims);
2417 ndata = matrix_double_data(nmat,1);
2418 /* create a new matrix with smaller size */
2419 for (i=0; i< size; i++) {
2420 indx[0]=i;
2421 ndata[i] = data[matrix_get_offset(mat, indx)];
2422 }
2423 }
2424 return YAP_Unify(YAP_ARG3, tf);
2425 }
2426
2427 /* given a matrix M and a set of dims, sum out one of the dimensions
2428 */
2429 static int
matrix_sum_out(void)2430 matrix_sum_out(void)
2431 {
2432 int ndims, i, j, *dims, newdims, prdim;
2433 int indx[MAX_DIMS], nindx[MAX_DIMS];
2434 YAP_Term tpdim, tf;
2435 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2436 if (!mat) {
2437 /* Error */
2438 return FALSE;
2439 }
2440 /* we now have our target matrix, let us grab our conversion arguments */
2441 tpdim = YAP_ARG2;
2442 ndims = mat[MAT_NDIMS];
2443 dims = mat+MAT_DIMS;
2444 if (!YAP_IsIntTerm(tpdim)) {
2445 return FALSE;
2446 }
2447 prdim = YAP_IntOfTerm(tpdim);
2448 newdims = ndims-1;
2449 for (i=0, j=0; i< ndims; i++) {
2450 if (i != prdim) {
2451 nindx[j]= (mat+MAT_DIMS)[i];
2452 j++;
2453 }
2454 }
2455 if (mat[MAT_TYPE] == INT_MATRIX) {
2456 long int *data, *ndata;
2457
2458 /* create a new matrix with the same size */
2459 tf = new_int_matrix(newdims,nindx,NULL);
2460 if (tf == YAP_TermNil())
2461 return FALSE;
2462 /* in case the matrix moved */
2463 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2464 nmat = (int *)YAP_BlobOfTerm(tf);
2465 data = matrix_long_data(mat,ndims);
2466 ndata = matrix_long_data(nmat,newdims);
2467 /* create a new matrix with smaller size */
2468 for (i=0;i<nmat[MAT_SIZE];i++)
2469 ndata[i] = 0;
2470 for (i=0; i< mat[MAT_SIZE]; i++) {
2471 int j, k;
2472 /*
2473 not very efficient, we could try to take advantage of the fact
2474 that we usually only change an index at a time
2475 */
2476 matrix_get_index(mat, i, indx);
2477 for (j = 0, k=0; j < ndims; j++) {
2478 if (j != prdim) {
2479 nindx[k++]= indx[j];
2480 }
2481 }
2482 ndata[matrix_get_offset(nmat, nindx)] += data[i];
2483 }
2484 } else {
2485 double *data, *ndata;
2486
2487 /* create a new matrix with the same size */
2488 tf = new_float_matrix(newdims,nindx,NULL);
2489 if (tf == YAP_TermNil())
2490 return FALSE;
2491 /* in case the matrix moved */
2492 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2493 nmat = (int *)YAP_BlobOfTerm(tf);
2494 data = matrix_double_data(mat,ndims);
2495 ndata = matrix_double_data(nmat,newdims);
2496 /* create a new matrix with smaller size */
2497 for (i=0;i<nmat[MAT_SIZE];i++)
2498 ndata[i] = 0.0;
2499 for (i=0; i< mat[MAT_SIZE]; i++) {
2500 int j, k;
2501 /*
2502 not very efficient, we could try to take advantage of the fact
2503 that we usually only change an index at a time
2504 */
2505 matrix_get_index(mat, i, indx);
2506 for (j = 0, k=0; j < ndims; j++) {
2507 if (j != prdim) {
2508 nindx[k++]= indx[j];
2509 }
2510 }
2511 ndata[matrix_get_offset(nmat, nindx)] += data[i];
2512 }
2513 }
2514 return YAP_Unify(YAP_ARG3, tf);
2515 }
2516
2517 /* given a matrix M and a set of dims, sum out one of the dimensions
2518 */
2519 static int
matrix_sum_out_several(void)2520 matrix_sum_out_several(void)
2521 {
2522 int ndims, i, *dims, newdims;
2523 int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
2524 YAP_Term tf, tconv;
2525 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2526 if (!mat) {
2527 /* Error */
2528 return FALSE;
2529 }
2530 ndims = mat[MAT_NDIMS];
2531 dims = mat+MAT_DIMS;
2532 /* we now have our target matrix, let us grab our conversion arguments */
2533 tconv = YAP_ARG2;
2534 for (i=0, newdims=0; i < ndims; i++) {
2535 YAP_Term th;
2536
2537 if (!YAP_IsPairTerm(tconv))
2538 return FALSE;
2539 th = YAP_HeadOfTerm(tconv);
2540 if (!YAP_IsIntTerm(th))
2541 return FALSE;
2542 conv[i] = YAP_IntOfTerm(th);
2543 if (!conv[i]) {
2544 nindx[newdims++] = dims[i];
2545 }
2546 tconv = YAP_TailOfTerm(tconv);
2547 }
2548 if (mat[MAT_TYPE] == INT_MATRIX) {
2549 long int *data, *ndata;
2550
2551 /* create a new matrix with the same size */
2552 tf = new_int_matrix(newdims,nindx,NULL);
2553 if (tf == YAP_TermNil())
2554 return FALSE;
2555 /* in case the matrix moved */
2556 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2557 nmat = (int *)YAP_BlobOfTerm(tf);
2558 data = matrix_long_data(mat,ndims);
2559 ndata = matrix_long_data(nmat,newdims);
2560 /* create a new matrix with smaller size */
2561 for (i=0;i<nmat[MAT_SIZE];i++)
2562 ndata[i] = 0;
2563 for (i=0; i< mat[MAT_SIZE]; i++) {
2564 int j, k;
2565 /*
2566 not very efficient, we could try to take advantage of the fact
2567 that we usually only change an index at a time
2568 */
2569 matrix_get_index(mat, i, indx);
2570 for (j = 0, k=0; j < ndims; j++) {
2571 if (!conv[j]) {
2572 nindx[k++]= indx[j];
2573 }
2574 }
2575 ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
2576 }
2577 } else {
2578 double *data, *ndata;
2579
2580 /* create a new matrix with the same size */
2581 tf = new_float_matrix(newdims,nindx,NULL);
2582 if (tf == YAP_TermNil())
2583 return FALSE;
2584 /* in case the matrix moved */
2585 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2586 nmat = (int *)YAP_BlobOfTerm(tf);
2587 data = matrix_double_data(mat,ndims);
2588 ndata = matrix_double_data(nmat,newdims);
2589 /* create a new matrix with smaller size */
2590 for (i=0;i<nmat[MAT_SIZE];i++)
2591 ndata[i] = 0.0;
2592 for (i=0; i< mat[MAT_SIZE]; i++) {
2593 int j, k;
2594 /*
2595 not very efficient, we could try to take advantage of the fact
2596 that we usually only change an index at a time
2597 */
2598 matrix_get_index(mat, i, indx);
2599 for (j = 0, k=0; j < ndims; j++) {
2600 if (!conv[j]) {
2601 nindx[k++]= indx[j];
2602 }
2603 }
2604 ndata[matrix_get_offset(nmat, nindx)] = log(exp(ndata[matrix_get_offset(nmat, nindx)]) + exp(data[i]));
2605 }
2606 }
2607 return YAP_Unify(YAP_ARG3, tf);
2608 }
2609
2610 /* given a matrix M and a set of dims, sum out one of the dimensions
2611 */
2612 static int
matrix_sum_out_logs(void)2613 matrix_sum_out_logs(void)
2614 {
2615 int ndims, i, j, *dims, newdims, prdim;
2616 int indx[MAX_DIMS], nindx[MAX_DIMS];
2617 YAP_Term tpdim, tf;
2618 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2619 if (!mat) {
2620 /* Error */
2621 return FALSE;
2622 }
2623 /* we now have our target matrix, let us grab our conversion arguments */
2624 tpdim = YAP_ARG2;
2625 ndims = mat[MAT_NDIMS];
2626 dims = mat+MAT_DIMS;
2627 if (!YAP_IsIntTerm(tpdim)) {
2628 return FALSE;
2629 }
2630 prdim = YAP_IntOfTerm(tpdim);
2631 newdims = ndims-1;
2632 for (i=0, j=0; i< ndims; i++) {
2633 if (i != prdim) {
2634 nindx[j]= (mat+MAT_DIMS)[i];
2635 j++;
2636 }
2637 }
2638 if (mat[MAT_TYPE] == INT_MATRIX) {
2639 long int *data, *ndata;
2640
2641 /* create a new matrix with the same size */
2642 tf = new_int_matrix(newdims,nindx,NULL);
2643 if (tf == YAP_TermNil())
2644 return FALSE;
2645 /* in case the matrix moved */
2646 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2647 nmat = (int *)YAP_BlobOfTerm(tf);
2648 data = matrix_long_data(mat,ndims);
2649 ndata = matrix_long_data(nmat,newdims);
2650 /* create a new matrix with smaller size */
2651 for (i=0;i<nmat[MAT_SIZE];i++)
2652 ndata[i] = 0;
2653 for (i=0; i< mat[MAT_SIZE]; i++) {
2654 int j, k;
2655 /*
2656 not very efficient, we could try to take advantage of the fact
2657 that we usually only change an index at a time
2658 */
2659 matrix_get_index(mat, i, indx);
2660 for (j = 0, k=0; j < ndims; j++) {
2661 if (j != prdim) {
2662 nindx[k++]= indx[j];
2663 }
2664 }
2665 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2666 }
2667 for (i=0; i< nmat[MAT_SIZE]; i++) {
2668 ndata[i] = log(ndata[i]);
2669 }
2670 } else {
2671 double *data, *ndata;
2672
2673 /* create a new matrix with the same size */
2674 tf = new_float_matrix(newdims,nindx,NULL);
2675 if (tf == YAP_TermNil())
2676 return FALSE;
2677 /* in case the matrix moved */
2678 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2679 nmat = (int *)YAP_BlobOfTerm(tf);
2680 data = matrix_double_data(mat,ndims);
2681 ndata = matrix_double_data(nmat,newdims);
2682 /* create a new matrix with smaller size */
2683 for (i=0;i<nmat[MAT_SIZE];i++)
2684 ndata[i] = 0.0;
2685 for (i=0; i< mat[MAT_SIZE]; i++) {
2686 int j, k;
2687 /*
2688 not very efficient, we could try to take advantage of the fact
2689 that we usually only change an index at a time
2690 */
2691 matrix_get_index(mat, i, indx);
2692 for (j = 0, k=0; j < ndims; j++) {
2693 if (j != prdim) {
2694 nindx[k++]= indx[j];
2695 }
2696 }
2697 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2698 }
2699 for (i=0; i< nmat[MAT_SIZE]; i++) {
2700 ndata[i] = log(ndata[i]);
2701 }
2702 }
2703 return YAP_Unify(YAP_ARG3, tf);
2704 }
2705
2706 /* given a matrix M and a set of dims, sum out one of the dimensions
2707 */
2708 static int
matrix_sum_out_logs_several(void)2709 matrix_sum_out_logs_several(void)
2710 {
2711 int ndims, i, *dims, newdims;
2712 int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
2713 YAP_Term tf, tconv;
2714 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2715 if (!mat) {
2716 /* Error */
2717 return FALSE;
2718 }
2719 ndims = mat[MAT_NDIMS];
2720 dims = mat+MAT_DIMS;
2721 /* we now have our target matrix, let us grab our conversion arguments */
2722 tconv = YAP_ARG2;
2723 for (i=0, newdims=0; i < ndims; i++) {
2724 YAP_Term th;
2725
2726 if (!YAP_IsPairTerm(tconv))
2727 return FALSE;
2728 th = YAP_HeadOfTerm(tconv);
2729 if (!YAP_IsIntTerm(th))
2730 return FALSE;
2731 conv[i] = YAP_IntOfTerm(th);
2732 if (!conv[i]) {
2733 nindx[newdims++] = dims[i];
2734 }
2735 tconv = YAP_TailOfTerm(tconv);
2736 }
2737 if (mat[MAT_TYPE] == INT_MATRIX) {
2738 long int *data, *ndata;
2739
2740 /* create a new matrix with the same size */
2741 tf = new_int_matrix(newdims,nindx,NULL);
2742 if (tf == YAP_TermNil())
2743 return FALSE;
2744 /* in case the matrix moved */
2745 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2746 nmat = (int *)YAP_BlobOfTerm(tf);
2747 data = matrix_long_data(mat,ndims);
2748 ndata = matrix_long_data(nmat,newdims);
2749 /* create a new matrix with smaller size */
2750 for (i=0;i<nmat[MAT_SIZE];i++)
2751 ndata[i] = 0;
2752 for (i=0; i< mat[MAT_SIZE]; i++) {
2753 int j, k;
2754 /*
2755 not very efficient, we could try to take advantage of the fact
2756 that we usually only change an index at a time
2757 */
2758 matrix_get_index(mat, i, indx);
2759 for (j = 0, k=0; j < ndims; j++) {
2760 if (!conv[j]) {
2761 nindx[k++]= indx[j];
2762 }
2763 }
2764 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2765 }
2766 } else {
2767 double *data, *ndata;
2768
2769 /* create a new matrix with the same size */
2770 tf = new_float_matrix(newdims,nindx,NULL);
2771 if (tf == YAP_TermNil())
2772 return FALSE;
2773 /* in case the matrix moved */
2774 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2775 nmat = (int *)YAP_BlobOfTerm(tf);
2776 data = matrix_double_data(mat,ndims);
2777 ndata = matrix_double_data(nmat,newdims);
2778 /* create a new matrix with smaller size */
2779 for (i=0;i<nmat[MAT_SIZE];i++)
2780 ndata[i] = 0.0;
2781 for (i=0; i< mat[MAT_SIZE]; i++) {
2782 int j, k;
2783 /*
2784 not very efficient, we could try to take advantage of the fact
2785 that we usually only change an index at a time
2786 */
2787 matrix_get_index(mat, i, indx);
2788 for (j = 0, k=0; j < ndims; j++) {
2789 if (!conv[j]) {
2790 nindx[k++]= indx[j];
2791 }
2792 }
2793 ndata[matrix_get_offset(nmat, nindx)] += exp(data[i]);
2794 }
2795 for (i=0; i< nmat[MAT_SIZE]; i++) {
2796 ndata[i] = log(ndata[i]);
2797 }
2798 }
2799 return YAP_Unify(YAP_ARG3, tf);
2800 }
2801
2802 /* given a matrix M and a set of dims, build contract a matrix to follow
2803 the new order
2804 */
2805 static int
matrix_expand(void)2806 matrix_expand(void)
2807 {
2808 int ndims, i, *dims, newdims=0, olddims = 0;
2809 int new[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
2810 YAP_Term tconv, tf;
2811 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2812 if (!mat) {
2813 /* Error */
2814 return FALSE;
2815 }
2816 /* we now have our target matrix, let us grab our conversion matrix */
2817 tconv = YAP_ARG2;
2818 ndims = mat[MAT_NDIMS];
2819 dims = mat+MAT_DIMS;
2820 for (i=0; i < MAX_DIMS; i++) {
2821 YAP_Term th;
2822 long int j;
2823
2824 if (!YAP_IsPairTerm(tconv)) {
2825 if (tconv != YAP_TermNil())
2826 return FALSE;
2827 break;
2828 }
2829 th = YAP_HeadOfTerm(tconv);
2830 if (!YAP_IsIntTerm(th))
2831 return FALSE;
2832 newdims++;
2833 j = YAP_IntOfTerm(th);
2834 if (j==0) {
2835 new[i] = 0;
2836 nindx[i] = dims[olddims];
2837 olddims++;
2838 } else {
2839 new[i] = 1;
2840 nindx[i] = j;
2841 }
2842 tconv = YAP_TailOfTerm(tconv);
2843 }
2844 if (olddims != ndims)
2845 return FALSE;
2846 if (mat[MAT_TYPE] == INT_MATRIX) {
2847 long int *data, *ndata;
2848
2849 /* create a new matrix with the same size */
2850 tf = new_int_matrix(newdims,nindx,NULL);
2851 if (tf == YAP_TermNil())
2852 return FALSE;
2853 /* in case the matrix moved */
2854 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2855 nmat = (int *)YAP_BlobOfTerm(tf);
2856 data = matrix_long_data(mat,ndims);
2857 ndata = matrix_long_data(nmat,newdims);
2858 /* create a new matrix with the same size */
2859 for (i=0; i< nmat[MAT_SIZE]; i++) {
2860 int j,k=0;
2861 /*
2862 not very efficient, we could try to take advantage of the fact
2863 that we usually only change an index at a time
2864 */
2865 matrix_get_index(nmat, i, indx);
2866 for (j = 0; j < newdims; j++) {
2867 if (!new[j])
2868 nindx[k++] = indx[j];
2869 }
2870 ndata[i] = data[matrix_get_offset(mat, nindx)];
2871 }
2872 } else {
2873 double *data, *ndata;
2874
2875 /* create a new matrix with the same size */
2876 tf = new_float_matrix(newdims,nindx,NULL);
2877 if (tf == YAP_TermNil())
2878 return FALSE;
2879 /* in case the matrix moved */
2880 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2881 nmat = (int *)YAP_BlobOfTerm(tf);
2882 data = matrix_double_data(mat,ndims);
2883 ndata = matrix_double_data(nmat,newdims);
2884 /* create a new matrix with the same size */
2885 for (i=0; i < newdims; i++)
2886 indx[i] = 0;
2887 for (i=0; i< nmat[MAT_SIZE]; i++) {
2888 int j,k=0;
2889 /*
2890 not very efficient, we could try to take advantage of the fact
2891 that we usually only change an index at a time
2892 */
2893 for (j = 0; j < newdims; j++) {
2894 if (!new[j])
2895 nindx[k++] = indx[j];
2896 }
2897 ndata[i] = data[matrix_get_offset(mat, nindx)];
2898 matrix_next_index(nmat+MAT_DIMS, newdims, indx);
2899 }
2900 }
2901 return YAP_Unify(YAP_ARG3, tf);
2902 }
2903
2904 /* given a matrix M and a set of dims, build contract a matrix to follow
2905 the new order
2906 */
2907 static int
matrix_set_all_that_disagree(void)2908 matrix_set_all_that_disagree(void)
2909 {
2910 int ndims, i, *dims;
2911 int indx[MAX_DIMS];
2912 YAP_Term tf;
2913 int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
2914 int dim = YAP_IntOfTerm(YAP_ARG2);
2915 int pos = YAP_IntOfTerm(YAP_ARG3);
2916
2917 if (!mat) {
2918 /* Error */
2919 return FALSE;
2920 }
2921 ndims = mat[MAT_NDIMS];
2922 dims = mat+MAT_DIMS;
2923 if (mat[MAT_TYPE] == INT_MATRIX) {
2924 long int *data, *ndata, val;
2925
2926 /* create a new matrix with the same size */
2927 tf = new_int_matrix(ndims,dims,NULL);
2928 if (tf == YAP_TermNil())
2929 return FALSE;
2930 /* in case the matrix moved */
2931 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2932 nmat = (int *)YAP_BlobOfTerm(tf);
2933 data = matrix_long_data(mat,ndims);
2934 ndata = matrix_long_data(nmat,ndims);
2935 if (!YAP_IsIntTerm(YAP_ARG4))
2936 return FALSE;
2937 val = YAP_IntOfTerm(YAP_ARG4);
2938 /* create a new matrix with the same size */
2939 for (i=0; i< nmat[MAT_SIZE]; i++) {
2940
2941 /*
2942 not very efficient, we could try to take advantage of the fact
2943 that we usually only change an index at a time
2944 */
2945 matrix_get_index(mat, i, indx);
2946 if (indx[dim] != pos)
2947 ndata[i] = val;
2948 else
2949 ndata[i] = data[i];
2950 }
2951 } else {
2952 double *data, *ndata, val;
2953
2954 /* create a new matrix with the same size */
2955 tf = new_float_matrix(ndims,dims,NULL);
2956 if (tf == YAP_TermNil())
2957 return FALSE;
2958 /* in case the matrix moved */
2959 mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
2960 nmat = (int *)YAP_BlobOfTerm(tf);
2961 data = matrix_double_data(mat,ndims);
2962 ndata = matrix_double_data(nmat,ndims);
2963 if (YAP_IsFloatTerm(YAP_ARG4))
2964 val = YAP_FloatOfTerm(YAP_ARG4);
2965 else if (YAP_IsIntTerm(YAP_ARG4))
2966 val = YAP_IntOfTerm(YAP_ARG4);
2967 else
2968 return FALSE;
2969 /* create a new matrix with the same size */
2970 for (i=0; i< nmat[MAT_SIZE]; i++) {
2971
2972 /*
2973 not very efficient, we could try to take advantage of the fact
2974 that we usually only change an index at a time
2975 */
2976 matrix_get_index(mat, i, indx);
2977 if (indx[dim] != pos)
2978 ndata[i] = val;
2979 else
2980 ndata[i] = data[i];
2981 }
2982 }
2983 return YAP_Unify(YAP_ARG5, tf);
2984 }
2985
2986 void PROTO(init_matrix, (void));
2987
2988 void
init_matrix(void)2989 init_matrix(void)
2990 {
2991 YAP_UserCPredicate("new_ints_matrix", new_ints_matrix, 4);
2992 YAP_UserCPredicate("new_ints_matrix_set", new_ints_matrix_set, 4);
2993 YAP_UserCPredicate("new_floats_matrix", new_floats_matrix, 4);
2994 YAP_UserCPredicate("new_floats_matrix_set", new_floats_matrix_set, 4);
2995 YAP_UserCPredicate("matrix_set", matrix_set, 3);
2996 YAP_UserCPredicate("matrix_set_all", matrix_set_all, 2);
2997 YAP_UserCPredicate("matrix_add", matrix_add, 3);
2998 YAP_UserCPredicate("matrix_get", do_matrix_access, 3);
2999 YAP_UserCPredicate("matrix_inc", do_matrix_inc, 2);
3000 YAP_UserCPredicate("matrix_dec", do_matrix_dec, 2);
3001 YAP_UserCPredicate("matrix_inc", do_matrix_inc2, 3);
3002 YAP_UserCPredicate("matrix_dec", do_matrix_dec2, 3);
3003 YAP_UserCPredicate("matrix_to_list", matrix_to_list, 2);
3004 YAP_UserCPredicate("matrix_dims", matrix_dims, 2);
3005 YAP_UserCPredicate("matrix_ndims", matrix_ndims, 2);
3006 YAP_UserCPredicate("matrix_size", matrix_size, 2);
3007 YAP_UserCPredicate("matrix_type_as_number", matrix_type, 2);
3008 YAP_UserCPredicate("matrix_arg_to_offset", matrix_arg_to_offset, 3);
3009 YAP_UserCPredicate("matrix_offset_to_arg", matrix_offset_to_arg, 3);
3010 YAP_UserCPredicate("matrix_max", matrix_max, 2);
3011 YAP_UserCPredicate("matrix_maxarg", matrix_maxarg, 2);
3012 YAP_UserCPredicate("matrix_min", matrix_min, 2);
3013 YAP_UserCPredicate("matrix_minarg", matrix_minarg, 2);
3014 YAP_UserCPredicate("matrix_sum", matrix_sum, 2);
3015 YAP_UserCPredicate("matrix_shuffle", matrix_transpose, 3);
3016 YAP_UserCPredicate("matrix_expand", matrix_expand, 3);
3017 YAP_UserCPredicate("matrix_select", matrix_select, 4);
3018 YAP_UserCPredicate("matrix_column", matrix_column, 3);
3019 YAP_UserCPredicate("matrix_to_logs", matrix_log_all,1);
3020 YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1);
3021 YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1);
3022 YAP_UserCPredicate("matrix_to_logs", matrix_log_all2,2);
3023 YAP_UserCPredicate("matrix_to_exps", matrix_exp_all2, 2);
3024 YAP_UserCPredicate("matrix_sum_out", matrix_sum_out, 3);
3025 YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3);
3026 YAP_UserCPredicate("matrix_sum_logs_out", matrix_sum_out_logs, 3);
3027 YAP_UserCPredicate("matrix_sum_logs_out_several", matrix_sum_out_logs_several, 3);
3028 YAP_UserCPredicate("matrix_set_all_that_disagree", matrix_set_all_that_disagree, 5);
3029 YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
3030 YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
3031 YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);
3032 YAP_UserCPredicate("do_matrix_op_to_all", matrix_op_to_all, 4);
3033 YAP_UserCPredicate("do_matrix_op_to_lines", matrix_op_to_lines, 4);
3034 YAP_UserCPredicate("do_matrix_op_to_cols", matrix_op_to_cols, 4);
3035 }
3036
3037 #ifdef _WIN32
3038
3039 int WINAPI PROTO(win_matrixs, (HANDLE, DWORD, LPVOID));
3040
win_matrixs(HANDLE hinst,DWORD reason,LPVOID reserved)3041 int WINAPI win_matrixs(HANDLE hinst, DWORD reason, LPVOID reserved)
3042 {
3043 switch (reason)
3044 {
3045 case DLL_PROCESS_ATTACH:
3046 break;
3047 case DLL_PROCESS_DETACH:
3048 break;
3049 case DLL_THREAD_ATTACH:
3050 break;
3051 case DLL_THREAD_DETACH:
3052 break;
3053 }
3054 return 1;
3055 }
3056 #endif
3057