diff --git a/docs/yap.tex b/docs/yap.tex index d9d4f73f4..41f96bf71 100644 --- a/docs/yap.tex +++ b/docs/yap.tex @@ -204,6 +204,7 @@ Subnodes of Library * DBUsage:: Information bout data base usage. * LineUtilities:: Line Manipulation Utilities * Lists:: List Manipulation +* MapArgs:: Apply on Arguments of Compound Terms. * MapList:: SWI-Compatible Apply library. * matrix:: Matrix Objects * MATLAB:: Matlab Interface @@ -8730,6 +8731,7 @@ Library, Extensions, Built-ins, Top * DBUsage:: Information bout data base usage. * Lists:: List Manipulation * LineUtilities:: Line Manipulation Utilities +* MapArgs:: Apply on Arguments of Compound Terms. * MapList:: SWI-Compatible Apply library. * matrix:: Matrix Objects * MATLAB:: Matlab Interface @@ -9460,7 +9462,7 @@ unification using @code{memberchk/2}. The complexity is See @code{ord_subtract/3}. @end table -@node LineUtilities, MapList, Lists, Library +@node LineUtilities, MapArgs, Lists, Library @section Line Manipulation Utilities @cindex Line Utilities Library @@ -9619,7 +9621,95 @@ Same as @code{file_filter/3}, but before starting the filter execute -@node MapList, matrix, LineUtilities, Library +@node MapArgs, MapList, LineUtilities, Library +@section Maplist +@cindex macros + +This library provides a set of utilities for applying a predicate to +all sub-terms of a term. They allow to +easily perform the most common do-loop constructs in Prolog. To avoid +performance degradation due to @code{apply/2}, each call creates an +equivalent Prolog program, without meta-calls, which is executed by +the Prolog engine instead. + +@table @code +@item mapargs(+@var{Pred}, +@var{TermIn}) +@findex mapargs/2 +@snindex mapargs/2 +@cnindex mapargs/2 + Applies the predicate @var{Pred} to all + arguments of @var{TermIn} + +@item mapargs(+@var{Pred}, +@var{TermIn}, ?@var{TermOut}) +@findex mapargs/3 +@snindex mapargs/3 +@cnindex mapargs/3 + Creates @var{TermOut} by applying the predicate @var{Pred} to all + arguments of @var{TermIn} + +@item mapargs(+@var{Pred}, +@var{TermIn}, ?@var{TermOut1}, ?@var{TermOut2}) +@findex mapargs/4 +@snindex mapargs/4 +@cnindex mapargs/4 + Creates @var{TermOut1} and @var{TermOut2} by applying the predicate @var{Pred} to all + arguments of @var{TermIn} + +@item mapargs(+@var{Pred}, +@var{TermIn}, ?@var{TermOut1}, + ?@var{TermOut2}, ?@var{TermOut3}) +@findex mapargs/5 +@snindex mapargs/5 +@cnindex mapargs/5 + Creates @var{TermOut1}, @var{TermOut2} and @var{TermOut3} by applying the predicate @var{Pred} to all + arguments of @var{TermIn} + +@item mapargs(+@var{Pred}, +@var{TermIn}, ?@var{TermOut1}, + ?@var{TermOut2}, ?@var{TermOut3}, ?@var{TermOut4}) +@findex mapargs/6 +@snindex mapargs/6 +@cnindex mapargs/6 + Creates @var{TermOut1}, @var{TermOut2}, @var{TermOut3} and @var{TermOut4} by applying the predicate @var{Pred} to all + arguments of @var{TermIn} + + +@item foldargs(+@var{Pred}, +@var{Term}, ?@var{AccIn}, ?@var{AccOut}) +@findex foldargs/4 +@snindex foldargs/4 +@cnindex foldargs/4 + Calls the predicate @var{Pred} on all arguments of @var{Term} and +collects a result in @var{Accumulator} + +@item foldargs(+@var{Pred}, +@var{Term}, +@var{Term1}, ?@var{AccIn}, ?@var{AccOut}) +@findex foldargs/5 +@snindex foldargs/5 +@cnindex foldargs/5 + Calls the predicate @var{Pred} on all arguments of @var{Term} and @var{Term1} and +collects a result in @var{Accumulator} + +@item foldargs(+@var{Pred}, +@var{Term}, +@var{Term1}, +@var{Term2}, ?@var{AccIn}, ?@var{AccOut}) +@findex foldargs/6 +@snindex foldargs/6 +@cnindex foldargs/6 + Calls the predicate @var{Pred} on all arguments of @var{Term}, +@var{Term1} and @var{Term2} and +collects a result in @var{Accumulator} + +@item foldargs(+@var{Pred}, +@var{Term}, +@var{Term1}, +@var{Term2}, +@var{Term3}, ?@var{AccIn}, ?@var{AccOut}) +@findex foldargs/7 +@snindex foldargs/7 +@cnindex foldargs/7 + Calls the predicate @var{Pred} on all arguments of @var{Term}, +@var{Term1}, +@var{Term2} and @var{Term3} and +collects a result in @var{Accumulator} + +@item sumargs(+@var{Pred}, +@var{Term}, ?@var{AccIn}, ?@var{AccOut}) +@findex sumargs/4 +@snindex sumargs/4 +@cnindex sumargs/4 + Calls the predicate @var{Pred} on all arguments of @var{Term} and +collects a result in @var{Accumulator} (uses reverse order than foldargs). + +@end table + + +@node MapList, matrix, MapArgs, Library @section Maplist @cindex macros @@ -9742,7 +9832,14 @@ applying the predicate @var{Pred} to all list elements on which @findex foldl2/7 @snindex foldl2/7 @cnindex foldl2/7 - Calls @var{Pred} on all elements of @code{List} and collects a result in + Calls @var{Pred} on all elements of @var{List} and @var{List1} and collects a result in +@var{X} and @var{Y}. + +@item foldl2(:@var{Pred}, +@var{List}, ?@var{List1}, ?@var{List2}, ?@var{X0}, ?@var{X}, ?@var{Y0}, ?@var{Y}) +@findex foldl2/8 +@snindex foldl2/8 +@cnindex foldl2/8 + Calls @var{Pred} on all elements of @var{List}, @var{List1} and @var{List2} and collects a result in @var{X} and @var{Y}. @item foldl3(:@var{Pred}, +@var{List1}, ?@var{List2}, ?@var{X0}, @@ -9796,20 +9893,6 @@ result in @var{X}, @var{Y}, @var{Z} and @var{W}. @cnindex scanl/7 Left scan of list. -@item mapargs(+@var{Pred}, ?@var{TermIn}, ?@var{TermOut}) -@findex mapargs/3 -@snindex mapargs/3 -@cnindex mapargs/3 - Creates @var{TermOut} by applying the predicate @var{Pred} to all - arguments of @var{TermIn} - -@item sumargs(+@var{Pred}, +@var{Term}, ?@var{AccIn}, ?@var{AccOut}) -@findex sumargs/4 -@snindex sumargs/4 -@cnindex sumargs/4 - Calls the predicate @var{Pred} on all arguments of @var{Term} and -collects a result in @var{Accumulator} - @item mapnodes(+@var{Pred}, +@var{TermIn}, ?@var{TermOut}) @findex mapnodes/3 @snindex mapnodes/3 @@ -9902,13 +9985,149 @@ This package provides a fast implementation of multi-dimensional matrices of integers and floats. In contrast to dynamic arrays, these matrices are multi-dimensional and compact. In contrast to static arrays. these arrays are allocated in the stack. Matrices are available -by loading the library @code{library(matrix)}. +by loading the library @code{library(matrix)}. They are multimensional +objects of type:@itemize + @item @t{terms}: Prolog terms +@item @t{ints}: bounded integers, represented as an opaque term. The +maximum integer depends on hardware, but should be obtained from the +natural size of the machine. +@item @t{floats}: floating-poiny numbers, represented as an opaque term. +@end itemize -Notice that the functionality in this library is only partial. Please +Matrix elements can be accessed through the @code{matrix_get/2} +predicate or through an @t{R}-inspired access notation (that uses the ciao +style extension to @code{[]}. Examples include: + +@table @code +@item @var{E} <== @var{X}[2,3] +Access the second row, third column of matrix @t{X}. Indices start from +@code{0}, +@item @var{L} <== @var{X}[2,_] +Access all the second row, the output is a list ofe elements. +@item @var{L} <== @var{X}[2..4,_] +Access all the second, thrd and fourth rows, the output is a list of elements. +@item @var{L} <== @var{X}[2..4+3,_] +Access all the fifth, sixth and eight rows, the output is a list of elements. +@end table + +The matrix library also supports a B-Prolog/ECliPSe inspired @code{for} operator to iterate over +elements of a matrix: + +@table @code +@item for(I in 0..N1, X[I] <== Y[I]) +Copies a vector, element by element. +@item for([I in 0..N1, J in I..N1], Z[I,J] <== X[I,J] - X[I,J]) +The lower-triangular matrix @var{Z} is the difference between the +lower-triangular and upper-triangular parts of @var{X}. +@item for([I in 0..N1, J in 0..N1], plus(X[I,J]), 0, Sum) +Add all elements of a matrix by using @var{Sum} as an accumulator. +@end table + +Notice that the library does not support all known matrix operations. Please contact the YAP maintainers if you require extra functionality. @table @code +@item @var{X} = array[@var{Dim1},...,@var{Dimn}] of @var{Objects} +@findex of/2 +@snindex of/2 +@cnindex of/2 +The @code{of/2} operator can be used to create a new array of +@var{Objects}. The objects supported are: + +@table @code +@item Unbound Variable +create an array of free variables +@item ints +create an array of integers +@item floats +create an array of floating-point numbers +@item @var{I}:@var{J} +create an array with integers from @var{I} to @var{J} +@item [..] +create an array from the values in a list + @end table + +The dimensions can be given as an integer, and the matrix will be +indexed @code{C}-style from @code{0..(@var{Max}-1)}, or can be given +as an interval @code{@var{Base}..@var{Limit}}. In the latter case, +matrices of integers and of floating-point numbers should have the same +@var{Base} on every dimension. + +@item ?@var{LHS} <== @var{RHS} +@findex <==/2 +@snindex <==/2 +@cnindex <==/2 +General matrix assignment operation. It evaluates the right-hand side +and then acts different according to the +left-hand side and to the matrix: +@itemize +@item if @var{LHS} is part of an integer or floating-point matrix, +perform non-backtrackable assignment. +@item other unify left-hand side and right-hand size. +@end itemize + +The right-hand side supports the following operators: +@table @code +@item []/2 +written as @var{M}[@var{Offset}]: obtain an element or list of elements +of matrix @var{M} at offset @var{Offset}. +@item matrix/1 +create a vector from a list +@item matrix/2 +create a matrix from a list. Oprions are: +@table @code +@item dim= +a list of dimensiona +@item type= +integers, floating-point or terms +@item base= +a list of base offsets per dimension (all must be the same for arrays of +integers and floating-points +@end table +@item matrix/3 +create matrix giving two options +@item dim/1 +list with matrix dimensions +@item nrow/1 +number of rows in bi-dimensional matrix +@item ncol/1 +number of columns in bi-dimensional matrix +@item length/1 +size of a matrix +@item size/1 +size of a matrix +@item max/1 +maximum element of a numeric matrix +@item maxarg/1 +argument of maximum element of a numeric matrix +@item min/1 +minimum element of a numeric matrix +@item minarg/1 +argument of minimum element of a numeric matrix +@item list/1 +represent matrix as a list +@item lists/2 +represent matrix as list of embedded lists +@item ../2 +@var{I}..@var{J} generates a list with all integers from @var{I} to +@var{J}, included. +@item +/2 +add two numbers, add two matrices element-by-element, or add a number to +all elements of a matrix or list +@item -/2 +subtract two numbers, subtract two matrices or lists element-by-element, or subtract a number from +all elements of a matrix or list +@item */2 +multiply two numbers, multiply two matrices or lists element-by-element, or multiply a number from +all elements of a matrix or list +@item log/1 +natural logarithm of a number, matrix or list +@item exp/1 +natural exponentiation of a number, matrix or list +@end table + + @item matrix_new(+@var{Type},+@var{Dims},-@var{Matrix}) @findex matrix_new/3 @snindex matrix_new/3 @@ -10143,7 +10362,7 @@ and @var{Matrix2}. Currently, only addition (@code{+}) is supported. @var{Result} is the result of applying @var{Op} to all elements of @var{Matrix1}, with @var{Operand} as the second argument. Currently, -only addition (@code{+}), multiplication (@code{+}), and division +only addition (@code{+}), multiplication (@code{*}), and division (@code{/}) are supported. @item matrix_op_to_lines(+@var{Matrix1},+@var{Lines},+@var{Op},-@var{Result}) diff --git a/library/Makefile.in b/library/Makefile.in index ee5b558b5..53433693a 100644 --- a/library/Makefile.in +++ b/library/Makefile.in @@ -50,7 +50,9 @@ PROGRAMS= \ $(srcdir)/lists.yap \ $(srcdir)/nb.yap \ $(srcdir)/ordsets.yap \ + $(srcdir)/mapargs.yap \ $(srcdir)/maplist.yap \ + $(srcdir)/maputils.yap \ $(srcdir)/matlab.yap \ $(srcdir)/matrix.yap \ $(srcdir)/prandom.yap \ diff --git a/library/matrix.yap b/library/matrix.yap index bbc10232e..78fc55ebe 100644 --- a/library/matrix.yap +++ b/library/matrix.yap @@ -20,7 +20,7 @@ may have a number of dimensions. These routines implement a number of routine manipulation procedures. - matrix(Type,D1,D2,...,Dn,data(......)) + '$matrix'(Type,D1,D2,...,Dn,data(......)) Type = int, float @@ -39,10 +39,11 @@ typedef enum { :- module( matrix, [op(100, yf, []), - (<==)/2, op(500, xfx, '<=='), + (<==)/2, op(600, xfx, '<=='), op(700, xfx, in), op(700, xfx, ins), op(450, xfx, ..), % should bind more tightly than \/ + op(710, xfx, of), of/2, matrix_new/3, matrix_new/4, matrix_new_set/4, @@ -99,11 +100,50 @@ typedef enum { :- load_foreign_files([matrix], [], init_matrix). +:- multifile rhs_opaque/1, array_extension/2. + :- meta_predicate for(+,0), for(+,2, +, -). :- use_module(library(maplist)). +:- use_module(library(mapargs)). :- use_module(library(lists)). +( X = '[]'(Dims0, array) of V ) :- + var(V), !, + foldl( norm_dim, Dims0, Dims, Bases, 1, Size ), + length( L, Size ), + X <== matrix( L, [dim=Dims,base=Bases] ). +( X = '[]'(Dims0, array) of ints ) :- !, + foldl( norm_dim, Dims0, Dims, Bases, 1, _Size ), + matrix_new( ints , Dims, X ), + matrix_base(X, Bases). +( X = '[]'(Dims0, array) of floats ) :- !, + foldl( norm_dim, Dims0, Dims, Bases, 1, _Size ), + matrix_new( floats , Dims, X ), + matrix_base(X, Bases). +( X = '[]'(Dims0, array) of (I:J) ) :- !, + foldl( norm_dim, Dims0, Dims, Bases, 1, Size ), + matrix_seq(I, J, Dims, X), + matrixn_size(X, Size), + matrix_base(X, Bases). +( X = '[]'(Dims0, array) of L ) :- + length( L, Size ), !, + foldl( norm_dim, Dims0, Dims, Bases, 1, Size ), + X <== matrix( L, [dim=Dims,base=Bases] ). +( X = '[]'(Dims0, array) of Pattern ) :- + array_extension(Pattern, Goal), + foldl( norm_dim, Dims0, Dims, Bases, 1, Size ), + call(Goal, Pattern, Dims, Size, L), + X <== matrix( L, [dim=Dims,base=Bases] ). + + +norm_dim( I..J, D, I, P0, P) :- !, + D is J+1-I, + P is P0*D. +norm_dim( I, I, 0, P0, P ) :- + P is P0*I. + + ( LHS <== RHS ) :- rhs(RHS, R), set_lhs( LHS, R). @@ -113,7 +153,7 @@ rhs(RHS, RHS) :- var(RHS), !. rhs(A, A) :- atom(A), !. rhs(RHS, RHS) :- number(RHS), !. rhs(RHS, RHS) :- opaque(RHS), !. -rhs(RHS, RHS) :- RHS = m(_, _, _, _), !. +rhs(RHS, RHS) :- RHS = '$matrix'(_, _, _, _, _), !. rhs(matrix(List), RHS) :- !, rhs( List, A1), new_matrix(A1, [], RHS). @@ -123,15 +163,6 @@ rhs(matrix(List, Opt1), RHS) :- !, rhs(matrix(List, Opt1, Opt2), RHS) :- !, rhs( List, A1), new_matrix(A1, [Opt1, Opt2], RHS). -rhs(matrix(List, Opt1, Opt2, Opt3), RHS) :- !, - rhs( List, A1), - new_matrix(A1, [Opt1, Opt2, Opt3], RHS). -rhs(matrix(List, Opt1, Opt2, Opt3, Opt4), RHS) :- !, - rhs( List, A1), - new_matrix(A1, [Opt1, Opt2, Opt3, Opt4], RHS). -rhs(matrix(List, Opt1, Opt2, Opt3, Opt4, Opt5), RHS) :- !, - rhs( List, A1), - new_matrix(A1, [Opt1, Opt2, Opt3, Opt4, Opt5], RHS). rhs(dim(RHS), Dims) :- !, rhs(RHS, X1), matrix_dims( X1, Dims ). @@ -156,19 +187,23 @@ rhs(max(RHS), Size) :- !, rhs(min(RHS), Size) :- !, rhs(RHS, X1), matrix_min( X1, Size ). +rhs(maxarg(RHS), Size) :- !, + rhs(RHS, X1), + matrix_maxarg( X1, Size ). +rhs(minarg(RHS), Size) :- !, + rhs(RHS, X1), + matrix_minarg( X1, Size ). rhs(list(RHS), List) :- !, rhs(RHS, X1), matrix_to_list( X1, List ). rhs(lists(RHS), List) :- !, rhs(RHS, X1), matrix_to_lists( X1, List ). -rhs(A=B, NA=NB) :- !, - rhs(A, NA), - rhs(B, NB). -rhs('[]'(Args, RHS), Val) :- !, +rhs('[]'(Args, RHS), Val) :- + !, rhs(RHS, X1), - matrix_dims( X1, Dims ), - maplist( index(Range), Args, Dims, NArgs), + matrix_dims( X1, Dims, Bases), + maplist( index(Range), Args, Dims, Bases, NArgs), ( var(Range) -> @@ -183,10 +218,23 @@ rhs('..'(I, J), [I1|Is]) :- !, rhs([H|T], [NH|NT]) :- !, rhs(H, NH), rhs(T, NT). -rhs(':'(I, J), [I1|Is]) :- !, - rhs(I, I1), - rhs(J, J1), - once( foldl(inc, Is, I1, J1) ). +rhs(log(RHS), Logs ) :- !, + rhs(RHS, X1), + matrix_to_logs( X1, Logs ). +rhs(exp(RHS), Logs ) :- !, + rhs(RHS, X1), + matrix_to_exps( X1, Logs ). +rhs(S, NS) :- + rhs_opaque( S ), !, + S = NS. +rhs(E1+E2, V) :- !, + rhs(E1, R1), + rhs(E2, R2), + mplus(R1, R2, V). +rhs(E1-E2, V) :- !, + rhs(E1, R1), + rhs(E2, R2), + msub(R1, R2, V). rhs(S, NS) :- S =.. [N|As], maplist(rhs, As, Bs), @@ -194,23 +242,32 @@ rhs(S, NS) :- set_lhs(V, R) :- var(V), !, V = R. set_lhs(V, R) :- number(V), !, V = R. -set_lhs(V, R) :- V = '[]'(Indx, M), !, - matrix_set( M, Indx, R). +set_lhs('[]'(Args, M), Val) :- + matrix_dims( M, Dims, Bases), + maplist( index(Range), Args, Dims, Bases, NArgs), + ( + var(Range) + -> + matrix_set( M, NArgs, Val ) + ; + matrix_set_range( M, NArgs, Val ) + ). % % ranges of arguments % -index(Range, V, M, Indx) :- var(V), !, - index(Range, 0..(M-1), M, Indx). -index(Range, '*', M, Indx) :- !, - index(Range, 0..(M-1), M, Indx). -index(Range, Exp, M, Indx) :- !, +index(Range, V, M, Base, Indx) :- var(V), !, + Max is (M-1)+Base, + index(Range, Base..Max, M, Base, Indx). +index(Range, '*', M, Base, Indx) :- !, + Max is (M-1)+Base, + index(Range, Base..Max, M, Base, Indx). +index(Range, Exp, M, _Base, Indx) :- !, index(Exp, M, Indx0), ( integer(Indx0) -> Indx = Indx0 ; Indx0 = [Indx] -> true ; Indx0 = Indx, Range = range ). - index(I, _M, I ) :- integer(I), !. index(I..J, _M, [I|O] ) :- !, I1 is I, J1 is J, @@ -218,23 +275,23 @@ index(I..J, _M, [I|O] ) :- !, index(I:J, _M, [I|O] ) :- !, I1 is I, J1 is J, once( foldl(inc, O, I1, J1) ). -index(I+J, _M, O ) :- +index(I+J, _M, O ) :- !, index(I, M, I1), index(J, M, J1), add_index(I1, J1, O). -index(I-J, _M, O ) :- +index(I-J, _M, O ) :- !, index(I, M, I1), index(J, M, J1), - add_index(I1, J1, O). -index(I*J, _M, O ) :- + sub_index(I1, J1, O). +index(I*J, _M, O ) :- !, index(I, M, I1), index(J, M, J1), O is I1*J1. -index(I div J, _M, O ) :- +index(I div J, _M, O ) :- !, index(I, M, I1), index(J, M, J1), O is I1 div J1. -index(I rem J, _M, O ) :- +index(I rem J, _M, O ) :- !, index(I, M, I1), index(J, M, J1), O is I1 rem J1. @@ -273,13 +330,73 @@ minus(X, Y, Z) :- Z is X-Y. rminus(X, Y, Z) :- Z is Y-X. +times(X, Y, Z) :- Z is Y*X. + +div(X, Y, Z) :- Z is X/Y. + +rdiv(X, Y, Z) :- Z is Y/X. + +zdiv(X, Y, Z) :- (X == 0 -> Z = 0 ; X == 0.0 -> Z = 0.0 ; Z is X / Y ). + +mplus(I1, I2, V) :- + number(I1) -> + ( number(I2) -> V is I1+I2 ; + '$matrix'(I2) -> matrix_op_to_all(I1, +, I2, V) ; + is_list(I2) -> maplist(plus(I1), I2, V) ; + V = I1+I2 ) ; + matrix(I1) -> + ( number(I2) -> matrix_op_to_all(I1, +, I2, V) ; + '$matrix'(I2) -> matrix_op(I1, I2, +, V) ; + V = I1+I2 ) ; + is_list(I1) -> + ( number(I2) -> maplist(plus(I2), I1, V) ; + is_list(I2) -> maplist(plus, I1, I2, V) ; + V = I1+I2 ) ; + V = I1 +I2. + +msub(I1, I2, V) :- + number(I1) -> + ( number(I2) -> V is I1-I2 ; + matrix(I2) -> matrix_op_to_all(I1, -, NI2, V) ; + is_list(I2) -> maplist(minus(I1), I2, V) ; + V = I1-I2 ) ; + matrix(I1) -> + ( number(I2) -> NI2 is -I2, matrix_op_to_all(I1, +, NI2, V) ; + matrix(I2) -> matrix_op(I1, I2, -, V) ; + V = I1-I2 ) ; + is_list(I1) -> + ( number(I2) -> NI2 is -I2, maplist(plus(NI2), I1, V) ; + is_list(I2) -> maplist(minus, I1, I2, V) ; + V = I1-I2 ) ; + V = I1-I2. + + +mtimes(I1, I2, V) :- + number(I1) -> + ( number(I2) -> V is I1*I2 ; + matrix(I2) -> matrix_op_to_all(I1, *, I2, V) ; + is_list(I2) -> maplist(times(I1), I2, V) ; + V = I1*I2 ) ; + matrix(I1) -> + ( number(I2) -> matrix_op_to_all(I1, *, I2, V) ; + matrix(I2) -> matrix_op(I1, I2, *, V) ; + V = I1*I2 ) ; + is_list(I1) -> + ( number(I2) -> maplist(times(I2), I1, V) ; + is_list(I2) -> maplist(times, I1, I2, V) ; + V = I1*I2 ) ; + V = I1 *I2. + + + % % three types of matrix: integers, floats and general terms. % -matrix_new(terms,Dims, m(Dims, NDims, Size, Matrix) ) :- +matrix_new(terms,Dims, '$matrix'(Dims, NDims, Size, Offsets, Matrix) ) :- length(Dims,NDims), foldl(size, Dims, 1, Size), + maplist(zero, Dims, Offsets), functor( Matrix, c, Size). matrix_new(ints,Dims,Matrix) :- length(Dims,NDims), @@ -289,9 +406,10 @@ matrix_new(floats,Dims,Matrix) :- new_floats_matrix_set(NDims, Dims, 0.0, Matrix). -matrix_new(terms, Dims, Data, m(Dims, NDims, Size, Matrix) ) :- +matrix_new(terms, Dims, Data, '$matrix'(Dims, NDims, Size, Offsets, Matrix) ) :- length(Dims,NDims), foldl(size, Dims, 1, Size), + maplist(zero, Dims, Offsets), functor( Matrix, c, Size), Matrix =.. [c|Data]. matrix_new(ints,Dims,Data,Matrix) :- @@ -304,19 +422,23 @@ matrix_new(floats,Dims,Data,Matrix) :- matrix_dims( Mat, Dims) :- ( opaque(Mat) -> matrixn_dims( Mat, Dims ) ; - Mat = m( Dims, _, _, _) ). + Mat = '$matrix'( Dims, _, _, _, _) ). + +matrix_dims( Mat, Dims, Bases) :- + ( opaque(Mat) -> matrixn_dims( Mat, Dims, Bases ) ; + Mat = '$matrix'( Dims, _, _, Bases, _) ). matrix_ndims( Mat, NDims) :- ( opaque(Mat) -> matrixn_ndims( Mat, NDims ) ; - Mat = m( _, NDims, _, _) ). + Mat = '$matrix'( _, NDims, _, _, _) ). matrix_size( Mat, Size) :- ( opaque(Mat) -> matrixn_size( Mat, Size ) ; - Mat = m( _, _, Size, _) ). + Mat = '$matrix'( _, _, Size, _, _) ). matrix_to_list( Mat, ToList) :- ( opaque(Mat) -> matrixn_to_list( Mat, ToList ) ; - Mat = m( _, _, _, M), M=.. [_|ToList] ). + Mat = '$matrix'( _, _, _, _, M), M=.. [_|ToList] ). matrix_to_lists( Mat, ToList) :- matrix_dims( Mat, [D|Dims] ), @@ -350,6 +472,10 @@ add_index_prefix( [L|Els0] , H ) --> [[H|L]], add_index_prefix( Els0 , H ). +matrix_set_range( Mat, Pos, Els) :- + slice(Pos, Keys), + maplist( matrix_set(Mat), Keys, Els). + matrix_set( Mat, Pos, El) :- ( opaque(Mat) -> matrixn_set( Mat, Pos, El ) ; m_set(Mat, Pos, El) ). @@ -367,30 +493,70 @@ matrix_type(Matrix,Type) :- opaque( Matrix ) -> Type = floats ; Type = terms ). +matrix_base(Matrix, Bases) :- + ( opaque( Matrix ) -> maplist('='(Base), Bases), matrixn_set_base( Matrix, Base ) ; + nb_setarg(4, Matrix, Bases ) ). + matrix_arg_to_offset(M, Index, Offset) :- ( opaque(M) -> matrixn_arg_to_offset( M, Index, Offset ) ; - M = m(Dims, _, Size, _) -> foldl2(indx, Index, Dims, Size, _, 0, Offset) ). + M = '$matrix'(Dims, _, Size, Bases, _) -> foldl2(indx, Index, Dims, Bases, Size, _, 0, Offset) ). matrix_offset_to_arg(M, Offset, Index) :- ( opaque(M) -> matrixn_offset_to_arg( M, Offset, Index ) ; - M = m(Dims, _, Size, _) -> foldl2(offset, Index, Dims, Size, _, Offset, _) ). + M = '$matrix'(Dims, _, Size, Bases, _) -> foldl2(offset, Index, Dims, Bases, Size, _, Offset, _) ). matrix_max(M, Max) :- ( opaque(M) -> matrixn_max( M, Max ) ; - M = m(_, _, _, M) -> fail ). + M = '$matrix'(_, _, _, _, C) -> + arg(1,C,V0), foldargs(max, M, V0, Max) ; + M = [V0|L], foldl(max, L, V0, Max) ). + +max(New, Old, Max) :- ( New >= Old -> New = Max ; Old = Max ). -matrix_maxarg(M, Max) :- - ( opaque(M) -> matrixn_maxarg( M, Max ) ; - M = m(_, _, _, _) -> fail ). +matrix_maxarg(M, MaxArg) :- + ( opaque(M) -> matrixn_maxarg( M, MaxArg ); + M = '$matrix'(_, _, _, _, C) -> + arg(1,C,V0), foldargs(maxarg, M, V0-0-0, _-Offset-_), matrix_offset_to_arg(M, Offset, MaxArg) ; + M = [V0|L], foldl(maxarg, L, V0-0-1, _Max-Off-_ ), MaxArg = [Off] ). + +maxarg(New, Old-OPos-I0, Max-MPos-I) :- I is I0+1, ( New > Old -> New = Max, MPos = I0 ; Old = Max, MPos = OPos ). matrix_min(M, Min) :- ( opaque(M) -> matrixn_min( M, Min ) ; - M = m(_, _, _, M) -> fail ). - -matrix_minarg(M, Min) :- - ( opaque(M) -> matrixn_minarg( M, Min ) ; - M = m(_Dims, _, _Size, _) -> fail ). - + M = '$matrix'(_, _, _, _, C) -> + arg(1,C,V0), foldargs(min, M, V0, Max) ; + M = [V0|L], foldl(min, L, V0, Max) ). + +min(New, Old, Max) :- ( New =< Old -> New = Max ; Old = Max ). + +matrix_minarg(M, MinArg) :- + ( opaque(M) -> matrixn_minarg( M, MinArg ); + M = '$matrix'(_, _, _, _, C) -> + arg(1,C,V0), foldargs(minarg, M, V0-0-0, _-Offset-_), matrix_offset_to_arg(M, Offset, MinArg) ; + M = [V0|L], foldl(minarg, L, V0-0-1, _Min-Off-_ ), MinArg = [Off] ). + +minarg(New, Old-OPos-I0, Min-MPos-I) :- I is I0+1, ( New < Old -> New = Min, MPos = I0 ; Old = Min, MPos = OPos ). + +matrix_to_logs(M, LogM) :- + ( opaque(M) -> matrixn_to_logs( M, LogM ) ; + M = '$matrix'(A, B, D, E, C) -> + LogM = '$matrix'(A, B, D, E, LogC), + mapargs(log, C, LogC) ; + M = [V0|L] -> maplist(log, [V0|L], LogM ) ; + LogM is log(M) ). + +log(X, Y) :- Y is log(X). + +matrix_to_exps(M, ExpM) :- + ( opaque(M) -> matrixn_to_exps( M, ExpM ) ; + M = '$matrix'(A, B, D, E, C) -> + ExpM = '$matrix'(A, B, D, E, ExpC), + mapargs(exp, C, ExpC) ; + M = [V0|L] -> maplist(exp, [V0|L], ExpM ) ; + ExpM is exp(M) ). + +exp(X, Y) :- Y is exp(X). + matrix_agg_lines(M1,+,NM) :- do_matrix_agg_lines(M1,0,NM). /* other operations: *, logprod */ @@ -400,24 +566,77 @@ matrix_agg_cols(M1,+,NM) :- /* other operations: *, logprod */ matrix_op(M1,M2,+,NM) :- - do_matrix_op(M1,M2,0,NM). + ( opaque(M1), opaque(M2) -> + do_matrix_op(M1,M2,0,NM) ; + matrix_m(M1, '$matrix'(A,B,D,E,C1)), + matrix_m(M2, '$matrix'(A,B,D,E,C2)), + mapargs(plus, C1, C2, C), + NM = '$matrix'(A,B,D,E,C) ). matrix_op(M1,M2,-,NM) :- - do_matrix_op(M1,M2,1,NM). + ( opaque(M1), opaque(M2) -> + do_matrix_op(M1,M2,1,NM) ; + matrix_m(M1, '$matrix'(A,B,D,E,C1)), + matrix_m(M2, '$matrix'(A,B,D,E,C2)), + mapargs(minus, C1, C2, C), + NM = '$matrix'(A,B,D,E,C) ). matrix_op(M1,M2,*,NM) :- - do_matrix_op(M1,M2,2,NM). + ( opaque(M1), opaque(M2) -> + do_matrix_op(M1,M2,2,NM) ; + matrix_m(M1, '$matrix'(A,B,D,E,C1)), + matrix_m(M2, '$matrix'(A,B,D,E,C2)), + mapargs(times, C1, C2, C), + NM = '$matrix'(A,B,D,E,C) ). matrix_op(M1,M2,/,NM) :- - do_matrix_op(M1,M2,3,NM). + ( opaque(M1), opaque(M2) -> + do_matrix_op(M1,M2,3,NM) ; + matrix_m(M1, '$matrix'(A,B,D,E,C1)), + matrix_m(M2, '$matrix'(A,B,D,E,C2)), + mapargs(div, C1, C2, C), + NM = '$matrix'(A,B,D,E,C) ). matrix_op(M1,M2,zdiv,NM) :- - do_matrix_op(M1,M2,5,NM). + ( opaque(M1), opaque(M2) -> + do_matrix_op(M1,M2,5,NM) ; + matrix_m(M1, '$matrix'(A,B,D,E,C1)), + matrix_m(M2, '$matrix'(A,B,D,E,C2)), + mapargs(zdiv, C1, C2, C), + NM = '$matrix'(A,B,D,E,C) ). + matrix_op_to_all(M1,+,Num,NM) :- - do_matrix_op_to_all(M1,0,Num,NM). + ( opaque(M1) -> + do_matrix_op_to_all(M1,0,Num,NM) + ; + M1 = '$matrix'(A,B,D,E,C), + mapargs(plus(Num), C, NC), + NM = '$matrix'(A,B,D,E,NC) + ). +matrix_op_to_all(M1,-,Num,NM) :- + ( opaque(M1) -> + do_matrix_op_to_all(M1,1,Num,NM) + ; + M1 = '$matrix'(A,B,D,E,C), + mapargs(minus(Num), C, NC), + NM = '$matrix'(A,B,D,E,NC) + ). matrix_op_to_all(M1,*,Num,NM) :- - do_matrix_op_to_all(M1,2,Num,NM). + ( opaque(M1) -> + do_matrix_op_to_all(M1,2,Num,NM) + ; + M1 = '$matrix'(A,B,D,E,C), + mapargs(times(Num), C, NC), + NM = '$matrix'(A,B,D,E,NC) + ). matrix_op_to_all(M1,/,Num,NM) :- % can only use floats. FNum is float(Num), - do_matrix_op_to_all(M1,3,FNum,NM). + ( opaque(M1) -> + do_matrix_op_to_all(M1,3,FNum,NM) + ; + M1 = '$matrix'(A,B,D,E,C), + mapargs(div(Num), C, NC), + NM = '$matrix'(A,B,D,E,NC) + ). + /* other operations: *, logprod */ matrix_op_to_lines(M1,M2,/,NM) :- @@ -436,21 +655,21 @@ size(N0, N1, N2) :- N2 is N0*N1. % use 1 to get access to matrix -m_get(m(Dims, _, Sz, M), Indx, V) :- - foldl2(indx, Indx, Dims, Sz, _, 1, Offset), +m_get('$matrix'(Dims, _, Sz, Bases, M), Indx, V) :- + foldl2(indx, Indx, Dims, Bases, Sz, _, 1, Offset), arg(Offset, M, V). -m_set(m(Dims, _, Sz, M), Indx, V) :- - foldl2(indx, Indx, Dims, Sz, _, 1, Offset), - nb_setarg(Offset, M, V). +m_set('$matrix'(Dims, _, Sz, Bases, M), Indx, V) :- + foldl2(indx, Indx, Dims, Bases, Sz, _, 1, Offset), + arg(Offset, M, V). -indx( I, Dim, BlkSz, NBlkSz, I0, IF) :- - NBlkSz is BlkSz div Dim, - IF is I*NBlkSz + I0. +indx( I, Dim, Base, BlkSz, NBlkSz, I0, IF) :- + NBlkSz is BlkSz div Dim , + IF is (I-Base)*NBlkSz + I0. -offset( I, Dim, BlkSz, NBlkSz, I0, IF) :- +offset( I, Dim, BlkSz, NBlkSz, Base, I0, IF) :- NBlkSz is BlkSz div Dim, - I is I0 div NBlkSz, + I is I0 div NBlkSz + Base, IF is I0 rem NBlkSz. inc(I1, I, I1) :- @@ -460,7 +679,7 @@ new_matrix(M0, Opts0, M) :- opaque(M), !, matrix_to_list(M0, L), new_matrix(L, Opts0, M). -new_matrix(m(_,_,_,C), Opts0, M) :- !, +new_matrix('$matrix'(_,_,_,_,C), Opts0, M) :- !, C =..[_|L], new_matrix(L, Opts0, M). new_matrix(C, Opts0, M) :- @@ -470,15 +689,17 @@ new_matrix(C, Opts0, M) :- new_matrix(List, Opts0, M) :- foldl2(el_list(MDims), List, Flat, [], 0, Dim), !, fix_opts(Opts0, Opts), - foldl2(process_new_opt, Opts, Type, TypeF, [Dim|MDims], Dims), + foldl2(process_new_opt, Opts, Type, TypeF, [Dim|MDims], Dims, Base), ( var(TypeF) -> guess_type( Flat, Type ) ; true ), - matrix_new( Type, Dims, Flat, M). + matrix_new( Type, Dims, Flat, M), + ( nonvar(Base) -> matrix_base(M, Base); true ). new_matrix([H|List], Opts0, M) :- length( [H|List], Size), fix_opts(Opts0, Opts), - foldl2(process_new_opt, Opts, Type, TypeF, [Size], Dims), + foldl2(process_new_opt(Base), Opts, Type, TypeF, [Size], Dims), ( var(TypeF) -> guess_type( [H|List], Type ) ; true ), - matrix_new( Type, Dims, [H|List], M). + matrix_new( Type, Dims, [H|List], M), + ( nonvar(Base) -> matrix_base(M, Base); true ). fix_opts(V, _) :- var(V), !, @@ -498,9 +719,10 @@ guess_type( List, Type ) :- Type = floats. guess_type( _List, terms ). -process_new_opt(dim=Dim, Type, Type, _, Dim) :- !. -process_new_opt(type=Type, _, Type, Dim, Dim) :- !. -process_new_opt(Opt, _, _Type, Dim, Dim) :- +process_new_opt(_Base, dim=Dim, Type, Type, _, Dim) :- !. +process_new_opt(_Base, type=Type, _, Type, Dim, Dim) :- !. +process_new_opt( Base, base=Base, Type, Type, Dim, Dim) :- !. +process_new_opt(_Base, Opt, Type, Type, Dim, Dim) :- throw(error(domain_error(opt=Opt), new_matrix)). el_list(_, V, _Els, _NEls, _I0, _I1) :- @@ -515,87 +737,106 @@ el_list([N], El, Els, NEls, I0, I1) :- append(El, NEls, Els), I1 is I0+1. -for( Domain, M:(Locals^Goal)) :- !, - global_variables( Domain, Locals, Goal, GlobalVars ), - iterate( Domain, [], GlobalVars, M:Goal, [], [] ). +for( Domain, Goal) :- + strip_module(Goal, M, Locals^NG), !, + term_variables(Domain+Locals, LocalVarsL), + LocalVars =.. [vs|LocalVarsL], + iterate( Domain, [], LocalVars, M:NG, [], [] ), + terms:reset_variables( LocalVars ). for( Domain, Goal ) :- - global_variables( Domain, [], Goal, GlobalVars ), - iterate( Domain, [], GlobalVars, Goal, [], [] ). + strip_module(Goal, M, NG), + term_variables(Domain, LocalVarsL), + LocalVars =.. [vs|LocalVarsL], + iterate( Domain, [], LocalVars, M:NG, [], [] ), + terms:reset_variables( LocalVars ). -for( Domain, M:(Locals^Goal), Inp, Out) :- !, - global_variables( Domain, Locals, Goal, GlobalVars ), - iterate( Domain, [], GlobalVars, M:Goal, [], [], Inp, Out). +for( Domain, Goal, Inp, Out) :- + strip_module(Goal, M, Locals^NG), !, + term_variables(Domain+Locals, LocalVarsL), + LocalVars =.. [vs|LocalVarsL], + iterate( Domain, [], LocalVars, M:NG, [], [], Inp, Out). for( Domain, Goal, Inp, Out ) :- - global_variables( Domain, [], Goal, GlobalVars ), - iterate( Domain, [], GlobalVars, Goal, [], [], Inp, Out ). + strip_module(Goal, M, NG), + term_variables(Domain, LocalVarsL), + LocalVars =.. [vs|LocalVarsL], + iterate( Domain, [], LocalVars, M:NG, [], [], Inp, Out ). -global_variables( Domain, Locals, Goal, GlobalVars ) :- - term_variables( Domain+Locals, Pars ), - term_variables( Goal, DGVs, Pars), - sort( DGVs, GVs ), - foldl( delv, Pars, GVs, GlobalVars ). - -delv( V, [V1|Vs], Vs) :- V == V1, !. -delv( V, [V1|Vs], [V1|NVs]) :- - delv( V, Vs, NVs). - -iterate( [], [], GlobalVars, Goal, Vs, Bs ) :- - copy_term(t(Vs, Goal, GlobalVars), t(Bs, G, GlobalVars) ), - strip_module(G, M, NG), - once( M:NG ). -iterate( [], [H|Cont], GlobalVars, Goal, Vs, Bs ) :- - iterate(H, Cont, GlobalVars, Goal, Vs, Bs ). -iterate( [H|L], Cont, GlobalVars, Goal, Vs, Bs ) :- !, +iterate( [], [], LocalVars, Goal, Vs, Bs ) :- + terms:freshen_variables(LocalVars), + Vs = Bs, + MG <== Goal, + once( MG ), + terms:reset_variables(LocalVars). +iterate( [], [H|Cont], LocalVars, Goal, Vs, Bs ) :- + iterate(H, Cont, LocalVars, Goal, Vs, Bs ). +iterate( [H|L], [], LocalVars, Goal, Vs, Bs ) :- !, + iterate(H, L, LocalVars, Goal, Vs, Bs ). +iterate( [H|L], Cont, LocalVars, Goal, Vs, Bs ) :- !, append(L, Cont, LCont), - iterate(H, LCont, GlobalVars, Goal, Vs, Bs ). -iterate( [] ins _A .. _B, Cont, GlobalVars, Goal, Vs, Bs ) :- !, - iterate(Cont, [], GlobalVars, Goal, Vs, Bs ). -iterate( [V|Ps] ins A..B, Cont, GlobalVars, Goal, Vs, Bs ) :- + iterate(H, LCont, LocalVars, Goal, Vs, Bs ). +iterate( [] ins _A .. _B, [H|L], LocalVars, Goal, Vs, Bs ) :- !, + iterate(H, L, LocalVars, Goal, Vs, Bs ). +iterate( [] ins _A .. _B, [], LocalVars, Goal, Vs, Bs ) :- !, + iterate([], [], LocalVars, Goal, Vs, Bs ). +iterate( [V|Ps] ins A..B, Cont, LocalVars, Goal, Vs, Bs ) :- eval(A, Vs, Bs, NA), eval(B, Vs, Bs, NB), ( NA > NB -> true ; A1 is NA+1, - iterate( Ps ins NA..NB, Cont, GlobalVars, Goal, [V|Vs], [NA|Bs] ), - iterate( [V|Ps] ins A1..NB, Cont, GlobalVars, Goal, Vs, Bs ) + iterate( Ps ins NA..NB, Cont, LocalVars, Goal, [V|Vs], [NA|Bs] ), + iterate( [V|Ps] ins A1..NB, Cont, LocalVars, Goal, Vs, Bs ) ). -iterate( V in A..B, Cont, GlobalVars, Goal, Vs, Bs) :- +iterate( V in A..B, Cont, LocalVars, Goal, Vs, Bs) :- var(V), eval(A, Vs, Bs, NA), eval(B, Vs, Bs, NB), ( NA > NB -> true ; A1 is NA+1, - iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs] ), - iterate( V in A1..NB, Cont, GlobalVars, Goal, Vs, Bs ) + (Cont = [H|L] -> + iterate( H, L, LocalVars, Goal, [V|Vs], [NA|Bs] ) + ; + iterate( [], [], LocalVars, Goal, [V|Vs], [NA|Bs] ) + ), + iterate( V in A1..NB, Cont, LocalVars, Goal, Vs, Bs ) ). -iterate( [], [], GlobalVars, Goal, Vs, Bs, Inp, Out ) :- - copy_term(t(Vs, Goal, GlobalVars), t(Bs, G, GlobalVars) ), - strip_module(G, M, NG), - MG <== NG, - once( call(M:MG, Inp, Out) ). -iterate( [], [H|Cont], GlobalVars, Goal, Vs, Bs, Inp, Out ) :- - iterate(H, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ). -iterate( [H|L], Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ) :- !, +iterate( [], [], LocalVars, Goal, Vs, Bs, Inp, Out ) :- + terms:freshen_variables(LocalVars), + Vs = Bs, + MG <== Goal, + once( call(MG, Inp, Out) ), + terms:reset_variables(LocalVars). +iterate( [], [H|Cont], LocalVars, Goal, Vs, Bs, Inp, Out ) :- + iterate(H, Cont, LocalVars, Goal, Vs, Bs, Inp, Out ). +iterate( [H|L], [], LocalVars, Goal, Vs, Bs, Inp, Out ) :- !, + iterate(H, L, LocalVars, Goal, Vs, Bs, Inp, Out ). +iterate( [H|L], Cont, LocalVars, Goal, Vs, Bs, Inp, Out ) :- !, append(L, Cont, LCont), - iterate(H, LCont, GlobalVars, Goal, Vs, Bs, Inp, Out ). -iterate( [] ins _A .. _B, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ) :- !, - iterate(Cont, [], GlobalVars, Goal, Vs, Bs, Inp, Out ). -iterate( [V|Ps] ins A..B, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ) :- + iterate(H, LCont, LocalVars, Goal, Vs, Bs, Inp, Out ). +iterate( [] ins _A .. _B, [], LocalVars, Goal, Vs, Bs, Inp, Out ) :- !, + iterate([], [], LocalVars, Goal, Vs, Bs, Inp, Out ). +iterate( [] ins _A .. _B, [H|L], LocalVars, Goal, Vs, Bs, Inp, Out ) :- !, + iterate(H, L, LocalVars, Goal, Vs, Bs, Inp, Out ). +iterate( [V|Ps] ins A..B, Cont, LocalVars, Goal, Vs, Bs, Inp, Out ) :- eval(A, Vs, Bs, NA), eval(B, Vs, Bs, NB), ( NA > NB -> Inp = Out ; A1 is NA+1, - iterate( Ps ins A..B, Cont, GlobalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ), - iterate( [V|Ps] ins A1..NB, Cont, GlobalVars, Goal, Vs, Bs, Mid, Out ) + iterate( Ps ins A..B, Cont, LocalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ), + iterate( [V|Ps] ins A1..NB, Cont, LocalVars, Goal, Vs, Bs, Mid, Out ) ). -iterate( V in A..B, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out) :- +iterate( V in A..B, Cont, LocalVars, Goal, Vs, Bs, Inp, Out) :- var(V), eval(A, Vs, Bs, NA), eval(B, Vs, Bs, NB), ( NA > NB -> Inp = Out ; A1 is NA+1, - iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ), - iterate( V in A1..NB, Cont, GlobalVars, Goal, Vs, Bs, Mid, Out ) + (Cont = [H|L] -> + iterate( H, L, LocalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ) + ; + iterate( [], [], LocalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ) + ), + iterate( V in A1..NB, Cont, LocalVars, Goal, Vs, Bs, Mid, Out ) ). @@ -604,3 +845,12 @@ eval(I, Vs, Bs, NI) :- copy_term(I+Vs, IA+Bs), NI <== IA. +matrix_seq(A, B, Dims, M) :- + ints(A, B, L), + matrix_new(ints, Dims, L, M). + +ints(A,B,O) :- + ( A > B -> O = [] ; O = [A|L], A1 is A+1, ints(A1,B,L) ). + +zero(_, 0). + diff --git a/library/matrix/matrix.c b/library/matrix/matrix.c index 3141b1df6..c8949ade3 100644 --- a/library/matrix/matrix.c +++ b/library/matrix/matrix.c @@ -29,6 +29,7 @@ A matrix is something of the form TYPE = {int,double} + BASE = integer #DIMS = an int DIM1 ... @@ -49,10 +50,11 @@ typedef enum { typedef enum { MAT_TYPE=0, - MAT_NDIMS=1, - MAT_SIZE=2, - MAT_ALIGN=3, - MAT_DIMS=4, + MAT_BASE=1, + MAT_NDIMS=2, + MAT_SIZE=3, + MAT_ALIGN=4, + MAT_DIMS=5, } mat_type; typedef enum { @@ -66,6 +68,9 @@ typedef enum { MAT_EXP=7 } op_type; +YAP_Functor FunctorM; +YAP_Atom AtomC; + static long int * matrix_long_data(int *mat, int ndims) { @@ -85,11 +90,13 @@ matrix_get_offset(int *mat, int* indx) /* find where we are */ for (i = 0; i < mat[MAT_NDIMS]; i++) { + int v; pos /= mat[MAT_DIMS+i]; - if (indx[i] >= mat[MAT_DIMS+i]) { + v = indx[i]-mat[MAT_BASE]; + if (v >= mat[MAT_DIMS+i]) { return off; } - off += pos*indx[i]; + off += pos*v; } return off; } @@ -144,6 +151,7 @@ new_int_matrix(int ndims, int dims[], long int data[]) } 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++) { @@ -177,6 +185,7 @@ new_float_matrix(int ndims, int dims[], double data[]) 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++) { @@ -436,6 +445,40 @@ mk_int_list(int nelems, int *data) 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) { @@ -869,6 +912,21 @@ matrix_to_list(void) return YAP_Unify(YAP_ARG2, tf); } +static int +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 int matrix_dims(void) { @@ -884,6 +942,22 @@ matrix_dims(void) return YAP_Unify(YAP_ARG2, tf); } +static int +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 int matrix_size(void) { @@ -967,7 +1041,7 @@ matrix_offset_to_arg(void) } off = YAP_IntOfTerm(ti); matrix_get_index(mat, off, indx); - tf = mk_int_list(mat[MAT_NDIMS], indx); + tf = mk_int_list2(mat[MAT_NDIMS], mat[MAT_BASE], indx); return YAP_Unify(YAP_ARG3, tf); } @@ -1133,7 +1207,27 @@ matrix_log_all2(void) return FALSE; } if (mat[MAT_TYPE] == INT_MATRIX) { - return FALSE; + 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; @@ -1221,7 +1315,27 @@ matrix_exp_all2(void) return FALSE; } if (mat[MAT_TYPE] == INT_MATRIX) { - return FALSE; + 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; @@ -2175,6 +2289,12 @@ matrix_op_to_all(void) 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; @@ -2228,6 +2348,15 @@ matrix_op_to_all(void) } } 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; @@ -3060,11 +3189,69 @@ matrix_set_all_that_disagree(void) return YAP_Unify(YAP_ARG5, tf); } +static int +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 int +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; +} + void PROTO(init_matrix, (void)); 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); @@ -3080,7 +3267,9 @@ init_matrix(void) 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); @@ -3098,8 +3287,8 @@ init_matrix(void) YAP_UserCPredicate("matrix_to_logs", matrix_log_all,1); YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1); YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1); - YAP_UserCPredicate("matrix_to_logs", matrix_log_all2,2); - YAP_UserCPredicate("matrix_to_exps", matrix_exp_all2, 2); + YAP_UserCPredicate("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); @@ -3111,6 +3300,8 @@ init_matrix(void) YAP_UserCPredicate("do_matrix_op_to_all", matrix_op_to_all, 4); YAP_UserCPredicate("do_matrix_op_to_lines", matrix_op_to_lines, 4); YAP_UserCPredicate("do_matrix_op_to_cols", matrix_op_to_cols, 4); + YAP_UserCPredicate("matrix_m", matrix_m, 2); + YAP_UserCPredicate("matrix", is_matrix, 1); } #ifdef _WIN32