diogo/yap-6.3
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

implement sophisticated operations with matrices.

Browse Source 
git-svn-id: https://yap.svn.sf.net/svnroot/yap/trunk@2022 b08c6af1-5177-4d33-ba66-4b1c6b8b522a
This commit is contained in:
vsc 2007-11-16 14:58:41 +00:00
parent 3fb3776842
commit d908c8633f
14 changed files with 764 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
C/adtdefs.c
View File

@ -599,6 +599,10 @@ Yap_NewPredPropByFunctor(FunctorEntry *fe, Term cur_mod)
WRITE_UNLOCK(fe->FRWLock);
return NULL;
}
if (cur_mod == TermProlog)
p->ModuleOfPred = 0L;
else
p->ModuleOfPred = cur_mod;
if (fe->PropsOfFE) {
UInt hsh = PRED_HASH(fe, cur_mod, PredHashTableSize);

27
C/c_interface.c
View File

@ -10,8 +10,11 @@
* File: c_interface.c *
* comments: c_interface primitives definition *
* *
* Last rev: $Date: 2007-11-01 20:50:31 $,$Author: vsc $ *
* Last rev: $Date: 2007-11-16 14:58:40 $,$Author: vsc $ *
* $Log: not supported by cvs2svn $
* Revision 1.102 2007/11/01 20:50:31 vsc
* fix YAP_LeaveGoal (again)
*
* Revision 1.101 2007/10/29 22:48:54 vsc
* small fixes
*
@ -581,11 +584,11 @@ YAP_MkBlobTerm(unsigned int sz)
MP_INT *dst;
BACKUP_H();
I = AbsAppl(H);
while (H+(sz+sizeof(MP_INT)/sizeof(CELL)+2) > ASP-1024) {
if (!doexpand((sz+sizeof(MP_INT)/sizeof(CELL)+2)*sizeof(CELL)))
return TermNil;
}
I = AbsAppl(H);
H[0] = (CELL)FunctorBigInt;
dst = (MP_INT *)(H+1);
dst->_mp_size = 0L;
@ -749,7 +752,10 @@ YAP_MkPairTerm(Term t1, Term t2)
Term t;
BACKUP_H();
t = MkPairTerm(t1, t2);
if (H > ASP-1024)
t = TermNil;
else
t = MkPairTerm(t1, t2);
RECOVER_H();
return t;
@ -761,7 +767,10 @@ YAP_MkNewPairTerm()
Term t;
BACKUP_H();
t = Yap_MkNewPairTerm();
if (H > ASP-1024)
t = TermNil;
else
t = Yap_MkNewPairTerm();
RECOVER_H();
return t;
@ -785,7 +794,10 @@ YAP_MkApplTerm(Functor f,unsigned long int arity, Term args[])
Term t;
BACKUP_H();
t = Yap_MkApplTerm(f, arity, args);
if (H+arity > ASP-1024)
t = TermNil;
else
t = Yap_MkApplTerm(f, arity, args);
RECOVER_H();
return t;
@ -797,7 +809,10 @@ YAP_MkNewApplTerm(Functor f,unsigned long int arity)
Term t;
BACKUP_H();
t = Yap_MkNewApplTerm(f, arity);
if (H+arity > ASP-1024)
t = TermNil;
else
t = Yap_MkNewApplTerm(f, arity);
RECOVER_H();
return t;

43
C/cdmgr.c
View File

@ -11,8 +11,11 @@
* File: cdmgr.c *
* comments: Code manager *
* *
* Last rev: $Date: 2007-11-08 09:53:01 $,$Author: vsc $ *
* Last rev: $Date: 2007-11-16 14:58:40 $,$Author: vsc $ *
* $Log: not supported by cvs2svn $
* Revision 1.211 2007/11/08 09:53:01 vsc
* YAP would always say the system has tabling!
*
* Revision 1.210 2007/11/07 09:25:27 vsc
* speedup meta-calls
*
@ -1851,6 +1854,26 @@ is_fact(Term t)
return FALSE;
}
static void
mark_preds_with_this_func(Functor f, Prop p0)
{
PredEntry *pe = RepPredProp(p0);
UInt i;
pe->PredFlags |= GoalExPredFlag;
for (i = 0; i < PredHashTableSize; i++) {
PredEntry *p = PredHash[i];
while (p) {
Prop nextp = p->NextOfPE;
if (p->FunctorOfPred == f)
p->PredFlags |= GoalExPredFlag;
p = RepPredProp(nextp);
}
}
}
static int
addclause(Term t, yamop *cp, int mode, Term mod, Term *t4ref)
/*
@ -1941,18 +1964,16 @@ addclause(Term t, yamop *cp, int mode, Term mod, Term *t4ref)
} else if (IsApplTerm(tg)) {
FunctorEntry *fe = (FunctorEntry *)FunctorOfTerm(tg);
Prop p0;
int found = FALSE;
p0 = fe->PropsOfFE;
while (p0) {
PredEntry *pe = RepPredProp(p0);
pe->PredFlags |= GoalExPredFlag;
p0 = pe->NextOfPE;
found = TRUE;
}
if (!found) {
PredEntry *npe = RepPredProp(PredPropByFunc(fe,IDB_MODULE));
if (p0) {
mark_preds_with_this_func(FunctorOfTerm(tg), p0);
} else {
Term mod = CurrentModule;
PredEntry *npe;
if (CurrentModule == PROLOG_MODULE)
mod = IDB_MODULE;
npe = RepPredProp(PredPropByFunc(fe,mod));
npe->PredFlags |= GoalExPredFlag;
}
}

80
C/cmppreds.c
View File

@ -120,7 +120,12 @@ static int compare_complex(register CELL *pt0, register CELL *pt0_end, register
out = IntOfTerm(d0) - LongIntOfTerm(d1);
#ifdef USE_GMP
} else if (IsBigIntTerm(d1)) {
out = -mpz_cmp_si(Yap_BigIntOfTerm(d1), IntOfTerm(d0));
MP_INT *b1 = Yap_BigIntOfTerm(d1);
if (!mpz_size(b1)) {
out = -1;
} else {
out = -mpz_cmp_si(b1, IntOfTerm(d0));
}
#endif
} else if (IsRefTerm(d1))
out = 1 ;
@ -146,7 +151,12 @@ static int compare_complex(register CELL *pt0, register CELL *pt0_end, register
out = LongIntOfTerm(d0) - LongIntOfTerm(d1);
#ifdef USE_GMP
} else if (IsBigIntTerm(d1)) {
out = -mpz_cmp_si(Yap_BigIntOfTerm(d1), LongIntOfTerm(d0));
MP_INT *b1 = Yap_BigIntOfTerm(d1);
if (!mpz_size(b1)) {
out = -1;
} else {
out = -mpz_cmp_si(b1, LongIntOfTerm(d0));
}
#endif
} else if (IsRefTerm(d1)) {
out = 1 ;
@ -158,15 +168,30 @@ static int compare_complex(register CELL *pt0, register CELL *pt0_end, register
}
#ifdef USE_GMP
else if (IsBigIntTerm(d0)) {
if (IsIntTerm(d1))
out = mpz_cmp_si(Yap_BigIntOfTerm(d0), IntOfTerm(d1));
MP_INT *b0 = Yap_BigIntOfTerm(d0);
if (!mpz_size(b0)) {
if (IsBigIntTerm(d1)) {
MP_INT *b1 = Yap_BigIntOfTerm(d1);
out = b0-b1;
} else {
out = 1;
}
} else if (IsIntTerm(d1))
out = mpz_cmp_si(b0, IntOfTerm(d1));
else if (IsFloatTerm(d1)) {
out = 1;
} else if (IsLongIntTerm(d1))
out = mpz_cmp_si(Yap_BigIntOfTerm(d0), LongIntOfTerm(d1));
else if (IsBigIntTerm(d1))
out = mpz_cmp(Yap_BigIntOfTerm(d0), Yap_BigIntOfTerm(d1));
else if (IsRefTerm(d1))
out = mpz_cmp_si(b0, LongIntOfTerm(d1));
else if (IsBigIntTerm(d1)) {
MP_INT *b1 = Yap_BigIntOfTerm(d1);
if (!mpz_size(b1)) {
out = -1;
} else {
out = mpz_cmp(b0, b1);
}
} else if (IsRefTerm(d1))
out = 1 ;
else out = -1;
if (out != 0)
@ -327,7 +352,12 @@ compare(Term t1, Term t2) /* compare terms t1 and t2 */
}
#ifdef USE_GMP
if (IsBigIntTerm(t2)) {
return -mpz_cmp_si(Yap_BigIntOfTerm(t2),IntOfTerm(t1));
MP_INT *b1 = Yap_BigIntOfTerm(t2);
if (!mpz_size(b1)) {
return -1;
} else {
return -mpz_cmp_si(b1,IntOfTerm(t1));
}
}
#endif
if (IsRefTerm(t2))
@ -377,8 +407,14 @@ compare(Term t1, Term t2) /* compare terms t1 and t2 */
if (IsLongIntTerm(t2))
return LongIntOfTerm(t1) - LongIntOfTerm(t2);
#ifdef USE_GMP
if (IsBigIntTerm(t2))
return -mpz_cmp_si(Yap_BigIntOfTerm(t2), LongIntOfTerm(t1));
if (IsBigIntTerm(t2)) {
MP_INT *b1 = Yap_BigIntOfTerm(t2);
if (!mpz_size(b1)) {
return -1;
} else {
return -mpz_cmp_si(b1, LongIntOfTerm(t1));
}
}
#endif
if (IsRefTerm(t2))
return 1;
@ -387,15 +423,31 @@ compare(Term t1, Term t2) /* compare terms t1 and t2 */
#ifdef USE_GMP
case big_int_e:
{
if (IsIntTerm(t2))
MP_INT *b0 = Yap_BigIntOfTerm(t1);
if (!mpz_size(b0)) {
if (IsBigIntTerm(t2)) {
MP_INT *b1 = Yap_BigIntOfTerm(t2);
return b0-b1;
} else {
return 1;
}
} else if (IsIntTerm(t2))
return mpz_cmp_si(Yap_BigIntOfTerm(t1), IntOfTerm(t2));
if (IsFloatTerm(t2)) {
return 1;
}
if (IsLongIntTerm(t2))
return mpz_cmp_si(Yap_BigIntOfTerm(t1), LongIntOfTerm(t2));
if (IsBigIntTerm(t2))
return mpz_cmp(Yap_BigIntOfTerm(t1), Yap_BigIntOfTerm(t2));
if (IsBigIntTerm(t2)) {
MP_INT *b1 = Yap_BigIntOfTerm(t2);
if (!mpz_size(b1)) {
return -1;
} else {
return mpz_cmp(b0, b1);
}
}
if (IsRefTerm(t2))
return 1;
return -1;

2
C/exec.c
View File

@ -210,7 +210,7 @@ do_execute(Term t, Term mod)
goto restart_exec;
}
} else {
return(CallMetaCall(mod));
return CallMetaCall(mod);
}
}
/* now let us do what we wanted to do from the beginning !! */

