improve and document matrix package

This commit is contained in:
Vítor Santos Costa 2013-09-28 11:09:32 +01:00
parent dd17f5a3aa
commit 3d863b1058
4 changed files with 826 additions and 164 deletions

View File

@ -204,6 +204,7 @@ Subnodes of Library
* DBUsage:: Information bout data base usage. * DBUsage:: Information bout data base usage.
* LineUtilities:: Line Manipulation Utilities * LineUtilities:: Line Manipulation Utilities
* Lists:: List Manipulation * Lists:: List Manipulation
* MapArgs:: Apply on Arguments of Compound Terms.
* MapList:: SWI-Compatible Apply library. * MapList:: SWI-Compatible Apply library.
* matrix:: Matrix Objects * matrix:: Matrix Objects
* MATLAB:: Matlab Interface * MATLAB:: Matlab Interface
@ -8730,6 +8731,7 @@ Library, Extensions, Built-ins, Top
* DBUsage:: Information bout data base usage. * DBUsage:: Information bout data base usage.
* Lists:: List Manipulation * Lists:: List Manipulation
* LineUtilities:: Line Manipulation Utilities * LineUtilities:: Line Manipulation Utilities
* MapArgs:: Apply on Arguments of Compound Terms.
* MapList:: SWI-Compatible Apply library. * MapList:: SWI-Compatible Apply library.
* matrix:: Matrix Objects * matrix:: Matrix Objects
* MATLAB:: Matlab Interface * MATLAB:: Matlab Interface
@ -9460,7 +9462,7 @@ unification using @code{memberchk/2}. The complexity is
See @code{ord_subtract/3}. See @code{ord_subtract/3}.
@end table @end table
@node LineUtilities, MapList, Lists, Library @node LineUtilities, MapArgs, Lists, Library
@section Line Manipulation Utilities @section Line Manipulation Utilities
@cindex Line Utilities Library @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 @section Maplist
@cindex macros @cindex macros
@ -9742,7 +9832,14 @@ applying the predicate @var{Pred} to all list elements on which
@findex foldl2/7 @findex foldl2/7
@snindex foldl2/7 @snindex foldl2/7
@cnindex 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}. @var{X} and @var{Y}.
@item foldl3(:@var{Pred}, +@var{List1}, ?@var{List2}, ?@var{X0}, @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 @cnindex scanl/7
Left scan of list. 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}) @item mapnodes(+@var{Pred}, +@var{TermIn}, ?@var{TermOut})
@findex mapnodes/3 @findex mapnodes/3
@snindex 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 of integers and floats. In contrast to dynamic arrays, these
matrices are multi-dimensional and compact. In contrast to static matrices are multi-dimensional and compact. In contrast to static
arrays. these arrays are allocated in the stack. Matrices are available 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. contact the YAP maintainers if you require extra functionality.
@table @code @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}) @item matrix_new(+@var{Type},+@var{Dims},-@var{Matrix})
@findex matrix_new/3 @findex matrix_new/3
@snindex 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{Result} is the result of applying @var{Op} to all elements of
@var{Matrix1}, with @var{Operand} as the second argument. Currently, @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. (@code{/}) are supported.
@item matrix_op_to_lines(+@var{Matrix1},+@var{Lines},+@var{Op},-@var{Result}) @item matrix_op_to_lines(+@var{Matrix1},+@var{Lines},+@var{Op},-@var{Result})

View File

@ -50,7 +50,9 @@ PROGRAMS= \
$(srcdir)/lists.yap \ $(srcdir)/lists.yap \
$(srcdir)/nb.yap \ $(srcdir)/nb.yap \
$(srcdir)/ordsets.yap \ $(srcdir)/ordsets.yap \
$(srcdir)/mapargs.yap \
$(srcdir)/maplist.yap \ $(srcdir)/maplist.yap \
$(srcdir)/maputils.yap \
$(srcdir)/matlab.yap \ $(srcdir)/matlab.yap \
$(srcdir)/matrix.yap \ $(srcdir)/matrix.yap \
$(srcdir)/prandom.yap \ $(srcdir)/prandom.yap \

View File

