junction tree algorithm

git-svn-id: https://yap.svn.sf.net/svnroot/yap/trunk@2031 b08c6af1-5177-4d33-ba66-4b1c6b8b522a
This commit is contained in:
vsc 2007-11-28 23:52:14 +00:00
parent 3beda27d14
commit 1bd96722de
15 changed files with 925 additions and 88 deletions

View File

@ -10,8 +10,11 @@
* * * *
* File: absmi.c * * File: absmi.c *
* comments: Portable abstract machine interpreter * * comments: Portable abstract machine interpreter *
* Last rev: $Date: 2007-11-26 23:43:07 $,$Author: vsc $ * * Last rev: $Date: 2007-11-28 23:52:14 $,$Author: vsc $ *
* $Log: not supported by cvs2svn $ * $Log: not supported by cvs2svn $
* Revision 1.231 2007/11/26 23:43:07 vsc
* fixes to support threads and assert correctly, even if inefficiently.
*
* Revision 1.230 2007/11/08 15:52:15 vsc * Revision 1.230 2007/11/08 15:52:15 vsc
* fix some bugs in new dbterm code. * fix some bugs in new dbterm code.
* *
@ -2688,11 +2691,6 @@ Yap_absmi(int inp)
#ifndef NO_CHECKING #ifndef NO_CHECKING
check_stack(NoStackCall, H); check_stack(NoStackCall, H);
#endif #endif
if (pt->PredFlags & LogUpdatePredFlag) {
if (pt->OpcodeOfPred != LOCKPRED_OPCODE &&
pt->ModuleOfPred != IDB_MODULE && pt->OpcodeOfPred != UNDEF_OPCODE)
fprintf(stderr,"OOPS\n");
}
ENV = ENV_YREG; ENV = ENV_YREG;
/* Try to preserve the environment */ /* Try to preserve the environment */
ENV_YREG = (CELL *) (((char *) ENV_YREG) + PREG->u.sla.s); ENV_YREG = (CELL *) (((char *) ENV_YREG) + PREG->u.sla.s);

View File

