diff --git a/C/adtdefs.c b/C/adtdefs.c
index d8c3b3701..b669e727b 100644
--- a/C/adtdefs.c
+++ b/C/adtdefs.c
@@ -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);
diff --git a/C/c_interface.c b/C/c_interface.c
index a7dd88abb..3acad9de1 100644
--- a/C/c_interface.c
+++ b/C/c_interface.c
@@ -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;
diff --git a/C/cdmgr.c b/C/cdmgr.c
index e1b5ba857..5c0874eff 100644
--- a/C/cdmgr.c
+++ b/C/cdmgr.c
@@ -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;
}
}
diff --git a/C/cmppreds.c b/C/cmppreds.c
index a56bd9671..e9302bbc6 100644
--- a/C/cmppreds.c
+++ b/C/cmppreds.c
@@ -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;
diff --git a/C/exec.c b/C/exec.c
index ef01a5458..066d2d6a2 100644
--- a/C/exec.c
+++ b/C/exec.c
@@ -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 !! */
diff --git a/C/tracer.c b/C/tracer.c
index 528932fc2..207bf7e3e 100644
--- a/C/tracer.c
+++ b/C/tracer.c
@@ -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);
diff --git a/CLPBN/Makefile.in b/CLPBN/Makefile.in
index 7984459a0..5d44adf2a 100644
--- a/CLPBN/Makefile.in
+++ b/CLPBN/Makefile.in
@@ -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 \
diff --git a/CLPBN/clpbn/discrete_utils.yap b/CLPBN/clpbn/discrete_utils.yap
index d3a3f2a19..7a1114f2b 100644
--- a/CLPBN/clpbn/discrete_utils.yap
+++ b/CLPBN/clpbn/discrete_utils.yap
@@ -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,
diff --git a/CLPBN/clpbn/dists.yap b/CLPBN/clpbn/dists.yap
index 967dfca5e..08ec2707f 100644
--- a/CLPBN/clpbn/dists.yap
+++ b/CLPBN/clpbn/dists.yap
@@ -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, _, _, _, _).
diff --git a/CLPBN/clpbn/vel.yap b/CLPBN/clpbn/vel.yap
index 040cb2716..c574dd11b 100644
--- a/CLPBN/clpbn/vel.yap
+++ b/CLPBN/clpbn/vel.yap
@@ -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) :-
diff --git a/changes-5.1.html b/changes-5.1.html
index b500f552a..362b20855 100644
--- a/changes-5.1.html
+++ b/changes-5.1.html
@@ -17,6 +17,8 @@
Yap-5.1.3:
+- FIXED: use matrices to implement variavel elimination, also fix
+ some overflow bugs with matrices.
- FIXED: Yap_shift_visit assumed we were using AUX DL_MALLOC (obs
from Bernd Gutmann).
- FIXED: auxiliary stack overflow in term_vars (obs from Bernd
diff --git a/docs/yap.tex b/docs/yap.tex
index 104b04294..b367765b7 100644
--- a/docs/yap.tex
+++ b/docs/yap.tex
@@ -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
diff --git a/library/matrix.yap b/library/matrix.yap
index 791fe4dba..82a996b16 100644
--- a/library/matrix.yap
+++ b/library/matrix.yap
@@ -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).
+
diff --git a/library/matrix/matrix.c b/library/matrix/matrix.c
index d127dd746..56c3819cc 100644
--- a/library/matrix/matrix.c
+++ b/library/matrix/matrix.c
@@ -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