@ -20,7 +20,7 @@
may have a number of dimensions. These routines implement a number of may have a number of dimensions. These routines implement a number of
routine manipulation procedures. routine manipulation procedures.
matrix(Type,D1,D2,...,Dn,data(......)) '$matrix'(Type,D1,D2,...,Dn,data(......))
Type = int, float Type = int, float
@ -39,10 +39,11 @@ typedef enum {
:- module( matrix, :- module( matrix,
[op(100, yf, []), [op(100, yf, []),
(<==)/2, op(500, xfx, '<=='), (<==)/2, op(600, xfx, '<=='),
op(700, xfx, in), op(700, xfx, in),
op(700, xfx, ins), op(700, xfx, ins),
op(450, xfx, ..), % should bind more tightly than \/ op(450, xfx, ..), % should bind more tightly than \/
op(710, xfx, of), of/2,
matrix_new/3, matrix_new/3,
matrix_new/4, matrix_new/4,
matrix_new_set/4, matrix_new_set/4,
@ -99,11 +100,50 @@ typedef enum {
:- load_foreign_files([matrix], [], init_matrix). :- load_foreign_files([matrix], [], init_matrix).
:- multifile rhs_opaque/1, array_extension/2.
:- meta_predicate for(+,0), for(+,2, +, -). :- meta_predicate for(+,0), for(+,2, +, -).
:- use_module(library(maplist)). :- use_module(library(maplist)).
:- use_module(library(mapargs)).
:- use_module(library(lists)). :- 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 ) :- ( LHS <== RHS ) :-
rhs(RHS, R), rhs(RHS, R),
set_lhs( LHS, R). set_lhs( LHS, R).
@ -113,7 +153,7 @@ rhs(RHS, RHS) :- var(RHS), !.
rhs(A, A) :- atom(A), !. rhs(A, A) :- atom(A), !.
rhs(RHS, RHS) :- number(RHS), !. rhs(RHS, RHS) :- number(RHS), !.
rhs(RHS, RHS) :- opaque(RHS), !. rhs(RHS, RHS) :- opaque(RHS), !.
rhs(RHS, RHS) :- RHS = m(_, _, _, _), !. rhs(RHS, RHS) :- RHS = '$matrix'(_, _, _, _, _), !.
rhs(matrix(List), RHS) :- !, rhs(matrix(List), RHS) :- !,
rhs( List, A1), rhs( List, A1),
new_matrix(A1, [], RHS). new_matrix(A1, [], RHS).
@ -123,15 +163,6 @@ rhs(matrix(List, Opt1), RHS) :- !,
rhs(matrix(List, Opt1, Opt2), RHS) :- !, rhs(matrix(List, Opt1, Opt2), RHS) :- !,
rhs( List, A1), rhs( List, A1),
new_matrix(A1, [Opt1, Opt2], RHS). 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(dim(RHS), Dims) :- !,
rhs(RHS, X1), rhs(RHS, X1),
matrix_dims( X1, Dims ). matrix_dims( X1, Dims ).
@ -156,19 +187,23 @@ rhs(max(RHS), Size) :- !,
rhs(min(RHS), Size) :- !, rhs(min(RHS), Size) :- !,
rhs(RHS, X1), rhs(RHS, X1),
matrix_min( X1, Size ). 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(list(RHS), List) :- !,
rhs(RHS, X1), rhs(RHS, X1),
matrix_to_list( X1, List ). matrix_to_list( X1, List ).
rhs(lists(RHS), List) :- !, rhs(lists(RHS), List) :- !,
rhs(RHS, X1), rhs(RHS, X1),
matrix_to_lists( X1, List ). matrix_to_lists( X1, List ).
rhs(A=B, NA=NB) :- !, rhs('[]'(Args, RHS), Val) :-
rhs(A, NA), !,
rhs(B, NB).
rhs('[]'(Args, RHS), Val) :- !,
rhs(RHS, X1), rhs(RHS, X1),
matrix_dims( X1, Dims ), matrix_dims( X1, Dims, Bases),
maplist( index(Range), Args, Dims, NArgs), maplist( index(Range), Args, Dims, Bases, NArgs),
( (
var(Range) var(Range)
-> ->
@ -183,10 +218,23 @@ rhs('..'(I, J), [I1|Is]) :- !,
rhs([H|T], [NH|NT]) :- !, rhs([H|T], [NH|NT]) :- !,
rhs(H, NH), rhs(H, NH),
rhs(T, NT). rhs(T, NT).
rhs(':'(I, J), [I1|Is]) :- !, rhs(log(RHS), Logs ) :- !,
rhs(I, I1), rhs(RHS, X1),
rhs(J, J1), matrix_to_logs( X1, Logs ).
once( foldl(inc, Is, I1, J1) ). 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) :- rhs(S, NS) :-
S =.. [N|As], S =.. [N|As],
maplist(rhs, As, Bs), maplist(rhs, As, Bs),
@ -194,23 +242,32 @@ rhs(S, NS) :-
set_lhs(V, R) :- var(V), !, V = R. set_lhs(V, R) :- var(V), !, V = R.
set_lhs(V, R) :- number(V), !, V = R. set_lhs(V, R) :- number(V), !, V = R.
set_lhs(V, R) :- V = '[]'(Indx, M), !, set_lhs('[]'(Args, M), Val) :-
matrix_set( M, Indx, R). 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 % ranges of arguments
% %
index(Range, V, M, Indx) :- var(V), !, index(Range, V, M, Base, Indx) :- var(V), !,
index(Range, 0..(M-1), M, Indx). Max is (M-1)+Base,
index(Range, '*', M, Indx) :- !, index(Range, Base..Max, M, Base, Indx).
index(Range, 0..(M-1), M, Indx). index(Range, '*', M, Base, Indx) :- !,
index(Range, Exp, M, Indx) :- !, Max is (M-1)+Base,
index(Range, Base..Max, M, Base, Indx).
index(Range, Exp, M, _Base, Indx) :- !,
index(Exp, M, Indx0), index(Exp, M, Indx0),
( integer(Indx0) -> Indx = Indx0 ; ( integer(Indx0) -> Indx = Indx0 ;
Indx0 = [Indx] -> true ; Indx0 = [Indx] -> true ;
Indx0 = Indx, Range = range ). Indx0 = Indx, Range = range ).
index(I, _M, I ) :- integer(I), !. index(I, _M, I ) :- integer(I), !.
index(I..J, _M, [I|O] ) :- !, index(I..J, _M, [I|O] ) :- !,
I1 is I, J1 is J, I1 is I, J1 is J,
@ -218,23 +275,23 @@ index(I..J, _M, [I|O] ) :- !,
index(I:J, _M, [I|O] ) :- !, index(I:J, _M, [I|O] ) :- !,
I1 is I, J1 is J, I1 is I, J1 is J,
once( foldl(inc, O, I1, J1) ). once( foldl(inc, O, I1, J1) ).
index(I+J, _M, O ) :- index(I+J, _M, O ) :- !,
index(I, M, I1), index(I, M, I1),
index(J, M, J1), index(J, M, J1),
add_index(I1, J1, O). add_index(I1, J1, O).
index(I-J, _M, O ) :- index(I-J, _M, O ) :- !,
index(I, M, I1), index(I, M, I1),
index(J, M, J1), index(J, M, J1),
add_index(I1, J1, O). sub_index(I1, J1, O).
index(I*J, _M, O ) :- index(I*J, _M, O ) :- !,
index(I, M, I1), index(I, M, I1),
index(J, M, J1), index(J, M, J1),
O is I1*J1. O is I1*J1.
index(I div J, _M, O ) :- index(I div J, _M, O ) :- !,
index(I, M, I1), index(I, M, I1),
index(J, M, J1), index(J, M, J1),
O is I1 div J1. O is I1 div J1.
index(I rem J, _M, O ) :- index(I rem J, _M, O ) :- !,
index(I, M, I1), index(I, M, I1),
index(J, M, J1), index(J, M, J1),
O is I1 rem 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. 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. % 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), length(Dims,NDims),
foldl(size, Dims, 1, Size), foldl(size, Dims, 1, Size),
maplist(zero, Dims, Offsets),
functor( Matrix, c, Size). functor( Matrix, c, Size).
matrix_new(ints,Dims,Matrix) :- matrix_new(ints,Dims,Matrix) :-
length(Dims,NDims), length(Dims,NDims),
@ -289,9 +406,10 @@ matrix_new(floats,Dims,Matrix) :-
new_floats_matrix_set(NDims, Dims, 0.0, 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), length(Dims,NDims),
foldl(size, Dims, 1, Size), foldl(size, Dims, 1, Size),
maplist(zero, Dims, Offsets),
functor( Matrix, c, Size), functor( Matrix, c, Size),
Matrix =.. [c|Data]. Matrix =.. [c|Data].
matrix_new(ints,Dims,Data,Matrix) :- matrix_new(ints,Dims,Data,Matrix) :-
@ -304,19 +422,23 @@ matrix_new(floats,Dims,Data,Matrix) :-
matrix_dims( Mat, Dims) :- matrix_dims( Mat, Dims) :-
( opaque(Mat) -> matrixn_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) :- matrix_ndims( Mat, NDims) :-
( opaque(Mat) -> matrixn_ndims( Mat, NDims ) ; ( opaque(Mat) -> matrixn_ndims( Mat, NDims ) ;
Mat = m( _, NDims, _, _) ). Mat = '$matrix'( _, NDims, _, _, _) ).
matrix_size( Mat, Size) :- matrix_size( Mat, Size) :-
( opaque(Mat) -> matrixn_size( Mat, Size ) ; ( opaque(Mat) -> matrixn_size( Mat, Size ) ;
Mat = m( _, _, Size, _) ). Mat = '$matrix'( _, _, Size, _, _) ).
matrix_to_list( Mat, ToList) :- matrix_to_list( Mat, ToList) :-
( opaque(Mat) -> matrixn_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_to_lists( Mat, ToList) :-
matrix_dims( Mat, [D|Dims] ), matrix_dims( Mat, [D|Dims] ),
@ -350,6 +472,10 @@ add_index_prefix( [L|Els0] , H ) --> [[H|L]],
add_index_prefix( Els0 , H ). 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) :- matrix_set( Mat, Pos, El) :-
( opaque(Mat) -> matrixn_set( Mat, Pos, El ) ; ( opaque(Mat) -> matrixn_set( Mat, Pos, El ) ;
m_set(Mat, Pos, El) ). m_set(Mat, Pos, El) ).
@ -367,30 +493,70 @@ matrix_type(Matrix,Type) :-
opaque( Matrix ) -> Type = floats ; opaque( Matrix ) -> Type = floats ;
Type = terms ). 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) :- matrix_arg_to_offset(M, Index, Offset) :-
( opaque(M) -> matrixn_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) :- matrix_offset_to_arg(M, Offset, Index) :-
( opaque(M) -> matrixn_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) :- matrix_max(M, Max) :-
( opaque(M) -> matrixn_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) :- matrix_maxarg(M, MaxArg) :-
( opaque(M) -> matrixn_maxarg( M, Max ) ; ( opaque(M) -> matrixn_maxarg( M, MaxArg );
M = m(_, _, _, _) -> fail ). 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) :- matrix_min(M, Min) :-
( opaque(M) -> matrixn_min( M, Min ) ; ( opaque(M) -> matrixn_min( M, Min ) ;
M = m(_, _, _, M) -> fail ). M = '$matrix'(_, _, _, _, C) ->
arg(1,C,V0), foldargs(min, M, V0, Max) ;
matrix_minarg(M, Min) :- M = [V0|L], foldl(min, L, V0, Max) ).
( opaque(M) -> matrixn_minarg( M, Min ) ;
M = m(_Dims, _, _Size, _) -> fail ). 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) :- matrix_agg_lines(M1,+,NM) :-
do_matrix_agg_lines(M1,0,NM). do_matrix_agg_lines(M1,0,NM).
/* other operations: *, logprod */ /* other operations: *, logprod */
@ -400,24 +566,77 @@ matrix_agg_cols(M1,+,NM) :-
/* other operations: *, logprod */ /* other operations: *, logprod */
matrix_op(M1,M2,+,NM) :- 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) :- 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) :- 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) :- 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) :- 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) :- 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) :- 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) :- matrix_op_to_all(M1,/,Num,NM) :-
% can only use floats. % can only use floats.
FNum is float(Num), 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 */ /* other operations: *, logprod */
matrix_op_to_lines(M1,M2,/,NM) :- matrix_op_to_lines(M1,M2,/,NM) :-
@ -436,21 +655,21 @@ size(N0, N1, N2) :-
N2 is N0*N1. N2 is N0*N1.
% use 1 to get access to matrix % use 1 to get access to matrix
m_get(m(Dims, _, Sz, M), Indx, V) :- m_get('$matrix'(Dims, _, Sz, Bases, M), Indx, V) :-
foldl2(indx, Indx, Dims, Sz, _, 1, Offset), foldl2(indx, Indx, Dims, Bases, Sz, _, 1, Offset),
arg(Offset, M, V). arg(Offset, M, V).
m_set(m(Dims, _, Sz, M), Indx, V) :- m_set('$matrix'(Dims, _, Sz, Bases, M), Indx, V) :-
foldl2(indx, Indx, Dims, Sz, _, 1, Offset), foldl2(indx, Indx, Dims, Bases, Sz, _, 1, Offset),
nb_setarg(Offset, M, V). arg(Offset, M, V).
indx( I, Dim, BlkSz, NBlkSz, I0, IF) :- indx( I, Dim, Base, BlkSz, NBlkSz, I0, IF) :-
NBlkSz is BlkSz div Dim, NBlkSz is BlkSz div Dim ,
IF is I*NBlkSz + I0. 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, NBlkSz is BlkSz div Dim,
I is I0 div NBlkSz, I is I0 div NBlkSz + Base,
IF is I0 rem NBlkSz. IF is I0 rem NBlkSz.
inc(I1, I, I1) :- inc(I1, I, I1) :-
@ -460,7 +679,7 @@ new_matrix(M0, Opts0, M) :-
opaque(M), !, opaque(M), !,
matrix_to_list(M0, L), matrix_to_list(M0, L),
new_matrix(L, Opts0, M). new_matrix(L, Opts0, M).
new_matrix(m(_,_,_,C), Opts0, M) :- !, new_matrix('$matrix'(_,_,_,_,C), Opts0, M) :- !,
C =..[_|L], C =..[_|L],
new_matrix(L, Opts0, M). new_matrix(L, Opts0, M).
new_matrix(C, Opts0, M) :- new_matrix(C, Opts0, M) :-
@ -470,15 +689,17 @@ new_matrix(C, Opts0, M) :-
new_matrix(List, Opts0, M) :- new_matrix(List, Opts0, M) :-
foldl2(el_list(MDims), List, Flat, [], 0, Dim), !, foldl2(el_list(MDims), List, Flat, [], 0, Dim), !,
fix_opts(Opts0, Opts), 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 ), ( 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) :- new_matrix([H|List], Opts0, M) :-
length( [H|List], Size), length( [H|List], Size),
fix_opts(Opts0, Opts), 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 ), ( 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, _) :- fix_opts(V, _) :-
var(V), !, var(V), !,
@ -498,9 +719,10 @@ guess_type( List, Type ) :-
Type = floats. Type = floats.
guess_type( _List, terms ). guess_type( _List, terms ).
process_new_opt(dim=Dim, Type, Type, _, Dim) :- !. process_new_opt(_Base, dim=Dim, Type, Type, _, Dim) :- !.
process_new_opt(type=Type, _, Type, Dim, Dim) :- !. process_new_opt(_Base, type=Type, _, Type, Dim, Dim) :- !.
process_new_opt(Opt, _, _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)). throw(error(domain_error(opt=Opt), new_matrix)).
el_list(_, V, _Els, _NEls, _I0, _I1) :- el_list(_, V, _Els, _NEls, _I0, _I1) :-
@ -515,87 +737,106 @@ el_list([N], El, Els, NEls, I0, I1) :-
append(El, NEls, Els), append(El, NEls, Els),
I1 is I0+1. I1 is I0+1.
for( Domain, M:(Locals^Goal)) :- !, for( Domain, Goal) :-
global_variables( Domain, Locals, Goal, GlobalVars ), strip_module(Goal, M, Locals^NG), !,
iterate( Domain, [], GlobalVars, M:Goal, [], [] ). term_variables(Domain+Locals, LocalVarsL),
LocalVars =.. [vs|LocalVarsL],
iterate( Domain, [], LocalVars, M:NG, [], [] ),
terms:reset_variables( LocalVars ).
for( Domain, Goal ) :- for( Domain, Goal ) :-
global_variables( Domain, [], Goal, GlobalVars ), strip_module(Goal, M, NG),
iterate( Domain, [], GlobalVars, Goal, [], [] ). term_variables(Domain, LocalVarsL),
LocalVars =.. [vs|LocalVarsL],
iterate( Domain, [], LocalVars, M:NG, [], [] ),
terms:reset_variables( LocalVars ).
for( Domain, M:(Locals^Goal), Inp, Out) :- !, for( Domain, Goal, Inp, Out) :-
global_variables( Domain, Locals, Goal, GlobalVars ), strip_module(Goal, M, Locals^NG), !,
iterate( Domain, [], GlobalVars, M:Goal, [], [], Inp, Out). term_variables(Domain+Locals, LocalVarsL),
LocalVars =.. [vs|LocalVarsL],
iterate( Domain, [], LocalVars, M:NG, [], [], Inp, Out).
for( Domain, Goal, Inp, Out ) :- for( Domain, Goal, Inp, Out ) :-
global_variables( Domain, [], Goal, GlobalVars ), strip_module(Goal, M, NG),
iterate( Domain, [], GlobalVars, Goal, [], [], Inp, Out ). term_variables(Domain, LocalVarsL),
LocalVars =.. [vs|LocalVarsL],
iterate( Domain, [], LocalVars, M:NG, [], [], Inp, Out ).
global_variables( Domain, Locals, Goal, GlobalVars ) :- iterate( [], [], LocalVars, Goal, Vs, Bs ) :-
term_variables( Domain+Locals, Pars ), terms:freshen_variables(LocalVars),
term_variables( Goal, DGVs, Pars), Vs = Bs,
sort( DGVs, GVs ), MG <== Goal,
foldl( delv, Pars, GVs, GlobalVars ). once( MG ),
terms:reset_variables(LocalVars).
delv( V, [V1|Vs], Vs) :- V == V1, !. iterate( [], [H|Cont], LocalVars, Goal, Vs, Bs ) :-
delv( V, [V1|Vs], [V1|NVs]) :- iterate(H, Cont, LocalVars, Goal, Vs, Bs ).
delv( V, Vs, NVs). iterate( [H|L], [], LocalVars, Goal, Vs, Bs ) :- !,
iterate(H, L, LocalVars, Goal, Vs, Bs ).
iterate( [], [], GlobalVars, Goal, Vs, Bs ) :- iterate( [H|L], Cont, LocalVars, 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 ) :- !,
append(L, Cont, LCont), append(L, Cont, LCont),
iterate(H, LCont, GlobalVars, Goal, Vs, Bs ). iterate(H, LCont, LocalVars, Goal, Vs, Bs ).
iterate( [] ins _A .. _B, Cont, GlobalVars, Goal, Vs, Bs ) :- !, iterate( [] ins _A .. _B, [H|L], LocalVars, Goal, Vs, Bs ) :- !,
iterate(Cont, [], GlobalVars, Goal, Vs, Bs ). iterate(H, L, LocalVars, Goal, Vs, Bs ).
iterate( [V|Ps] ins A..B, Cont, GlobalVars, 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(A, Vs, Bs, NA),
eval(B, Vs, Bs, NB), eval(B, Vs, Bs, NB),
( NA > NB -> true ; ( NA > NB -> true ;
A1 is NA+1, A1 is NA+1,
iterate( Ps ins NA..NB, Cont, GlobalVars, Goal, [V|Vs], [NA|Bs] ), iterate( Ps ins NA..NB, Cont, LocalVars, Goal, [V|Vs], [NA|Bs] ),
iterate( [V|Ps] ins A1..NB, Cont, GlobalVars, Goal, Vs, 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), var(V),
eval(A, Vs, Bs, NA), eval(A, Vs, Bs, NA),
eval(B, Vs, Bs, NB), eval(B, Vs, Bs, NB),
( NA > NB -> true ; ( NA > NB -> true ;
A1 is NA+1, A1 is NA+1,
iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs] ), (Cont = [H|L] ->
iterate( V in A1..NB, Cont, GlobalVars, Goal, Vs, Bs ) 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 ) :- iterate( [], [], LocalVars, Goal, Vs, Bs, Inp, Out ) :-
copy_term(t(Vs, Goal, GlobalVars), t(Bs, G, GlobalVars) ), terms:freshen_variables(LocalVars),
strip_module(G, M, NG), Vs = Bs,
MG <== NG, MG <== Goal,
once( call(M:MG, Inp, Out) ). once( call(MG, Inp, Out) ),
iterate( [], [H|Cont], GlobalVars, Goal, Vs, Bs, Inp, Out ) :- terms:reset_variables(LocalVars).
iterate(H, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ). iterate( [], [H|Cont], LocalVars, Goal, Vs, Bs, Inp, Out ) :-
iterate( [H|L], Cont, GlobalVars, 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), append(L, Cont, LCont),
iterate(H, LCont, GlobalVars, Goal, Vs, Bs, Inp, Out ). iterate(H, LCont, LocalVars, Goal, Vs, Bs, Inp, Out ).
iterate( [] ins _A .. _B, Cont, GlobalVars, Goal, Vs, Bs, Inp, Out ) :- !, iterate( [] ins _A .. _B, [], LocalVars, Goal, Vs, Bs, Inp, Out ) :- !,
iterate(Cont, [], GlobalVars, Goal, Vs, Bs, Inp, Out ). iterate([], [], LocalVars, Goal, Vs, Bs, Inp, Out ).
iterate( [V|Ps] ins A..B, Cont, GlobalVars, 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(A, Vs, Bs, NA),
eval(B, Vs, Bs, NB), eval(B, Vs, Bs, NB),
( NA > NB -> Inp = Out ; ( NA > NB -> Inp = Out ;
A1 is NA+1, A1 is NA+1,
iterate( Ps ins A..B, Cont, GlobalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ), iterate( Ps ins A..B, Cont, LocalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ),
iterate( [V|Ps] ins A1..NB, Cont, GlobalVars, Goal, Vs, Bs, Mid, Out ) 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), var(V),
eval(A, Vs, Bs, NA), eval(A, Vs, Bs, NA),
eval(B, Vs, Bs, NB), eval(B, Vs, Bs, NB),
( NA > NB -> Inp = Out ; ( NA > NB -> Inp = Out ;
A1 is NA+1, A1 is NA+1,
iterate( Cont, [], GlobalVars, Goal, [V|Vs], [NA|Bs], Inp, Mid ), (Cont = [H|L] ->
iterate( V in A1..NB, Cont, GlobalVars, Goal, Vs, Bs, Mid, Out ) 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), copy_term(I+Vs, IA+Bs),
NI <== IA. 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).

View File

@ -29,6 +29,7 @@
A matrix is something of the form A matrix is something of the form
TYPE = {int,double} TYPE = {int,double}
BASE = integer
#DIMS = an int #DIMS = an int
DIM1 DIM1
... ...
@ -49,10 +50,11 @@ typedef enum {
typedef enum { typedef enum {
MAT_TYPE=0, MAT_TYPE=0,
MAT_NDIMS=1, MAT_BASE=1,
MAT_SIZE=2, MAT_NDIMS=2,
MAT_ALIGN=3, MAT_SIZE=3,
MAT_DIMS=4, MAT_ALIGN=4,
MAT_DIMS=5,
} mat_type; } mat_type;
typedef enum { typedef enum {
@ -66,6 +68,9 @@ typedef enum {
MAT_EXP=7 MAT_EXP=7
} op_type; } op_type;
YAP_Functor FunctorM;
YAP_Atom AtomC;
static long int * static long int *
matrix_long_data(int *mat, int ndims) matrix_long_data(int *mat, int ndims)
{ {
@ -85,11 +90,13 @@ matrix_get_offset(int *mat, int* indx)
/* find where we are */ /* find where we are */
for (i = 0; i < mat[MAT_NDIMS]; i++) { for (i = 0; i < mat[MAT_NDIMS]; i++) {
int v;
pos /= mat[MAT_DIMS+i]; 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; return off;
} }
off += pos*indx[i]; off += pos*v;
} }
return off; return off;
} }
@ -144,6 +151,7 @@ new_int_matrix(int ndims, int dims[], long int data[])
} }
mat = (int *)YAP_BlobOfTerm(blob); mat = (int *)YAP_BlobOfTerm(blob);
mat[MAT_TYPE] = INT_MATRIX; mat[MAT_TYPE] = INT_MATRIX;
mat[MAT_BASE] = 0;
mat[MAT_NDIMS] = ndims; mat[MAT_NDIMS] = ndims;
mat[MAT_SIZE] = nelems; mat[MAT_SIZE] = nelems;
for (i=0;i< ndims;i++) { for (i=0;i< ndims;i++) {
@ -177,6 +185,7 @@ new_float_matrix(int ndims, int dims[], double data[])
return blob; return blob;
mat = YAP_BlobOfTerm(blob); mat = YAP_BlobOfTerm(blob);
mat[MAT_TYPE] = FLOAT_MATRIX; mat[MAT_TYPE] = FLOAT_MATRIX;
mat[MAT_BASE] = 0;
mat[MAT_NDIMS] = ndims; mat[MAT_NDIMS] = ndims;
mat[MAT_SIZE] = nelems; mat[MAT_SIZE] = nelems;
for (i=0;i< ndims;i++) { for (i=0;i< ndims;i++) {
@ -436,6 +445,40 @@ mk_int_list(int nelems, int *data)
return tf; 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 static YAP_Term
mk_long_list(int nelems, long int *data) mk_long_list(int nelems, long int *data)
{ {
@ -869,6 +912,21 @@ matrix_to_list(void)
return YAP_Unify(YAP_ARG2, tf); 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 static int
matrix_dims(void) matrix_dims(void)
{ {
@ -884,6 +942,22 @@ matrix_dims(void)
return YAP_Unify(YAP_ARG2, tf); 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 static int
matrix_size(void) matrix_size(void)
{ {
@ -967,7 +1041,7 @@ matrix_offset_to_arg(void)
} }
off = YAP_IntOfTerm(ti); off = YAP_IntOfTerm(ti);
matrix_get_index(mat, off, indx); 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); return YAP_Unify(YAP_ARG3, tf);
} }
@ -1133,7 +1207,27 @@ matrix_log_all2(void)
return FALSE; return FALSE;
} }
if (mat[MAT_TYPE] == INT_MATRIX) { 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 { } else {
YAP_Term out; YAP_Term out;
double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata; double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata;
@ -1221,7 +1315,27 @@ matrix_exp_all2(void)
return FALSE; return FALSE;
} }
if (mat[MAT_TYPE] == INT_MATRIX) { 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 { } else {
YAP_Term out; YAP_Term out;
double *data = matrix_double_data(mat, mat[MAT_NDIMS]), *ndata; 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++) { for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] + num; 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) { } else if (op == MAT_TIMES) {
int i; int i;
@ -2228,6 +2348,15 @@ matrix_op_to_all(void)
} }
} }
break; break;
case MAT_SUB:
{
int i;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = num-data[i];
}
}
break;
case MAT_TIMES: case MAT_TIMES:
{ {
int i; int i;
@ -3060,11 +3189,69 @@ matrix_set_all_that_disagree(void)
return YAP_Unify(YAP_ARG5, tf); 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 PROTO(init_matrix, (void));
void void
init_matrix(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", new_ints_matrix, 4);
YAP_UserCPredicate("new_ints_matrix_set", new_ints_matrix_set, 4); YAP_UserCPredicate("new_ints_matrix_set", new_ints_matrix_set, 4);
YAP_UserCPredicate("new_floats_matrix", new_floats_matrix, 4); YAP_UserCPredicate("new_floats_matrix", new_floats_matrix, 4);
@ -3080,7 +3267,9 @@ init_matrix(void)
YAP_UserCPredicate("matrix_inc", do_matrix_inc2, 3); YAP_UserCPredicate("matrix_inc", do_matrix_inc2, 3);
YAP_UserCPredicate("matrix_dec", do_matrix_dec2, 3); YAP_UserCPredicate("matrix_dec", do_matrix_dec2, 3);
YAP_UserCPredicate("matrixn_to_list", matrix_to_list, 2); 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_dims, 2);
YAP_UserCPredicate("matrixn_dims", matrix_dims3, 3);
YAP_UserCPredicate("matrixn_ndims", matrix_ndims, 2); YAP_UserCPredicate("matrixn_ndims", matrix_ndims, 2);
YAP_UserCPredicate("matrixn_size", matrix_size, 2); YAP_UserCPredicate("matrixn_size", matrix_size, 2);
YAP_UserCPredicate("matrix_type_as_number", matrix_type, 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_logs", matrix_log_all,1);
YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1); YAP_UserCPredicate("matrix_to_exps", matrix_exp_all, 1);
YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1); YAP_UserCPredicate("matrix_to_exps2", matrix_exp2_all, 1);
YAP_UserCPredicate("matrix_to_logs", matrix_log_all2,2); YAP_UserCPredicate("matrixn_to_logs", matrix_log_all2,2);
YAP_UserCPredicate("matrix_to_exps", matrix_exp_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", matrix_sum_out, 3);
YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3); YAP_UserCPredicate("matrix_sum_out_several", matrix_sum_out_several, 3);
YAP_UserCPredicate("matrix_sum_logs_out", matrix_sum_out_logs, 3); YAP_UserCPredicate("matrix_sum_logs_out", 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_all", matrix_op_to_all, 4);
YAP_UserCPredicate("do_matrix_op_to_lines", matrix_op_to_lines, 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("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 #ifdef _WIN32