4
C/tracer.c
View File

@ -161,6 +161,10 @@ low_level_trace(yap_low_level_port port, PredEntry *pred, CELL *args)
LOCK(Yap_heap_regs->low_level_trace_lock);
sc = Yap_heap_regs;
vsc_count++;
if (vsc_count == 5723445)
jmp_deb(1);
if (vsc_count < 5723000)
return;
#ifdef COMMENTED
if (vsc_count == 40650191LL)
jmp_deb(1);

1
CLPBN/Makefile.in
View File

@ -40,6 +40,7 @@ CLPBN_PROGRAMS= \
$(CLPBN_SRCDIR)/graphs.yap \
$(CLPBN_SRCDIR)/graphviz.yap \
$(CLPBN_SRCDIR)/hmm.yap \
$(CLPBN_SRCDIR)/matrix_cpt_utils.yap \
$(CLPBN_SRCDIR)/topsort.yap \
$(CLPBN_SRCDIR)/utils.yap \
$(CLPBN_SRCDIR)/vel.yap \

1
CLPBN/clpbn/discrete_utils.yap
View File

@ -69,7 +69,6 @@ project_inner_loop(I,Sz,[ok|Evs],NBase,F,Table,Ent0,Ent) :- !,
project_inner_loop(I,Sz,[_|Evs],NBase,F,Table,Ent0,Ent) :- !,
I1 is I+1,
project_inner_loop(I1,Sz,Evs,NBase,F,Table,Ent0,Ent).
%
% Given a set of variables Vs0 and a discrete CPT T0,