@ -11,8 +11,11 @@
* File: cdmgr.c * * File: cdmgr.c *
* comments: Code manager * * comments: Code manager *
* * * *
* Last rev: $Date: 2007-11-26 23:43:07 $,$Author: vsc $ * * Last rev: $Date: 2007-11-28 23:52:14 $,$Author: vsc $ *
* $Log: not supported by cvs2svn $ * $Log: not supported by cvs2svn $
* Revision 1.213 2007/11/26 23:43:07 vsc
* fixes to support threads and assert correctly, even if inefficiently.
*
* Revision 1.212 2007/11/16 14:58:40 vsc * Revision 1.212 2007/11/16 14:58:40 vsc
* implement sophisticated operations with matrices. * implement sophisticated operations with matrices.
* *
@ -2966,7 +2969,7 @@ p_undefined(void)
{ /* '$undefined'(P,Mod) */ { /* '$undefined'(P,Mod) */
PredEntry *pe; PredEntry *pe;
pe = get_pred(Deref(ARG1), CurrentModule, "undefined/1"); pe = get_pred(Deref(ARG1), Deref(ARG2), "undefined/1");
if (EndOfPAEntr(pe)) if (EndOfPAEntr(pe))
return TRUE; return TRUE;
LOCK(pe->PELock); LOCK(pe->PELock);
@ -2992,7 +2995,7 @@ p_kill_dynamic(void)
{ /* '$kill_dynamic'(P,M) */ { /* '$kill_dynamic'(P,M) */
PredEntry *pe; PredEntry *pe;
pe = get_pred(Deref(ARG1), CurrentModule, "kill_dynamic/1"); pe = get_pred(Deref(ARG1), Deref(ARG2), "kill_dynamic/1");
if (EndOfPAEntr(pe)) if (EndOfPAEntr(pe))
return TRUE; return TRUE;
LOCK(pe->PELock); LOCK(pe->PELock);

View File

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

View File

@ -29,6 +29,9 @@
check_if_vel_done/1 check_if_vel_done/1
]). ]).
:- use_module('clpbn/jt', [jt/3
]).
:- use_module('clpbn/bnt', [do_bnt/3, :- use_module('clpbn/bnt', [do_bnt/3,
check_if_bnt_done/1 check_if_bnt_done/1
]). ]).
@ -57,7 +60,7 @@
:- dynamic solver/1,output/1,use/1. :- dynamic solver/1,output/1,use/1.
solver(vel). solver(jt).
%output(xbif(user_error)). %output(xbif(user_error)).
%output(gviz(user_error)). %output(gviz(user_error)).
@ -142,6 +145,8 @@ get_clpbn_vars([_|GVars],CLPBNGVars) :-
write_out(vel, GVars, AVars, DiffVars) :- write_out(vel, GVars, AVars, DiffVars) :-
vel(GVars, AVars, DiffVars). vel(GVars, AVars, DiffVars).
write_out(jt, GVars, AVars, DiffVars) :-
jt(GVars, AVars, DiffVars).
write_out(gibbs, GVars, AVars, DiffVars) :- write_out(gibbs, GVars, AVars, DiffVars) :-
gibbs(GVars, AVars, DiffVars). gibbs(GVars, AVars, DiffVars).
write_out(bnt, GVars, AVars, DiffVars) :- write_out(bnt, GVars, AVars, DiffVars) :-
@ -224,6 +229,9 @@ bind_clpbn(_, Var, _, _, _, _) :-
bind_clpbn(_, Var, _, _, _, _) :- bind_clpbn(_, Var, _, _, _, _) :-
use(vel), use(vel),
check_if_vel_done(Var), !. check_if_vel_done(Var), !.
bind_clpbn(_, Var, _, _, _, _) :-
use(jt),
check_if_vel_done(Var), !.
bind_clpbn(T, Var, Key0, _, _, _) :- bind_clpbn(T, Var, Key0, _, _, _) :-
get_atts(Var, [key(Key)]), !, get_atts(Var, [key(Key)]), !,
( (

469
CLPBN/clpbn/jt.yap Normal file
View File

@ -0,0 +1,469 @@
:- use_module(library(dgraphs),
[dgraph_new/1,
dgraph_add_edges/3,
dgraph_add_vertex/3,
dgraph_add_vertices/3,
dgraph_edges/2,
dgraph_vertices/2,
dgraph_transpose/2,
dgraph_to_ugraph/2,
ugraph_to_dgraph/2,
dgraph_neighbors/3
]).
:- use_module(library(undgraphs),
[undgraph_new/1,
undgraph_add_edge/4,
undgraph_add_edges/3,
undgraph_del_vertex/3,
undgraph_del_vertices/3,
undgraph_vertices/2,
undgraph_edges/2,
undgraph_neighbors/3,
undgraph_edge/3,
dgraph_to_undgraph/2
]).
:- use_module(library(wundgraphs),
[wundgraph_new/1,
wundgraph_max_tree/3,
wundgraph_add_edges/3,
wundgraph_add_vertices/3,
wundgraph_to_undgraph/2
]).
:- use_module(library(rbtrees),
[rb_new/1,
rb_insert/4,
rb_lookup/3]).
:- use_module(library(ordsets),
[ord_subset/2,
ord_insert/3,
ord_intersection/3,
ord_del_element/3,
ord_memberchk/2]).
:- use_module(library(lists),
[reverse/2]).
:- use_module(library('clpbn/dists'),
[get_dist_domain_size/2,
get_dist_domain/2,
get_dist_matrix/5]).
:- use_module(library('clpbn/matrix_cpt_utils'),
[project_from_CPT/3,
reorder_CPT/5,
unit_CPT/2,
multiply_CPTs/4,
divide_CPTs/3,
normalise_CPT/2,
expand_CPT/4,
get_CPT_sizes/2,
reset_CPT_that_disagrees/5,
sum_out_from_CPT/4,
list_from_CPT/2]).
:- use_module(library('clpbn/display'), [
clpbn_bind_vals/3]).
jt(LVs,Vs0,AllDiffs) :-
get_graph(Vs0, BayesNet, CPTs, Evidence),
build_jt(BayesNet, CPTs, JTree),
% JTree is a dgraph
% now our tree has cpts
fill_with_cpts(JTree, NewTree),
propagate_evidence(Evidence, NewTree, EvTree),
message_passing(EvTree, MTree),
get_margin(MTree, LVs, LPs),
clpbn_bind_vals(LVs,LPs,AllDiffs).
get_graph(LVs, BayesNet, CPTs, Evidence) :-
run_vars(LVs, Edges, Vertices, CPTs, Evidence),
dgraph_new(V0),
dgraph_add_edges(Edges, V0, V1),
dgraph_add_vertices(Vertices, V1, V2),
dgraph_to_ugraph(V2, BayesNet).
run_vars([], [], [], [], []).
run_vars([V|LVs], Edges, [V|Vs], [CPTVars-dist([V|Parents],Id)|CPTs], Ev) :-
clpbn:get_atts(V, [dist(Id,Parents)]),
add_evidence(V, Id, Ev, Ev0),
sort([V|Parents],CPTVars),
add_edges(Parents, V, Edges, Edges0),
run_vars(LVs, Edges0, Vs, CPTs, Ev0).
add_evidence(V, Id, [e(V,P)|Evs], Evs) :-
clpbn:get_atts(V, [evidence(Ev)]), !,
get_dist_domain(Id, D),
find_nth0(D, Ev, 0, P).
add_evidence(_, _, Evs, Evs).
find_nth0([Id|_], Id, P, P) :- !.
find_nth0([_|D], Id, P0, P) :-
P1 is P0+1,
find_nth0(D, Id, P1, P).
add_edges([], _, Edges, Edges).
add_edges([P|Parents], V, [V-P|Edges], Edges0) :-
add_edges(Parents, V, Edges, Edges0).
build_jt(BayesNet, CPTs, Tree) :-
init_undgraph(BayesNet, Moral0),
moralised(BayesNet, Moral0, Markov),
undgraph_vertices(Markov, Vertices),
triangulate(Vertices, Markov, Markov, _, Cliques0),
cliques(Cliques0, EndCliques),
wundgraph_max_tree(EndCliques, J0Tree, _),
root(J0Tree, JTree),
populate(CPTs, JTree, Tree).
initial_graph(_,Parents, CPTs) :-
test_graph(0, Graph0, CPTs),
dgraph_new(V0),
dgraph_add_edges(Graph0, V0, V1),
% OK, this is a bit silly, I could have written the transposed graph
% from the very beginning.
dgraph_transpose(V1, V2),
dgraph_to_ugraph(V2, Parents).
problem_graph([], []).
problem_graph([V|BNet], GraphF) :-
clpbn:get_atts(V, [dist(_,_,Parents)]),
add_parents(Parents, V, Graph0, GraphF),
problem_graph(BNet, Graph0).
add_parents([], _, Graph, Graph).
add_parents([P|Parents], V, Graph0, [P-V|GraphF]) :-
add_parents(Parents, V, Graph0, GraphF).
% From David Page's lectures
test_graph(0,
[1-3,2-3,2-4,5-4,5-7,10-7,10-9,11-9,3-6,4-6,7-8,9-8,6-12,8-12],
[[1]-a,
[2]-b,
[1,2,3]-c,
[2,4,5]-d,
[5]-e,
[3,4,6]-f,
[5,7,10]-g,
[7,8,9]-h,
[9]-i,
[10]-j,
[11]-k,
[6,8,12]-l
]).
test_graph(1,[a-b,a-c,b-d,c-e,d-f,e-f],
[]).
init_undgraph(Parents, UndGraph) :-
ugraph_to_dgraph(Parents, DGraph),
dgraph_to_undgraph(DGraph, UndGraph).
get_par_keys([], []).
get_par_keys([P|Parents],[K|KPars]) :-
clpbn:get_atts(P, [key(K)]),
get_par_kets(Parents,KPars).
moralised([],Moral,Moral).
moralised([_-KPars|Ks],Moral0,MoralF) :-
add_moral_edges(KPars, Moral0, MoralI),
moralised(Ks,MoralI,MoralF).
add_moral_edges([], Moral, Moral).
add_moral_edges([_], Moral, Moral).
add_moral_edges([K1,K2|KPars], Moral0, MoralF) :-
undgraph_add_edge(K1,K2,Moral0,MoralI),
add_moral_edges([K1|KPars], MoralI, MoralJ),
add_moral_edges([K2|KPars],MoralJ,MoralF).
triangulate([], _, Triangulated, Triangulated, []) :- !.
triangulate(Vertices, S0, T0, Tf, Cliques) :-
choose(Vertices, S0, +inf, [], -1, BestVertex, _, Cliques0, Cliques, Edges),
ord_del_element(Vertices, BestVertex, NextVertices),
undgraph_add_edges(Edges, T0, T1),
undgraph_del_vertex(BestVertex, S0, Si),
undgraph_add_edges(Edges, Si, Si2),
triangulate(NextVertices, Si2, T1, Tf, Cliques0).
choose([], _, _, NewEdges, Best, Best, Clique, Cliques0, [Clique|Cliques0], NewEdges).
choose([V|Vertices], Graph, Score0, _, _, Best, _, Cliques0, Cliques, EdgesF) :-
undgraph_neighbors(V, Graph, Neighbors),
ord_insert(Neighbors, V, PossibleClique),
new_edges(Neighbors, Graph, NewEdges),
(
% simplicial edge
NewEdges == []
->
!,
Best = V,
NewEdges = EdgesF,
length(PossibleClique,L),
Cliques = [L-PossibleClique|Cliques0]
;
length(PossibleClique,CL),
CL < Score0, !,
choose(Vertices,Graph,CL,NewEdges, V, Best, CL-PossibleClique, Cliques0,Cliques,EdgesF)
).
choose([_|Vertices], Graph, Score0, Edges0, BestSoFar, Best, Clique, Cliques0, Cliques, EdgesF) :-
choose(Vertices,Graph,Score0,Edges0, BestSoFar, Best, Clique, Cliques0,Cliques,EdgesF).
new_edges([], _, []).
new_edges([N|Neighbors], Graph, NewEdgesF) :-
new_edges(Neighbors,N,Graph,NewEdges0, NewEdgesF),
new_edges(Neighbors, Graph, NewEdges0).
new_edges([],_,_,NewEdges, NewEdges).
new_edges([N1|Neighbors],N,Graph,NewEdges0, NewEdgesF) :-
undgraph_edge(N, N1, Graph), !,
new_edges(Neighbors,N,Graph,NewEdges0, NewEdgesF).
new_edges([N1|Neighbors],N,Graph,NewEdges0, [N-N1|NewEdgesF]) :-
new_edges(Neighbors,N,Graph,NewEdges0, NewEdgesF).
%
% This is simple stuff, I just have to remove cliques that
% are subset of the others.
%
cliques(CliqueList, CliquesF) :-
wundgraph_new(Cliques0),
% first step, order by size,
keysort(CliqueList,Sort),
reverse(Sort, Rev),
get_links(Rev, [], Vertices, [], Edges),
wundgraph_add_vertices(Vertices, Cliques0, CliquesI),
wundgraph_add_edges(Edges, CliquesI, CliquesF).
% stupid quadratic algorithm, needs to be improved.
get_links([], Vertices, Vertices, Edges, Edges).
get_links([Sz-Clique|Cliques], SoFar, Vertices, Edges0, Edges) :-
add_clique_edges(SoFar, Clique, Sz, Edges0, EdgesI), !,
get_links(Cliques, [Clique|SoFar], Vertices, EdgesI, Edges).
get_links([_|Cliques], SoFar, Vertices, Edges0, Edges) :-
get_links(Cliques, SoFar, Vertices, Edges0, Edges).
add_clique_edges([], _, _, Edges, Edges).
add_clique_edges([Clique1|Cliques], Clique, Sz, Edges0, EdgesF) :-
ord_intersection(Clique1, Clique, Int),
Int \== Clique,
(
Int = [] ->
add_clique_edges(Cliques, Clique, Sz, Edges0, EdgesF)
;
% we connect
length(Int, LSz),
add_clique_edges(Cliques, Clique, Sz, [Clique-(Clique1-LSz)|Edges0], EdgesF)
).
root(WTree, JTree) :-
wundgraph_to_undgraph(WTree, Tree),
remove_leaves(Tree, SmallerTree),
undgraph_vertices(SmallerTree, InnerVs),
pick_root(InnerVs, Root),
rb_new(M0),
build_tree(Root, M0, Tree, JTree, _).
remove_leaves(Tree, SmallerTree) :-
undgraph_vertices(Tree, Vertices),
Vertices = [_,_,_|_],
get_leaves(Vertices, Tree, Leaves),
Leaves = [_|_], !,
undgraph_del_vertices(Leaves, Tree, NTree),
remove_leaves(NTree, SmallerTree).
remove_leaves(Tree, Tree).
get_leaves([], _, []).
get_leaves([V|Vertices], Tree, [V|Leaves]) :-
undgraph_neighbors(V, Tree, [_]), !,
get_leaves(Vertices, Tree, Leaves).
get_leaves([_|Vertices], Tree, Leaves) :-
get_leaves(Vertices, Tree, Leaves).
pick_root([V|_],V).
direct_edges([], _, [], []) :- !.
direct_edges([], NewVs, RemEdges, Directed) :-
direct_edges(RemEdges, NewVs, [], Directed).
direct_edges([V1-V2|Edges], NewVs0, RemEdges, [V1-V2|Directed]) :-
ord_memberchk(V1, NewVs0), !,
ord_insert(NewVs0, V2, NewVs),
direct_edges(Edges, NewVs, RemEdges, Directed).
direct_edges([V1-V2|Edges], NewVs0, RemEdges, [V2-V1|Directed]) :-
ord_memberchk(V2, NewVs0), !,
ord_insert(NewVs0, V1, NewVs),
direct_edges(Edges, NewVs, RemEdges, Directed).
direct_edges([Edge|Edges], NewVs, RemEdges, Directed) :-
direct_edges(Edges, NewVs, [Edge|RemEdges], Directed).
populate(CPTs, JTree, NewJTree) :-
keysort(CPTs, KCPTs),
populate_cliques(JTree, KCPTs, NewJTree, []).
populate_cliques(tree(Clique,Kids), CPTs, tree(Clique-MyCPTs,NewKids), RemCPTs) :-
get_cpts(CPTs, Clique, MyCPTs, MoreCPTs),
populate_trees_with_cliques(Kids, MoreCPTs, NewKids, RemCPTs).
populate_trees_with_cliques([], MoreCPTs, [], MoreCPTs).
populate_trees_with_cliques([Node|Kids], MoreCPTs, [NewNode|NewKids], RemCPts) :-
populate_cliques(Node, MoreCPTs, NewNode, ExtraCPTs),
populate_trees_with_cliques(Kids, ExtraCPTs, NewKids, RemCPts).
get_cpts([], _, [], []).
get_cpts([CPT|CPts], [], [], [CPT|CPts]) :- !.
get_cpts([[I|MCPT]-Info|CPTs], [J|Clique], MyCPTs, MoreCPTs) :-
compare(C,I,J),
( C == < ->
% our CPT cannot be a part of the clique.
MoreCPTs = [[I|MCPT]-Info|LeftoverCPTs],
get_cpts(CPTs, [J|Clique], MyCPTs, LeftoverCPTs)
;
C == = ->
% our CPT cannot be a part of the clique.
get_cpt(MCPT, Clique, I, Info, MyCPTs, MyCPTs0, MoreCPTs, MoreCPTs0),
get_cpts(CPTs, [J|Clique], MyCPTs0, MoreCPTs0)
;
% the first element in our CPT may not be in a clique
get_cpts([[I|MCPT]-Info|CPTs], Clique, MyCPTs, MoreCPTs)
).
get_cpt(MCPT, Clique, I, Info, [[I|MCPT]-Info|MyCPTs], MyCPTs, MoreCPTs, MoreCPTs) :-
ord_subset(MCPT, Clique), !.
get_cpt(MCPT, _, I, Info, MyCPTs, MyCPTs, [[I|MCPT]-Info|MoreCPTs], MoreCPTs).
translate_edges([], [], []).
translate_edges([E1-E2|Edges], [(E1-A)-(E2-B)|NEdges], [E1-A,E2-B|Vs]) :-
translate_edges(Edges, NEdges, Vs).
match_vs(_,[]).
match_vs([K-A|Cls],[K1-B|KVs]) :-
compare(C, K, K1),
(C == = ->
A = B,
match_vs([K-A|Cls], KVs)
;
C = < ->
match_vs(Cls,[K1-B|KVs])
;
match_vs([K-A|Cls],KVs)
).
fill_with_cpts(tree(Clique-Dists,Leafs), tree(Clique-NewDists,NewLeafs)) :-
compile_cpts(Dists, Clique, NewDists),
fill_tree_with_cpts(Leafs, NewLeafs).
fill_tree_with_cpts([], []).
fill_tree_with_cpts([L|Leafs], [NL|NewLeafs]) :-
fill_with_cpts(L, NL),
fill_tree_with_cpts(Leafs, NewLeafs).
transform([], []).
transform([Clique-Dists|Nodes],[Clique-NewDist|NewNodes]) :-
compile_cpts(Dists, Clique, NewDist),
transform(Nodes, NewNodes).
compile_cpts([Vs-dist(OVs,Id)|Dists], Clique, TAB) :-
OVs = [_|Ps], !,
get_dist_matrix(Id, Ps, _, _, TAB0),
reorder_CPT(OVs, TAB0, Vs, TAB1, Sz1),
multiply_dists(Dists,Vs,TAB1,Sz1,Vars2,ITAB),
expand_CPT(ITAB,Vars2,Clique,TAB).
compile_cpts([], [V|Clique], TAB) :-
unit_CPT(V, CPT0),
expand_CPT(CPT0, [V], [V|Clique], TAB).
multiply_dists([],Vs,TAB,_,Vs,TAB).
multiply_dists([Vs-dist(OVs,Id)|Dists],MVs,TAB2,Sz2,FVars,FTAB) :-
OVs = [_|Ps],
get_dist_matrix(Id, Ps, _, _, TAB0),
reorder_CPT(OVs, TAB0, Vs, TAB1, Sz1),
multiply_CPTs(tab(TAB1,Vs,Sz1),tab(TAB2,MVs,Sz2),tab(TAB3,NVs,Sz),_),
multiply_dists(Dists,NVs,TAB3,Sz,FVars,FTAB).
build_tree(Root, Leafs, WTree, tree(Root,Leaves), NewLeafs) :-
rb_insert(Leafs, Root, [], Leafs0),
undgraph_neighbors(Root, WTree, Children),
build_trees(Children, Leafs0, WTree, Leaves, NewLeafs).
build_trees( [], Leafs, _, [], Leafs).
build_trees([V|Children], Leafs, WTree, NLeaves, NewLeafs) :-
% back pointer
rb_lookup(V, _, Leafs), !,
build_trees(Children, Leafs, WTree, NLeaves, NewLeafs).
build_trees([V|Children], Leafs, WTree, [VT|NLeaves], NewLeafs) :-
build_tree(V, Leafs, WTree, VT, Leafs1),
build_trees(Children, Leafs1, WTree, NLeaves, NewLeafs).
propagate_evidence([], NewTree, NewTree).
propagate_evidence([e(V,P)|Evs], Tree0, NewTree) :-
add_evidence_to_matrix(Tree0, V, P, Tree1), !,
propagate_evidence(Evs, Tree1, NewTree).
add_evidence_to_matrix(tree(Clique-Dist,Kids), V, P, tree(Clique-NDist,Kids)) :-
ord_memberchk(V, Clique), !,
reset_CPT_that_disagrees(Dist, Clique, V, P, NDist).
add_evidence_to_matrix(tree(C,Kids), V, P, tree(C,NKids)) :-
add_evidence_to_kids(Kids, V, P, NKids).
add_evidence_to_kids([K|Kids], V, P, [NK|Kids]) :-
add_evidence_to_matrix(K, V, P, NK), !.
add_evidence_to_kids([K|Kids], V, P, [K|NNKids]) :-
add_evidence_to_kids(Kids, V, P, NNKids).
message_passing(tree(Clique-Dist,Kids), tree(Clique-NDist,NKids)) :-
get_CPT_sizes(Dist, Sizes),
upward(Kids, Clique, tab(Dist, Clique, Sizes), IKids, ITab),
ITab = tab(NDist, _, _),
downward(IKids, Clique, ITab, NKids).
upward([], _, Dist, [], Dist).
upward([tree(Clique1-Dist1,DistKids)|Kids], Clique, Tab, [tree(Clique1-(NewDist1,EDist1),NDistKids)|Kids], NewTab) :-
get_CPT_sizes(Dist1, Sizes1),
upward(DistKids, Clique1, tab(Dist1,Clique1,Sizes1), NDistKids, NewTab1),
NewTab1 = tab(NewDist1,_,_),
ord_intersection(Clique1, Clique, Int),
sum_out_from_CPT(Int, NewDist1, Clique1, Tab1),
multiply_CPTs(Tab, Tab1, NewTab, EDist1).
downward([], _, _, []).
downward([tree(Clique1-(Dist1,Msg1),DistKids)|Kids], Clique, Tab, [tree(Clique1-NDist1,NDistKids)|Kids]) :-
get_CPT_sizes(Dist1, Sizes1),
ord_intersection(Clique1, Clique, Int),
Tab = tab(Dist,_,_),
divide_CPTs(Dist, Msg1, Div),
sum_out_from_CPT(Int, Div, Clique, STab),
multiply_CPTs(STab, tab(Dist1, Clique1, Sizes1), NewTab, _),
NewTab = tab(NDist1,_,_),
downward(DistKids, Clique1, NewTab, NDistKids).
get_margin(NewTree, LVs0, LPs) :-
sort(LVs0, LVs),
find_clique(NewTree, LVs, Clique, Dist),
sum_out_from_CPT(LVs, Dist, Clique, tab(TAB,_,_)),
reorder_CPT(LVs, TAB, LVs0, NTAB, _),
normalise_CPT(NTAB, Ps),
list_from_CPT(Ps, LPs).
find_clique(tree(Clique-Dist,_), LVs, Clique, Dist) :-
ord_subset(LVs, Clique), !.
find_clique(tree(_,Kids), LVs, Clique, Dist) :-
find_clique_from_kids(Kids, LVs, Clique, Dist).
find_clique_from_kids([K|_], LVs, Clique, Dist) :-
find_clique(K, LVs, Clique, Dist), !.
find_clique_from_kids([_|Kids], LVs, Clique, Dist) :-
find_clique_from_kids(Kids, LVs, Clique, Dist).

View File

@ -1,15 +1,24 @@
:- module(clpbn_matrix_utils, [init_CPT/2, :- module(clpbn_matrix_utils,
[init_CPT/2,
project_from_CPT/3, project_from_CPT/3,
reorder_CPT/5, reorder_CPT/5,
get_dist_size/2, get_CPT_sizes/2,
normalise_CPT/2, normalise_CPT/2,
multiply_CPTs/3, multiply_CPTs/4,
divide_CPTs/3,
expand_CPT/4,
reset_CPT_that_disagrees/5,
unit_CPT/2,
sum_out_from_CPT/4,
list_from_CPT/2]). list_from_CPT/2]).
:- use_module(dists, [get_dist_domain_size/2, :- use_module(dists,
[get_dist_domain_size/2,
get_dist_domain/2]). get_dist_domain/2]).
:- use_module(library(matrix), [matrix_new/4, :- use_module(library(matrix),
[matrix_new/4,
matrix_new_set/4,
matrix_select/4, matrix_select/4,
matrix_dims/2, matrix_dims/2,
matrix_shuffle/3, matrix_shuffle/3,
@ -18,7 +27,9 @@
matrix_dims/2, matrix_dims/2,
matrix_sum/2, matrix_sum/2,
matrix_sum_out/3, matrix_sum_out/3,
matrix_sum_out_several/3,
matrix_op_to_all/4, matrix_op_to_all/4,
matrix_set_all_that_disagree/5,
matrix_to_list/2]). matrix_to_list/2]).
:- use_module(library(lists), [nth0/3]). :- use_module(library(lists), [nth0/3]).
@ -51,21 +62,21 @@ reorder_CPT(Vs0,T0,Vs,TF,Sizes) :-
var(Vs), !, var(Vs), !,
order_vec(Vs0,Vs,Map), order_vec(Vs0,Vs,Map),
( (
Vs == V0 Vs == Vs0
-> ->
matrix_shuffle(T0,Map,TF)
;
TF = T0 TF = T0
;
matrix_shuffle(T0,Map,TF)
), ),
matrix_dims(TF, Sizes). matrix_dims(TF, Sizes).
reorder_CPT(Vs0,T0,Vs,TF,Sizes) :- reorder_CPT(Vs0,T0,Vs,TF,Sizes) :-
mapping(Vs0,Vs,Map), mapping(Vs0,Vs,Map),
( (
Vs == V0 Vs == Vs0
-> ->
matrix_shuffle(T0,Map,TF)
;
TF = T0 TF = T0
;
matrix_shuffle(T0,Map,TF)
), ),
matrix_dims(TF, Sizes). matrix_dims(TF, Sizes).
@ -98,7 +109,11 @@ split_map([], []).
split_map([_-M|Is], [M|Map]) :- split_map([_-M|Is], [M|Map]) :-
split_map(Is, Map). split_map(Is, Map).
multiply_CPTs(tab(Tab1, Deps1, Sz1), tab(Tab2, Deps2, Sz2), tab(OT, NDeps, NSz)) :- divide_CPTs(Tab1, Tab2, OT) :-
matrix_op(Tab1,Tab2,/,OT).
multiply_CPTs(tab(Tab1, Deps1, Sz1), tab(Tab2, Deps2, Sz2), tab(OT, NDeps, NSz), NTab2) :-
expand_tabs(Deps1, Sz1, Deps2, Sz2, Map1, Map2, NDeps), expand_tabs(Deps1, Sz1, Deps2, Sz2, Map1, Map2, NDeps),
matrix_expand(Tab1, Map1, NTab1), matrix_expand(Tab1, Map1, NTab1),
matrix_expand(Tab2, Map2, NTab2), matrix_expand(Tab2, Map2, NTab2),
@ -140,4 +155,39 @@ normalise_CPT(MAT,NMAT) :-
list_from_CPT(MAT, List) :- list_from_CPT(MAT, List) :-
matrix_to_list(MAT, List). matrix_to_list(MAT, List).
expand_CPT(MAT0, Dims0, DimsNew, MAT) :-
generate_map(DimsNew, Dims0, Map),
matrix_expand(MAT0, Map, MAT).
generate_map([], [], []).
generate_map([V|DimsNew], [V0|Dims0], [0|Map]) :- V == V0, !,
generate_map(DimsNew, Dims0, Map).
generate_map([V|DimsNew], Dims0, [Sz|Map]) :-
clpbn:get_atts(V, [dist(Id,_)]),
get_dist_domain_size(Id, Sz),
generate_map(DimsNew, Dims0, Map).
unit_CPT(V,CPT) :-
clpbn:get_atts(V, [dist(Id,_)]),
get_dist_domain_size(Id, Sz),
matrix_new_set(floats,[Sz],1.0,CPT).
reset_CPT_that_disagrees(CPT, Vars, V, Pos, NCPT) :-
vnth(Vars, 0, V, Dim, _),
matrix_set_all_that_disagree(CPT, Dim, Pos, 0.0, NCPT).
sum_out_from_CPT(Vs,Table,Deps,tab(NewTable,Vs,Sz)) :-
conversion_matrix(Vs, Deps, Conv),
matrix_sum_out_several(Table, Conv, NewTable),
matrix_dims(NewTable, Sz).
conversion_matrix([], [], []).
conversion_matrix([], [_|Deps], [1|Conv]) :-
conversion_matrix([], Deps, Conv).
conversion_matrix([V|Vs], [V1|Deps], [0|Conv]) :- V==V1, !,
conversion_matrix(Vs, Deps, Conv).
conversion_matrix([V|Vs], [_|Deps], [1|Conv]) :-
conversion_matrix([V|Vs], Deps, Conv).
get_CPT_sizes(CPT, Sizes) :-
matrix_dims(CPT, Sizes).

View File

@ -25,7 +25,8 @@
:- use_module(library('clpbn/graphviz'), [clpbn2gviz/4]). :- use_module(library('clpbn/graphviz'), [clpbn2gviz/4]).
:- use_module(library('clpbn/dists'), [ :- use_module(library('clpbn/dists'),
[
get_dist_domain_size/2, get_dist_domain_size/2,
get_dist_matrix/5]). get_dist_matrix/5]).
@ -39,8 +40,7 @@
:- use_module(library('clpbn/matrix_cpt_utils'), :- use_module(library('clpbn/matrix_cpt_utils'),
[project_from_CPT/3, [project_from_CPT/3,
reorder_CPT/5, reorder_CPT/5,
get_dist_size/2, multiply_CPTs/4,
multiply_CPTs/3,
normalise_CPT/2, normalise_CPT/2,
list_from_CPT/2]). list_from_CPT/2]).
@ -178,7 +178,7 @@ find_best([V|LV], V0, Threshold, VF, WorkTables, [V|LVF], Inputs) :-
multiply_tables([Table], Table) :- !. multiply_tables([Table], Table) :- !.
multiply_tables([TAB1, TAB2| Tables], Out) :- multiply_tables([TAB1, TAB2| Tables], Out) :-
multiply_CPTs(TAB1, TAB2, TAB), multiply_CPTs(TAB1, TAB2, TAB, _),
multiply_tables([TAB| Tables], Out). multiply_tables([TAB| Tables], Out).

View File

@ -17,6 +17,7 @@
<h2>Yap-5.1.3:</h2> <h2>Yap-5.1.3:</h2>
<ul> <ul>
<li> FIXED: implement JT for CLP(BN).</li>
<li> FIXED: use safe locking to ensure that dynamic predicates <li> FIXED: use safe locking to ensure that dynamic predicates
run correctly.</li> run correctly.</li>
<li> FIXED: use matrices to implement variavel elimination, also fix <li> FIXED: use matrices to implement variavel elimination, also fix

View File

@ -31,6 +31,7 @@
dgraph_min_path/5, dgraph_min_path/5,
dgraph_max_path/5, dgraph_max_path/5,
dgraph_min_paths/3, dgraph_min_paths/3,
dgraph_isomorphic/4,
dgraph_path/3]). dgraph_path/3]).
:- use_module(library(rbtrees), :- use_module(library(rbtrees),
@ -387,3 +388,24 @@ do_children([_|Children], G, SoFar, Path) :-
do_children(Children, G, SoFar, Path). do_children(Children, G, SoFar, Path).
dgraph_isomorphic(Vs, Vs2, G1, G2) :-
rb_new(Map0),
mapping(Vs,Vs2,Map0,Map),
dgraph_edges(G1,Edges),
translate_edges(Edges,Map,TEdges),
dgraph_new(G20),
dgraph_add_vertices(Vs2,G20,G21),
dgraph_add_edges(TEdges,G21,G2).
mapping([],[],Map,Map).
mapping([V1|Vs],[V2|Vs2],Map0,Map) :-
rb_insert(Map0,V1,V2,MapI),
mapping(Vs,Vs2,MapI,Map).
translate_edges([],_,[]).
translate_edges([V1-V2|Edges],Map,[NV1-NV2|TEdges]) :-
rb_lookup(V1,NV1,Map),
rb_lookup(V2,NV2,Map),
translate_edges(Edges,Map,TEdges).

View File

@ -62,6 +62,7 @@ typedef enum {
matrix_minarg/2, matrix_minarg/2,
matrix_sum/2, matrix_sum/2,
matrix_sum_out/3, matrix_sum_out/3,
matrix_sum_out_several/3,
matrix_add_to_all/2, matrix_add_to_all/2,
matrix_agg_lines/3, matrix_agg_lines/3,
matrix_agg_cols/3, matrix_agg_cols/3,
@ -71,6 +72,7 @@ typedef enum {
matrix_op_to_cols/4, matrix_op_to_cols/4,
matrix_shuffle/3, matrix_shuffle/3,
matrix_transpose/2, matrix_transpose/2,
matrix_set_all_that_disagree/5,
matrix_expand/3, matrix_expand/3,
matrix_select/4 matrix_select/4
]). ]).
@ -121,6 +123,8 @@ matrix_op(M1,M2,-,NM) :-
do_matrix_op(M1,M2,1,NM). do_matrix_op(M1,M2,1,NM).
matrix_op(M1,M2,*,NM) :- matrix_op(M1,M2,*,NM) :-
do_matrix_op(M1,M2,2,NM). do_matrix_op(M1,M2,2,NM).
matrix_op(M1,M2,/,NM) :-
do_matrix_op(M1,M2,3,NM).
matrix_op_to_all(M1,+,Num,NM) :- matrix_op_to_all(M1,+,Num,NM) :-
do_matrix_op_to_all(M1,0,Num,NM). do_matrix_op_to_all(M1,0,Num,NM).