15
CLPBN/clpbn/dists.yap
View File

@ -7,6 +7,7 @@
dist/3,
dists/1,
get_dist/4,
get_dist_matrix/5,
get_dist_domain/2,
get_dist_params/2,
get_dist_domain_size/2,
@ -16,6 +17,8 @@
:- use_module(library(lists),[is_list/1]).
:- use_module(library(matrix),[matrix_new/4]).
/*
:- mode dist(+, -).
@ -146,6 +149,18 @@ dist(Id) :-
get_dist(Id, Type, Domain, Tab) :-
db(Id, Tab, Type, Domain, _, _).
get_dist_matrix(Id, Parents, Type, Domain, Mat) :-
db(Id, Tab, Type, Domain, _, DomainSize),
get_dsizes(Parents, Sizes, []),
matrix_new(floats, [DomainSize|Sizes], Tab, Mat).
get_dsizes([], Sizes, Sizes).
get_dsizes([P|Parents], [Sz|Sizes], Sizes0) :-
clpbn:get_atts(P,dist(Dist,_)),
get_dist_domain_size(Dist, Sz),
get_dsizes(Parents, Sizes, Sizes0).
get_dist_params(Id, Parms) :-
db(Id, Parms, _, _, _, _).

98
CLPBN/clpbn/vel.yap
View File

@ -27,7 +27,7 @@
:- use_module(library('clpbn/dists'), [
get_dist_domain_size/2,
get_dist/4]).
get_dist_matrix/5]).
:- use_module(library('clpbn/utils'), [
clpbn_not_var_member/2,
@ -36,9 +36,13 @@
:- use_module(library('clpbn/display'), [
clpbn_bind_vals/3]).
:- use_module(library('clpbn/discrete_utils'), [
project_from_CPT/3,
reorder_CPT/5]).
:- use_module(library('clpbn/matrix_cpt_utils'),
[project_from_CPT/3,
reorder_CPT/5,
get_dist_size/2,
multiply_CPTs/3,
normalise_CPT/2,
list_from_CPT/2]).
:- use_module(library(lists),
[
@ -62,9 +66,9 @@ do_vel(LVs,Vs0,AllDiffs) :-
(clpbn:output(xbif(XBifStream)) -> clpbn2xbif(XBifStream,vel,Vs) ; true),
(clpbn:output(gviz(XBifStream)) -> clpbn2gviz(XBifStream,vel,Vs,LVs) ; true),
process(LVi, LVs, tab(Dist,_,_)),
Dist =.. [_|Ps0],
normalise(Ps0,Ps),
clpbn_bind_vals(LVs,Ps,AllDiffs).
normalise_CPT(Dist,Ps),
list_from_CPT(Ps, LPs),
clpbn_bind_vals(LVs,LPs,AllDiffs).
%
% some variables might already have evidence in the data-base.
@ -93,10 +97,10 @@ find_all_clpbn_vars([V|Vs], [Var|LV], ProcessedVars, [table(I,Table,Parents,Size
var_with_deps(V, Table, Deps, Sizes, Ev, Vals) :-
clpbn:get_atts(V, [dist(Id,Parents)]),
get_dist(Id,_,Vals,OTable),
get_dist_matrix(Id,Parents,_,Vals,TAB0),
( clpbn:get_atts(V, [evidence(Ev)]) -> true ; true),
reorder_CPT([V|Parents],OTable,Deps0,Table0,Sizes0),
simplify_evidence(Deps0, Table0, Deps0, Sizes0, Table, Deps, Sizes).
reorder_CPT([V|Parents],TAB0,Deps0,TAB1,Sizes1),
simplify_evidence(Deps0, TAB1, Deps0, Sizes1, Table, Deps, Sizes).
find_all_table_deps(Tables0, LV) :-
find_dep_graph(Tables0, DepGraph0),
@ -143,8 +147,7 @@ process(LV0, InputVs, Out) :-
V \== V0, !,
%format('1 ~w: ~w~n',[V,WorkTables]),
multiply_tables(WorkTables, tab(Tab0,Deps0,_)),
Tab0 =.. [_|LTab0],
reorder_CPT(Deps0,LTab0,Deps,Tab,Sizes),
reorder_CPT(Deps0,Tab0,Deps,Tab,Sizes),
Table = tab(Tab,Deps,Sizes),
%format('2 ~w: ~w~n',[V,Table]),
project_from_CPT(V,Table,NewTable),
@ -174,9 +177,9 @@ find_best([V|LV], V0, Threshold, VF, WorkTables, [V|LVF], Inputs) :-
find_best(LV, V0, Threshold, VF, WorkTables, LVF, Inputs).
multiply_tables([Table], Table) :- !.
multiply_tables([tab(Tab1,Deps1,Szs1), tab(Tab2,Deps2,Sz2)| Tables], Out) :-
multiply_table(Tab1, Deps1, Szs1, Tab2, Deps2, Sz2, NTab, NDeps, NSz),
multiply_tables([tab(NTab,NDeps,NSz)| Tables], Out).
multiply_tables([TAB1, TAB2| Tables], Out) :-
multiply_CPTs(TAB1, TAB2, TAB),
multiply_tables([TAB| Tables], Out).
simplify_evidence([], Table, Deps, Sizes, Table, Deps, Sizes).
@ -191,56 +194,6 @@ fetch_tables([], []).
fetch_tables([var(_,_,_,_,_,_,Deps,_)|LV0], Tables) :-
append(Deps,Tables0,Tables),
fetch_tables(LV0, Tables0).
multiply_table(Tab1, Deps1, Szs1, Tab2, Deps2, Szs2, NTab, NDeps, NSzs) :-
deps_union(Deps1,Szs1,Fs10,Deps2,Szs2,Fs20,NDeps,NSzs),
factors(NSzs, Fs, Total),
factors(Fs10, Fs1, _),
factors(Fs20, Fs2, _),
elements(0, Total, Fs, Fs1, Fs2, Tab1, Tab2, LTab),
NTab =.. [t|LTab].
deps_union([],[],[],[],[],[],[],[]) :- !.
deps_union([],[],Fs1,[V2|Deps2],[Sz|Szs2],[Sz|Szs2],[V2|Deps2],[Sz|Szs2]) :- !,
mk_zeros([Sz|Szs2],Fs1).
deps_union([V1|Deps1],[Sz|Szs1],[Sz|Szs1],[],[],Fs2,[V1|Deps1],[Sz|Szs1]) :- !,
mk_zeros([Sz|Szs1],Fs2).
deps_union([V1|Deps1],[Sz|Szs1],[Sz|Fs1],[V2|Deps2],[Sz|Szs2],[Sz|Fs2],[V1|NDeps],[Sz|NSzs]) :- V1 == V2, !,
deps_union(Deps1,Szs1,Fs1,Deps2,Szs2,Fs2,NDeps,NSzs).
deps_union([V1|Deps1],[Sz1|Szs1],[Sz1|Fs1],[V2|Deps2],Szs2,[0|Fs2],[V1|NDeps],[Sz1|NSzs]) :- V1 @< V2, !,
deps_union(Deps1,Szs1,Fs1,[V2|Deps2],Szs2,Fs2,NDeps,NSzs).
deps_union([V1|Deps1],Szs1,[0|Fs1],[V2|Deps2],[Sz|Szs2],[Sz|Fs2],[V2|NDeps],[Sz|NSzs]) :-
deps_union([V1|Deps1],Szs1,Fs1,Deps2,Szs2,Fs2,NDeps,NSzs).
mk_zeros([],[]).
mk_zeros([_|Szs],[0|Fs]) :-
mk_zeros(Szs,Fs).
factors([], [], 1).
factors([0|Ls], [0|NLs], Prod) :- !,
factors(Ls, NLs, Prod).
factors([N|Ls], [Prod0|NLs], Prod) :-
factors(Ls, NLs, Prod0),
Prod is Prod0*N.
elements(Total, Total, _, _, _, _, _, []) :- !.
elements(I, Total, Fs, Fs1, Fs2, Tab1, Tab2, [El|Els]) :-
element(Fs, I, 1, Fs1, 1, Fs2, Tab1, Tab2, El),
I1 is I+1,
elements(I1, Total, Fs, Fs1, Fs2, Tab1, Tab2, Els).
element([], _, P1, [], P2, [], Tab1, Tab2, El) :-
arg(P1, Tab1, El1),
arg(P2, Tab2, El2),
El is El1*El2.
element([F|Fs], I, P1, [F1|Fs1], P2, [F2|Fs2], Tab1, Tab2, El) :-
FF is I // F,
NP1 is P1+F1*FF,
NP2 is P2+F2*FF,
NI is I mod F,
element(Fs, NI, NP1, Fs1, NP2, Fs2, Tab1, Tab2, El).
%
include([],_,_,[]).
include([var(V,P,VSz,D,Parents,Ev,Tabs,Est)|LV],tab(T,Vs,Sz),V1,[var(V,P,VSz,D,Parents,Ev,Tabs,Est)|NLV]) :-
@ -258,21 +211,6 @@ update_tables([tab(Tab0,Vs,Sz)|Tabs],[tab(Tab0,Vs,Sz)|NTabs],Table,V) :-
update_tables([_|Tabs],NTabs,Table,V) :-
update_tables(Tabs,NTabs,Table,V).
normalise(Ps0,Ps) :-
add_all(Ps0,0.0,Sum),
divide_by_sum(Ps0,Sum,Ps).
add_all([],Sum,Sum).
add_all([P|Ps0],Sum0,Sum) :-
SumI is Sum0+P,
add_all(Ps0,SumI,Sum).
divide_by_sum([],_,[]).
divide_by_sum([P|Ps0],Sum,[PN|Ps]) :-
PN is P/Sum,
divide_by_sum(Ps0,Sum,Ps).
vel_get_dist_size(V,Sz) :-
get_atts(V, [size(Sz)]), !.
vel_get_dist_size(V,Sz) :-

2
changes-5.1.html
View File

@ -17,6 +17,8 @@
<h2>Yap-5.1.3:</h2>
<ul>
<li> FIXED: use matrices to implement variavel elimination, also fix
some overflow bugs with matrices.</li>
<li> FIXED: Yap_shift_visit assumed we were using AUX DL_MALLOC (obs
from Bernd Gutmann).</li>
<li> FIXED: auxiliary stack overflow in term_vars (obs from Bernd

37
docs/yap.tex
View File

@ -8360,6 +8360,43 @@ second argument. Currently, only division (@code{/}) is supported.
second argument. Currently, only addition (@code{+}) is
supported. Notice that @var{Cols} will have n-1 dimensions.
@item matrix_shuffle(+@var{Matrix},+@var{NewOrder},-@var{Shuffle})
@findex matrix_shuffle/3
@snindex matrix_shuffle/3
@cnindex matrix_shuffle/3
Shuffle the dimensions of matrix @var{Matrix} according to
@var{NewOrder}. The list @var{NewOrder} must have all the dimensions of
@var{Matrix}, starting from 0.
@item matrix_transpose(+@var{Matrix},-@var{Transpose})
@findex matrix_reorder/3
@snindex matrix_reorder/3
@cnindex matrix_reorder/3
Transpose matrix @var{Matrix} to @var{Transpose}. Equivalent to:
@example
matrix_transpose(Matrix,Transpose) :-
matrix_shuffle(Matrix,[1,0],Transpose).
@end example
@item matrix_expand(+@var{Matrix},+@var{NewDimensions},-@var{New})
@findex matrix_expand/3
@snindex matrix_expand/3
@cnindex matrix_expand/3
Expand @var{Matrix} to occupy new dimensions. The elements in
@var{NewDimensions} are either 0, for an existing dimension, or a
positive integer with the size of the new dimension.
@item matrix_select(+@var{Matrix},+@var{Dimension},+@var{Index},-@var{New})
@findex matrix_select/3
@snindex matrix_select/3
@cnindex matrix_select/3
Select from @var{Matrix} the elements who have @var{Index} at
@var{Dimension}.
@end table
@node MATLAB, Non-Backtrackable Data Structures, matrix, Library

18
library/matrix.yap
View File

@ -61,13 +61,18 @@ typedef enum {
matrix_min/2,
matrix_minarg/2,
matrix_sum/2,
matrix_sum_out/3,
matrix_add_to_all/2,
matrix_agg_lines/3,
matrix_agg_cols/3,
matrix_op/4,
matrix_op_to_all/4,
matrix_op_to_lines/4,
matrix_op_to_cols/4
matrix_op_to_cols/4,
matrix_shuffle/3,
matrix_transpose/2,
matrix_expand/3,
matrix_select/4
]).
:- load_foreign_files([matrix], [], init_matrix).
@ -77,7 +82,7 @@ matrix_new(ints,Dims,Matrix) :-
new_ints_matrix_set(NDims, Dims, 0, Matrix).
matrix_new(floats,Dims,Matrix) :-
length(Dims,NDims),
new_float_matrix_set(NDims, Dims, 0.0, Matrix).
new_floats_matrix_set(NDims, Dims, 0.0, Matrix).
matrix_new(ints,Dims,Data,Matrix) :-
@ -85,7 +90,7 @@ matrix_new(ints,Dims,Data,Matrix) :-
new_ints_matrix(NDims, Dims, Data, Matrix).
matrix_new(floats,Dims,Data,Matrix) :-
length(Dims,NDims),
new_float_matrix(NDims, Dims, Data, Matrix).
new_floats_matrix(NDims, Dims, Data, Matrix).
matrix_new_set(ints,Dims,Elem,Matrix) :-
@ -112,6 +117,10 @@ matrix_agg_cols(M1,+,NM) :-
matrix_op(M1,M2,+,NM) :-
do_matrix_op(M1,M2,0,NM).
matrix_op(M1,M2,-,NM) :-
do_matrix_op(M1,M2,1,NM).
matrix_op(M1,M2,*,NM) :-
do_matrix_op(M1,M2,2,NM).
matrix_op_to_all(M1,+,Num,NM) :-
do_matrix_op_to_all(M1,0,Num,NM).
@ -132,4 +141,7 @@ matrix_op_to_cols(M1,M2,+,NM) :-
/* other operations: *, logprod */
matrix_transpose(M1,M2) :-
matrix_shuffle(M1,[1,0],M2).

597
library/matrix/matrix.c
View File

@ -112,8 +112,12 @@ new_int_matrix(int ndims, int dims[], long int data[])
YAP_Term blob;
int *mat;
long int *bdata;
int idims[MAX_DIMS];
/* in case we don't have enough room and need to shift the stack, we can't
really afford to keep a pointer to the global */
for (i=0;i< ndims;i++) {
idims[i] = dims[i];
nelems *= dims[i];
}
sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+nelems*sizeof(long int))/sizeof(YAP_CELL);
@ -125,7 +129,7 @@ new_int_matrix(int ndims, int dims[], long int data[])
mat[MAT_NDIMS] = ndims;
mat[MAT_SIZE] = nelems;
for (i=0;i< ndims;i++) {
mat[MAT_DIMS+i] = dims[i];
mat[MAT_DIMS+i] = idims[i];
}
bdata = matrix_long_data(mat,ndims);
if (data)
@ -141,8 +145,12 @@ new_float_matrix(int ndims, int dims[], double data[])
YAP_Term blob;
int *mat;
double *bdata;
int idims[MAX_DIMS];
/* in case we don't have enough room and need to shift the stack, we can't
really afford to keep a pointer to the global */
for (i=0;i< ndims;i++) {
idims[i] = dims[i];
nelems *= dims[i];
}
sz = ((MAT_DIMS+1)*sizeof(int)+ndims*sizeof(int)+(nelems+1)*sizeof(double)+(sizeof(YAP_CELL)-1))/sizeof(YAP_CELL);
@ -154,7 +162,7 @@ new_float_matrix(int ndims, int dims[], double data[])
mat[MAT_NDIMS] = ndims;
mat[MAT_SIZE] = nelems;
for (i=0;i< ndims;i++) {
mat[MAT_DIMS+i] = dims[i];
mat[MAT_DIMS+i] = idims[i];
}
bdata = matrix_double_data(mat,ndims);
if (data)
@ -1360,6 +1368,76 @@ matrix_double_add_data(double *nmat, int siz, double mat1[], double mat2[])
}
}
static void
matrix_long_sub_data(long int *nmat, int siz, long int mat1[], long int mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]-mat2[i];
}
}
static void
matrix_long_double_sub_data(double *nmat, int siz, long int mat1[], double mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]-mat2[i];
}
}
static void
matrix_long_double_rsub_data(double *nmat, int siz, double mat1[], long int mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat2[i]-mat1[i];
}
}
static void
matrix_double_sub_data(double *nmat, int siz, double mat1[], double mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]-mat2[i];
}
}
static void
matrix_long_mult_data(long int *nmat, int siz, long int mat1[], long int mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]*mat2[i];
}
}
static void
matrix_long_double_mult_data(double *nmat, int siz, long int mat1[], double mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]*mat2[i];
}
}
static void
matrix_double_mult_data(double *nmat, int siz, double mat1[], double mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]*mat2[i];
}
}
static int
matrix_op(void)
{
@ -1389,31 +1467,47 @@ matrix_op(void)
data1 = matrix_long_data(mat1, dims);
if (mat2[MAT_TYPE] == INT_MATRIX) {
long int *data2 = matrix_long_data(mat2, 1);
if (op == MAT_PLUS) {
long int *ndata;
long int *data2 = matrix_long_data(mat2, dims);
long int *ndata;
tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_long_data(nmat, dims);
tf = new_int_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_long_data(nmat, dims);
switch (op) {
case MAT_PLUS:
matrix_long_add_data(ndata, mat1[MAT_SIZE], data1, data2);
} else {
break;
case MAT_SUB:
matrix_long_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
case MAT_TIMES:
matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default:
return FALSE;
}
} else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
double *data2 = matrix_double_data(mat2, 1);
if (op == MAT_PLUS) {
double *ndata;
double *data2 = matrix_double_data(mat2, dims);
double *ndata;
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
switch (op) {
case MAT_PLUS:
matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
} else {
break;
case MAT_SUB:
matrix_long_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
case MAT_TIMES:
matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default:
return FALSE;
}
} else {
@ -1426,31 +1520,47 @@ matrix_op(void)
data1 = matrix_double_data(mat1, dims);
if (mat2[MAT_TYPE] == INT_MATRIX) {
long int *data2 = matrix_long_data(mat2, 1);
if (op == MAT_PLUS) {
double *ndata;
long int *data2 = matrix_long_data(mat2, dims);
double *ndata;
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
switch (op) {
case MAT_PLUS:
matrix_long_double_add_data(ndata, mat1[MAT_SIZE], data2, data1);
} else {
break;
case MAT_SUB:
matrix_long_double_rsub_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
case MAT_TIMES:
matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
break;
default:
return FALSE;
}
} else if (mat2[MAT_TYPE] == FLOAT_MATRIX) {
double *data2 = matrix_double_data(mat2, 1);
if (op == MAT_PLUS) {
double *ndata;
double *data2 = matrix_double_data(mat2, dims);
double *ndata;
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
tf = new_float_matrix(dims,mat1+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
nmat = YAP_BlobOfTerm(tf);
ndata = matrix_double_data(nmat, dims);
switch (op) {
case MAT_PLUS:
matrix_double_add_data(ndata, mat1[MAT_SIZE], data1, data2);
} else {
break;
case MAT_SUB:
matrix_double_sub_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
case MAT_TIMES:
matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default:
return FALSE;
}
} else {
@ -1673,30 +1783,415 @@ matrix_op_to_all(void)
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_double_data(mat, dims);
ndata = matrix_double_data(nmat, dims);
if (op == MAT_PLUS) {
int i;
switch(op) {
case MAT_PLUS:
{
int i;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] + num;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] + num;
}
}
} else if (op == MAT_TIMES) {
break;
case MAT_TIMES:
{
int i;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] * num;
}
} else if (op == MAT_DIV) {
int i;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] / num;
}
} else
break;
case MAT_DIV:
{
int i;
for (i = 0; i < mat[MAT_SIZE]; i++) {
ndata[i] = data[i] / num;
}
}
break;
default:
return FALSE;
}
}
return YAP_Unify(YAP_ARG4,tf);
}
/* given a matrix M and a set of dims, build a new reordered matrix to follow
the new order
*/
static int
matrix_transpose(void)
{
int ndims, i, *dims, *dimsn;
int conv[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
YAP_Term tconv, tf;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
if (!mat) {
/* Error */
return FALSE;
}
ndims = mat[MAT_NDIMS];
if (mat[MAT_TYPE] == INT_MATRIX) {
/* create a new matrix with the same size */
tf = new_int_matrix(ndims,mat+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
} else {
/* create a new matrix with the same size */
tf = new_float_matrix(ndims,mat+MAT_DIMS,NULL);
if (tf == YAP_TermNil())
return FALSE;
}
/* just in case there was an overflow */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
dims = mat+MAT_DIMS;
dimsn = nmat+MAT_DIMS;
/* we now have our target matrix, let us grab our conversion matrix */
tconv = YAP_ARG2;
for (i=0; i < ndims; i++) {
YAP_Term th;
long int j;
if (!YAP_IsPairTerm(tconv))
return FALSE;
th = YAP_HeadOfTerm(tconv);
if (!YAP_IsIntTerm(th))
return FALSE;
conv[i] = j = YAP_IntOfTerm(th);
dimsn[i] = dims[j];
tconv = YAP_TailOfTerm(tconv);
}
/*
we now got all the dimensions set up, so what we need to do
next is to copy the elements to the new matrix.
*/
if (mat[MAT_TYPE] == INT_MATRIX) {
long int *data = matrix_long_data(mat,ndims);
/* create a new matrix with the same size */
for (i=0; i< mat[MAT_SIZE]; i++) {
long int x = data[i];
int j;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(mat, i, indx);
for (j = 0; j < ndims; j++) {
nindx[j] = indx[conv[j]];
}
matrix_long_set(nmat, nindx, x);
}
} else {
double *data = matrix_double_data(mat,ndims);
/* create a new matrix with the same size */
for (i=0; i< mat[MAT_SIZE]; i++) {
double x = data[i];
long j;
matrix_get_index(mat, i, indx);
for (j = 0; j < ndims; j++)
nindx[j] = indx[conv[j]];
matrix_float_set(nmat, nindx, x);
}
}
return YAP_Unify(YAP_ARG3, tf);
}
/* given a matrix M and a set of dims, fold one of the dimensions of the
matrix on one of the elements
*/
static int
matrix_select(void)
{
int ndims, i, j, *dims, newdims, prdim, leftarg;
int indx[MAX_DIMS], nindx[MAX_DIMS];
YAP_Term tpdim, tdimarg, tf;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
if (!mat) {
/* Error */
return FALSE;
}
/* we now have our target matrix, let us grab our conversion arguments */
tpdim = YAP_ARG2;
ndims = mat[MAT_NDIMS];
dims = mat+MAT_DIMS;
if (!YAP_IsIntTerm(tpdim)) {
return FALSE;
}
prdim = YAP_IntOfTerm(tpdim);
tdimarg = YAP_ARG3;
if (!YAP_IsIntTerm(tdimarg)) {
return FALSE;
}
leftarg = YAP_IntOfTerm(tdimarg);
for (i=0, j=0; i< ndims; i++) {
if (i != prdim) {
nindx[j]= (mat+MAT_DIMS)[i];
j++;
}
}
newdims = ndims-1;
if (mat[MAT_TYPE] == INT_MATRIX) {
long int *data, *ndata;
/* create a new matrix with the same size */
tf = new_int_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_long_data(mat,ndims);
ndata = matrix_long_data(nmat,newdims);
/* create a new matrix with smaller size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
int j,k;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(nmat, i, indx);
for (j = 0, k=0; j < newdims; j++,k++) {
if (j == prdim) {
nindx[k] = leftarg;
k++;
}
nindx[k]= indx[j];
}
if (k == prdim) {
nindx[k] = leftarg;
}
ndata[i] = data[matrix_get_offset(mat, nindx)];
}
} else {
double *data, *ndata;
/* create a new matrix with the same size */
tf = new_float_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_double_data(mat,ndims);
ndata = matrix_double_data(nmat,newdims);
/* create a new matrix with the same size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
int j,k;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(nmat, i, indx);
for (j = 0, k=0; j < newdims; j++,k++) {
if (j == prdim) {
nindx[k] = leftarg;
k++;
}
nindx[k]= indx[j];
}
if (k == prdim) {
nindx[k] = leftarg;
}
ndata[i] = data[matrix_get_offset(mat, nindx)];
}
}
return YAP_Unify(YAP_ARG4, tf);
}
/* given a matrix M and a set of dims, sum out one of the dimensions
*/
static int
matrix_sum_out(void)
{
int ndims, i, j, *dims, newdims, prdim;
int indx[MAX_DIMS], nindx[MAX_DIMS];
YAP_Term tpdim, tf;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
if (!mat) {
/* Error */
return FALSE;
}
/* we now have our target matrix, let us grab our conversion arguments */
tpdim = YAP_ARG2;
ndims = mat[MAT_NDIMS];
dims = mat+MAT_DIMS;
if (!YAP_IsIntTerm(tpdim)) {
return FALSE;
}
prdim = YAP_IntOfTerm(tpdim);
newdims = ndims-1;
for (i=0, j=0; i< ndims; i++) {
if (i != prdim) {
nindx[j]= (mat+MAT_DIMS)[i];
j++;
}
}
if (mat[MAT_TYPE] == INT_MATRIX) {
long int *data, *ndata;
/* create a new matrix with the same size */
tf = new_int_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_long_data(mat,ndims);
ndata = matrix_long_data(nmat,newdims);
/* create a new matrix with smaller size */
for (i=0;i<nmat[MAT_SIZE];i++)
ndata[i] = 0;
for (i=0; i< mat[MAT_SIZE]; i++) {
int j, k;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(mat, i, indx);
for (j = 0, k=0; j < ndims; j++) {
if (j != prdim) {
nindx[k++]= indx[j];
}
}
ndata[matrix_get_offset(nmat, nindx)] += data[i];
}
} else {
double *data, *ndata;
/* create a new matrix with the same size */
tf = new_float_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_double_data(mat,ndims);
ndata = matrix_double_data(nmat,newdims);
/* create a new matrix with smaller size */
for (i=0;i<nmat[MAT_SIZE];i++)
ndata[i] = 0.0;
for (i=0; i< mat[MAT_SIZE]; i++) {
int j, k;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(mat, i, indx);
for (j = 0, k=0; j < ndims; j++) {
if (j != prdim) {
nindx[k++]= indx[j];
}
}
ndata[matrix_get_offset(nmat, nindx)] += data[i];
}
}
return YAP_Unify(YAP_ARG3, tf);
}
/* given a matrix M and a set of dims, build contract a matrix to follow
the new order
*/
static int
matrix_expand(void)
{
int ndims, i, *dims, newdims=0, olddims = 0;
int new[MAX_DIMS], indx[MAX_DIMS], nindx[MAX_DIMS];
YAP_Term tconv, tf;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
if (!mat) {
/* Error */
return FALSE;
}
/* we now have our target matrix, let us grab our conversion matrix */
tconv = YAP_ARG2;
ndims = mat[MAT_NDIMS];
dims = mat+MAT_DIMS;
for (i=0; i < MAX_DIMS; i++) {
YAP_Term th;
long int j;
if (!YAP_IsPairTerm(tconv)) {
if (tconv != YAP_TermNil())
return FALSE;
break;
}
th = YAP_HeadOfTerm(tconv);
if (!YAP_IsIntTerm(th))
return FALSE;
newdims++;
j = YAP_IntOfTerm(th);
if (j==0) {
new[i] = 0;
nindx[i] = dims[olddims];
olddims++;
} else {
new[i] = 1;
nindx[i] = j;
}
tconv = YAP_TailOfTerm(tconv);
}
if (olddims != ndims)
return FALSE;
if (mat[MAT_TYPE] == INT_MATRIX) {
long int *data, *ndata;
/* create a new matrix with the same size */
tf = new_int_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_long_data(mat,ndims);
ndata = matrix_long_data(nmat,newdims);
/* create a new matrix with the same size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
int j,k=0;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(nmat, i, indx);
for (j = 0; j < newdims; j++) {
if (!new[j])
nindx[k++] = indx[j];
}
ndata[i] = data[matrix_get_offset(mat, nindx)];
}
} else {
double *data, *ndata;
/* create a new matrix with the same size */
tf = new_float_matrix(newdims,nindx,NULL);
if (tf == YAP_TermNil())
return FALSE;
/* in case the matrix moved */
mat = (int *)YAP_BlobOfTerm(YAP_ARG1);
nmat = (int *)YAP_BlobOfTerm(tf);
data = matrix_double_data(mat,ndims);
ndata = matrix_double_data(nmat,newdims);
/* create a new matrix with the same size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
int j,k=0;
/*
not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time
*/
matrix_get_index(nmat, i, indx);
for (j = 0; j < newdims; j++) {
if (!new[j])
nindx[k++] = indx[j];
}
ndata[i] = data[matrix_get_offset(mat, nindx)];
}
}
return YAP_Unify(YAP_ARG3, tf);
}
void PROTO(init_matrix, (void));
void
@ -1726,7 +2221,11 @@ init_matrix(void)
YAP_UserCPredicate("matrix_min", matrix_min, 2);
YAP_UserCPredicate("matrix_minarg", matrix_minarg, 2);
YAP_UserCPredicate("matrix_sum", matrix_sum, 2);
YAP_UserCPredicate("matrix_shuffle", matrix_transpose, 3);
YAP_UserCPredicate("matrix_expand", matrix_expand, 3);
YAP_UserCPredicate("matrix_select", matrix_select, 4);
YAP_UserCPredicate("matrix_add_to_all", matrix_sum, 2);
YAP_UserCPredicate("matrix_sum_out", matrix_sum_out, 3);
YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);