View File

@ -104,6 +104,20 @@ matrix_get_index(int *mat, unsigned int offset, int* indx)
} }
} }
static void
matrix_next_index(int *dims, int ndims, int* indx)
{
unsigned int i;
/* find where we are */
for (i = ndims; i >0; ) {
i--;
indx[i]++;
if (indx[i]!=dims[i]) return;
indx[i] = 0;
}
}
static YAP_Term static YAP_Term
new_int_matrix(int ndims, int dims[], long int data[]) new_int_matrix(int ndims, int dims[], long int data[])
{ {
@ -1438,6 +1452,46 @@ matrix_double_mult_data(double *nmat, int siz, double mat1[], double mat2[])
} }
} }
static void
matrix_long_div_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_div_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_div2_data(double *nmat, int siz, double mat1[], long int mat2[])
{
int i;
for (i=0; i< siz; i++) {
nmat[i] = mat1[i]/mat2[i];
}
}
static void
matrix_double_div_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 static int
matrix_op(void) matrix_op(void)
{ {
@ -1492,6 +1546,9 @@ matrix_op(void)
case MAT_TIMES: case MAT_TIMES:
matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2); matrix_long_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break; break;
case MAT_DIV:
matrix_long_div_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default: default:
return FALSE; return FALSE;
} }
@ -1521,6 +1578,9 @@ matrix_op(void)
case MAT_TIMES: case MAT_TIMES:
matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2); matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break; break;
case MAT_DIV:
matrix_long_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default: default:
return FALSE; return FALSE;
} }
@ -1559,6 +1619,9 @@ matrix_op(void)
case MAT_TIMES: case MAT_TIMES:
matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1); matrix_long_double_mult_data(ndata, mat1[MAT_SIZE], data2, data1);
break; break;
case MAT_DIV:
matrix_long_double_div2_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default: default:
return FALSE; return FALSE;
} }
@ -1588,6 +1651,9 @@ matrix_op(void)
case MAT_TIMES: case MAT_TIMES:
matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2); matrix_double_mult_data(ndata, mat1[MAT_SIZE], data1, data2);
break; break;
case MAT_DIV:
matrix_double_div_data(ndata, mat1[MAT_SIZE], data1, data2);
break;
default: default:
return FALSE; return FALSE;
} }
@ -2120,6 +2186,99 @@ matrix_sum_out(void)
return YAP_Unify(YAP_ARG3, tf); return YAP_Unify(YAP_ARG3, tf);
} }
/* given a matrix M and a set of dims, sum out one of the dimensions
*/
static int
matrix_sum_out_several(void)
{
int ndims, i, *dims, newdims;
int indx[MAX_DIMS], nindx[MAX_DIMS], conv[MAX_DIMS];
YAP_Term tf, tconv;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
if (!mat) {
/* Error */
return FALSE;
}
ndims = mat[MAT_NDIMS];
dims = mat+MAT_DIMS;
/* we now have our target matrix, let us grab our conversion arguments */
tconv = YAP_ARG2;
for (i=0, newdims=0; i < ndims; i++) {
YAP_Term th;
if (!YAP_IsPairTerm(tconv))
return FALSE;
th = YAP_HeadOfTerm(tconv);
if (!YAP_IsIntTerm(th))
return FALSE;
conv[i] = YAP_IntOfTerm(th);
if (!conv[i]) {
nindx[newdims++] = dims[i];
}
tconv = YAP_TailOfTerm(tconv);
}
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 (!conv[j]) {
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 (!conv[j]) {
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 /* given a matrix M and a set of dims, build contract a matrix to follow
the new order the new order
*/ */
@ -2203,13 +2362,15 @@ matrix_expand(void)
data = matrix_double_data(mat,ndims); data = matrix_double_data(mat,ndims);
ndata = matrix_double_data(nmat,newdims); ndata = matrix_double_data(nmat,newdims);
/* create a new matrix with the same size */ /* create a new matrix with the same size */
for (i=0; i < newdims; i++)
indx[i] = 0;
for (i=0; i< nmat[MAT_SIZE]; i++) { for (i=0; i< nmat[MAT_SIZE]; i++) {
int j,k=0; int j,k=0;
/* /*
not very efficient, we could try to take advantage of the fact not very efficient, we could try to take advantage of the fact
that we usually only change an index at a time that we usually only change an index at a time
*/ */
matrix_get_index(nmat, i, indx); matrix_next_index(nmat+MAT_DIMS, newdims, indx);
for (j = 0; j < newdims; j++) { for (j = 0; j < newdims; j++) {
if (!new[j]) if (!new[j])
nindx[k++] = indx[j]; nindx[k++] = indx[j];
@ -2220,6 +2381,88 @@ matrix_expand(void)
return YAP_Unify(YAP_ARG3, tf); 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_set_all_that_disagree(void)
{
int ndims, i, *dims;
int indx[MAX_DIMS];
YAP_Term tf;
int *mat = (int *)YAP_BlobOfTerm(YAP_ARG1), *nmat;
int dim = YAP_IntOfTerm(YAP_ARG2);
int pos = YAP_IntOfTerm(YAP_ARG3);
if (!mat) {
/* Error */
return FALSE;
}
ndims = mat[MAT_NDIMS];
dims = mat+MAT_DIMS;
if (mat[MAT_TYPE] == INT_MATRIX) {
long int *data, *ndata, val;
/* create a new matrix with the same size */
tf = new_int_matrix(ndims,dims,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,ndims);
if (!YAP_IsIntTerm(YAP_ARG4))
return FALSE;
val = YAP_IntOfTerm(YAP_ARG4);
/* create a new matrix with the same size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
/*
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);
if (indx[dim] != pos)
ndata[i] = val;
else
ndata[i] = data[i];
}
} else {
double *data, *ndata, val;
/* create a new matrix with the same size */
tf = new_float_matrix(ndims,dims,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,ndims);
if (YAP_IsFloatTerm(YAP_ARG4))
val = YAP_FloatOfTerm(YAP_ARG4);
else if (YAP_IsIntTerm(YAP_ARG4))
val = YAP_IntOfTerm(YAP_ARG4);
else
return FALSE;
/* create a new matrix with the same size */
for (i=0; i< nmat[MAT_SIZE]; i++) {
/*
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);
if (indx[dim] != pos)
ndata[i] = val;
else
ndata[i] = data[i];
}
}
return YAP_Unify(YAP_ARG5, tf);
}
void PROTO(init_matrix, (void)); void PROTO(init_matrix, (void));
void void
@ -2254,6 +2497,8 @@ init_matrix(void)
YAP_UserCPredicate("matrix_select", matrix_select, 4); YAP_UserCPredicate("matrix_select", matrix_select, 4);
YAP_UserCPredicate("matrix_add_to_all", matrix_sum, 2); YAP_UserCPredicate("matrix_add_to_all", matrix_sum, 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_set_all_that_disagree", matrix_set_all_that_disagree, 5);
YAP_UserCPredicate("do_matrix_op", matrix_op, 4); YAP_UserCPredicate("do_matrix_op", matrix_op, 4);
YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3); YAP_UserCPredicate("do_matrix_agg_lines", matrix_agg_lines, 3);
YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3); YAP_UserCPredicate("do_matrix_agg_cols", matrix_agg_cols, 3);

View File

@ -41,9 +41,6 @@
ord_memberchk/2 % Element X Set ord_memberchk/2 % Element X Set
]). ]).
:- use_module(library(lists),
[memberchk/2]).
/* /*
:- mode :- mode
list_to_ord_set(+, ?), list_to_ord_set(+, ?),
@ -354,6 +351,7 @@ ord_union_all(N,Sets0,Union,Sets) :-
ord_empty([]). ord_empty([]).
ord_memberchk(Element, Set) :- ord_memberchk(Element, [E|_]) :- E == Element, !.
memberchk(Element, Set). ord_memberchk(Element, [_|Set]) :-
ord_memberchk(Element, Set).

View File

@ -467,7 +467,7 @@ del_max(black(L,K0,V0,R), K, V, Nil, NT, Flag) :-
delete_red_node(L,L,L,done) :- !. delete_red_node(L1,L2,L1,done) :- L1 == L2, !.
delete_red_node(black([],[],[],[]),R,R,done) :- !. delete_red_node(black([],[],[],[]),R,R,done) :- !.
delete_red_node(L,black([],[],[],[]),L,done) :- !. delete_red_node(L,black([],[],[],[]),L,done) :- !.
delete_red_node(L,R,OUT,Done) :- delete_red_node(L,R,OUT,Done) :-
@ -475,7 +475,7 @@ delete_red_node(L,R,OUT,Done) :-
fixup_right(Done0,red(L,NK,NV,NR),OUT,Done). fixup_right(Done0,red(L,NK,NV,NR),OUT,Done).
delete_black_node(L,L,L,not_done) :- !. delete_black_node(L1,L2,L1,not_done) :- L1 == L2, !.
delete_black_node(black([],[],[],[]),red(L,K,V,R),black(L,K,V,R),done) :- !. delete_black_node(black([],[],[],[]),red(L,K,V,R),black(L,K,V,R),done) :- !.
delete_black_node(black([],[],[],[]),R,R,not_done) :- !. delete_black_node(black([],[],[],[]),R,R,not_done) :- !.
delete_black_node(red(L,K,V,R),black([],[],[],[]),black(L,K,V,R),done) :- !. delete_black_node(red(L,K,V,R),black([],[],[],[]),black(L,K,V,R),done) :- !.
@ -543,6 +543,7 @@ fixup_right(not_done,T,NT,Done) :-
fixup3(T,NT,Done). fixup3(T,NT,Done).
% %
% case 1: x moves down, so we have to try to fix it again. % case 1: x moves down, so we have to try to fix it again.
% case 1 -> 2,3,4 -> done % case 1 -> 2,3,4 -> done
@ -659,33 +660,37 @@ partial_map(red(L,K,V,R),Map,MapF,Nil,Goal,red(NL,K,NV,NR)) :-
partial_map(L,Map,MapI,Nil,Goal,NL), partial_map(L,Map,MapI,Nil,Goal,NL),
( (
MapI == [] -> MapI == [] ->
NR = R, NV = V NR = R, NV = V, MapF = []
; ;
MapI = [K1|MapR], MapI = [K1|MapR],
( (
K == K1 -> K == K1
once(call(Goal,V,NV)), ->
Map2 = MapR ( call(Goal,V,NV) -> true ; NV = V ),
MapN = MapR
; ;
Map2 = MapI, NV = V NV = V,
MapN = MapI
), ),
partial_map(R,Map2,MapF,Nil,Goal,NR) partial_map(R,MapN,MapF,Nil,Goal,NR)
). ).
partial_map(black(L,K,V,R),Map,MapF,Nil,Goal,black(NL,K,NV,NR)) :- partial_map(black(L,K,V,R),Map,MapF,Nil,Goal,black(NL,K,NV,NR)) :-
partial_map(L,Map,MapI,Nil,Goal,NL), partial_map(L,Map,MapI,Nil,Goal,NL),
( (
MapI == [] -> MapI == [] ->
NR = R, NV = V NR = R, NV = V, MapF = []
; ;
MapI = [K1|MapR], MapI = [K1|MapR],
( (
K == K1 -> K == K1
once(call(Goal,V,NV)), ->
Map2 = MapR ( call(Goal,V,NV) -> true ; NV = V ),
MapN = MapR
; ;
Map2 = MapI, NV = V NV = V,
MapN = MapI
), ),
partial_map(R,Map2,MapF,Nil,Goal,NR) partial_map(R,MapN,MapF,Nil,Goal,NR)
). ).
@ -706,6 +711,22 @@ keys(black(L,K,_,R),L0,Lf) :-
keys(L,[K|L1],Lf), keys(L,[K|L1],Lf),
keys(R,L0,L1). keys(R,L0,L1).
ord_list_to_rbtree(List,Tree) :-
list_to_rbtree(List,Tree).
list_to_rbtree(List,Tree) :-
rb_new(T0),
list_to_rbtree(List,T0,Tree).
list_to_rbtree([],Tree,Tree).
list_to_rbtree([K-V|List],T0,Tree) :-
rb_insert(T0, K, V, T1),
list_to_rbtree(List,T1,Tree).
/*
list_to_rbtree(List,t(Nil,Tree)) :- list_to_rbtree(List,t(Nil,Tree)) :-
Nil = black([], [], [], []), Nil = black([], [], [], []),
sort(List,Sorted), sort(List,Sorted),
@ -718,6 +739,7 @@ ord_list_to_rbtree(List,t(Nil,Tree)) :-
Ar =.. [seq|List], Ar =.. [seq|List],
functor(Ar,_,L), functor(Ar,_,L),
construct_rbtree(1, L, Ar, black, Nil, Tree). construct_rbtree(1, L, Ar, black, Nil, Tree).
*/
construct_rbtree(L, M, _, _, Nil, Nil) :- M < L, !. construct_rbtree(L, M, _, _, Nil, Nil) :- M < L, !.
construct_rbtree(L, L, Ar, Color, Nil, Node) :- !, construct_rbtree(L, L, Ar, Color, Nil, Node) :- !,
@ -762,8 +784,8 @@ rbtree(T) :-
rbtree1(black(L,K,_,R)) :- rbtree1(black(L,K,_,R)) :-
find_path_blacks(L, 0, Bls), find_path_blacks(L, 0, Bls),
check_rbtree(L,-1000000,K,Bls), check_rbtree(L,-inf,K,Bls),
check_rbtree(R,K,1000000,Bls). check_rbtree(R,K,+inf,Bls).
rbtree1(red(_,_,_,_)) :- rbtree1(red(_,_,_,_)) :-
throw(msg("root should be black",[])). throw(msg("root should be black",[])).
@ -793,7 +815,7 @@ check_height(0,_,_) :- !.
check_height(Bls0,Min,Max) :- check_height(Bls0,Min,Max) :-
throw(msg("Unbalance ~d between ~w and ~w~n",[Bls0,Min,Max])). throw(msg("Unbalance ~d between ~w and ~w~n",[Bls0,Min,Max])).
check_val(K, Min, Max) :- K > Min, K < Max, !. check_val(K, Min, Max) :- ( K @> Min ; Min == -inf), (K @< Max ; Max == +inf), !.
check_val(K, Min, Max) :- check_val(K, Min, Max) :-
throw(msg("not ordered: ~w not between ~w and ~w~n",[K,Min,Max])). throw(msg("not ordered: ~w not between ~w and ~w~n",[K,Min,Max])).

View File

@ -101,16 +101,16 @@ edges2wgraphl([V|Vertices], SortedEdges, [V-[]|GraphL]) :-
wdgraph_add_edges([],[]) --> []. wdgraph_add_edges([],[]) --> [].
wdgraph_add_edges([V|Vs],[V-(V1-W)|Es]) --> !, wdgraph_add_edges([VA|Vs],[VB-(V1-W)|Es]) --> { VA == VB }, !,
{ get_extra_children(Es,V,Children,REs) }, { get_extra_children(Es,VA,Children,REs) },
wdgraph_update_vertex(V,[V1-W|Children]), wdgraph_update_vertex(VA,[V1-W|Children]),
wdgraph_add_edges(Vs,REs). wdgraph_add_edges(Vs,REs).
wdgraph_add_edges([V|Vs],Es) --> !, wdgraph_add_edges([V|Vs],Es) --> !,
wdgraph_update_vertex(V,[]), wdgraph_update_vertex(V,[]),
wdgraph_add_edges(Vs,Es). wdgraph_add_edges(Vs,Es).
get_extra_children([V-(C-W)|Es],V,[C-W|Children],REs) :- !, get_extra_children([VA-(C-W)|Es],VB,[C-W|Children],REs) :- VA == VB, !,
get_extra_children(Es,V,Children,REs). get_extra_children(Es,VB,Children,REs).
get_extra_children(Es,_,[],Es). get_extra_children(Es,_,[],Es).
@ -120,9 +120,9 @@ wdgraph_update_vertex(V,Edges,WG0,WGF) :-
wdgraph_update_vertex(V,Edges,WG0,WGF) :- wdgraph_update_vertex(V,Edges,WG0,WGF) :-
rb_insert(WG0, V, Edges, WGF). rb_insert(WG0, V, Edges, WGF).
key_union([], [], []). key_union([], [], []) :- !.
key_union([], [C|Children], [C|Children]). key_union([], [C|Children], [C|Children]).
key_union([C|Children], [], [C|Children]). key_union([C|Children], [], [C|Children]) :- !.
key_union([K-W|ToAdd], [K1-W1|Children0], NewUnion) :- key_union([K-W|ToAdd], [K1-W1|Children0], NewUnion) :-
( K == K1 -> ( K == K1 ->
NewUnion = [K-W|NewChildren], NewUnion = [K-W|NewChildren],

View File

@ -219,6 +219,14 @@ wundgraph_to_undgraph(G1, G2) :-
wundgraph_min_tree(G, T, C) :- wundgraph_min_tree(G, T, C) :-
rb_visit(G, Els0), rb_visit(G, Els0),
generate_min_tree(Els0, T, C).
generate_min_tree([], T, 0) :- !,
wundgraph_new(T).
generate_min_tree([El-_], T, 0) :- !,
wundgraph_new(T0),
wundgraph_add_vertex(El,T0,T).
generate_min_tree(Els0, T, C) :-
mk_list_of_edges(Els0, Edges), mk_list_of_edges(Els0, Edges),
keysort(Edges, SortedEdges), keysort(Edges, SortedEdges),
rb_new(V0), rb_new(V0),
@ -228,6 +236,14 @@ wundgraph_min_tree(G, T, C) :-
wundgraph_max_tree(G, T, C) :- wundgraph_max_tree(G, T, C) :-
rb_visit(G, Els0), rb_visit(G, Els0),
generate_max_tree(Els0, T, C).
generate_max_tree([], T, 0) :- !,
wundgraph_new(T).
generate_max_tree([El-_], T, 0) :- !,
wundgraph_new(T0),
wundgraph_add_vertex(El,T0,T).
generate_max_tree(Els0, T, C) :-
mk_list_of_edges(Els0, Edges), mk_list_of_edges(Els0, Edges),
keysort(Edges, SortedEdges), keysort(Edges, SortedEdges),
reverse(SortedEdges, ReversedEdges), reverse(SortedEdges, ReversedEdges),