Merge branch 'master' of ssh://yap.git.sourceforge.net/gitroot/yap/yap-6.3

C/c_interface.c

@@ -334,6 +334,7 @@
 #include "clause.h"
 #include "yapio.h"
 #include "attvar.h"
+#include "SWI-Stream.h"
 #if HAVE_STDARG_H
 #include <stdarg.h>
 #endif
@@ -458,6 +459,7 @@ X_API int     STD_PROTO(YAP_GoalHasException,(Term *));
 X_API void    STD_PROTO(YAP_ClearExceptions,(void));
 X_API int     STD_PROTO(YAP_ContinueGoal,(void));
 X_API void    STD_PROTO(YAP_PruneGoal,(void));
+X_API IOSTREAM   *STD_PROTO(YAP_TermToStream,(Term));
 X_API IOSTREAM   *STD_PROTO(YAP_InitConsult,(int, char *));
 X_API void    STD_PROTO(YAP_EndConsult,(IOSTREAM *));
 X_API Term    STD_PROTO(YAP_Read, (IOSTREAM *));
@@ -2488,6 +2490,22 @@ YAP_InitConsult(int mode, char *filename)
   return st;
 }
 
+X_API IOSTREAM *
+YAP_TermToStream(Term t)
+{
+  CACHE_REGS
+  IOSTREAM *s;
+  int rc;
+  BACKUP_MACHINE_REGS();
+
+  if ( (rc=PL_get_stream_handle(Yap_InitSlot(t PASS_REGS), &s)) ) {
+    RECOVER_MACHINE_REGS();
+    return s;
+  }
+  RECOVER_MACHINE_REGS();
+  return NULL;
+}
+
 X_API void
 YAP_EndConsult(IOSTREAM *s)
 {

C/scanner.c

@@ -854,6 +854,15 @@ Yap_tokenizer(IOSTREAM *inp_stream, Term *tposp)
 	  while (chtype(ch) == BS) {
 	    ch = getchr(inp_stream);
 	  }
+	  if (ASP-H < 1024) {
+	    Yap_ErrorMessage = "Stack Overflow";
+	    Yap_Error_TYPE = OUT_OF_STACK_ERROR;	            
+	    Yap_Error_Size = 0L;	            
+	    if (p)
+	      p->Tok = Ord(kind = eot_tok);
+	    /* serious error now */
+	    return l;
+	  }
 	  *tposp = Yap_StreamPosition(inp_stream);
 	}
 	goto restart;
@@ -1161,6 +1170,15 @@ Yap_tokenizer(IOSTREAM *inp_stream, Term *tposp)
 	  while (chtype(ch) == BS) {
 	    ch = getchr(inp_stream);
 	  }
+	  if (ASP-H < 1024) {
+	    Yap_ErrorMessage = "Stack Overflow";
+	    Yap_Error_TYPE = OUT_OF_STACK_ERROR;	            
+	    Yap_Error_Size = 0L;	            
+	    if (p)
+	      p->Tok = Ord(kind = eot_tok);
+	    /* serious error now */
+	    return l;
+	  }
 	  *tposp = Yap_StreamPosition(inp_stream);
 	}
 	goto restart;

console/LGPL/pl-ntmain.c

@@ -1035,6 +1035,7 @@ win32main(rlc_console c, int argc, TCHAR **argv)
   if ( !PL_initialise(argc, av) )
     PL_halt(1);
 
+  rlc_bind_terminal(c);
   PL_halt(PL_toplevel() ? 0 : 1);
 
   return 0;
@@ -1050,7 +1051,6 @@ int PASCAL
 WinMain(HINSTANCE hInstance, HINSTANCE hPrevInstance,
 	LPSTR lpszCmdLine, int nCmdShow)
 { LPTSTR cmdline;
-  fprintf(stderr,"Hello\n");
 
   InitializeCriticalSection(&mutex);
 

include/YapInterface.h

@@ -325,6 +325,9 @@ extern X_API int PROTO(YAP_Init,(YAP_init_args *));
 /*  int YAP_FastInit(const char *) */
 extern X_API int PROTO(YAP_FastInit,(CONST char *));
 
+/*  void * YAP_TermToStream(YAP_Term) */
+extern X_API void * PROTO(YAP_TermToStream,(YAP_Term));
+
 /*  void * YAP_InitConsult(int, const char *) */
 extern X_API void * PROTO(YAP_InitConsult,(int, CONST char *));
 

packages/CLPBN/Makefile.in

@@ -51,7 +51,7 @@ CLPBN_PROGRAMS= \
 	$(CLPBN_SRCDIR)/table.yap \
 	$(CLPBN_SRCDIR)/topsort.yap \
 	$(CLPBN_SRCDIR)/utils.yap \
-	$(CLPBN_SRCDIR)/vel.yap \
+	$(CLPBN_SRCDIR)/ve.yap \
 	$(CLPBN_SRCDIR)/viterbi.yap \
 	$(CLPBN_SRCDIR)/xbif.yap
 

packages/CLPBN/clpbn.yap

@@ -7,6 +7,7 @@
 		  clpbn_key/2,
 		  clpbn_init_solver/4,
 		  clpbn_run_solver/3,
+		  clpbn_finalize_solver/1,
 		  clpbn_init_solver/5,
 		  clpbn_run_solver/4,
 		  clpbn_init_graph/1,
@@ -28,13 +29,22 @@
 :- attribute key/1, dist/2, evidence/1, starter/0.
 
 
-:- use_module('clpbn/vel',
+:- use_module('clpbn/ve',
-	      [vel/3,
+	      [ve/3,
-	       check_if_vel_done/1,
+	       check_if_ve_done/1,
-	       init_vel_solver/4,
+	       init_ve_solver/4,
-	       run_vel_solver/3
+	       run_ve_solver/3
 	      ]).
 
+:- use_module('clpbn/bp',
+	      [bp/3,
+               check_if_bp_done/1,
+	       init_bp_solver/4,
+	       run_bp_solver/3,
+	       finalize_bp_solver/1
+	      ]).
+
+
 :- use_module('clpbn/jt',
 	      [jt/3,
 	       init_jt_solver/4,
@@ -53,6 +63,14 @@
 	       run_gibbs_solver/3
 	      ]).
 
+:- use_module('clpbn/bp',
+	      [bp/3,
+	       check_if_bp_done/1,
+	       init_bp_solver/4,
+	       run_bp_solver/3,
+	       finalize_bp_solver/1
+	      ]).
+
 :- use_module('clpbn/pgrammar',
 	      [init_pcg_solver/4,
 	       run_pcg_solver/3,
@@ -92,8 +110,8 @@
 
 :- dynamic solver/1,output/1,use/1,suppress_attribute_display/1, parameter_softening/1, em_solver/1.
 
-solver(vel).
+solver(ve).
-em_solver(vel).
+em_solver(ve).
 
 %output(xbif(user_error)).
 %output(gviz(user_error)).
@@ -142,6 +160,18 @@ clpbn_flag(parameter_softening,Before,After) :-
 %	,writeln({Var = Key with Dist})
 .
 
+%
+% make sure a query variable is reachable by the garbage collector.
+%
+store_var(El) :- 
+	catch(b_getval(clpbn_qvars,Q.Tail), _, init_clpbn_vars(El, Q, Tail)),
+	Tail = [El|NewTail],
+	b_setval(clpbn_qvars, [Q|NewTail]).
+	
+init_clpbn_vars(El, Q, Tail) :-
+	Q = [El|Tail],
+	b_setval(clpbn_qvars, [Q|Tail]).
+
 check_constraint(Constraint, _, _, Constraint) :- var(Constraint), !.
 check_constraint((A->D), _, _, (A->D)) :- var(A), !.
 check_constraint((([A|B].L)->D), Vars, NVars, (([A|B].NL)->D)) :- !,
@@ -162,8 +192,10 @@ add_evidence(V,Key,Distinfo,NV) :-
 	nonvar(V), !,
 	get_evidence_position(V, Distinfo, Pos),
 	check_stored_evidence(Key, Pos),
+	store_var(NV),
 	clpbn:put_atts(NV,evidence(Pos)).
 add_evidence(V,K,_,V) :-
+	store_var(V),
 	add_evidence(K,V).
 
 clpbn_marginalise(V, Dist) :-
@@ -183,7 +215,7 @@ project_attributes(GVars, AVars) :-
 	clpbn_vars(AVars, DiffVars, AllVars),
 	get_clpbn_vars(GVars,CLPBNGVars0),
 	simplify_query_vars(CLPBNGVars0, CLPBNGVars),
-	(output(xbif(XBifStream)) -> clpbn2xbif(XBifStream,vel,AllVars) ; true),
+	(output(xbif(XBifStream)) -> clpbn2xbif(XBifStream,ve,AllVars) ; true),
 	(output(gviz(XBifStream)) -> clpbn2gviz(XBifStream,sort,AllVars,GVars) ; true),
 	(
 	    Solver = graphs
@@ -225,10 +257,12 @@ get_rid_of_ev_vars([V|LVs0],[V|LVs]) :-
 % do nothing if we don't have query variables to compute.
 write_out(graphs, _, AVars, _) :-
 	clpbn2graph(AVars).
-write_out(vel, GVars, AVars, DiffVars) :-
+write_out(ve, GVars, AVars, DiffVars) :-
-	vel(GVars, AVars, DiffVars).
+	ve(GVars, AVars, DiffVars).
 write_out(jt, GVars, AVars, DiffVars) :-
 	jt(GVars, AVars, DiffVars).
+write_out(bp, GVars, AVars, DiffVars) :-
+	bp(GVars, AVars, DiffVars).
 write_out(gibbs, GVars, AVars, DiffVars) :-
 	gibbs(GVars, AVars, DiffVars).
 write_out(bnt, GVars, AVars, DiffVars) :-
@@ -315,11 +349,14 @@ bind_clpbn(_, Var, _, _, _, _, []) :-
 	use(bnt),
 	check_if_bnt_done(Var), !.
 bind_clpbn(_, Var, _, _, _, _, []) :-
-	use(vel),
+	use(ve),
-	check_if_vel_done(Var), !.
+	check_if_ve_done(Var), !.
+bind_clpbn(_, Var, _, _, _, _, []) :-
+	use(bp),
+	check_if_bp_done(Var), !.
 bind_clpbn(_, Var, _, _, _, _, []) :-
 	use(jt),
-	check_if_vel_done(Var), !.
+	check_if_ve_done(Var), !.
 bind_clpbn(T, Var, Key0, _, _, _, []) :-
 	get_atts(Var, [key(Key)]), !,
 	(
@@ -384,13 +421,16 @@ clpbn_init_solver(LVs, Vs0, VarsWithUnboundKeys, State) :-
 
 clpbn_init_solver(gibbs, LVs, Vs0, VarsWithUnboundKeys, State) :-
 	init_gibbs_solver(LVs, Vs0, VarsWithUnboundKeys, State).
-clpbn_init_solver(vel, LVs, Vs0, VarsWithUnboundKeys, State) :-
+clpbn_init_solver(ve, LVs, Vs0, VarsWithUnboundKeys, State) :-
-	init_vel_solver(LVs, Vs0, VarsWithUnboundKeys, State).
+	init_ve_solver(LVs, Vs0, VarsWithUnboundKeys, State).
+clpbn_init_solver(bp, LVs, Vs0, VarsWithUnboundKeys, State) :-
+	init_bp_solver(LVs, Vs0, VarsWithUnboundKeys, State).
 clpbn_init_solver(jt, LVs, Vs0, VarsWithUnboundKeys, State) :-
 	init_jt_solver(LVs, Vs0, VarsWithUnboundKeys, State).
 clpbn_init_solver(pcg, LVs, Vs0, VarsWithUnboundKeys, State) :-
 	init_pcg_solver(LVs, Vs0, VarsWithUnboundKeys, State).
 
+
 %
 % LVs is the list of lists of variables to marginalise
 % Vs is the full graph
@@ -403,10 +443,16 @@ clpbn_run_solver(LVs, LPs, State) :-
 
 clpbn_run_solver(gibbs, LVs, LPs, State) :-
 	run_gibbs_solver(LVs, LPs, State).
-clpbn_run_solver(vel, LVs, LPs, State) :-
+
-	run_vel_solver(LVs, LPs, State).
+clpbn_run_solver(ve, LVs, LPs, State) :-
+	run_ve_solver(LVs, LPs, State).
+
+clpbn_run_solver(bp, LVs, LPs, State) :-
+	run_bp_solver(LVs, LPs, State).
+
 clpbn_run_solver(jt, LVs, LPs, State) :-
 	run_jt_solver(LVs, LPs, State).
+
 clpbn_run_solver(pcg, LVs, LPs, State) :-
 	run_pcg_solver(LVs, LPs, State).
 
@@ -415,3 +461,10 @@ add_keys(Key1+V1,_Key2,Key1+V1).
 clpbn_init_graph(pcg) :- !,
 	pcg_init_graph.
 clpbn_init_graph(_).
+
+clpbn_finalize_solver(State) :-
+	solver(bp), !,
+	functor(State, _, Last),
+	arg(Last, State, Info),
+	finalize_bp_solver(Info).
+clpbn_finalize_solver(_State).

packages/CLPBN/clpbn/aggregates.yap

@@ -52,7 +52,7 @@ cpt_average(AllVars, Key, Els0, Tab, Vs, NewVs) :-
 	cpt_average(AllVars, Key, Els0, 1.0, Tab, Vs, NewVs).
 
 % support variables with evidence from domain. This should make everyone's life easier.
-cpt_average([Ev|Vars], Key, Els0, Softness, p(Els0, CPT, NewParents), Vs, NewVs) :-
+cpt_average([Ev|Vars], Key, Els0, Softness, pf(Els0, MAT, NewParents), Vs, NewVs) :-
 	find_evidence(Vars, 0, TotEvidence, RVars),
 	build_avg_table(RVars, Vars, Els0, Key, TotEvidence, Softness, MAT0, NewParents0, Vs, IVs),
 	include_qevidence(Ev, MAT0, MAT, NewParents0, NewParents, Vs, IVs, NewVs).

packages/CLPBN/clpbn/bp.yap

@@ -1,157 +1,152 @@
 
-/***********************************
+/************************************************
 
   Belief Propagation in CLP(BN)
  
-  This should connect to C-code.
+**************************************************/
 
-*********************************/
+:- module(clpbn_bp,
-
+          [bp/3,
-:- module(clpbn_bp, [
-          bp/3,
            check_if_bp_done/1,
            init_bp_solver/4,
-          run_bp_solver/3]).
+           run_bp_solver/3,
-
+           finalize_bp_solver/1
-
-:- use_module(library('clpbn/aggregates'),
-          [check_for_agg_vars/2]).
-
-:- use_module(library('clpbn/connected'),
-          [init_influences/3,
-           influences/5
           ]).
 
+
 :- use_module(library('clpbn/dists'),
           [dist/4,
            get_dist_domain/2,
+           get_dist_domain_size/2,
            get_dist_params/2
           ]).
 
+
 :- use_module(library('clpbn/display'),
          [clpbn_bind_vals/3]).
 
-:-use_module(library(lists),
+
-          [append/3,
+:- use_module(library(atts)).
-           memberchk/2
+
-          ]).
+:- use_module(library(charsio)).
 
 :- load_foreign_files(['horus'], [], init_predicates).
 
-:- attribute all_diffs/1.
+:- attribute id/1.
+
+:- dynamic num_bayes_nets/1.
 
 
 check_if_bp_done(_Var).
 
-%
+num_bayes_nets(0).
-% implementation of belief propagation
+
-%
+
-% A1 = +QueryVars -> sets of independent marginalization variables
-% A2 = *AllVars   -> list of all variables
-% A3 = -Output    -> output probabilities
-%
-% Other important variables:
-%
-% State0 initialized graph, is used to pass data from initialization 
-% to query solving (eg, State might be the JT and be used to run
-% different queries). 
-%
 bp([[]],_,_) :- !.
-bp([QueryVars],AllVars,Output) :-
+bp([QueryVars], AllVars, Output) :-
-    writeln(queryVars:QueryVars),
+    init_bp_solver(_, AllVars, _, BayesNet),
-    writeln(allVars:AllVars),
+    run_bp_solver([QueryVars], LPs, BayesNet),
-    % init_bp_solver([QueryVars], AllVars, Output, State),
+    finalize_bp_solver(BayesNet),
-    run_bp_solver([QueryVars], [AllVars], _State),
+    clpbn_bind_vals([QueryVars], LPs, Output).
-    % bind probs back to variables so that they can be output.
-    clpbn_bind_vals([QueryVars],[LPs],Output).
-
-% initialise necessary data for query solver
-init_bp_solver(Qs, AllVars, _, graph(LVis)) :-
-    % replace average, max, min and friends by binary nodes.
-    check_for_agg_vars(AllVars, UnFoldedVars),
-    % replace the variables reachable from G
-    init_influences(UnfoldedVars, G, RG),
-    init_bp_solver_for_questions(Qs, G, RG, _, LVis).
-
-init_bp_solver_for_questions([], _, _, [], []).
-init_bp_solver_for_questions([Vs|MVs], G, RG, [NVs|MNVs0], [NVs|LVis]) :-
-    % find variables connectd to Vs
-    influences(Vs, _, NVs0, G, RG),
-    sort(NVs0, NVs),
-    init_bp_solver_for_questions(MVs, G, RG, MNVs0, LVis).
 
 
-% use a findall to recover space without needing for GC
+init_bp_solver(_, AllVars, _, (BayesNet, DistIds)) :-
-run_bp_solver(LVs, LPs, _) :-
+    %inc_num_bayes_nets,
-    findall(Ps, solve_bp(LVs, LPs, Ps), LPs).
+    %(showprofres(50) -> true ; true),
+    process_ids(AllVars, 0, DistIds0),
+    get_vars_info(AllVars, VarsInfo),
+    sort(DistIds0, DistIds),
+    %(num_bayes_nets(0) -> writeln(vars:VarsInfo) ; true),
+    %(num_bayes_nets(0) -> writeln(dists:DistsInfo) ; true),
+    create_network(VarsInfo, BayesNet).
+    %get_extra_vars_info(AllVars, ExtraVarsInfo),
+    %set_extra_vars_info(BayesNet, ExtraVarsInfo).
     
 
-solve_bp([LVs|_], [NVs0|_], Ps) :-
+process_ids([], _, []).
-    get_vars_info(NVs0, LVi),
+process_ids([V|Vs], VarId0, [DistId|DistIds]) :-
-    get_dists_info(NVs0, Dists),
+    clpbn:get_atts(V, [dist(DistId, _)]), !,
-    process(LVi, Dists, LVs, Ps).
+    put_atts(V, [id(VarId0)]),
-solve_bp([_|MoreLVs], [_|MoreLVis], Ps) :-
+    VarId is VarId0 + 1,
-    solve_bp(MoreLVs, MoreLVis, Ps).
+    process_ids(Vs, VarId, DistIds).
+process_ids([_|Vs], VarId, DistIds) :-
+    process_ids(Vs, VarId, DistIds).
 
 
 get_vars_info([], []).
-get_vars_info([V|Vs], [var(V, Id, Parents, NParents, Ev)|LV]) :-
+get_vars_info([V|Vs], [var(VarId, DSize, Ev, ParentIds, DistId)|VarsInfo]) :-
-    clpbn:get_atts(V, [dist(Id, Parents)]), !,
+    clpbn:get_atts(V, [dist(DistId, Parents)]), !,
-    length(Parents, NParents),
+    get_atts(V, [id(VarId)]),
+    get_dist_domain_size(DistId, DSize),
     get_evidence(V, Ev),
-    get_vars_info(Vs, LV).
+    vars2ids(Parents, ParentIds),
-get_vars_info([_|Vs], LV) :-
+    get_vars_info(Vs, VarsInfo).
-    get_vars_info(Vs, LV).
+get_vars_info([_|Vs], VarsInfo) :-
+    get_vars_info(Vs, VarsInfo).
+
+
+vars2ids([], []).
+vars2ids([V|QueryVars], [VarId|Ids]) :-
+    get_atts(V, [id(VarId)]),
+    vars2ids(QueryVars, Ids).
 
 
 get_evidence(V, Ev) :-
     clpbn:get_atts(V, [evidence(Ev)]), !.
-get_evidence(V, -1). % no evidence !!!
+get_evidence(_V, -1). % no evidence !!!
 
 
-get_dists_info([],[]).
+get_extra_vars_info([], []).
-get_dists_info([V|Vs], [dist(Id, Domain, DSize, Params, NParams) | Dists]) :-
+get_extra_vars_info([V|Vs], [v(VarId, Label, Domain)|VarsInfo]) :-
-    clpbn:get_atts(V, [dist(Id, _)]), !,
+    get_atts(V, [id(VarId)]), !,
-    get_dist_domain(Id, Domain),
+    clpbn:get_atts(V, [key(Key),dist(DistId, _)]),
-    length(Domain, DSize),
+    term_to_atom(Key, Label),
+    get_dist_domain(DistId, Domain0),
+    numbers2atoms(Domain0, Domain),
+    get_extra_vars_info(Vs, VarsInfo).
+get_extra_vars_info([_|Vs], VarsInfo) :-
+    get_extra_vars_info(Vs, VarsInfo).
+
+
+numbers2atoms([], []).
+numbers2atoms([Atom|L0], [Atom|L]) :-
+     atom(Atom), !,
+     numbers2atoms(L0, L).
+numbers2atoms([Number|L0], [Atom|L]) :-
+     number_atom(Number, Atom),
+     numbers2atoms(L0, L).
+
+
+run_bp_solver(QVsL0, LPs, (BayesNet, DistIds)) :-
+    get_dists_parameters(DistIds, DistsParams),
+    set_parameters(BayesNet, DistsParams),
+    process_query_list(QVsL0, QVsL),
+    %writeln(' qvs':QVsL),
+    %(num_bayes_nets(1506) -> writeln(qvs:QVsL) ; true),
+    run_solver(BayesNet, QVsL, LPs).
+
+
+process_query_list([], []).
+process_query_list([[V]|QueryVars], [VarId|Ids]) :- !,
+    get_atts(V, [id(VarId)]),
+    process_query_list(QueryVars, Ids).
+process_query_list([Vs|QueryVars], [VarIds|Ids]) :-
+    vars2ids(Vs, VarIds),
+    process_query_list(QueryVars, Ids).
+
+
+get_dists_parameters([],[]).
+get_dists_parameters([Id|Ids], [dist(Id, Params)|DistsInfo]) :-
     get_dist_params(Id, Params),
-    length(Params, NParams),
+    get_dists_parameters(Ids, DistsInfo).
-    get_dists_info(Vs, Dists).
-get_dists_info([_|Vs], Dists) :-
-    get_dists_info(Vs, Dists).
 
 
-% +Vars  is a list containing info about every clpbn variables
+finalize_bp_solver((BayesNet, _)) :-
-% +Dists is a list containing info about distributions
+    delete_bayes_net(BayesNet).
-% +QVs   is a list containing the query variables
-% -Out   is some output term stating the probabilities
-process(Vars, Dists, QVs, Out) :-
-    write('vars  = '), writeln(Vars),
-    order_vars(Vars, [], OrderedVars),
-    write('ovars = '), writeln(OrderedVars),
-    write('dists = '), writeln(Dists),
-    write('qvs   = '), writeln(QVs),
-    length(OrderedVars, NVars),
-    length(Dists, NDists),
-    %create_network(OrderedVars, NVars, Dists, NDists, BNet),
-    length(QVs, NQVs),
-    run_solver(BNet, QVs, NQVs, Out),
-    free_memory(BNet).
 
 
-order_vars([V], _, [V]) :- !.
+inc_num_bayes_nets :-
-order_vars([var(V, Id, Parents, NParents, Ev)|Vs], ParsedVars, [var(V, Id, Parents, NParents, Ev)|OrderedVars]) :-
+    retract(num_bayes_nets(Count0)),
-    \+ memberchk(V, ParsedVars), 
+    Count is Count0 + 1,
-    parents_defined(Parents, ParsedVars), !,
+    assert(num_bayes_nets(Count)).
-    order_vars(Vs, [V|ParsedVars], OrderedVars).
-order_vars([var(V, Id, Parents, NParents, Ev)|Vs], ParsedVars, OrderedVars) :- 
-    append(Vs, [var(V, Id, Parents, NParents, Ev)], NVs),
-    order_vars(NVs, ParsedVars, OrderedVars).
 
-
-parents_defined([], _) :- !.
-parents_defined([Parent|Parents], ParsedVars) :- 
-    memberchk(Parent, ParsedVars),
-    parents_defined(Parents, ParsedVars).
-
-% f(F), b(B). ----> FAIL

packages/CLPBN/clpbn/bp/BPSolver.cpp

@@ -0,0 +1,891 @@
+#include <cstdlib>
+#include <time.h>
+#include <algorithm>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+
+#include "BPSolver.h"
+#include "BpNode.h"
+
+
+BPSolver* Edge::klass = 0;
+StatisticMap Statistics::stats_;
+unsigned Statistics::numCreatedNets     = 0;
+unsigned Statistics::numSolvedPolyTrees = 0;
+unsigned Statistics::numSolvedLoopyNets = 0;
+unsigned Statistics::numUnconvergedRuns = 0;
+unsigned Statistics::maxIterations      = 0;
+unsigned Statistics::totalOfIterations  = 0;
+
+
+BPSolver::BPSolver (const BayesNet& bn) : Solver (&bn)
+{
+  bn_ = &bn;
+  forceGenericSolver_ = false;
+  //forceGenericSolver_ = true;
+  schedule_           = S_SEQ_FIXED;
+  //schedule_           = S_SEQ_RANDOM;
+  //schedule_           = S_PARALLEL;
+  //schedule_           = S_MAX_RESIDUAL;
+  maxIter_            = 205;
+  accuracy_           = 0.000001;
+}
+
+
+
+BPSolver::~BPSolver (void)
+{
+  for (unsigned i = 0; i < msgs_.size(); i++) {
+    delete msgs_[i];
+  }
+}
+
+
+
+void
+BPSolver::runSolver (void)
+{
+  if (DL >= 1) {
+    //bn_->printNetwork();
+  }
+
+  clock_t start_ = clock();
+  if (bn_->isSingleConnected() && !forceGenericSolver_) {
+    runPolyTreeSolver();
+    Statistics::numSolvedPolyTrees ++;
+  } else {
+    runGenericSolver();
+    Statistics::numSolvedLoopyNets ++;
+    if (nIter_ >= maxIter_) {
+      Statistics::numUnconvergedRuns ++;
+    } else {
+      Statistics::updateIterations (nIter_);
+    }
+    if (DL >= 1) {
+      cout << endl;
+      if (nIter_ < maxIter_) {
+        cout << "Belief propagation converged in " ; 
+        cout << nIter_ << " iterations" << endl;
+      } else {
+        cout << "The maximum number of iterations was hit, terminating..." ;
+        cout << endl;
+      }
+    }
+  }
+  double time = (double (clock() - start_)) / CLOCKS_PER_SEC;
+  unsigned size = bn_->getNumberOfNodes();
+  Statistics::updateStats (size, time);
+  //if (size > 30) {
+  //   stringstream ss;
+  //  ss << size << "." << Statistics::getCounting (size) << ".dot" ;
+  //  bn_->exportToDotFile (ss.str().c_str());
+  //}
+}
+
+
+
+ParamSet
+BPSolver::getPosterioriOf (const Variable* var) const
+{
+  assert (var);
+  assert (var == bn_->getNode (var->getVarId()));
+  assert (var->getIndex() < msgs_.size());
+  return msgs_[var->getIndex()]->getBeliefs();
+}
+
+
+
+
+ParamSet
+BPSolver::getJointDistribution (const NodeSet& jointVars) const
+{
+  if (DL >= 1) {
+    cout << "calculating joint distribuition on: " ;
+    for (unsigned i = 0; i < jointVars.size(); i++) {
+      cout << jointVars[i]->getLabel() << " " ;
+    }
+    cout << endl;
+  }
+
+  //BayesNet* workingNet = bn_->pruneNetwork (bn_->getNodes());
+  //FIXME see if this works:
+  BayesNet* workingNet = bn_->pruneNetwork (jointVars);
+  BayesNode* node      = workingNet->getNode (jointVars[0]->getVarId());
+
+  BayesNet* tempNet    = workingNet->pruneNetwork (node);
+  BPSolver solver (*tempNet);
+  solver.runSolver();
+
+  NodeSet observedVars = { jointVars[0] };
+
+  node = tempNet->getNode (jointVars[0]->getVarId());
+  ParamSet prevBeliefs = solver.getPosterioriOf (node);
+
+  delete tempNet;
+
+  for (unsigned i = 1; i < jointVars.size(); i++) {
+    node = workingNet->getNode (observedVars[i - 1]->getVarId());
+    if (!node->hasEvidence()) {
+      node->setEvidence (0);
+    }
+    node = workingNet->getNode (jointVars[i]->getVarId());
+    tempNet = workingNet->pruneNetwork (node);
+
+    ParamSet allBeliefs;
+    vector<DomainConf> confs = 
+        BayesNet::getDomainConfigurationsOf (observedVars);
+    for (unsigned j = 0; j < confs.size(); j++) {
+      for (unsigned k = 0; k < observedVars.size(); k++) {
+        node = tempNet->getNode (observedVars[k]->getVarId());
+        if (!observedVars[k]->hasEvidence()) {
+          if (node) {
+            node->setEvidence (confs[j][k]);
+          } else {
+            // FIXME try optimize
+            //assert (false);
+            cout << observedVars[k]->getLabel();
+            cout << " is not in temporary net!" ;
+            cout << endl;
+          }
+        } else {
+          cout << observedVars[k]->getLabel();
+          cout << " already has evidence in original net!" ;
+          cout << endl;
+        }
+      }
+      BPSolver solver (*tempNet);
+      node = tempNet->getNode (jointVars[i]->getVarId());
+      solver.runSolver();
+      ParamSet beliefs = solver.getPosterioriOf (node);
+      for (unsigned k = 0; k < beliefs.size(); k++) {
+        allBeliefs.push_back (beliefs[k]);
+      }
+    }
+
+    int count = -1;
+    for (unsigned j = 0; j < allBeliefs.size(); j++) {
+      if (j % jointVars[i]->getDomainSize() == 0) {
+        count ++;
+      }
+      allBeliefs[j] *= prevBeliefs[count];
+    }
+    prevBeliefs = allBeliefs;
+    observedVars.push_back (jointVars[i]);
+    delete tempNet;
+  }
+  delete workingNet;
+  return prevBeliefs;
+}
+
+
+
+void 
+BPSolver::initializeSolver (void)
+{
+  if (DL >= 1) {
+    cout << "Initializing solver" << endl;
+    cout << "-> schedule        = ";
+    if (forceGenericSolver_) {
+      switch (schedule_) {
+        case S_SEQ_FIXED:    cout << "sequential fixed" ;  break;
+        case S_SEQ_RANDOM:   cout << "sequential random" ; break;
+        case S_PARALLEL:     cout << "parallel" ;          break;
+        case S_MAX_RESIDUAL: cout << "max residual" ;      break;
+      }
+    } else {
+      cout << "polytree solver" ;
+    }
+    cout << endl;
+    cout << "-> max iters       = " << maxIter_ << endl;
+    cout << "-> accuracy        = " << accuracy_ << endl;
+    cout << endl;
+  }
+
+  const NodeSet& nodes = bn_->getNodes();
+  for (unsigned i = 0; i < msgs_.size(); i++) {
+    delete msgs_[i];
+  }
+  msgs_.clear();
+  msgs_.reserve (nodes.size());
+  updateOrder_.clear();
+  sortedOrder_.clear();
+  edgeMap_.clear();
+
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    msgs_.push_back (new BpNode (nodes[i]));
+  }
+
+  NodeSet roots = bn_->getRootNodes();
+  for (unsigned i = 0; i < roots.size(); i++) {
+    const ParamSet& params = roots[i]->getParameters();
+    ParamSet& piVals = M(roots[i])->getPiValues();
+    for (int ri = 0; ri < roots[i]->getDomainSize(); ri++) {
+      piVals[ri] = params[ri];
+    }
+  }
+}
+
+
+
+void 
+BPSolver::incorporateEvidence (BayesNode* x)
+{
+  ParamSet& piVals = M(x)->getPiValues();
+  ParamSet& ldVals = M(x)->getLambdaValues();
+  for (int xi = 0; xi < x->getDomainSize(); xi++) {
+    piVals[xi] = 0.0;
+    ldVals[xi] = 0.0;
+  }
+  piVals[x->getEvidence()] = 1.0;
+  ldVals[x->getEvidence()] = 1.0;
+}
+
+
+
+void
+BPSolver::runPolyTreeSolver (void)
+{
+  initializeSolver();
+  const NodeSet& nodes = bn_->getNodes();
+
+  // Hack: I need this else this can happen with bayes ball
+  // Variable: 174
+  // Id:       174
+  // Domain:   -1, 0, 1
+  // Evidence: 1
+  // Parents:  
+  // Childs:   176
+  // cpt
+  // ----------------------------------------------------
+  // -1             0           0           0           0 ...
+  // 0       0.857143    0.857143    0.857143    0.857143 ...
+  // 1       0.142857    0.142857    0.142857    0.142857 ...
+  // the cpt for this node would be 0,0,0
+
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    if (nodes[i]->hasEvidence()) {
+      incorporateEvidence (nodes[i]);
+    }
+  } 
+
+  // first compute all node marginals ...
+  NodeSet roots = bn_->getRootNodes();
+  for (unsigned i = 0; i < roots.size(); i++) {
+    const NodeSet& childs = roots[i]->getChilds();
+    for (unsigned j = 0; j < childs.size(); j++) {
+      polyTreePiMessage (roots[i], childs[j]);
+    }
+  }
+  // then propagate the evidence
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    if (nodes[i]->hasEvidence()) {
+      incorporateEvidence (nodes[i]); 
+      const NodeSet& parents = nodes[i]->getParents();
+      for (unsigned j = 0; j < parents.size(); j++) {
+        if (!parents[j]->hasEvidence()) {
+          polyTreeLambdaMessage (nodes[i], parents[j]);
+        }
+      }
+      const NodeSet& childs = nodes[i]->getChilds();
+      for (unsigned j = 0; j < childs.size(); j++) {
+        polyTreePiMessage (nodes[i], childs[j]);
+      }
+    }
+  }
+}
+
+
+
+void
+BPSolver::polyTreePiMessage (BayesNode* z, BayesNode* x)
+{
+  if (DL >= 1) {
+    cout << PI << " (" << z->getLabel();
+    cout << " --> " << x->getLabel();
+    cout << ")" << endl;
+  }
+  calculateNextPiMessage (z, x);
+  updatePiMessage (z, x);
+
+  if (!x->hasEvidence()) {
+    updatePiValues (x);
+    const NodeSet& xChilds = x->getChilds();
+    for (unsigned i = 0; i < xChilds.size(); i++) {
+      polyTreePiMessage (x, xChilds[i]);
+    }
+  }
+
+  if (M(x)->hasReceivedChildInfluence()) {
+    const NodeSet& xParents = x->getParents();
+    for (unsigned i = 0; i < xParents.size(); i++) {
+      if (xParents[i] != z && !xParents[i]->hasEvidence()) {
+        polyTreeLambdaMessage (x, xParents[i]);
+      }
+    }
+  }
+}
+
+
+
+void 
+BPSolver::polyTreeLambdaMessage (BayesNode* y, BayesNode* x)
+{
+  if (DL >= 1) {
+    cout << LD << " (" << y->getLabel();
+    cout << " --> " << x->getLabel();
+    cout << ")" << endl;
+  }
+  calculateNextLambdaMessage (y, x);
+  updateLambdaMessage (y, x);
+  updateLambdaValues (x);
+
+  const NodeSet& xParents = x->getParents();
+  for (unsigned i = 0; i < xParents.size(); i++) {
+    if (!xParents[i]->hasEvidence()) {
+      polyTreeLambdaMessage (x, xParents[i]);
+    }
+  }
+
+  const NodeSet& xChilds = x->getChilds();
+  for (unsigned i = 0; i < xChilds.size(); i++) {
+    if (xChilds[i] != y) {
+      polyTreePiMessage (x, xChilds[i]);
+    }
+  }
+}
+
+
+
+void
+BPSolver::runGenericSolver()
+{ 
+  initializeSolver();
+  const NodeSet& nodes = bn_->getNodes();
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    if (nodes[i]->hasEvidence()) {
+      incorporateEvidence (nodes[i]);
+    }
+  }
+ 
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    // pi messages
+    const NodeSet& childs = nodes[i]->getChilds();
+    for (unsigned j = 0; j < childs.size(); j++) {
+      updateOrder_.push_back (Edge (nodes[i], childs[j], PI_MSG));
+    }
+    // lambda messages
+    const NodeSet& parents = nodes[i]->getParents();
+    for (unsigned j = 0; j < parents.size(); j++) {
+      if (!parents[j]->hasEvidence()) {
+        updateOrder_.push_back (Edge (nodes[i], parents[j], LAMBDA_MSG));
+      }
+    }
+  }
+
+  nIter_ = 0;
+  while (!converged() && nIter_ < maxIter_) {
+
+    nIter_++;
+    if (DL >= 1) {
+      cout << endl;
+      cout << "****************************************" ;
+      cout << "****************************************" ;
+      cout << endl;
+      cout << " Iteration " << nIter_ << endl;
+      cout << "****************************************" ;
+      cout << "****************************************" ;
+      cout << endl;
+    }
+ 
+    switch (schedule_) {
+
+      case S_SEQ_RANDOM:
+        random_shuffle (updateOrder_.begin(), updateOrder_.end());
+        // no break
+  
+      case S_SEQ_FIXED:
+        for (unsigned i = 0; i < updateOrder_.size(); i++) {
+          calculateNextMessage (updateOrder_[i]);
+          updateMessage (updateOrder_[i]);
+          updateValues (updateOrder_[i]);
+        }
+        break;
+
+      case S_PARALLEL:
+        for (unsigned i = 0; i < updateOrder_.size(); i++) {
+          calculateNextMessage (updateOrder_[i]);
+        }
+        for (unsigned i = 0; i < updateOrder_.size(); i++) {
+          updateMessage (updateOrder_[i]);
+          updateValues (updateOrder_[i]);
+        }
+        break;
+ 
+      case S_MAX_RESIDUAL:
+        maxResidualSchedule();
+        break;
+
+    }
+  }
+}
+
+
+
+void
+BPSolver::maxResidualSchedule (void)
+{
+  if (nIter_ == 1) {
+    Edge::klass = this;
+    for (unsigned i = 0; i < updateOrder_.size(); i++) {
+      calculateNextMessage (updateOrder_[i]);
+      updateResidual (updateOrder_[i]);
+      SortedOrder::iterator it = sortedOrder_.insert (updateOrder_[i]);
+      edgeMap_.insert (make_pair (updateOrder_[i].getId(), it));
+    }
+    return;
+  }
+
+  for (unsigned c = 0; c < sortedOrder_.size(); c++) {
+    if (DL >= 1) {
+      for (set<Edge, compare>::iterator it = sortedOrder_.begin();
+          it != sortedOrder_.end(); it ++) {
+        cout << it->toString() << " residual = " ;
+        cout << getResidual (*it) << endl;
+      }
+    }
+ 
+    set<Edge, compare>::iterator it = sortedOrder_.begin();
+    Edge e = *it;
+    if (getResidual (e) < accuracy_) {
+      return;
+    }
+    updateMessage (e);
+    updateValues (e);
+    clearResidual (e);
+    sortedOrder_.erase (it);
+    assert (edgeMap_.find (e.getId()) != edgeMap_.end());
+    edgeMap_.find (e.getId())->second = sortedOrder_.insert (e);
+
+     // update the messages that depend on message source --> destination
+    const NodeSet& childs = e.destination->getChilds();
+    for (unsigned i = 0; i < childs.size(); i++) {
+      if (childs[i] != e.source) {
+        Edge neighbor (e.destination, childs[i], PI_MSG);
+        calculateNextMessage (neighbor);
+        updateResidual (neighbor);
+        assert (edgeMap_.find (neighbor.getId()) != edgeMap_.end());
+        EdgeMap::iterator iter = edgeMap_.find (neighbor.getId());
+        sortedOrder_.erase (iter->second);
+        iter->second = sortedOrder_.insert (neighbor);
+      }
+    }
+    const NodeSet& parents = e.destination->getParents();
+    for (unsigned i = 0; i < parents.size(); i++) {
+      if (parents[i] != e.source && !parents[i]->hasEvidence()) {
+        Edge neighbor (e.destination, parents[i], LAMBDA_MSG);
+        calculateNextMessage (neighbor);
+        updateResidual (neighbor);
+        assert (edgeMap_.find (neighbor.getId()) != edgeMap_.end());
+        EdgeMap::iterator iter = edgeMap_.find (neighbor.getId());
+        sortedOrder_.erase (iter->second);
+        iter->second = sortedOrder_.insert (neighbor);
+      }
+    }
+  }
+}
+
+
+
+bool
+BPSolver::converged (void) const
+{
+  bool converged = true;
+  if (schedule_ == S_MAX_RESIDUAL) {
+    if (nIter_ <= 2) {
+      return false;
+    }
+    // this can happen if every node does not have neighbors
+    if (sortedOrder_.size() == 0) {
+      return true;
+    }
+    Param maxResidual = getResidual (*(sortedOrder_.begin()));
+    if (maxResidual > accuracy_) {
+      return false;
+    }
+  } else {
+    if (nIter_ == 0) {
+      return false;
+    }
+    const NodeSet& nodes = bn_->getNodes();
+    for (unsigned i = 0; i < nodes.size(); i++) {
+      if (!nodes[i]->hasEvidence()) {
+        double change = M(nodes[i])->getBeliefChange();
+        if (DL >= 1) {
+          cout << nodes[i]->getLabel() + " belief change = " ;
+          cout << change << endl;
+        }
+        if (change > accuracy_) {
+          converged = false;
+          if (DL == 0) break;
+        }
+      }
+    }
+  }
+
+  return converged;
+}
+
+
+
+void
+BPSolver::updatePiValues (BayesNode* x)
+{
+  // π(Xi)
+  const NodeSet& parents = x->getParents();
+  const vector<CptEntry>& entries = x->getCptEntries();
+  assert (parents.size() != 0);
+  stringstream* calcs1;
+  stringstream* calcs2;
+
+  ParamSet messageProducts (entries.size());
+  for (unsigned k = 0; k < entries.size(); k++) {
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+    }
+    double messageProduct = 1.0;
+    const DomainConf& conf = entries[k].getParentConfigurations();
+    for (unsigned i = 0; i < parents.size(); i++) {
+      messageProduct *= M(parents[i])->getPiMessageValue(x, conf[i]);
+      if (DL >= 5) {
+        if (i != 0) *calcs1 << "." ;
+        if (i != 0) *calcs2 << "*" ;
+        *calcs1 << PI << "(" << x->getLabel() << ")" ;
+        *calcs1 << "[" << parents[i]->getDomain()[conf[i]] << "]";
+        *calcs2 << M(parents[i])->getPiMessageValue(x, conf[i]);
+      }
+    }
+    messageProducts[k] = messageProduct;
+    if (DL >= 5) {
+      cout << "    mp" << k;
+      cout << " = " << (*calcs1).str();
+      if (parents.size() == 1) {
+        cout << " = " << messageProduct << endl;
+      } else {
+        cout << " = " << (*calcs2).str();
+        cout << " = " << messageProduct << endl;
+      }
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+
+  for (int xi = 0; xi < x->getDomainSize(); xi++) {
+    double sum = 0.0;
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+    }
+    for (unsigned k = 0; k < entries.size(); k++) {
+      sum += x->getProbability (xi, entries[k]) * messageProducts[k];
+      if (DL >= 5) {
+        if (k != 0) *calcs1 << " + " ;
+        if (k != 0) *calcs2 << " + " ;
+        *calcs1 << x->cptEntryToString (xi, entries[k]); 
+        *calcs1 << ".mp" << k;
+        *calcs2 << x->getProbability (xi, entries[k]);
+        *calcs2 << "*" << messageProducts[k];
+      }
+    }
+    M(x)->setPiValue (xi, sum);
+    if (DL >= 5) {
+      cout << "    " << PI << "(" << x->getLabel() << ")" ;
+      cout << "[" << x->getDomain()[xi] << "]" ;
+      cout << " = " << (*calcs1).str();
+      cout << " = " << (*calcs2).str();
+      cout << " = " << sum << endl;
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+}
+
+
+
+void
+BPSolver::updateLambdaValues (BayesNode* x)
+{
+  // λ(Xi)
+  const NodeSet& childs = x->getChilds();
+  assert (childs.size() != 0);
+  stringstream* calcs1;
+  stringstream* calcs2;
+
+  for (int xi = 0; xi < x->getDomainSize(); xi++) {
+    double product = 1.0;
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+    }
+    for (unsigned i = 0; i < childs.size(); i++) {
+      product *= M(x)->getLambdaMessageValue(childs[i], xi);
+      if (DL >= 5) {
+        if (i != 0) *calcs1 << "." ;
+        if (i != 0) *calcs2 << "*" ;
+        *calcs1 << LD << "(" << childs[i]->getLabel();
+        *calcs1 << "-->" << x->getLabel() << ")" ;
+        *calcs1 << "[" << x->getDomain()[xi] << "]" ;
+        *calcs2 << M(x)->getLambdaMessageValue(childs[i], xi);
+      }
+    }
+    M(x)->setLambdaValue (xi, product);
+    if (DL >= 5) {
+      cout << "    " << LD << "(" << x->getLabel() << ")" ;
+      cout << "[" << x->getDomain()[xi] << "]" ;
+      cout << " = " << (*calcs1).str();
+      if (childs.size() == 1) {
+        cout << " = " << product << endl;
+      } else {
+        cout << " = " << (*calcs2).str();
+        cout << " = " << product << endl;
+      }
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+}
+
+
+
+void
+BPSolver::calculateNextPiMessage (BayesNode* z, BayesNode* x)
+{
+  // πX(Zi)
+  ParamSet& zxPiNextMessage = M(z)->piNextMessageReference (x);
+  const NodeSet& zChilds = z->getChilds();
+  stringstream* calcs1;
+  stringstream* calcs2;
+
+  for (int zi = 0; zi < z->getDomainSize(); zi++) {
+    double product = M(z)->getPiValue (zi);
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+      *calcs1 << PI << "(" << z->getLabel() << ")";
+      *calcs1 << "[" << z->getDomain()[zi] << "]" ;
+      *calcs2 << product;
+    }
+    for (unsigned i = 0; i < zChilds.size(); i++) {
+      if (zChilds[i] != x) {
+        product *= M(z)->getLambdaMessageValue(zChilds[i], zi);
+        if (DL >= 5) {
+          *calcs1 << "." << LD << "(" << zChilds[i]->getLabel();
+          *calcs1 << "-->" << z->getLabel() << ")";
+          *calcs1 << "[" << z->getDomain()[zi] + "]" ;
+          *calcs2 << " * " << M(z)->getLambdaMessageValue(zChilds[i], zi);
+        }
+      }
+    }
+    zxPiNextMessage[zi] = product;
+    if (DL >= 5) {
+      cout << "    " << PI << "(" << z->getLabel();
+      cout << "-->" << x->getLabel() << ")" ;
+      cout << "["  << z->getDomain()[zi] << "]" ;
+      cout << " = " << (*calcs1).str();
+      if (zChilds.size() == 1) {
+        cout << " = " << product << endl;
+      } else {
+        cout << " = " << (*calcs2).str();
+        cout << " = " << product << endl;
+      }
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+}
+
+
+
+void
+BPSolver::calculateNextLambdaMessage (BayesNode* y, BayesNode* x)
+{
+  // λY(Xi)
+  //if (!y->hasEvidence() && !M(y)->hasReceivedChildInfluence()) {
+  //  if (DL >= 5) {
+  //    cout << "unnecessary calculation" << endl;
+  //  }
+  //  return;
+  //}
+  ParamSet& yxLambdaNextMessage       = M(x)->lambdaNextMessageReference (y);
+  const NodeSet& yParents             = y->getParents();
+  const vector<CptEntry>& allEntries  = y->getCptEntries();
+  int parentIndex                     = y->getIndexOfParent (x);
+  stringstream* calcs1;
+  stringstream* calcs2;
+ 
+  vector<CptEntry> entries;
+  DomainConstr constr = make_pair (parentIndex, 0);
+  for (unsigned i = 0; i < allEntries.size(); i++) {
+    if (allEntries[i].matchConstraints(constr)) {
+      entries.push_back (allEntries[i]);
+    }
+  }
+
+  ParamSet messageProducts (entries.size());
+  for (unsigned k = 0; k < entries.size(); k++) {
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+    }
+    double messageProduct = 1.0;
+    const DomainConf& conf = entries[k].getParentConfigurations();
+    for (unsigned i = 0; i < yParents.size(); i++) {
+      if (yParents[i] != x) {
+        if (DL >= 5) {
+          if (messageProduct != 1.0) *calcs1 << "*" ;
+          if (messageProduct != 1.0) *calcs2 << "*" ;
+          *calcs1 << PI << "(" << yParents[i]->getLabel();
+          *calcs1 << "-->" << y->getLabel() << ")" ;
+          *calcs1 << "[" << yParents[i]->getDomain()[conf[i]] << "]" ;
+          *calcs2 << M(yParents[i])->getPiMessageValue(y, conf[i]);
+        }
+        messageProduct *= M(yParents[i])->getPiMessageValue(y, conf[i]);
+      }
+    }
+    messageProducts[k] = messageProduct;
+    if (DL >= 5) {
+      cout << "    mp" << k;
+      cout << " = " << (*calcs1).str();
+      if  (yParents.size() == 1) {
+        cout << 1 << endl;
+      } else if (yParents.size() == 2) {
+        cout << " = " << messageProduct << endl;
+      } else {
+        cout << " = " << (*calcs2).str();
+        cout << " = " << messageProduct << endl;
+      }
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+
+  for (int xi = 0; xi < x->getDomainSize(); xi++) {
+    if (DL >= 5) {
+      calcs1 = new stringstream;
+      calcs2 = new stringstream;
+    }
+    vector<CptEntry> entries;
+    DomainConstr constr = make_pair (parentIndex, xi);
+    for (unsigned i = 0; i < allEntries.size(); i++) {
+      if (allEntries[i].matchConstraints(constr)) {
+        entries.push_back (allEntries[i]);
+      }
+    }
+    double outerSum = 0.0;
+    for (int yi = 0; yi < y->getDomainSize(); yi++) {
+      if (DL >= 5) {
+        (yi != 0) ? *calcs1 << " + {" : *calcs1 << "{" ;
+        (yi != 0) ? *calcs2 << " + {" : *calcs2 << "{" ;
+      } 
+      double innerSum = 0.0;
+      for (unsigned k = 0; k < entries.size(); k++) {
+        if (DL >= 5) {
+          if (k != 0) *calcs1 << " + " ;
+          if (k != 0) *calcs2 << " + " ;
+          *calcs1 << y->cptEntryToString (yi, entries[k]);
+          *calcs1 << ".mp" << k;
+          *calcs2 << y->getProbability (yi, entries[k]);
+          *calcs2 << "*" << messageProducts[k];
+        }
+        innerSum += y->getProbability (yi, entries[k]) * messageProducts[k];
+      }
+      outerSum += innerSum * M(y)->getLambdaValue (yi);
+      if (DL >= 5) {
+        *calcs1 << "}." << LD << "(" << y->getLabel() << ")" ;
+        *calcs1 << "["  << y->getDomain()[yi] << "]";
+        *calcs2 << "}*" << M(y)->getLambdaValue (yi);
+      }
+    }
+    yxLambdaNextMessage[xi] = outerSum;
+    if (DL >= 5) {
+      cout << "    " << LD << "(" << y->getLabel();
+      cout << "-->" << x->getLabel() << ")" ;
+      cout << "[" << x->getDomain()[xi] << "]" ;
+      cout << " = " << (*calcs1).str();
+      cout << " = " << (*calcs2).str();
+      cout << " = " << outerSum << endl;
+      delete calcs1;
+      delete calcs2;
+    }
+  }
+}
+
+
+
+void
+BPSolver::printMessageStatusOf (const BayesNode* var) const
+{
+  cout << left;
+  cout << setw (10) << "domain" ;
+  cout << setw (20) << PI << "(" + var->getLabel() + ")" ;
+  cout << setw (20) << LD << "(" + var->getLabel() + ")" ;
+  cout << setw (16) << "belief" ;
+  cout << endl;
+  cout << "--------------------------------" ;
+  cout << "--------------------------------" ;
+  cout << endl;
+
+  BpNode*         x        = M(var);
+  ParamSet&       piVals   = x->getPiValues();
+  ParamSet&       ldVals   = x->getLambdaValues();
+  ParamSet        beliefs  = x->getBeliefs();
+  const Domain&   domain   = var->getDomain();
+  const NodeSet&  childs   = var->getChilds();
+
+  for (int xi = 0; xi < var->getDomainSize(); xi++) {
+    cout << setw (10) << domain[xi];
+    cout << setw (19) << piVals[xi];
+    cout << setw (19) << ldVals[xi];
+    cout.precision (PRECISION);
+    cout << setw (16) << beliefs[xi];
+    cout << endl;
+  }
+  cout << endl;
+  if (childs.size() > 0) {
+    string s = "(" + var->getLabel() + ")" ;
+    for (unsigned j = 0; j < childs.size(); j++) {
+      cout << setw (10) << "domain" ;
+      cout << setw (28) << PI + childs[j]->getLabel() + s;
+      cout << setw (28) << LD + childs[j]->getLabel() + s;
+      cout << endl;
+      cout << "--------------------------------" ;
+      cout << "--------------------------------" ;
+      cout << endl;
+      const ParamSet& piMessage     = x->getPiMessage (childs[j]);
+      const ParamSet& lambdaMessage = x->getLambdaMessage (childs[j]);
+      for (int xi = 0; xi < var->getDomainSize(); xi++) {
+        cout << setw (10) << domain[xi];
+        cout.precision (PRECISION);
+        cout << setw (27) << piMessage[xi];
+        cout.precision (PRECISION);
+        cout << setw (27) << lambdaMessage[xi];
+        cout << endl;
+      }
+      cout << endl;
+    }
+  }
+}
+
+
+
+void 
+BPSolver::printAllMessageStatus (void) const
+{ 
+  const NodeSet& nodes = bn_->getNodes();
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    printMessageStatusOf (nodes[i]);
+  }
+}
+
+

packages/CLPBN/clpbn/bp/BPSolver.h

@@ -0,0 +1,450 @@
+#ifndef BP_BPSOLVER_H
+#define BP_BPSOLVER_H
+
+#include <vector>
+#include <string>
+#include <set>
+
+#include "Solver.h"
+#include "BayesNet.h"
+#include "BpNode.h"
+#include "Shared.h"
+
+using namespace std;
+
+class BPSolver;
+
+static const string PI = "pi" ;
+static const string LD = "ld" ;
+
+enum MessageType {PI_MSG, LAMBDA_MSG};
+
+class BPSolver;
+struct Edge
+{
+  Edge (BayesNode* s, BayesNode* d, MessageType t)
+  { 
+    source       = s;
+    destination  = d;
+    type         = t;
+  }
+  string getId (void) const 
+  {
+    stringstream ss;
+    type == PI_MSG ? ss << PI : ss << LD;
+    ss << source->getVarId() << "." << destination->getVarId();
+    return ss.str();
+  }
+  string toString (void) const
+  {
+    stringstream ss;
+    type == PI_MSG ? ss << PI << "(" : ss << LD << "(" ;
+    ss << source->getLabel() << " --> " ;
+    ss << destination->getLabel();
+    ss << ")" ;
+    return ss.str();
+  }
+  BayesNode*        source;
+  BayesNode*        destination;
+  MessageType       type;
+  static BPSolver*  klass;
+};
+
+
+
+/*
+class BPMessage
+{
+  BPMessage (BayesNode* parent, BayesNode* child) 
+  {
+    parent_ = parent;
+    child_  = child;
+    currPiMsg_.resize (child->getDomainSize(),  1);
+    currLdMsg_.resize (parent->getDomainSize(), 1);
+    nextLdMsg_.resize (parent->getDomainSize(), 1);
+    nextPiMsg_.resize (child->getDomainSize(),  1);
+    piResidual_ = 1.0;
+    ldResidual_ = 1.0;
+  }
+
+  Param getPiMessageValue (int idx) const
+  {
+    assert (idx >=0 && idx < child->getDomainSize());
+    return currPiMsg_[idx];
+  }
+
+  Param getLambdaMessageValue (int idx) const
+  {
+    assert (idx >=0 && idx < parent->getDomainSize());
+    return currLdMsg_[idx];
+  }
+
+  const ParamSet& getPiMessage (void) const
+  {
+    return currPiMsg_;
+  }
+
+  const ParamSet& getLambdaMessage (void) const
+  {
+    return currLdMsg_;
+  }
+
+  ParamSet& piNextMessageReference (void)
+  {
+    return nextPiMsg_;
+  }
+
+  ParamSet& lambdaNextMessageReference (const BayesNode* source)
+  {
+    return nextLdMsg_;
+  }
+
+  void updatePiMessage (void)
+  {
+    currPiMsg_ = nextPiMsg_;
+    Util::normalize (currPiMsg_);
+  }
+
+  void updateLambdaMessage (void)
+  {
+    currLdMsg_ = nextLdMsg_;
+    Util::normalize (currLdMsg_);
+  }
+
+  double getPiResidual (void)
+  {
+    return piResidual_;
+  }
+
+  double getLambdaResidual (void)
+  {
+    return ldResidual_;
+  }
+
+  void updatePiResidual (void)
+  {
+    piResidual_ = Util::getL1dist (currPiMsg_, nextPiMsg_);
+  }
+
+  void updateLambdaResidual (void)
+  {
+    ldResidual_ = Util::getL1dist (currLdMsg_, nextLdMsg_);
+  }
+
+  void clearPiResidual (void)
+  {
+    piResidual_ = 0.0;
+  }
+
+  void clearLambdaResidual (void)
+  {
+    ldResidual_ = 0.0;
+  }
+
+  BayesNode*        parent_;
+  BayesNode*        child_;
+  ParamSet          currPiMsg_; // current pi messages
+  ParamSet          currLdMsg_; // current lambda messages
+  ParamSet          nextPiMsg_;
+  ParamSet          nextLdMsg_;
+  Param             piResidual_;
+  Param             ldResidual_;
+};
+
+
+
+class NodeInfo
+{
+  NodeInfo (BayesNode* node)
+  {
+    node_ = node;
+    piVals_.resize (node->getDomainSize(), 1);
+    ldVals_.resize (node->getDomainSize(), 1);
+  }
+
+  ParamSet getBeliefs (void) const
+  {
+    double sum = 0.0;
+    ParamSet beliefs (node_->getDomainSize());
+    for (int xi = 0; xi < node_->getDomainSize(); xi++) {
+      double prod = piVals_[xi] * ldVals_[xi];
+      beliefs[xi] = prod;
+      sum += prod;
+    }
+    assert (sum);
+    //normalize the beliefs
+    for (int xi = 0; xi < node_->getDomainSize(); xi++) {
+      beliefs[xi] /= sum;
+    }
+    return beliefs;
+  }
+
+  double getPiValue (int idx) const
+  {
+    assert (idx >=0 && idx < node_->getDomainSize());
+    return piVals_[idx];
+  }
+
+  void setPiValue (int idx, double value)
+  {
+    assert (idx >=0 && idx < node_->getDomainSize());
+    piVals_[idx] = value;
+  }
+
+  double getLambdaValue (int idx) const
+  {
+    assert (idx >=0 && idx < node_->getDomainSize());
+    return ldVals_[idx];
+  }
+
+  void setLambdaValue (int idx, double value)
+  {
+    assert (idx >=0 && idx < node_->getDomainSize());
+    ldVals_[idx] = value;
+  }
+
+  ParamSet& getPiValues (void)
+  {
+    return piVals_;
+  }
+
+  ParamSet& getLambdaValues (void)
+  {
+    return ldVals_;
+  }
+
+  double getBeliefChange (void)
+  {
+    double change = 0.0;
+    if (oldBeliefs_.size() == 0) {
+      oldBeliefs_ = getBeliefs();
+      change = MAX_CHANGE_;
+    } else {
+      ParamSet currentBeliefs = getBeliefs();
+      for (int xi = 0; xi < node_->getDomainSize(); xi++) {
+        change += abs (currentBeliefs[xi] - oldBeliefs_[xi]);
+      }
+      oldBeliefs_ = currentBeliefs;
+    }
+    return change;
+  }
+
+  bool hasReceivedChildInfluence (void) const
+  {
+    // if all lambda values are equal, then neither
+    // this node neither its descendents have evidence,
+    // we can use this to don't send lambda messages his parents
+    bool childInfluenced = false;
+    for (int xi = 1; xi < node_->getDomainSize(); xi++) {
+      if (ldVals_[xi] != ldVals_[0]) {
+        childInfluenced = true;
+        break;
+      }
+    }
+    return childInfluenced;
+  }
+
+  BayesNode*             node_;
+  ParamSet               piVals_;     // pi values
+  ParamSet               ldVals_;     // lambda values
+  ParamSet               oldBeliefs_;
+};
+*/
+
+
+bool compareResidual (const Edge&, const Edge&);
+
+class BPSolver : public Solver
+{
+  public:
+    BPSolver (const BayesNet&);
+   ~BPSolver (void);
+
+    void                runSolver (void);
+    ParamSet            getPosterioriOf (const Variable* var) const;
+    ParamSet            getJointDistribution (const NodeSet&) const;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (BPSolver);
+
+    void                initializeSolver (void);
+    void                incorporateEvidence (BayesNode*);
+    void                runPolyTreeSolver (void);
+    void                polyTreePiMessage (BayesNode*, BayesNode*);
+    void                polyTreeLambdaMessage (BayesNode*, BayesNode*);
+    void                runGenericSolver (void);
+    void                maxResidualSchedule (void);
+    bool                converged (void) const;
+    void                updatePiValues (BayesNode*);
+    void                updateLambdaValues (BayesNode*);
+    void                calculateNextPiMessage (BayesNode*, BayesNode*);
+    void                calculateNextLambdaMessage (BayesNode*, BayesNode*);
+    void                printMessageStatusOf (const BayesNode*) const;
+    void                printAllMessageStatus (void) const;
+    // inlines
+    void                updatePiMessage (BayesNode*, BayesNode*);
+    void                updateLambdaMessage (BayesNode*, BayesNode*);
+    void                calculateNextMessage (const Edge&);
+    void                updateMessage (const Edge&);
+    void                updateValues (const Edge&);
+    double              getResidual (const Edge&) const;
+    void                updateResidual (const Edge&);
+    void                clearResidual (const Edge&);
+    BpNode*             M (const BayesNode*) const;
+    friend bool         compareResidual (const Edge&, const Edge&);
+
+    const BayesNet*     bn_;
+    vector<BpNode*>     msgs_;
+    Schedule            schedule_;
+    int                 nIter_;
+    int                 maxIter_;
+    double              accuracy_;
+    vector<Edge>        updateOrder_;
+    bool                forceGenericSolver_;
+
+    struct compare
+    {
+      inline bool operator() (const Edge& e1, const Edge& e2)
+      {
+        return compareResidual (e1, e2);
+      }
+    };
+
+    typedef multiset<Edge, compare> SortedOrder;
+    SortedOrder sortedOrder_;
+
+    typedef unordered_map<string, SortedOrder::iterator> EdgeMap;
+    EdgeMap edgeMap_;
+
+};
+
+
+
+inline void
+BPSolver::updatePiMessage (BayesNode* source, BayesNode* destination)
+{
+  M(source)->updatePiMessage(destination);
+}
+
+
+
+inline void
+BPSolver::updateLambdaMessage (BayesNode* source, BayesNode* destination)
+{
+  M(destination)->updateLambdaMessage(source);
+}
+
+
+
+inline void
+BPSolver::calculateNextMessage (const Edge& e)
+{
+  if (DL >= 1) {
+    cout << "calculating " << e.toString() << endl;
+  }
+  if (e.type == PI_MSG) {
+    calculateNextPiMessage (e.source, e.destination);
+  } else {
+    calculateNextLambdaMessage (e.source, e.destination);
+  }
+}
+
+
+
+inline void
+BPSolver::updateMessage (const Edge& e)
+{
+  if (DL >= 1) {
+      cout << "updating " << e.toString() << endl;
+    }
+  if (e.type == PI_MSG) {
+    M(e.source)->updatePiMessage(e.destination);
+  } else {
+    M(e.destination)->updateLambdaMessage(e.source);
+  }
+}
+
+
+
+inline void
+BPSolver::updateValues (const Edge& e)
+{
+  if (!e.destination->hasEvidence()) {
+    if (e.type == PI_MSG) {
+      updatePiValues (e.destination);
+    } else {
+      updateLambdaValues (e.destination);
+    }
+  }
+}
+
+
+
+inline double
+BPSolver::getResidual (const Edge& e) const
+{
+  if (e.type == PI_MSG) {
+    return M(e.source)->getPiResidual(e.destination);
+  } else {
+    return M(e.destination)->getLambdaResidual(e.source);
+  }
+}
+
+
+
+inline void
+BPSolver::updateResidual (const Edge& e)
+{
+  if (e.type == PI_MSG) {
+    M(e.source)->updatePiResidual(e.destination);
+  } else {
+    M(e.destination)->updateLambdaResidual(e.source);
+  }
+}
+
+
+
+inline void
+BPSolver::clearResidual (const Edge& e)
+{
+  if (e.type == PI_MSG) {
+    M(e.source)->clearPiResidual(e.destination);
+  } else {
+    M(e.destination)->clearLambdaResidual(e.source);
+  }
+}
+
+
+
+inline bool
+compareResidual (const Edge& e1, const Edge& e2)
+{
+  double residual1;
+  double residual2;
+  if (e1.type == PI_MSG) {
+    residual1 = Edge::klass->M(e1.source)->getPiResidual(e1.destination);
+  } else {
+    residual1 = Edge::klass->M(e1.destination)->getLambdaResidual(e1.source);
+  }
+  if (e2.type == PI_MSG) {
+    residual2 = Edge::klass->M(e2.source)->getPiResidual(e2.destination);
+  } else {
+    residual2 = Edge::klass->M(e2.destination)->getLambdaResidual(e2.source);
+  }
+  return residual1 > residual2;
+}
+
+
+
+inline BpNode*
+BPSolver::M (const BayesNode* node) const
+{
+  assert (node);
+  assert (node == bn_->getNode (node->getVarId()));
+  assert (node->getIndex() < msgs_.size());
+  return msgs_[node->getIndex()];
+}
+
+
+#endif
+

packages/CLPBN/clpbn/bp/BayesNet.cpp

@@ -0,0 +1,792 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <iomanip>
+#include <cassert>
+#include <cstdlib>
+#include <map>
+
+#include "xmlParser/xmlParser.h"
+
+#include "BayesNet.h"
+
+
+BayesNet::BayesNet (void)
+{
+}
+
+
+
+BayesNet::BayesNet (const char* fileName)
+{
+  map<string, Domain> domains;
+  XMLNode xMainNode = XMLNode::openFileHelper (fileName, "BIF");
+  // only the first network is parsed, others are ignored
+  XMLNode xNode = xMainNode.getChildNode ("NETWORK");
+  int nVars = xNode.nChildNode ("VARIABLE");
+  for (int i = 0; i < nVars; i++) {
+    XMLNode var = xNode.getChildNode ("VARIABLE", i);
+    string type = var.getAttribute ("TYPE");
+    if (type != "nature") {
+      cerr << "error: only \"nature\" variables are supported" << endl;
+      abort();
+    }
+    Domain domain;
+    string label = var.getChildNode("NAME").getText();
+    int domainSize = var.nChildNode ("OUTCOME");
+    for (int j = 0; j < domainSize; j++) {
+      if (var.getChildNode("OUTCOME", j).getText() == 0) {
+        stringstream ss;
+        ss << j + 1;
+        domain.push_back (ss.str());
+      } else {
+        domain.push_back (var.getChildNode("OUTCOME", j).getText());
+      }
+    }
+    domains.insert (make_pair (label, domain));
+  }
+
+  int nDefs = xNode.nChildNode ("DEFINITION");
+  if (nVars != nDefs) {
+    cerr << "error: different number of variables and definitions";
+    cerr << endl;
+  }
+
+  queue<int> indexes;
+  for (int i = 0; i < nDefs; i++) {
+    indexes.push (i);
+  }
+
+  while (!indexes.empty()) {
+    int index = indexes.front();
+    indexes.pop();
+    XMLNode def = xNode.getChildNode ("DEFINITION", index);
+    string label = def.getChildNode("FOR").getText();
+    map<string, Domain>::const_iterator iter;
+    iter = domains.find (label);
+    if (iter == domains.end()) {
+      cerr << "error: unknow variable `" << label << "'" << endl;
+      abort();
+    }
+    bool processItLatter = false;
+    NodeSet parents;
+    int nParams = iter->second.size();
+    for (int j = 0; j < def.nChildNode ("GIVEN"); j++) {
+      string parentLabel = def.getChildNode("GIVEN", j).getText();
+      BayesNode* parentNode = getNode (parentLabel);
+      if (parentNode) {
+        nParams *= parentNode->getDomainSize();
+        parents.push_back (parentNode);
+      }
+      else {
+        iter = domains.find (parentLabel);
+        if (iter == domains.end()) {
+          cerr << "error: unknow parent `" << parentLabel << "'" << endl;
+          abort();
+        } else {
+          // this definition contains a parent that doesn't
+          // have a corresponding bayesian node instance yet,
+          // so process this definition latter
+          indexes.push (index);
+          processItLatter = true;
+          break;
+        }
+      }
+    }
+
+    if (!processItLatter) {
+      int count = 0;
+      ParamSet params (nParams);
+      stringstream s (def.getChildNode("TABLE").getText());
+      while (!s.eof() && count < nParams) {
+        s >> params[count];
+        count ++;
+      }
+       if (count != nParams) {
+        cerr << "error: invalid number of parameters " ;
+        cerr << "for variable `" << label << "'" << endl;
+        abort();
+      } 
+      params = reorderParameters (params, iter->second.size());
+      addNode (label, iter->second, parents, params);
+    }
+  }
+  setIndexes();
+}
+
+
+
+BayesNet::~BayesNet (void)
+{
+  Statistics::writeStats();
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    delete nodes_[i];
+  }
+}
+
+
+
+BayesNode*
+BayesNet::addNode (unsigned varId)
+{
+  indexMap_.insert (make_pair (varId, nodes_.size()));
+  nodes_.push_back (new BayesNode (varId));
+  return nodes_.back();
+}
+
+
+
+BayesNode*
+BayesNet::addNode (unsigned varId,
+                   unsigned dsize,
+                   int evidence,
+                   NodeSet& parents,
+                   Distribution* dist)
+{
+  indexMap_.insert (make_pair (varId, nodes_.size()));
+  nodes_.push_back (new BayesNode (
+      varId, dsize, evidence, parents, dist));
+  return nodes_.back();
+}
+
+
+
+BayesNode*
+BayesNet::addNode (string label,
+                   Domain domain,
+                   NodeSet& parents, 
+                   ParamSet& params)
+{
+  indexMap_.insert (make_pair (nodes_.size(), nodes_.size()));
+  Distribution* dist = new Distribution (params);
+  BayesNode* node = new BayesNode (
+      nodes_.size(), label, domain, parents, dist);
+  dists_.push_back (dist);
+  nodes_.push_back (node);
+  return node;
+}
+
+
+
+BayesNode*
+BayesNet::getNode (unsigned varId) const
+{
+  IndexMap::const_iterator it = indexMap_.find(varId);
+  if (it == indexMap_.end()) {
+    return 0;
+  } else {
+    return nodes_[it->second];
+  }
+}
+
+
+
+BayesNode*
+BayesNet::getNode (string label) const
+{
+  BayesNode* node = 0;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    if (nodes_[i]->getLabel() == label) {
+      node = nodes_[i];
+      break;
+    }
+  }
+  return node;
+}
+
+
+
+void
+BayesNet::addDistribution (Distribution* dist)
+{
+  dists_.push_back (dist);
+}
+
+
+
+Distribution*
+BayesNet::getDistribution (unsigned distId) const
+{
+  Distribution* dist = 0;
+  for (unsigned i = 0; i < dists_.size(); i++) {
+    if (dists_[i]->id == distId) {
+      dist = dists_[i];
+      break;
+    }
+  }
+  return dist;
+}
+
+
+
+const NodeSet&
+BayesNet::getNodes (void) const
+{
+  return nodes_;
+}
+
+
+
+int
+BayesNet::getNumberOfNodes (void) const
+{
+  return nodes_.size();
+}
+
+
+
+NodeSet
+BayesNet::getRootNodes (void) const
+{
+  NodeSet roots;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    if (nodes_[i]->isRoot()) {
+      roots.push_back (nodes_[i]);
+    }
+  }
+  return roots;
+}
+
+
+
+NodeSet
+BayesNet::getLeafNodes (void) const
+{
+  NodeSet leafs;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    if (nodes_[i]->isLeaf()) {
+      leafs.push_back (nodes_[i]);
+    }
+  }
+  return leafs;
+}
+
+
+
+VarSet
+BayesNet::getVariables (void) const
+{
+  VarSet vars;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    vars.push_back (nodes_[i]);
+  }
+  return vars;
+}
+
+
+
+BayesNet*
+BayesNet::pruneNetwork (BayesNode* queryNode) const
+{
+  NodeSet queryNodes;
+  queryNodes.push_back (queryNode);
+  return pruneNetwork (queryNodes);
+}
+
+
+
+BayesNet*
+BayesNet::pruneNetwork (const NodeSet& interestedVars) const
+{
+  /*
+  cout << "interested vars: " ;
+  for (unsigned i = 0; i < interestedVars.size(); i++) {
+    cout << interestedVars[i]->getLabel() << " " ;
+  }
+  cout << endl;
+  */
+  vector<StateInfo*> states (nodes_.size(), 0);
+
+  Scheduling scheduling;
+  for (NodeSet::const_iterator it = interestedVars.begin();
+      it != interestedVars.end(); it++) {
+    scheduling.push (ScheduleInfo (*it, false, true));
+  }
+
+  while (!scheduling.empty()) {
+    ScheduleInfo& sch = scheduling.front();
+    StateInfo* state  = states[sch.node->getIndex()];
+    if (!state) {
+      state = new StateInfo();
+      states[sch.node->getIndex()] = state;
+    } else {
+      state->visited = true;
+    }
+    if (!sch.node->hasEvidence() && sch.visitedFromChild) {
+      if (!state->markedOnTop) {
+        state->markedOnTop = true;
+        scheduleParents (sch.node, scheduling);
+      }
+      if (!state->markedOnBottom) {
+        state->markedOnBottom = true;
+        scheduleChilds (sch.node, scheduling);
+      }
+    }
+    if (sch.visitedFromParent) {
+      if (sch.node->hasEvidence() && !state->markedOnTop) {
+        state->markedOnTop = true;
+        scheduleParents (sch.node, scheduling);
+      }
+      if (!sch.node->hasEvidence() && !state->markedOnBottom) {
+        state->markedOnBottom = true;
+        scheduleChilds (sch.node, scheduling);
+      }
+    }
+    scheduling.pop();
+  }
+  /*
+  cout << "\t\ttop\tbottom" << endl;
+  cout << "variable\t\tmarked\tmarked\tvisited\tobserved" << endl;
+  cout << "----------------------------------------------------------" ;
+  cout << endl;
+  for (unsigned i = 0; i < states.size(); i++) {
+    cout << nodes_[i]->getLabel() << ":\t\t" ;
+    if (states[i]) {
+      states[i]->markedOnTop    ? cout << "yes\t" : cout << "no\t" ;
+      states[i]->markedOnBottom ? cout << "yes\t" : cout << "no\t" ;
+      states[i]->visited        ? cout << "yes\t" : cout << "no\t" ;
+      nodes_[i]->hasEvidence()  ? cout << "yes"   : cout << "no" ;
+      cout << endl;
+    } else {
+      cout << "no\tno\tno\t" ;
+      nodes_[i]->hasEvidence() ? cout << "yes" : cout << "no" ;
+      cout << endl;
+    }
+  }
+  cout << endl;
+  */
+  BayesNet* bn = new BayesNet();
+  constructGraph (bn, states);
+
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    delete states[i];
+  }
+  return bn;
+}
+
+
+
+void
+BayesNet::constructGraph (BayesNet* bn,
+                          const vector<StateInfo*>& states) const
+{
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    bool isRequired = false;
+    if (states[i]) {
+      isRequired = (nodes_[i]->hasEvidence() && states[i]->visited)
+                   || 
+                   states[i]->markedOnTop;
+    }
+    if (isRequired) {
+      NodeSet parents;
+      if (states[i]->markedOnTop) {
+        const NodeSet& ps = nodes_[i]->getParents();
+        for (unsigned j = 0; j < ps.size(); j++) {
+          BayesNode* parent = bn->getNode (ps[j]->getVarId());
+          if (!parent) {
+            parent = bn->addNode (ps[j]->getVarId());
+          }
+          parents.push_back (parent);
+        }
+      }
+      BayesNode* node = bn->getNode (nodes_[i]->getVarId());
+      if (node) {
+        node->setData      (nodes_[i]->getDomainSize(),
+                            nodes_[i]->getEvidence(), parents,
+                            nodes_[i]->getDistribution());
+      } else {
+        node = bn->addNode (nodes_[i]->getVarId(),
+                            nodes_[i]->getDomainSize(),
+                            nodes_[i]->getEvidence(), parents, 
+                            nodes_[i]->getDistribution());
+      }
+      if (nodes_[i]->hasDomain()) {
+        node->setDomain (nodes_[i]->getDomain());
+      }
+      if (nodes_[i]->hasLabel()) {
+        node->setLabel (nodes_[i]->getLabel());
+      }
+    }
+  }
+  bn->setIndexes();
+}
+
+/*
+void
+BayesNet::constructGraph (BayesNet* bn,
+                          const vector<StateInfo*>& states) const
+{
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    if (states[i]) {
+      if (nodes_[i]->hasEvidence() && states[i]->visited) {
+        NodeSet parents;
+        if (states[i]->markedOnTop) {
+          const NodeSet& ps = nodes_[i]->getParents();
+          for (unsigned j = 0; j < ps.size(); j++) {
+            BayesNode* parent = bn->getNode (ps[j]->getVarId());
+            if (parent == 0) {
+              parent = bn->addNode (ps[j]->getVarId());
+            }
+            parents.push_back (parent);
+          }
+        }
+
+        BayesNode* n = bn->getNode (nodes_[i]->getVarId());
+        if (n) {
+          n->setData  (nodes_[i]->getDomainSize(),
+                       nodes_[i]->getEvidence(), parents,
+                       nodes_[i]->getDistribution());
+        } else {
+          bn->addNode (nodes_[i]->getVarId(),
+                       nodes_[i]->getDomainSize(),
+                       nodes_[i]->getEvidence(), parents, 
+                       nodes_[i]->getDistribution());
+        }
+        
+      } else if (states[i]->markedOnTop) {
+        NodeSet parents;
+        const NodeSet& ps = nodes_[i]->getParents();
+        for (unsigned j = 0; j < ps.size(); j++) {
+          BayesNode* parent = bn->getNode (ps[j]->getVarId());
+          if (parent == 0) {
+            parent = bn->addNode (ps[j]->getVarId());
+          }
+          parents.push_back (parent);
+        }
+
+        BayesNode* n = bn->getNode (nodes_[i]->getVarId());
+        if (n) {
+          n->setData  (nodes_[i]->getDomainSize(),
+                       nodes_[i]->getEvidence(), parents,
+                       nodes_[i]->getDistribution());
+        } else {
+          bn->addNode (nodes_[i]->getVarId(),
+                       nodes_[i]->getDomainSize(),
+                       nodes_[i]->getEvidence(), parents,
+                       nodes_[i]->getDistribution());
+        }
+      }
+    }
+  }
+}*/
+
+
+
+bool
+BayesNet::isSingleConnected (void) const
+{
+  return !containsUndirectedCycle();
+}
+
+
+
+vector<DomainConf>
+BayesNet::getDomainConfigurationsOf (const NodeSet& nodes)
+{
+  int nConfs = 1;
+  for (unsigned i = 0; i < nodes.size(); i++) {
+    nConfs *= nodes[i]->getDomainSize();
+  }
+
+  vector<DomainConf> confs (nConfs);
+  for (int i = 0; i < nConfs; i++) {
+    confs[i].resize (nodes.size());
+  }
+
+  int nReps = 1;
+  for (int i = nodes.size() - 1; i >= 0; i--) {
+    int index = 0;
+    while (index < nConfs) {
+      for (int j = 0; j < nodes[i]->getDomainSize(); j++) {
+        for (int r = 0; r < nReps; r++) {
+          confs[index][i] = j; 
+          index++;
+        }
+      }
+    }
+    nReps *= nodes[i]->getDomainSize();
+  }
+
+  return confs;
+}
+
+
+
+vector<string>
+BayesNet::getInstantiations (const NodeSet& parents_)
+{
+  int nParents = parents_.size();
+  int rowSize = 1;
+  for (unsigned i = 0; i < parents_.size(); i++) {
+    rowSize *= parents_[i]->getDomainSize();
+  }
+  int nReps     = 1;
+  vector<string> headers (rowSize);
+  for (int i = nParents - 1; i >= 0; i--) {
+    Domain domain = parents_[i]->getDomain();
+    int index = 0;
+    while (index < rowSize) {
+      for (int j = 0; j < parents_[i]->getDomainSize(); j++) {
+        for (int r = 0; r < nReps; r++) {
+          if (headers[index] != "") {
+            headers[index] = domain[j] + "," + headers[index];
+          } else {
+            headers[index] = domain[j];
+          }
+          index++;
+        }
+      }
+    }
+    nReps *= parents_[i]->getDomainSize();
+  }
+  return headers;
+}
+
+
+
+void
+BayesNet::setIndexes (void)
+{
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    nodes_[i]->setIndex (i);
+  }
+}
+
+
+
+void
+BayesNet::freeDistributions (void)
+{
+  for (unsigned i = 0; i < dists_.size(); i++) {
+    delete dists_[i];
+  }
+}
+
+
+
+void
+BayesNet::printNetwork (void) const
+{
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    cout << *nodes_[i];
+  }
+}
+
+
+
+void
+BayesNet::printNetworkToFile (const char* fileName) const
+{
+  string s = "../../" ;
+  s += fileName;
+  ofstream out (s.c_str());
+  if (!out.is_open()) {
+    cerr << "error: cannot open file to write at " ;
+    cerr << "BayesNet::printToFile()" << endl;
+    abort();
+  }
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    out << *nodes_[i];
+  }
+  out.close();
+}
+ 
+
+
+void
+BayesNet::exportToDotFile (const char* fileName,
+                           bool showNeighborless,
+                           const NodeSet& highlightNodes) const
+{
+  string s = "../../" ;
+  s+= fileName;
+  ofstream out (s.c_str());
+  if (!out.is_open()) {
+    cerr << "error: cannot open file to write at " ;
+    cerr << "BayesNet::exportToDotFile()" << endl;
+    abort();
+  }
+
+  out << "digraph \"" << fileName << "\" {" << endl;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    const NodeSet& childs = nodes_[i]->getChilds();
+    for (unsigned j = 0; j < childs.size(); j++) {
+      out << '"' << nodes_[i]->getLabel() << '"' << " -> " ;
+      out << '"' << childs[j]->getLabel() << '"' << endl;
+    }
+  }
+
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    if (showNeighborless || nodes_[i]->hasNeighbors()) {
+      out << '"' << nodes_[i]->getLabel() << '"' ;
+      if (nodes_[i]->hasEvidence()) {
+        out << " [style=filled, fillcolor=yellow]" << endl;
+      } else {
+        out << endl;
+      }
+    }
+  }
+
+  for (unsigned i = 0; i < highlightNodes.size(); i++) {
+    out << '"' << highlightNodes[i]->getLabel() << '"' ;
+    out << " [shape=box]" << endl;
+  }
+
+  out << "}" << endl;
+  out.close();
+}
+
+
+
+void
+BayesNet::exportToBifFile (const char* fileName) const
+{
+  string s = "../../" ;
+  s += fileName;
+  ofstream out (s.c_str());
+  if(!out.is_open()) {
+    cerr << "error: cannot open file to write at " ;
+    cerr << "BayesNet::exportToBifFile()" << endl;
+    abort();
+  }
+  out << "<?xml version=\"1.0\" encoding=\"US-ASCII\"?>" << endl;
+  out << "<BIF VERSION=\"0.3\">" << endl;
+  out << "<NETWORK>" << endl;
+  out << "<NAME>" << fileName << "</NAME>" << endl << endl;
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    out << "<VARIABLE TYPE=\"nature\">" << endl;
+    out << "\t<NAME>" << nodes_[i]->getLabel() << "</NAME>" << endl;
+    const Domain& domain = nodes_[i]->getDomain();
+    for (unsigned j = 0; j < domain.size(); j++) {
+      out << "\t<OUTCOME>" << domain[j] << "</OUTCOME>" << endl;
+    }
+    out << "</VARIABLE>" << endl << endl;
+  }
+
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    out << "<DEFINITION>" << endl;
+    out << "\t<FOR>" << nodes_[i]->getLabel() << "</FOR>" << endl;
+    const NodeSet& parents = nodes_[i]->getParents();
+    for (unsigned j = 0; j < parents.size(); j++) {
+      out << "\t<GIVEN>" << parents[j]->getLabel();
+      out << "</GIVEN>" << endl;
+    }
+    ParamSet params = revertParameterReorder (nodes_[i]->getParameters(),
+                                              nodes_[i]->getDomainSize());
+    out << "\t<TABLE>" ;
+    for (unsigned j = 0; j < params.size(); j++) {
+       out << " " << params[j];
+    }
+    out << " </TABLE>" << endl;
+    out << "</DEFINITION>" << endl << endl;
+  }
+  out << "</NETWORK>" << endl;
+  out << "</BIF>" << endl << endl;
+  out.close();  
+}
+
+
+
+bool
+BayesNet::containsUndirectedCycle (void) const
+{
+  vector<bool> visited (nodes_.size(), false);
+  for (unsigned i = 0; i < nodes_.size(); i++) {
+    int v = nodes_[i]->getIndex();
+    if (!visited[v]) {
+      if (containsUndirectedCycle (v, -1, visited)) {
+        return true;
+      }
+    }
+  }
+  return false;
+}
+
+
+
+bool
+BayesNet::containsUndirectedCycle (int v,
+                                   int p,
+                                   vector<bool>& visited) const
+{
+  visited[v] = true;
+  vector<int> adjacencies = getAdjacentNodes (v);
+  for (unsigned i = 0; i < adjacencies.size(); i++) {
+    int w = adjacencies[i];
+    if (!visited[w]) {
+      if (containsUndirectedCycle (w, v, visited)) {
+        return true;
+      }
+    }
+    else if (visited[w] && w != p) {
+      return true;
+    }
+  }
+  return false; // no cycle detected in this component
+}
+
+
+
+vector<int>
+BayesNet::getAdjacentNodes (int v) const
+{
+  vector<int> adjacencies;
+  const NodeSet& parents = nodes_[v]->getParents();
+  const NodeSet& childs  = nodes_[v]->getChilds();
+  for (unsigned i = 0; i < parents.size(); i++) {
+    adjacencies.push_back (parents[i]->getIndex());
+  }
+  for (unsigned i = 0; i < childs.size(); i++) {
+    adjacencies.push_back (childs[i]->getIndex());
+  }
+  return adjacencies;
+}
+
+
+
+ParamSet
+BayesNet::reorderParameters (const ParamSet& params,
+                             int domainSize) const
+{
+  // the interchange format for bayesian networks keeps the probabilities 
+  // in the following order:
+  // p(a1|b1,c1) p(a2|b1,c1) p(a1|b1,c2) p(a2|b1,c2) p(a1|b2,c1) p(a2|b2,c1)
+  // p(a1|b2,c2) p(a2|b2,c2).
+  //
+  // however, in clpbn we keep the probabilities in this order:
+  // p(a1|b1,c1) p(a1|b1,c2) p(a1|b2,c1) p(a1|b2,c2) p(a2|b1,c1) p(a2|b1,c2)
+  // p(a2|b2,c1) p(a2|b2,c2).
+  unsigned count    = 0;
+  unsigned rowSize  = params.size() / domainSize;
+  ParamSet reordered;
+  while (reordered.size() < params.size()) {
+    unsigned idx = count;
+    for (unsigned i = 0; i < rowSize; i++) {
+      reordered.push_back (params[idx]);
+      idx += domainSize;
+    }
+    count++;
+  }
+  return reordered;
+}
+
+
+
+ParamSet
+BayesNet::revertParameterReorder (const ParamSet& params,
+                                  int domainSize) const
+{
+  unsigned count    = 0;
+  unsigned rowSize  = params.size() / domainSize;
+  ParamSet reordered;
+  while (reordered.size() < params.size()) {
+    unsigned idx = count;
+    for (int i = 0; i < domainSize; i++) {
+      reordered.push_back (params[idx]);
+      idx += rowSize;
+    }
+    count ++;
+  }
+  return reordered;
+}
+

packages/CLPBN/clpbn/bp/BayesNet.h

@@ -0,0 +1,129 @@
+#ifndef BP_BAYES_NET_H
+#define BP_BAYES_NET_H
+
+#include <vector>
+#include <queue>
+#include <list>
+#include <string>
+#include <unordered_map>
+#include <map>
+
+#include "GraphicalModel.h"
+#include "BayesNode.h"
+#include "Shared.h"
+
+
+using namespace std;
+
+class Distribution;
+
+struct ScheduleInfo
+{
+  ScheduleInfo (BayesNode* n, bool vfp, bool vfc)
+  {
+    node               = n;
+    visitedFromParent  = vfp;
+    visitedFromChild   = vfc;
+  }
+  BayesNode*  node;
+  bool        visitedFromParent;
+  bool        visitedFromChild;
+};
+
+
+struct StateInfo
+{
+  StateInfo (void)
+  {
+    visited         = true;
+    markedOnTop     = false;
+    markedOnBottom  = false;
+  }
+  bool visited;
+  bool markedOnTop;
+  bool markedOnBottom;
+};
+
+typedef vector<Distribution*>  DistSet;
+typedef queue<ScheduleInfo, list<ScheduleInfo> > Scheduling;
+typedef unordered_map<unsigned, unsigned> Histogram;
+typedef unordered_map<unsigned, double>   Times;
+
+
+class BayesNet : public GraphicalModel
+{
+  public:
+    BayesNet (void);
+    BayesNet (const char*);
+   ~BayesNet (void);
+
+    BayesNode*           addNode (unsigned);
+    BayesNode*           addNode (unsigned, unsigned, int, NodeSet&, Distribution*);
+    BayesNode*           addNode (string, Domain, NodeSet&, ParamSet&);
+    BayesNode*           getNode (unsigned) const;
+    BayesNode*           getNode (string) const;
+    void                 addDistribution (Distribution*);
+    Distribution*        getDistribution (unsigned) const;
+    const NodeSet&       getNodes (void) const;
+    int                  getNumberOfNodes (void) const;
+    NodeSet              getRootNodes (void) const;
+    NodeSet              getLeafNodes (void) const;
+    VarSet               getVariables (void) const;
+    BayesNet*            pruneNetwork (BayesNode*) const;
+    BayesNet*            pruneNetwork (const NodeSet& queryNodes) const;
+    void                 constructGraph (BayesNet*, const vector<StateInfo*>&) const;
+    bool                 isSingleConnected (void) const;
+    static vector<DomainConf> getDomainConfigurationsOf (const NodeSet&);
+    static vector<string>     getInstantiations (const NodeSet& nodes);
+    void                 setIndexes (void);
+    void                 freeDistributions (void);
+    void                 printNetwork (void) const;
+    void                 printNetworkToFile (const char*) const;
+    void                 exportToDotFile (const char*, bool = true,
+                                          const NodeSet& = NodeSet()) const;
+    void                 exportToBifFile (const char*) const;
+
+    static Histogram     histogram_;
+    static Times         times_;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (BayesNet);
+
+    bool                 containsUndirectedCycle (void) const;
+    bool                 containsUndirectedCycle (int, int,
+                                                  vector<bool>&)const;
+    vector<int>          getAdjacentNodes (int) const ;
+    ParamSet             reorderParameters (const ParamSet&, int) const;
+    ParamSet             revertParameterReorder (const ParamSet&, int) const;
+    void                 scheduleParents (const BayesNode*, Scheduling&) const;
+    void                 scheduleChilds (const BayesNode*, Scheduling&) const;
+
+    NodeSet              nodes_;
+    DistSet              dists_;
+    IndexMap             indexMap_;
+};
+
+
+
+inline void
+BayesNet::scheduleParents (const BayesNode* n, Scheduling& sch) const
+{
+  const NodeSet& ps = n->getParents();
+  for (NodeSet::const_iterator it = ps.begin(); it != ps.end(); it++) {
+    sch.push (ScheduleInfo (*it, false, true));
+  }
+}
+
+
+
+inline void
+BayesNet::scheduleChilds (const BayesNode* n, Scheduling& sch) const
+{
+  const NodeSet& cs = n->getChilds();
+  for (NodeSet::const_iterator it = cs.begin(); it != cs.end(); it++) {
+    sch.push (ScheduleInfo (*it, true, false));
+  }
+}
+
+#endif
+

packages/CLPBN/clpbn/bp/BayesNode.cpp

@@ -0,0 +1,355 @@
+#include <iostream>
+#include <sstream>
+#include <iomanip>
+#include <cassert>
+#include <cstdlib>
+
+#include "BayesNode.h"
+
+
+BayesNode::BayesNode (unsigned varId) : Variable (varId)
+{
+}
+
+
+
+BayesNode::BayesNode (unsigned varId,
+                      unsigned dsize,
+                      int evidence,
+                      const NodeSet& parents,
+                      Distribution* dist) : Variable(varId, dsize, evidence)
+{
+  parents_   = parents;
+  dist_      = dist;
+  for (unsigned int i = 0; i < parents.size(); i++) {
+    parents[i]->addChild (this);
+  }
+}
+
+
+
+BayesNode::BayesNode (unsigned varId,
+                      string label,
+                      const Domain& domain,
+                      const NodeSet& parents,
+                      Distribution* dist) : Variable(varId, domain)
+{
+  label_     = new string (label);
+  parents_   = parents;
+  dist_      = dist;
+  for (unsigned int i = 0; i < parents.size(); i++) {
+    parents[i]->addChild (this);
+  }
+}
+
+
+
+void
+BayesNode::setData (unsigned dsize,
+                    int evidence,
+                    const NodeSet& parents,
+                    Distribution* dist)
+{
+  setDomainSize (dsize);
+  evidence_  = evidence;
+  parents_   = parents;
+  dist_      = dist;
+  for (unsigned int i = 0; i < parents.size(); i++) {
+    parents[i]->addChild (this);
+  } 
+}
+
+
+
+void 
+BayesNode::addChild (BayesNode* node)
+{
+  childs_.push_back (node);
+}
+
+
+
+Distribution*
+BayesNode::getDistribution (void)
+{
+  return dist_;
+}
+
+
+
+const ParamSet&
+BayesNode::getParameters (void)
+{
+  return dist_->params;
+}
+
+
+
+ParamSet
+BayesNode::getRow (int rowIndex) const
+{
+  int rowSize = getRowSize();
+  int offset  = rowSize * rowIndex;
+  ParamSet row (rowSize);
+  for (int i = 0; i < rowSize; i++) {
+    row[i] = dist_->params[offset + i] ;
+  }
+  return row;
+}
+
+
+
+bool 
+BayesNode::isRoot (void)
+{
+  return getParents().empty();
+}
+
+
+
+bool 
+BayesNode::isLeaf (void)
+{
+  return getChilds().empty();
+}
+
+
+
+bool
+BayesNode::hasNeighbors (void) const
+{
+  return childs_.size() != 0 || parents_.size() != 0;
+}
+
+
+int 
+BayesNode::getCptSize (void)
+{
+  return dist_->params.size();
+}
+
+
+
+const vector<CptEntry>&
+BayesNode::getCptEntries (void)
+{
+  if (dist_->entries.size() == 0) {
+    unsigned rowSize  = getRowSize();
+    unsigned nParents = parents_.size();
+    vector<DomainConf> confs (rowSize);
+
+    for (unsigned i = 0; i < rowSize; i++) {  
+      confs[i].resize (nParents);
+    }
+
+    int nReps = 1;
+    for (int i = nParents - 1; i >= 0; i--) {
+      unsigned index = 0;
+      while (index < rowSize) {
+        for (int j = 0; j < parents_[i]->getDomainSize(); j++) {
+          for (int r = 0; r < nReps; r++) {
+            confs[index][i] = j;
+            index++;
+          }
+        }
+      }
+      nReps *= parents_[i]->getDomainSize();
+    }
+  
+    dist_->entries.reserve (rowSize);
+    for (unsigned i = 0; i < rowSize; i++) {
+      dist_->entries.push_back (CptEntry (i, confs[i]));
+    }
+  }
+  return dist_->entries;
+}
+
+
+
+int
+BayesNode::getIndexOfParent (const BayesNode* parent) const
+{
+  for (unsigned int i = 0; i < parents_.size(); i++) {
+    if (parents_[i] == parent) {
+      return i;
+    }
+  }
+  return -1;
+}
+
+
+
+string
+BayesNode::cptEntryToString (const CptEntry& entry) const
+{
+  stringstream ss;
+  ss << "p(" ;
+  const DomainConf& conf = entry.getParentConfigurations();
+  int row = entry.getParameterIndex() / getRowSize();
+  ss << getDomain()[row]; 
+  if (parents_.size() > 0) {
+    ss << "|" ;
+    for (unsigned int i = 0; i < conf.size(); i++) {
+      if (i != 0) {
+        ss << ",";
+      }
+      ss << parents_[i]->getDomain()[conf[i]];
+    }
+  }
+  ss << ")" ;
+  return ss.str();
+}
+
+
+
+string
+BayesNode::cptEntryToString (int row, const CptEntry& entry) const
+{
+  stringstream ss;
+  ss << "p(" ;
+  const DomainConf& conf = entry.getParentConfigurations();
+  ss << getDomain()[row]; 
+  if (parents_.size() > 0) {
+    ss << "|" ;
+    for (unsigned int i = 0; i < conf.size(); i++) {
+      if (i != 0) {
+        ss << ",";
+      }
+      ss << parents_[i]->getDomain()[conf[i]];
+    }
+  }
+  ss << ")" ;
+  return ss.str();
+}
+
+
+
+vector<string> 
+BayesNode::getDomainHeaders (void) const
+{
+  int nParents = parents_.size();
+  int rowSize  = getRowSize(); 
+  int nReps     = 1;
+  vector<string> headers (rowSize);
+  for (int i = nParents - 1; i >= 0; i--) {
+    Domain domain = parents_[i]->getDomain();
+    int index = 0;
+    while (index < rowSize) {
+      for (int j = 0; j < parents_[i]->getDomainSize(); j++) {
+        for (int r = 0; r < nReps; r++) {
+          if (headers[index] != "") {
+            headers[index] = domain[j] + "," + headers[index];
+          } else {
+            headers[index] = domain[j];
+          }
+          index++;
+        }
+      }
+    }
+    nReps *= parents_[i]->getDomainSize();
+  }
+  return headers;
+}
+
+
+
+ostream& 
+operator << (ostream& o, const BayesNode& node)
+{
+  o << "variable " << node.getIndex() << endl;
+  o << "Var Id:   " << node.getVarId() << endl;
+  o << "Label:    " << node.getLabel() << endl;
+  
+  o << "Evidence: " ;
+  if (node.hasEvidence()) {
+    o << node.getEvidence();
+  }
+  else {
+    o << "no" ;
+  }
+  o << endl;
+
+  o << "Parents:  " ;
+  const NodeSet& parents = node.getParents();
+  if (parents.size() != 0) {
+    for (unsigned int i = 0; i < parents.size() - 1; i++) {
+      o << parents[i]->getLabel() << ", " ;
+    }
+    o << parents[parents.size() - 1]->getLabel();
+  }
+  o << endl;
+  
+  o << "Childs:   " ;
+  const NodeSet& childs = node.getChilds();
+  if (childs.size() != 0) {
+    for (unsigned int i = 0; i < childs.size() - 1; i++) {
+      o << childs[i]->getLabel() << ", " ;
+    }
+    o << childs[childs.size() - 1]->getLabel();
+  }
+  o << endl;
+
+  o << "Domain:   " ;
+  Domain domain = node.getDomain();
+  for (unsigned int i = 0; i < domain.size() - 1; i++) {
+    o << domain[i] << ", " ;
+  }
+  if (domain.size() != 0) {
+    o << domain[domain.size() - 1];
+  }
+  o << endl;
+
+  // min width of first column
+  const unsigned int MIN_DOMAIN_WIDTH = 4;
+  // min width of following columns
+  const unsigned int MIN_COMBO_WIDTH = 12;
+
+  unsigned int domainWidth = domain[0].length();
+  for (unsigned int i = 1; i < domain.size(); i++) {
+    if (domain[i].length() > domainWidth) {
+      domainWidth = domain[i].length();
+    }
+  }
+  domainWidth = (domainWidth < MIN_DOMAIN_WIDTH) 
+              ?  MIN_DOMAIN_WIDTH 
+              :  domainWidth;
+
+  o << left << setw (domainWidth) << "cpt" << right;
+
+  vector<int> widths;
+  int lineWidth = domainWidth;
+  vector<string> headers = node.getDomainHeaders();
+
+  if (!headers.empty()) {
+    for (unsigned int i = 0; i < headers.size(); i++) {
+      unsigned int len = headers[i].length();
+      int w = (len < MIN_COMBO_WIDTH) ? MIN_COMBO_WIDTH : len;
+      widths.push_back (w);
+      o << setw (w) << headers[i];
+      lineWidth += w;
+    }
+    o << endl;
+  } else {
+    cout << endl;
+    widths.push_back (domainWidth);
+    lineWidth += MIN_COMBO_WIDTH;
+  }
+
+  for (int i = 0; i < lineWidth; i++) {
+    o << "-" ;
+  }
+  o << endl;
+
+  for (unsigned int i = 0; i < domain.size(); i++) {   
+    ParamSet row = node.getRow (i);
+    o << left << setw (domainWidth) << domain[i] << right;
+    for (unsigned j = 0; j < node.getRowSize(); j++) {
+      o << setw (widths[j]) << row[j];
+    }
+    o << endl;
+  }
+  o << endl;
+
+  return o;
+}
+

packages/CLPBN/clpbn/bp/BayesNode.h

@@ -0,0 +1,91 @@
+#ifndef BP_BAYESNODE_H
+#define BP_BAYESNODE_H
+
+#include <vector>
+#include <string>
+#include <sstream>
+
+#include "Variable.h"
+#include "CptEntry.h"
+#include "Distribution.h"
+#include "Shared.h"
+
+using namespace std;
+
+
+class BayesNode : public Variable
+{
+  public:
+    BayesNode (unsigned);
+    BayesNode (unsigned, unsigned, int, const NodeSet&, Distribution*);
+    BayesNode (unsigned, string, const Domain&, const NodeSet&, Distribution*);
+
+    void                     setData (unsigned, int, const NodeSet&, Distribution*);
+    void                     addChild (BayesNode*);
+    Distribution*            getDistribution (void);
+    const ParamSet&          getParameters (void);
+    ParamSet                 getRow (int) const;
+    void                     setProbability (int, const CptEntry&, double);
+    bool                     isRoot (void);
+    bool                     isLeaf (void);
+    bool                     hasNeighbors (void) const;
+    int                      getCptSize (void);
+    const vector<CptEntry>&  getCptEntries (void);
+    int                      getIndexOfParent (const BayesNode*) const;
+    string                   cptEntryToString (const CptEntry&) const;
+    string                   cptEntryToString (int, const CptEntry&) const;
+    // inlines
+    const NodeSet&           getParents (void) const;
+    const NodeSet&           getChilds (void) const;
+    double                   getProbability (int, const CptEntry& entry);
+    unsigned                 getRowSize (void) const;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (BayesNode);
+
+    Domain                   getDomainHeaders (void) const;
+    friend ostream&          operator << (ostream&, const BayesNode&);
+
+    NodeSet                  parents_;
+    NodeSet                  childs_;
+    Distribution*            dist_;
+};
+
+ostream& operator << (ostream&, const BayesNode&);
+
+
+
+inline const NodeSet&
+BayesNode::getParents (void) const
+{
+  return parents_;
+}
+
+
+
+inline const NodeSet&
+BayesNode::getChilds (void) const
+{
+  return childs_;
+}
+
+
+
+inline double
+BayesNode::getProbability (int row, const CptEntry& entry)
+{
+  int col = entry.getParameterIndex();
+  int idx = (row * getRowSize()) + col;
+  return dist_->params[idx];
+}
+
+
+
+inline unsigned 
+BayesNode::getRowSize (void) const
+{ 
+  return dist_->params.size() / getDomainSize();
+}
+
+#endif
+

packages/CLPBN/clpbn/bp/BpNode.cpp

@@ -4,226 +4,44 @@
 
 #include "BpNode.h"
 
+bool BpNode::calculateMessageResidual_ = true;
 
-bool BpNode::parallelSchedule_ = false;
 
-BpNode::BpNode (string varName, 
+BpNode::BpNode (BayesNode* node)
-                vector<BayesianNode*> parents,
-                Distribution* dist,
-                int evidence) : BayesianNode (varName, parents, dist, evidence)
 {
-
+  ds_ = node->getDomainSize();
-}
+  const NodeSet& childs = node->getChilds();
-
+  piVals_.resize (ds_, 1);
-
+  ldVals_.resize (ds_, 1);
-
+  if (calculateMessageResidual_) {
-BpNode::~BpNode (void)
+    piResiduals_.resize (childs.size(), 0.0);
-{
+    ldResiduals_.resize (childs.size(), 0.0);
-  delete [] piValues_;
-  delete [] lambdaValues_;
-  delete [] oldBeliefs_;
-  map<BpNode*, double*>::iterator iter;
-  for (iter = lambdaMessages_.begin(); iter != lambdaMessages_.end(); ++iter) {
-    delete [] iter->second;
-  }
-  for (iter = piMessages_.begin(); iter != piMessages_.end(); ++iter) {
-    delete [] iter->second;
-  }
-  // FIXME delete new messages
-}
-
-
-
-void
-BpNode::enableParallelSchedule (void)
-{
-  parallelSchedule_ = true;
-}
-
-
-
-void 
-BpNode::allocateMemory (void)
-{
-  // FIXME do i need this !?
-  int domainSize = getDomainSize();
-  piValues_      = new double [domainSize];
-  lambdaValues_  = new double [domainSize];
-  if (parallelSchedule_) {
-    newPiMessages_     = new map<BpNode*, double*>;
-    newLambdaMessages_ = new map<BpNode*, double*>;
-  }
-  oldBeliefs_ = 0;
-  vector <BayesianNode*> childs = getChilds();
-  for (unsigned int i = 0; i < childs.size(); i++) {
-    BpNode* child = static_cast<BpNode*> (childs[i]);
-    piMessages_.insert     (make_pair (child, new double [domainSize]));
-    lambdaMessages_.insert (make_pair (child, new double [domainSize]));
-    if (parallelSchedule_) {
-      newPiMessages_->insert (make_pair (child, new double [domainSize]));
-      newLambdaMessages_->insert (make_pair (child, new double [domainSize]));
   }
+  childs_ = &childs;
+  for (unsigned i = 0; i < childs.size(); i++) {
+    //indexMap_.insert (make_pair (childs[i]->getVarId(), i));
+    currPiMsgs_.push_back (ParamSet (ds_, 1));
+    currLdMsgs_.push_back (ParamSet (ds_, 1));
+    nextPiMsgs_.push_back (ParamSet (ds_, 1));
+    nextLdMsgs_.push_back (ParamSet (ds_, 1));
   }
 }
 
 
 
-double* 
+ParamSet
-BpNode::getPiValues (void) const
+BpNode::getBeliefs (void) const
-{
-  return piValues_;
-}
-
-
-
-double 
-BpNode::getPiValue (int index) const
-{
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  return piValues_[index];
-}
-
-
-
-void 
-BpNode::setPiValue (int index, double value)
-{
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  piValues_[index] = value;
-}
-
-
-
-double* 
-BpNode::getLambdaValues (void) const
-{
-  return lambdaValues_;
-}
-
-
-
-double 
-BpNode::getLambdaValue (int index) const
-{
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  return lambdaValues_[index];
-}
-
-
-
-void 
-BpNode::setLambdaValue (int index, double value)
-{
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  lambdaValues_[index] = value;
-}
-
-
-
-double*
-BpNode::getPiMessages (BpNode* node) const
-{
-  assert (node);
-  map<BpNode*, double*>::const_iterator iter = piMessages_.find (node);
-  assert (iter != piMessages_.end());
-  return iter->second;
-}
-
-
-
-double 
-BpNode::getPiMessage (BpNode* node, int index) const
-{
-  assert (node);
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  map<BpNode*, double*>::const_iterator iter = piMessages_.find (node);
-  assert (iter != piMessages_.end());
-  return iter->second[index];
-}
-
-
-
-void 
-BpNode::setPiMessage (BpNode* node, int index, double probability)
-{
-  assert (node);
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  map<BpNode*, double*>::const_iterator iter;
-  if (parallelSchedule_) {
-    // cerr << "set_pi_message" << endl;
-    iter = newPiMessages_->find (node); 
-    assert (iter != newPiMessages_->end());
-  } else {
-    iter = piMessages_.find (node);
-    assert (iter != piMessages_.end());
-  }
-  iter->second[index] = probability;
-}
-
-
-
-double*
-BpNode::getLambdaMessages (BpNode* node) const
-{
-  assert (node);
-  map<BpNode*, double*>::const_iterator iter = lambdaMessages_.find (node);
-  assert (iter != piMessages_.end());
-  return iter->second;
-}
-
-
-
-double 
-BpNode::getLambdaMessage (BpNode* node, int index) const
-{
-  assert (node);
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  map<BpNode*, double*>::const_iterator iter = lambdaMessages_.find (node);
-  assert (iter != piMessages_.end());
-  return iter->second[index];
-}
-
-
-
-void 
-BpNode::setLambdaMessage (BpNode* node, int index, double probability)
-{
-  assert (node);
-  const int c = getDomainSize(); 
-  assert (index >=0 && index < c);
-  map<BpNode*, double*>::const_iterator iter;
-  if (parallelSchedule_) {
-    //cerr << "set_lambda_message" << endl;
-    iter = newLambdaMessages_->find (node);
-    assert (iter != newLambdaMessages_->end());
-  } else {
-    iter = lambdaMessages_.find (node);
-    assert (iter != lambdaMessages_.end());
-  }
-  iter->second[index] = probability;
-}
-
-
-
-double*
-BpNode::getBeliefs (void)
 {
   double sum = 0.0;
-  double* beliefs = new double [getDomainSize()];
+  ParamSet beliefs (ds_);
-  for (int xi = 0; xi < getDomainSize(); xi++) {
+  for (int xi = 0; xi < ds_; xi++) {
-    double prod = piValues_[xi] * lambdaValues_[xi];
+    double prod = piVals_[xi] * ldVals_[xi];
     beliefs[xi] = prod;
     sum += prod;
   }
-  // normalize the beliefs
+  assert (sum);
-  for (int xi = 0; xi < getDomainSize(); xi++) {
+  //normalize the beliefs
+  for (int xi = 0; xi < ds_; xi++) {
     beliefs[xi] /= sum;
   }
   return beliefs;
@@ -231,91 +49,202 @@ BpNode::getBeliefs (void)
 
 
 
+double
+BpNode::getPiValue (int idx) const
+{
+  assert (idx >=0 && idx < ds_);
+  return piVals_[idx];
+}
+
+
+
+void
+BpNode::setPiValue (int idx, double value)
+{
+  assert (idx >=0 && idx < ds_);
+  piVals_[idx] = value;
+}
+
+
+
+double 
+BpNode::getLambdaValue (int idx) const
+{
+  assert (idx >=0 && idx < ds_);
+  return ldVals_[idx];
+}
+
+
+
+void 
+BpNode::setLambdaValue (int idx, double value)
+{
+  assert (idx >=0 && idx < ds_);
+  ldVals_[idx] = value;
+}
+
+
+
+ParamSet&
+BpNode::getPiValues (void)
+{
+  return piVals_;
+}
+
+
+
+ParamSet&
+BpNode::getLambdaValues (void)
+{
+  return ldVals_;
+}
+
+
+
+double 
+BpNode::getPiMessageValue (const BayesNode* destination, int idx) const
+{
+  assert (idx >=0 && idx < ds_);
+  return currPiMsgs_[getIndex(destination)][idx];
+}
+
+
+
+double 
+BpNode::getLambdaMessageValue (const BayesNode* source, int idx) const
+{
+  assert (idx >=0 && idx < ds_);
+  return currLdMsgs_[getIndex(source)][idx];
+}
+
+
+
+const ParamSet&
+BpNode::getPiMessage (const BayesNode* destination) const
+{
+  return currPiMsgs_[getIndex(destination)];
+}
+
+
+
+const ParamSet&
+BpNode::getLambdaMessage (const BayesNode* source) const
+{
+  return currLdMsgs_[getIndex(source)];
+}
+
+
+
+ParamSet&
+BpNode::piNextMessageReference (const BayesNode* destination)
+{
+  return nextPiMsgs_[getIndex(destination)];
+}
+
+
+
+ParamSet&
+BpNode::lambdaNextMessageReference (const BayesNode* source)
+{
+  return nextLdMsgs_[getIndex(source)];
+}
+
+
+
+void
+BpNode::updatePiMessage (const BayesNode* destination)
+{
+  int idx = getIndex (destination);
+  currPiMsgs_[idx] = nextPiMsgs_[idx];
+  Util::normalize (currPiMsgs_[idx]);
+}
+
+
+
+void
+BpNode::updateLambdaMessage (const BayesNode* source)
+{
+  int idx = getIndex (source);
+  currLdMsgs_[idx] = nextLdMsgs_[idx];
+  Util::normalize (currLdMsgs_[idx]);
+}
+
+
+
 double
 BpNode::getBeliefChange (void)
 {
   double change = 0.0;
-  if (!oldBeliefs_) {
+  if (oldBeliefs_.size() == 0) {
     oldBeliefs_ = getBeliefs();
-    change = MAX_CHANGE_;
+    change = 9999999999.0;
   } else {
-    double* currentBeliefs = getBeliefs();
+    ParamSet currentBeliefs = getBeliefs();
-    for (int xi = 0; xi < getDomainSize(); xi++) {
+    for (int xi = 0; xi < ds_; xi++) {
       change += abs (currentBeliefs[xi] - oldBeliefs_[xi]);
     }
     oldBeliefs_ = currentBeliefs;
   }
-  //FIXME memory leak
   return change;
 }
 
 
 
 void
-BpNode::normalizeMessages (void)
+BpNode::updatePiResidual (const BayesNode* destination)
 {
-  map<BpNode*, double*>::iterator iter;
+  int idx = getIndex (destination);
-
+  Util::normalize (nextPiMsgs_[idx]);
-  iter = lambdaMessages_.begin();
+  //piResiduals_[idx] = Util::getL1dist (
-  while (iter != lambdaMessages_.end()) {
+  //    currPiMsgs_[idx], nextPiMsgs_[idx]);
-    double* v = iter->second;
+  piResiduals_[idx] = Util::getMaxNorm (
-    double sum = 0.0;
+      currPiMsgs_[idx], nextPiMsgs_[idx]);
-    for (int xi = 0; xi < getDomainSize(); xi++) {
-      sum += v[xi];
-    }
-    for (int xi = 0; xi < getDomainSize(); xi++) {
-      v[xi] /= sum;
-    }
-    iter ++;
-  }
-  
-  iter = piMessages_.begin();
-  while (iter != piMessages_.end()) {
-    double* v = iter->second;
-    double sum = 0.0;
-    for (int xi = 0; xi < getDomainSize(); xi++) {
-      sum += v[xi];
-    }
-    for (int xi = 0; xi < getDomainSize(); xi++) {
-      v[xi] /= sum;
-    }
-    iter ++;
-  }
 }
 
 
 
 void
-BpNode::swapMessages (void)
+BpNode::updateLambdaResidual (const BayesNode* source)
 {
-  //FIXME fast way to do this
+  int idx = getIndex (source);
-  map<BpNode*, double*>::iterator iter1;
+  Util::normalize (nextLdMsgs_[idx]);
-  map<BpNode*, double*>::iterator iter2;
+  //ldResiduals_[idx] = Util::getL1dist (
-
+  //    currLdMsgs_[idx], nextLdMsgs_[idx]);
-  iter1 = lambdaMessages_.begin();
+  ldResiduals_[idx] = Util::getMaxNorm (
-  iter2 = newLambdaMessages_->begin();
+      currLdMsgs_[idx], nextLdMsgs_[idx]);
-  while (iter1 != lambdaMessages_.end()) {
+}
-    double* v1 = iter1->second;
+
-    double* v2 = iter2->second;
+
-    for (int xi = 0; xi < getDomainSize(); xi++) {
+
-      //v1[xi] = v2[xi];
+void
-      v1[xi] = (v1[xi] + v2[xi]) / 2;
+BpNode::clearPiResidual (const BayesNode* destination)
-    }
+{
-    iter1 ++;
+  piResiduals_[getIndex(destination)] = 0;
-    iter2 ++;
+}
-  }
+
-  
+
-  iter1 = piMessages_.begin();
+
-  iter2 = newPiMessages_->begin();
+void
-  while (iter1 != piMessages_.end()) {
+BpNode::clearLambdaResidual (const BayesNode* source)
-    double* v1 = iter1->second;
+{
-    double* v2 = iter2->second;
+  ldResiduals_[getIndex(source)] = 0;
-    for (int xi = 0; xi < getDomainSize(); xi++) {
+}
-      //v1[xi] = v2[xi];
+
-      v1[xi] = (v1[xi] + v2[xi]) / 2;
+
-    }
+
-    iter1 ++;
+bool
-    iter2 ++;
+BpNode::hasReceivedChildInfluence (void) const
-  }
+{
+  // if all lambda values are equal, then neither
+  // this node neither its descendents have evidence,
+  // we can use this to don't send lambda messages his parents
+  bool childInfluenced = false;
+  for (int xi = 1; xi < ds_; xi++) {
+    if (ldVals_[xi] != ldVals_[0]) {
+      childInfluenced = true;
+      break;
+    }
+  }
+  return childInfluenced;
 }
 

packages/CLPBN/clpbn/bp/BpNode.h

@@ -1,56 +1,99 @@
-#ifndef BP_BP_NODE_H
+#ifndef BP_BPNODE_H
-#define BP_BP_NODE_H
+#define BP_BPNODE_H
 
 #include <vector>
 #include <map>
-#include <deque>
 #include <string>
+#include <unordered_map>
 
-#include "BayesianNode.h"
+#include "BayesNode.h"
+#include "Shared.h"
 
 using namespace std;
 
-class BpNode : public BayesianNode 
+class BpNode
 {
   public:
-    // constructs
+    BpNode (int);
-    BpNode (string, vector<BayesianNode*>, Distribution* dist, int = -1);
+    BpNode (BayesNode*);
-    // destruct
+
-    ~BpNode (void);
+    ParamSet               getBeliefs (void) const;
-    // methods
-    static void                enableParallelSchedule (void);
-    void                       allocateMemory (void);
-    double*                    getPiValues (void) const;
     double                 getPiValue (int) const;
     void                   setPiValue (int, double);
-    double*                    getLambdaValues (void) const;
     double                 getLambdaValue (int) const;
     void                   setLambdaValue (int, double);
-    double*                    getPiMessages (BpNode*) const;
+    ParamSet&              getPiValues (void);
-    double                     getPiMessage (BpNode*, int) const;
+    ParamSet&              getLambdaValues (void);
-    void                       setPiMessage (BpNode*, int, double);
+    double                 getPiMessageValue (const BayesNode*, int) const;
-    double*                    getLambdaMessages (BpNode*) const;
+    double                 getLambdaMessageValue (const BayesNode*, int) const;
-    double                     getLambdaMessage (BpNode*, int) const;
+    const ParamSet&        getPiMessage (const BayesNode*) const;
-    void                       setLambdaMessage (BpNode*, int, double);
+    const ParamSet&        getLambdaMessage (const BayesNode*) const;
-    double*                    getBeliefs (void);
+    ParamSet&              piNextMessageReference (const BayesNode*);
+    ParamSet&              lambdaNextMessageReference (const BayesNode*);
+    void                   updatePiMessage (const BayesNode*);
+    void                   updateLambdaMessage (const BayesNode*);
     double                 getBeliefChange (void);
-    void                       normalizeMessages (void);
+    void                   updatePiResidual (const BayesNode*);
-    void                       swapMessages (void);
+    void                   updateLambdaResidual (const BayesNode*);
+    void                   clearPiResidual (const BayesNode*);
+    void                   clearLambdaResidual (const BayesNode*);
+    bool                   hasReceivedChildInfluence (void) const;
+    // inlines
+    double                 getPiResidual (const BayesNode*);
+    double                 getLambdaResidual (const BayesNode*);
+    int                    getIndex (const BayesNode*) const;
  
   private:
-    BpNode (const BpNode&);          // disallow copy
+    DISALLOW_COPY_AND_ASSIGN (BpNode);
-    void operator= (const BpNode&);  // disallow assign
+
-    // members
+    IndexMap               indexMap_;
-    double*                    lambdaValues_;
+    ParamSet               piVals_;     // pi values
-    double*                    piValues_;
+    ParamSet               ldVals_;     // lambda values
-    map<BpNode*, double*>      piMessages_;
+    vector<ParamSet>       currPiMsgs_; // current pi messages
-    map<BpNode*, double*>      lambdaMessages_;
+    vector<ParamSet>       currLdMsgs_; // current lambda messages
-    map<BpNode*, double*>*     newPiMessages_;
+    vector<ParamSet>       nextPiMsgs_;
-    map<BpNode*, double*>*     newLambdaMessages_;
+    vector<ParamSet>       nextLdMsgs_;
-    double*                    oldBeliefs_;
+    ParamSet               oldBeliefs_;
-    static bool                parallelSchedule_;
+    ParamSet               piResiduals_;
-    static const double        MAX_CHANGE_ = 1.0;
+    ParamSet               ldResiduals_;
+    int                    ds_;
+    const NodeSet*         childs_;
+    static bool            calculateMessageResidual_;
+//    static const double    MAX_CHANGE_ = 10000000.0;
 };
 
-#endif // BP_BP_NODE_H
+
+
+inline double
+BpNode::getPiResidual (const BayesNode* destination)
+{
+  return piResiduals_[getIndex(destination)];
+}
+
+
+inline double
+BpNode::getLambdaResidual (const BayesNode* source)
+{
+  return ldResiduals_[getIndex(source)];
+}
+
+
+
+inline int
+BpNode::getIndex (const BayesNode* node) const
+{
+  assert (node);
+  //assert (indexMap_.find(node->getVarId()) != indexMap_.end());
+  //return indexMap_.find(node->getVarId())->second;
+  for (unsigned i = 0; childs_->size(); i++) {
+    if ((*childs_)[i]->getVarId() == node->getVarId()) {
+      return i;
+    }
+  }
+  assert (false);
+  return -1;
+}
+
+
+#endif
 

packages/CLPBN/clpbn/bp/CptEntry.h

@@ -1,23 +1,71 @@
-#ifndef CPT_ENTRY_H
+#ifndef BP_CPTENTRY_H
-#define CPT_ENTRY_H
+#define BP_CPTENTRY_H
 
 #include <vector>
 
+#include "Shared.h"
+
 using namespace std;
 
 class CptEntry
 {
   public:
-    // constructs
+    CptEntry (unsigned, const vector<unsigned>&);
-    CptEntry (int, vector<int>);
+
-    // methods
+    unsigned                getParameterIndex (void) const;
-    int                getCptIndex (void) const;
+    const vector<unsigned>& getParentConfigurations (void) const;
-    vector<int>        getDomainInstantiations (void) const;
+    bool                    matchConstraints (const DomainConstr&) const;
-    bool               matchConstraints (const vector<pair<int,int> >&) const;
+    bool                    matchConstraints (const vector<DomainConstr>&) const;
+
   private:
-    // members
+    unsigned                index_;
-    int                cptIndex_;
+    vector<unsigned>        confs_;
-    vector<int>        instantiations_;
 };
 
-#endif // CPT_ENTRY_H
+
+
+inline
+CptEntry::CptEntry (unsigned index, const vector<unsigned>& confs)
+{
+  index_ = index;
+  confs_ = confs;
+}
+
+
+
+inline unsigned
+CptEntry::getParameterIndex (void) const
+{
+  return index_;
+}
+
+
+
+inline const vector<unsigned>&
+CptEntry::getParentConfigurations (void) const
+{
+  return confs_;
+}
+
+
+
+inline bool
+CptEntry::matchConstraints (const DomainConstr& constr) const
+{
+  return confs_[constr.first] == constr.second;
+}
+
+
+
+inline bool
+CptEntry::matchConstraints (const vector<DomainConstr>& constrs) const
+{
+  for (unsigned j = 0; j < constrs.size(); j++) {
+    if (confs_[constrs[j].first] != constrs[j].second) {
+      return false;
+    }
+  }
+  return true;
+}
+
+#endif

packages/CLPBN/clpbn/bp/Distribution.h

@@ -1,24 +1,40 @@
-#ifndef DISTRIBUTION_H
+#ifndef BP_DISTRIBUTION_H
-#define DISTRIBUTION_H
+#define BP_DISTRIBUTION_H
 
 #include <vector>
 #include <string>
 
+#include "Shared.h"
+
 using namespace std;
 
-class CptEntry;
+struct Distribution
-
-class Distribution
 {
   public:
-    Distribution (int, double*, int, vector<string>);
+    Distribution (unsigned id)
-    Distribution (double*, int, vector<string>);
+    {
-    int                      id;
+      this->id     = id;
-    double*                  params;
+      this->params = params;
-    int                      nParams;
+    }
-    vector<string>           domain;
+
-    int*                     offsets;
+    Distribution (const ParamSet& params)
+    {
+      this->id     = -1;
+      this->params = params;
+    }
+
+    void updateParameters (const ParamSet& params)
+    {
+      this->params = params;
+    }
+
+    unsigned          id;
+    ParamSet          params;
+    vector<CptEntry>  entries;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (Distribution);
 };
 
-#endif // DISTRIBUTION
+#endif
 

packages/CLPBN/clpbn/bp/Factor.cpp

@@ -0,0 +1,346 @@
+#include <iostream>
+#include <sstream>
+#include <cstdlib>
+#include <cassert>
+
+#include "Factor.h"
+#include "FgVarNode.h"
+
+
+int Factor::indexCount_ = 0;
+
+Factor::Factor (FgVarNode* var) {
+  vs_.push_back (var);
+  int nParams = var->getDomainSize();
+  // create a uniform distribution
+  double val = 1.0 / nParams;
+  ps_ = ParamSet (nParams, val);
+  id_ = indexCount_;
+  indexCount_ ++;
+}
+
+
+
+Factor::Factor (const FgVarSet& vars) {
+  vs_ = vars;
+  int nParams = 1;
+  for (unsigned i = 0; i < vs_.size(); i++) {
+    nParams *= vs_[i]->getDomainSize();
+  }
+  // create a uniform distribution
+  double val = 1.0 / nParams;
+  ps_ = ParamSet (nParams, val);
+  id_ = indexCount_;
+  indexCount_ ++;
+}
+
+
+
+Factor::Factor (FgVarNode* var,
+                const ParamSet& params)
+{
+  vs_.push_back (var);
+  ps_ = params;
+  id_ = indexCount_;
+  indexCount_ ++;
+}
+
+
+
+Factor::Factor (const FgVarSet& vars,
+                const ParamSet& params)
+{
+  vs_ = vars;
+  ps_ = params;
+  id_ = indexCount_;
+  indexCount_ ++;
+}
+
+
+
+const FgVarSet&
+Factor::getFgVarNodes (void) const
+{
+  return vs_;
+}
+
+
+
+FgVarSet&
+Factor::getFgVarNodes (void)
+{
+  return vs_;
+}
+
+
+
+const ParamSet&
+Factor::getParameters (void) const
+{
+  return ps_;
+}
+
+
+
+ParamSet&
+Factor::getParameters (void)
+{
+  return ps_;
+}
+
+
+
+void
+Factor::setParameters (const ParamSet& params)
+{
+  //cout << "ps size:     " << ps_.size() << endl;
+  //cout << "params size: " << params.size() << endl;
+  assert (ps_.size() == params.size());
+  ps_ = params;
+}
+
+
+
+Factor&
+Factor::operator= (const Factor& g)
+{
+  FgVarSet vars = g.getFgVarNodes();
+  ParamSet params = g.getParameters();
+  return *this;
+}
+
+
+
+Factor&
+Factor::operator*= (const Factor& g)
+{
+  FgVarSet gVs        = g.getFgVarNodes();
+  const ParamSet& gPs = g.getParameters();
+
+  bool hasCommonVars = false;
+  vector<int> varIndexes;
+  for (unsigned i = 0; i < gVs.size(); i++) {
+    int idx = getIndexOf (gVs[i]);
+    if (idx == -1) {
+      insertVariable (gVs[i]);
+      varIndexes.push_back (vs_.size() - 1);
+    } else {
+      hasCommonVars = true;
+      varIndexes.push_back (idx);
+    }
+  }
+ 
+  if (hasCommonVars) {
+   vector<int> offsets (gVs.size());
+    offsets[gVs.size() - 1] = 1;
+    for (int i = gVs.size() - 2; i >= 0; i--) {
+      offsets[i] = offsets[i + 1] * gVs[i + 1]->getDomainSize();
+    }
+    vector<CptEntry> entries = getCptEntries();
+    for (unsigned i = 0; i < entries.size(); i++) {
+      int idx = 0;
+      const DomainConf conf = entries[i].getParentConfigurations();
+      for (unsigned j = 0; j < varIndexes.size(); j++) {
+        idx += offsets[j] * conf[varIndexes[j]];
+      }
+      //cout << "ps_[" << i << "] = " << ps_[i] << " * " ;
+      //cout << gPs[idx] << " , idx = " << idx << endl;
+      ps_[i] = ps_[i] * gPs[idx];
+    }
+  } else {
+    // if the originally factors doesn't have common factors.
+    // we don't have to make domain comparations
+    unsigned idx = 0;
+    for (unsigned i = 0; i < ps_.size(); i++) {
+      //cout << "ps_[" << i << "] = " << ps_[i] << " * " ;
+      //cout << gPs[idx] << " , idx = " << idx << endl;
+      ps_[i] = ps_[i] * gPs[idx];
+      idx ++;
+      if (idx >= gPs.size()) {
+        idx = 0;
+      }
+    }
+  }
+  return *this;
+}
+
+
+
+void
+Factor::insertVariable (FgVarNode* var)
+{
+  int c = 0;
+  ParamSet newPs (ps_.size() * var->getDomainSize());
+  for (unsigned i = 0; i < ps_.size(); i++) {
+    for (int j = 0; j < var->getDomainSize(); j++) {
+      newPs[c] = ps_[i];
+      c ++;
+    }
+  }
+  vs_.push_back (var);
+  ps_ = newPs;
+}
+
+
+
+void
+Factor::marginalizeVariable (const FgVarNode* var) {
+  int varIndex = getIndexOf (var);
+  marginalizeVariable (varIndex);
+}
+
+
+
+void
+Factor::marginalizeVariable (unsigned varIndex)
+{
+  assert (varIndex >= 0 && varIndex < vs_.size());
+  int distOffset    = 1;
+  int leftVarOffset = 1;
+  for (unsigned i = vs_.size() - 1; i > varIndex; i--) {
+    distOffset    *= vs_[i]->getDomainSize();
+    leftVarOffset *= vs_[i]->getDomainSize();
+  }
+  leftVarOffset *= vs_[varIndex]->getDomainSize();
+
+  int ds             = vs_[varIndex]->getDomainSize();
+  int count          = 0;
+  int offset         = 0;
+  int startIndex     = 0;
+  int currDomainIdx  = 0;
+  unsigned newPsSize = ps_.size() / ds;
+  ParamSet newPs;
+  newPs.reserve (newPsSize);
+
+  stringstream ss;
+  ss << "marginalizing " << vs_[varIndex]->getLabel();
+  ss << " from factor " << getLabel() << endl;
+  while (newPs.size() < newPsSize) {
+    ss << "  sum = ";
+    double sum = 0.0;
+    for (int j = 0; j < ds; j++) {
+      if (j != 0) ss << " + ";
+      ss << ps_[offset];
+      sum    = sum    + ps_[offset];
+      offset = offset + distOffset;
+    }
+    newPs.push_back (sum);
+    count ++;
+    if (varIndex == vs_.size() - 1) {
+      offset = count * ds;
+    } else {
+      offset = offset - distOffset + 1;
+      if ((offset % leftVarOffset) == 0) {
+        currDomainIdx ++;
+        startIndex = leftVarOffset * currDomainIdx;
+        offset = startIndex;
+        count = 0;
+      } else {
+        offset = startIndex + count;
+      }
+    }
+    ss << " = " << sum << endl;
+  }
+  //cout << ss.str() << endl;
+  ps_ = newPs;
+  vs_.erase (vs_.begin() + varIndex);
+}
+
+
+
+string
+Factor::getLabel (void) const
+{
+  stringstream ss;
+  ss << "f(" ;
+  // ss << "Φ(" ;
+  for (unsigned i = 0; i < vs_.size(); i++) {
+    if (i != 0) ss << ", " ;
+    ss << "v" << vs_[i]->getVarId();
+  }
+  ss << ")" ;
+  return ss.str();
+}
+
+
+
+string
+Factor::toString (void) const
+{
+  stringstream ss;
+  ss << "vars: " ;
+  for (unsigned i = 0; i < vs_.size(); i++) {
+    if (i != 0) ss << ", " ;
+    ss << "v" << vs_[i]->getVarId();
+  }
+  ss << endl;
+  vector<CptEntry> entries = getCptEntries();
+  for (unsigned i = 0; i < entries.size(); i++) {
+    ss << "Φ(" ;
+    char s = 'a' ;
+    const DomainConf& conf = entries[i].getParentConfigurations();
+    for (unsigned j = 0; j < conf.size(); j++) {
+      if (j != 0) ss << "," ;
+      ss << s << conf[j] + 1;
+      s++;
+    }
+    ss << ") = " << ps_[entries[i].getParameterIndex()] << endl;
+  }
+  return ss.str();
+}
+
+
+
+vector<CptEntry>
+Factor::getCptEntries (void) const
+{
+  vector<DomainConf> confs (ps_.size());
+  for (unsigned i = 0; i < ps_.size(); i++) {
+    confs[i].resize (vs_.size());
+  }
+
+  int nReps = 1;
+  for (int i = vs_.size() - 1; i >= 0; i--) {
+    unsigned index = 0;
+    while (index < ps_.size()) {
+      for (int j = 0; j < vs_[i]->getDomainSize(); j++) {
+        for (int r = 0; r < nReps; r++) {
+          confs[index][i] = j; 
+          index++;
+        }
+      }
+    }
+    nReps *= vs_[i]->getDomainSize();
+  }
+  
+  vector<CptEntry> entries;
+  for (unsigned i = 0; i < ps_.size(); i++) {
+    for (unsigned j = 0; j < vs_.size(); j++) {
+    }
+    entries.push_back (CptEntry (i, confs[i]));
+  }
+  return entries;
+}
+
+
+
+int
+Factor::getIndexOf (const FgVarNode* var) const
+{
+  for (unsigned i = 0; i < vs_.size(); i++) {
+    if (vs_[i] == var) {
+      return i;
+    }
+  }
+  return -1;
+}
+
+
+
+Factor operator* (const Factor& f, const Factor& g)
+{
+  Factor r = f;
+  r *= g;
+  return r;
+}
+

packages/CLPBN/clpbn/bp/Factor.h

@@ -0,0 +1,45 @@
+#ifndef BP_FACTOR_H
+#define BP_FACTOR_H
+
+#include <vector>
+
+#include "CptEntry.h"
+
+using namespace std;
+
+class FgVarNode;
+
+class Factor
+{
+  public:
+    Factor (FgVarNode*);
+    Factor (const FgVarSet&);
+    Factor (FgVarNode*, const ParamSet&);
+    Factor (const FgVarSet&, const ParamSet&);
+
+    const FgVarSet&            getFgVarNodes (void) const;
+    FgVarSet&                  getFgVarNodes (void);
+    const ParamSet&            getParameters (void) const;
+    ParamSet&                  getParameters (void);
+    void                       setParameters (const ParamSet&);
+    Factor&                    operator= (const Factor& f);
+    Factor&                    operator*= (const Factor& f);
+    void                       insertVariable (FgVarNode* index);
+    void                       marginalizeVariable (const FgVarNode* var);
+    void                       marginalizeVariable (unsigned);
+    string                     getLabel (void) const;
+    string                     toString (void) const;
+
+  private:
+    vector<CptEntry>           getCptEntries() const;
+    int                        getIndexOf (const FgVarNode*) const;
+
+    FgVarSet                   vs_;
+    ParamSet                   ps_;
+    int                        id_;
+    static int                 indexCount_;
+};
+
+Factor operator* (const Factor&, const Factor&);
+
+#endif

packages/CLPBN/clpbn/bp/FactorGraph.cpp

@@ -0,0 +1,173 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <cstdlib>
+
+#include "FactorGraph.h"
+#include "FgVarNode.h"
+#include "Factor.h"
+
+
+FactorGraph::FactorGraph (const char* fileName)
+{
+  string line;
+  ifstream is (fileName);
+  if (!is.is_open()) {
+    cerr << "error: cannot read from file " + std::string (fileName) << endl;
+    abort();
+  }
+
+  while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+  getline (is, line);
+  if (line != "MARKOV") {
+    cerr << "error: the network must be a MARKOV network " << endl;
+    abort();
+  }
+
+  while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+  int nVars;
+  is >> nVars;
+
+  while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+  vector<int> domainSizes (nVars);
+  for (int i = 0; i < nVars; i++) {
+    int ds;
+    is >> ds;
+    domainSizes[i] = ds;
+  }
+
+  while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+  for (int i = 0; i < nVars; i++) {
+    varNodes_.push_back (new FgVarNode (i, domainSizes[i]));
+  }
+
+  int nFactors;
+  is >> nFactors;
+  for (int i = 0; i < nFactors; i++) {
+    while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+    int nFactorVars;
+    is >> nFactorVars;
+    FgVarSet factorVars;
+    for (int j = 0; j < nFactorVars; j++) {
+      int varId;
+      is >> varId;
+      FgVarNode* var = getVariableById (varId);
+      if (var == 0) {
+        cerr << "error: invalid variable identifier (" << varId << ")" << endl;
+        abort();
+      }
+      factorVars.push_back (var);
+    }
+    Factor* f = new Factor (factorVars);
+    factors_.push_back (f);
+    for (unsigned j = 0; j < factorVars.size(); j++) {
+      factorVars[j]->addFactor (f);
+    }
+  }
+
+  for (int i = 0; i < nFactors; i++) {
+    while (is.peek() == '#' || is.peek() == '\n') getline (is, line);
+    int nParams;
+    is >> nParams;
+    ParamSet params (nParams);
+    for (int j = 0; j < nParams; j++) {
+      double param;
+      is >> param;
+      params[j] = param;
+    }
+    factors_[i]->setParameters (params);
+  }
+  is.close();
+
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    varNodes_[i]->setIndex (i);
+  }
+}
+
+
+
+FactorGraph::~FactorGraph (void)
+{
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    delete varNodes_[i];
+  }
+  for (unsigned i = 0; i < factors_.size(); i++) {
+    delete factors_[i];
+  }
+}
+
+
+
+FgVarSet
+FactorGraph::getFgVarNodes (void) const
+{
+  return varNodes_;
+}
+
+
+
+vector<Factor*>
+FactorGraph::getFactors (void) const
+{
+  return factors_;
+}
+
+
+
+VarSet
+FactorGraph::getVariables (void) const
+{
+  VarSet vars;
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    vars.push_back (varNodes_[i]);
+  }
+  return vars;
+}
+
+
+
+FgVarNode*
+FactorGraph::getVariableById (unsigned id) const
+{
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    if (varNodes_[i]->getVarId() == id) {
+      return varNodes_[i];
+    }
+  }
+  return 0;
+}
+
+
+
+FgVarNode*
+FactorGraph::getVariableByLabel (string label) const
+{
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    stringstream ss;
+    ss << "v" << varNodes_[i]->getVarId();
+    if (ss.str() == label) {
+      return varNodes_[i];
+    }
+  }
+  return 0;
+}
+
+
+
+void
+FactorGraph::printFactorGraph (void) const
+{
+  for (unsigned i = 0; i < varNodes_.size(); i++) {
+    cout << "variable number " << varNodes_[i]->getIndex() << endl;
+    cout << "Id          = "   << varNodes_[i]->getVarId() << endl;
+    cout << "Domain size = "   << varNodes_[i]->getDomainSize() << endl;
+    cout << "Evidence    = "   << varNodes_[i]->getEvidence() << endl;
+    cout << endl;
+  }
+  cout << endl;
+  for (unsigned i = 0; i < factors_.size(); i++) {
+    cout << factors_[i]->toString() << endl;
+  }
+}
+

packages/CLPBN/clpbn/bp/FactorGraph.h

@@ -0,0 +1,35 @@
+#ifndef BP_FACTORGRAPH_H
+#define BP_FACTORGRAPH_H
+
+#include <vector>
+#include <string>
+
+#include "GraphicalModel.h"
+#include "Shared.h"
+
+using namespace std;
+
+class FgVarNode;
+class Factor;
+
+class FactorGraph : public GraphicalModel
+{
+  public:
+    FactorGraph (const char* fileName);
+   ~FactorGraph (void);
+
+    FgVarSet            getFgVarNodes (void) const;
+    vector<Factor*>     getFactors (void) const;
+    VarSet              getVariables (void) const;
+    FgVarNode*          getVariableById (unsigned) const;
+    FgVarNode*          getVariableByLabel (string) const;
+    void                printFactorGraph (void) const;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (FactorGraph);
+
+    FgVarSet varNodes_;
+    vector<Factor*>     factors_;
+};
+
+#endif

packages/CLPBN/clpbn/bp/FgVarNode.h

@@ -0,0 +1,28 @@
+#ifndef BP_VARIABLE_H
+#define BP_VARIABLE_H
+
+#include <vector>
+#include <string>
+
+#include "Variable.h"
+#include "Shared.h"
+
+using namespace std;
+
+class Factor;
+
+class FgVarNode : public Variable
+{
+  public:
+    FgVarNode (int varId, int dsize) : Variable (varId, dsize) { }
+
+    void                addFactor (Factor* f)      { factors_.push_back (f); }
+    vector<Factor*>     getFactors (void) const    { return factors_; }
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (FgVarNode);
+    // members
+    vector<Factor*>     factors_;
+};
+
+#endif // BP_VARIABLE_H

packages/CLPBN/clpbn/bp/GraphicalModel.h

@@ -0,0 +1,17 @@
+#ifndef BP_GRAPHICALMODEL_H
+#define BP_GRAPHICALMODEL_H
+
+#include "Variable.h"
+#include "Shared.h"
+
+using namespace std;
+
+class GraphicalModel
+{
+  public:
+    virtual VarSet getVariables (void) const = 0;
+
+  private:
+};
+
+#endif

packages/CLPBN/clpbn/bp/HorusCli.cpp

@@ -0,0 +1,214 @@
+#include <iostream>
+#include <cstdlib>
+#include <sstream>
+
+#include "BayesNet.h"
+#include "BPSolver.h"
+
+#include "FactorGraph.h"
+#include "SPSolver.h"
+
+using namespace std;
+
+void BayesianNetwork (int, const char* []);
+void markovNetwork (int, const char* []);
+
+const string USAGE = "usage: \
+./hcli  FILE  [VARIABLE | OBSERVED_VARIABLE=EVIDENCE]..." ;
+
+
+int
+main (int argc, const char* argv[])
+{
+  if (!argv[1]) {
+    cerr << "error: no graphical model specified" << endl;
+    cerr << USAGE << endl;
+    exit (0);
+  }
+  string fileName  = argv[1];
+  string extension = fileName.substr (fileName.find_last_of ('.') + 1);
+  if (extension == "xml") {
+    BayesianNetwork (argc, argv);
+  } else if (extension == "uai") {
+    markovNetwork (argc, argv);
+  } else {
+    cerr << "error: the graphical model must be defined either " ; 
+    cerr << "in a xml file or uai file" << endl;
+    exit (0);
+  }
+  return 0;
+}
+
+
+
+void
+BayesianNetwork (int argc, const char* argv[])
+{
+  BayesNet bn (argv[1]);
+  //bn.printNetwork();
+
+  NodeSet queryVars;
+  for (int i = 2; i < argc; i++) {
+    string arg = argv[i];
+    if (arg.find ('=') == std::string::npos) {
+      BayesNode* queryVar = bn.getNode (arg);
+      if (queryVar) {
+        queryVars.push_back (queryVar);
+      } else {
+        cerr << "error: there isn't a variable labeled of " ;
+        cerr << "`" << arg << "'" ;
+        cerr << endl;
+        exit (0);
+      }
+    } else {
+      size_t pos   = arg.find ('=');
+      string label = arg.substr (0, pos);
+      string state = arg.substr (pos + 1);
+      if (label.empty()) {
+        cerr << "error: missing left argument" << endl;
+        cerr << USAGE << endl;
+        exit (0);
+      }
+      if (state.empty()) {
+        cerr << "error: missing right argument" << endl;
+        cerr << USAGE << endl;
+        exit (0);
+      }
+      BayesNode* node = bn.getNode (label);
+      if (node) {
+        if (node->isValidState (state)) {
+          node->setEvidence (state);
+        } else {
+          cerr << "error: `" << state << "' " ;
+          cerr << "is not a valid state for " ;
+          cerr << "`" << node->getLabel() << "'" ;
+          cerr << endl;
+          exit (0);
+        }
+      } else {
+        cerr << "error: there isn't a variable labeled of " ;
+        cerr << "`" << label << "'" ;
+        cerr << endl;
+        exit (0);
+      }
+    }
+  }
+
+  BPSolver solver (bn);
+  if (queryVars.size() == 0) {
+    solver.runSolver();
+    solver.printAllPosterioris();
+  } else if (queryVars.size() == 1) {
+    solver.runSolver();
+    solver.printPosterioriOf (queryVars[0]);
+  } else {
+    Domain domain   = BayesNet::getInstantiations(queryVars);
+    ParamSet params = solver.getJointDistribution (queryVars);
+    for (unsigned i = 0; i < params.size(); i++) {
+      cout << domain[i] << "\t" << params[i] << endl;
+    }
+  }
+  bn.freeDistributions();
+}
+
+
+
+void
+markovNetwork (int argc, const char* argv[])
+{
+  FactorGraph fg (argv[1]);
+  //fg.printFactorGraph();
+
+  VarSet queryVars;
+  for (int i = 2; i < argc; i++) {
+    string arg = argv[i];
+    if (arg.find ('=') == std::string::npos) {
+      if (!Util::isInteger (arg)) {
+        cerr << "error: `" << arg << "' " ;
+        cerr << "is not a valid variable id" ;
+        cerr << endl;
+        exit (0);
+      }
+      unsigned varId;
+      stringstream ss;
+      ss << arg;
+      ss >> varId;
+      Variable* queryVar = fg.getVariableById (varId);
+      if (queryVar) {
+        queryVars.push_back (queryVar);
+      } else {
+        cerr << "error: there isn't a variable with " ;
+        cerr << "`" << varId << "' as id" ;
+        cerr << endl;
+        exit (0);
+      }
+    } else {
+      size_t pos = arg.find ('=');
+      if (arg.substr (0, pos).empty()) {
+        cerr << "error: missing left argument" << endl;
+        cerr << USAGE << endl;
+        exit (0);
+      }
+      if (arg.substr (pos + 1).empty()) {
+        cerr << "error: missing right argument" << endl;
+        cerr << USAGE << endl;
+        exit (0);
+      }
+      if (!Util::isInteger (arg.substr (0, pos))) {
+        cerr << "error: `" << arg.substr (0, pos) << "' " ;
+        cerr << "is not a variable id" ;
+        cerr << endl;
+        exit (0);
+      }
+      unsigned varId;
+      stringstream ss; 
+      ss << arg.substr (0, pos);
+      ss >> varId;
+      Variable* var = fg.getVariableById (varId);
+      if (var) {
+        if (!Util::isInteger (arg.substr (pos + 1))) {
+          cerr << "error: `" << arg.substr (pos + 1) << "' " ;
+          cerr << "is not a state index" ;
+          cerr << endl;
+          exit (0);
+        }
+        int stateIndex;
+        stringstream ss; 
+        ss << arg.substr (pos + 1);
+        ss >> stateIndex;
+        cout << "si: " << stateIndex << endl;
+        if (var->isValidStateIndex (stateIndex)) {
+          var->setEvidence (stateIndex);
+        } else {
+          cerr << "error: `" << stateIndex << "' " ;
+          cerr << "is not a valid state index for variable " ;
+          cerr << "`" << var->getVarId() << "'" ;
+          cerr << endl;
+          exit (0);
+        }
+      } else {
+        cerr << "error: there isn't a variable with " ;
+        cerr << "`" << varId << "' as id" ;
+        cerr << endl;
+        exit (0);
+      }
+    }
+  }
+
+  SPSolver solver (fg);
+  if (queryVars.size() == 0) {
+    solver.runSolver();
+    solver.printAllPosterioris();
+  } else if (queryVars.size() == 1) {
+    solver.runSolver();
+    solver.printPosterioriOf (queryVars[0]);
+  } else {
+    assert (false); //FIXME
+    //Domain domain   = BayesNet::getInstantiations(queryVars);
+    //ParamSet params = solver.getJointDistribution (queryVars);
+    //for (unsigned i = 0; i < params.size(); i++) {
+    //  cout << domain[i] << "\t" << params[i] << endl;
+    //}
+  }
+}
+

packages/CLPBN/clpbn/bp/HorusYap.cpp

@@ -0,0 +1,229 @@
+#include <cstdlib>
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <string>
+
+#include <YapInterface.h>
+
+#include "callgrind.h"
+
+#include "BayesNet.h"
+#include "BayesNode.h"
+#include "BPSolver.h"
+
+using namespace std;
+
+int
+createNetwork (void)
+{
+  Statistics::numCreatedNets ++;
+  cout << "creating network number " << Statistics::numCreatedNets << endl;
+  if (Statistics::numCreatedNets == 1) {
+    //CALLGRIND_START_INSTRUMENTATION;
+  }
+  BayesNet* bn = new BayesNet();
+
+  YAP_Term varList = YAP_ARG1;
+  while (varList != YAP_TermNil()) {
+    YAP_Term var     = YAP_HeadOfTerm (varList);
+    unsigned varId   = (unsigned) YAP_IntOfTerm (YAP_ArgOfTerm (1, var));
+    unsigned dsize   = (unsigned) YAP_IntOfTerm (YAP_ArgOfTerm (2, var));
+    int evidence     = (int)      YAP_IntOfTerm (YAP_ArgOfTerm (3, var));
+    YAP_Term parentL =                           YAP_ArgOfTerm (4, var);
+    unsigned distId  = (unsigned) YAP_IntOfTerm (YAP_ArgOfTerm (5, var));
+    NodeSet parents;
+    while (parentL != YAP_TermNil()) {
+      unsigned parentId = (unsigned) YAP_IntOfTerm (YAP_HeadOfTerm (parentL));
+      BayesNode* parent = bn->getNode (parentId);
+      if (!parent) {
+        parent = bn->addNode (parentId);
+      }
+      parents.push_back (parent);
+      parentL = YAP_TailOfTerm (parentL);
+    }
+    Distribution* dist = bn->getDistribution (distId);
+    if (!dist) {
+      dist = new Distribution (distId);
+      bn->addDistribution (dist);
+    }
+    BayesNode* node = bn->getNode (varId);
+    if (node) {
+      node->setData (dsize, evidence, parents, dist);
+    } else {
+      bn->addNode (varId, dsize, evidence, parents, dist);
+    }
+    varList = YAP_TailOfTerm (varList);
+  }
+  bn->setIndexes();
+
+  if (Statistics::numCreatedNets == 1688) {
+    Statistics::writeStats();
+    //Statistics::writeStats();
+    //CALLGRIND_STOP_INSTRUMENTATION;
+    //CALLGRIND_DUMP_STATS;
+    //exit (0);
+  }
+  YAP_Int p = (YAP_Int) (bn);
+  return YAP_Unify (YAP_MkIntTerm (p), YAP_ARG2);
+}
+
+
+
+int
+setExtraVarsInfo (void)
+{
+  BayesNet* bn     = (BayesNet*) YAP_IntOfTerm (YAP_ARG1);
+  YAP_Term varsInfoL =                            YAP_ARG2;
+  while (varsInfoL != YAP_TermNil()) {
+    YAP_Term head    = YAP_HeadOfTerm (varsInfoL);
+    unsigned varId   = YAP_IntOfTerm  (YAP_ArgOfTerm (1, head));
+    YAP_Atom label   = YAP_AtomOfTerm (YAP_ArgOfTerm (2, head));
+    YAP_Term domainL =                 YAP_ArgOfTerm (3, head);
+    Domain domain;
+    while (domainL != YAP_TermNil()) {
+      YAP_Atom atom = YAP_AtomOfTerm (YAP_HeadOfTerm (domainL));
+      domain.push_back ((char*) YAP_AtomName (atom));
+      domainL = YAP_TailOfTerm (domainL);
+    }
+    BayesNode* node = bn->getNode (varId);
+    assert (node);
+    node->setLabel ((char*) YAP_AtomName (label));
+    node->setDomain (domain);
+    varsInfoL = YAP_TailOfTerm (varsInfoL);
+  }
+  return TRUE;
+}
+
+
+
+int
+setParameters (void)
+{
+  BayesNet* bn    = (BayesNet*) YAP_IntOfTerm (YAP_ARG1);
+  YAP_Term distList =                            YAP_ARG2;
+  while (distList != YAP_TermNil()) {
+    YAP_Term dist    = YAP_HeadOfTerm (distList);
+    unsigned distId  = (unsigned) YAP_IntOfTerm (YAP_ArgOfTerm (1, dist));
+    YAP_Term paramL  =                           YAP_ArgOfTerm (2, dist);
+    ParamSet params;
+    while (paramL!= YAP_TermNil()) {
+      params.push_back ((double) YAP_FloatOfTerm (YAP_HeadOfTerm (paramL)));
+      paramL = YAP_TailOfTerm (paramL);
+    }
+    bn->getDistribution(distId)->updateParameters(params);
+    distList = YAP_TailOfTerm (distList);
+  }
+  return TRUE;
+}
+
+
+
+int
+runSolver (void)
+{
+  BayesNet* bn      = (BayesNet*) YAP_IntOfTerm (YAP_ARG1);
+  YAP_Term taskList =                            YAP_ARG2;
+
+  vector<NodeSet> tasks;
+  NodeSet marginalVars;
+
+  while (taskList != YAP_TermNil()) {
+    if (YAP_IsPairTerm (YAP_HeadOfTerm (taskList))) {
+      NodeSet jointVars;
+      YAP_Term jointList = YAP_HeadOfTerm (taskList);      
+      while (jointList != YAP_TermNil()) {
+        unsigned varId = (unsigned) YAP_IntOfTerm (YAP_HeadOfTerm (jointList));
+        assert (bn->getNode (varId));
+        jointVars.push_back (bn->getNode (varId));
+        jointList = YAP_TailOfTerm (jointList);
+      }
+      tasks.push_back (jointVars);
+    } else {
+      unsigned varId = (unsigned) YAP_IntOfTerm (YAP_HeadOfTerm (taskList));
+      BayesNode* node = bn->getNode (varId);
+      assert (node);
+      tasks.push_back (NodeSet() = {node});
+      marginalVars.push_back (node);
+    }
+    taskList = YAP_TailOfTerm (taskList);
+  }
+  /*
+  cout << "tasks to resolve:" << endl;
+  for (unsigned i = 0; i < tasks.size(); i++) {
+    cout << "i" << ": " ;
+    if (tasks[i].size() == 1) {
+      cout << tasks[i][0]->getVarId() << endl;
+    } else {
+      for (unsigned j = 0; j < tasks[i].size(); j++) {
+        cout << tasks[i][j]->getVarId() << " " ;
+      }
+      cout << endl;
+    }
+  }
+  */
+  
+  cerr << "prunning now..." << endl;
+  BayesNet* prunedNet = bn->pruneNetwork (marginalVars);
+  bn->printNetworkToFile ("net.txt");
+  BPSolver solver (*prunedNet);
+  cerr << "solving marginals now..." << endl;
+  solver.runSolver();
+  cerr << "calculating joints now ..." << endl;
+
+  vector<ParamSet> results;
+  results.reserve (tasks.size());  
+  for (unsigned i = 0; i < tasks.size(); i++) {
+    if (tasks[i].size() == 1) {
+      BayesNode* node = prunedNet->getNode (tasks[i][0]->getVarId());
+      results.push_back (solver.getPosterioriOf (node));
+    } else {
+      BPSolver solver2 (*bn);
+      cout << "calculating an join dist on: " ;
+      for (unsigned j = 0; j < tasks[i].size(); j++) {
+        cout << tasks[i][j]->getVarId() << " " ;
+      }
+      cout << "..." << endl;
+      results.push_back (solver2.getJointDistribution (tasks[i]));
+    }
+  }
+
+  delete prunedNet;
+
+  YAP_Term list = YAP_TermNil();
+  for (int i = results.size() - 1; i >= 0; i--) {
+    const ParamSet& beliefs = results[i];
+    YAP_Term queryBeliefsL = YAP_TermNil();
+    for (int j = beliefs.size() - 1; j >= 0; j--) {
+      YAP_Term belief = YAP_MkFloatTerm (beliefs[j]);
+      queryBeliefsL = YAP_MkPairTerm (belief, queryBeliefsL);
+    }
+    list = YAP_MkPairTerm (queryBeliefsL, list);
+  }
+ 
+  return YAP_Unify (list, YAP_ARG3);
+}
+
+
+
+int
+deleteBayesNet (void)
+{
+  BayesNet* bn = (BayesNet*) YAP_IntOfTerm (YAP_ARG1);
+  bn->freeDistributions();
+  delete bn;
+  return TRUE;
+}
+
+
+
+extern "C" void
+init_predicates (void)
+{
+  YAP_UserCPredicate ("create_network",      createNetwork,    2);
+  YAP_UserCPredicate ("set_extra_vars_info", setExtraVarsInfo, 2);
+  YAP_UserCPredicate ("set_parameters",      setParameters,    2);
+  YAP_UserCPredicate ("run_solver",          runSolver,        3);
+  YAP_UserCPredicate ("delete_bayes_net",    deleteBayesNet,   1);
+}
+

packages/CLPBN/clpbn/bp/Makefile.in

@@ -21,7 +21,17 @@ YAPLIBDIR=@libdir@/Yap
 #
 CC=@CC@
 CXX=@CXX@
-CXXFLAGS= @SHLIB_CXXFLAGS@ $(YAP_EXTRAS) $(DEFS) -D_YAP_NOT_INSTALLED_=1 -I$(srcdir) -I../../../.. -I$(srcdir)/../../../../include  @CPPFLAGS@
+
+# normal
+CXXFLAGS= -std=c++0x @SHLIB_CXXFLAGS@ $(YAP_EXTRAS) $(DEFS) -D_YAP_NOT_INSTALLED_=1 -I$(srcdir) -I../../../.. -I$(srcdir)/../../../../include @CPPFLAGS@ -DNDEBUG
+
+# debug 
+#CXXFLAGS= -std=c++0x @SHLIB_CXXFLAGS@ $(YAP_EXTRAS) $(DEFS) -D_YAP_NOT_INSTALLED_=1 -I$(srcdir) -I../../../.. -I$(srcdir)/../../../../include @CPPFLAGS@ -g -O0
+
+# profiling (callgrind)
+#CXXFLAGS= -std=c++0x @SHLIB_CXXFLAGS@ $(YAP_EXTRAS) $(DEFS) -D_YAP_NOT_INSTALLED_=1 -I$(srcdir) -I../../../.. -I$(srcdir)/../../../../include @CPPFLAGS@ -g -DNDEBUG
+
+
 #
 #
 # You shouldn't need to change what follows.
@@ -38,64 +48,75 @@ CWD=$(PWD)
 
 
 HEADERS = \
-	$(srcdir)/BayesianNetwork.h \
+	$(srcdir)/GraphicalModel.h \
-	$(srcdir)/BayesianNode.h \
+	$(srcdir)/Variable.h \
-	$(srcdir)/BpNetwork.h \
+	$(srcdir)/BayesNet.h \
-	$(srcdir)/BpNode.h \
+	$(srcdir)/BayesNode.h \
 	$(srcdir)/Distribution.h \
 	$(srcdir)/CptEntry.h \
-	$(srcdir)/BifInterface.h \
+	$(srcdir)/FactorGraph.h \
+	$(srcdir)/FgVarNode.h \
+	$(srcdir)/Factor.h \
+	$(srcdir)/Solver.h \
+	$(srcdir)/BPSolver.h \
+	$(srcdir)/BpNode.h \
+	$(srcdir)/SPSolver.h \
+	$(srcdir)/Shared.h \
 	$(srcdir)/xmlParser/xmlParser.h
 	
 CPP_SOURCES = \
-	$(srcdir)/BayesianNetwork.cpp \
+	$(srcdir)/BayesNet.cpp \
-	$(srcdir)/BayesianNode.cpp \
+	$(srcdir)/BayesNode.cpp \
-	$(srcdir)/BpNetwork.cpp \
+	$(srcdir)/FactorGraph.cpp \
+	$(srcdir)/Factor.cpp \
+	$(srcdir)/BPSolver.cpp \
 	$(srcdir)/BpNode.cpp \
-	$(srcdir)/Distribution.cpp \
+	$(srcdir)/SPSolver.cpp \
-	$(srcdir)/CptEntry.cpp \
+	$(srcdir)/HorusYap.cpp \
-	$(srcdir)/Horus.cpp \
+	$(srcdir)/HorusCli.cpp \
-	$(srcdir)/BifInterface.cpp \
-	$(srcdir)/BifTest.cpp \
 	$(srcdir)/xmlParser/xmlParser.cpp
 
 OBJS = \
-	BayesianNetwork.o \
+	BayesNet.o \
-	BayesianNode.o \
+	BayesNode.o \
-	BpNetwork.o \
+	FactorGraph.o \
+	Factor.o \
+	BPSolver.o \
 	BpNode.o \
-    Distribution.o \
+	SPSolver.o \
-	CptEntry.o \
+	HorusYap.o
-	Horus.o
 
-BIF_OBJS = \
+HCLI_OBJS = \
-	BayesianNetwork.o \
+	BayesNet.o \
-	BayesianNode.o \
+	BayesNode.o \
-    BpNetwork.o \
+	FactorGraph.o \
+	Factor.o \
+	BPSolver.o \
 	BpNode.o \
-    Distribution.o \
+	SPSolver.o \
-	CptEntry.o \
+	HorusCli.o \
-    BifInterface.o \
-    BifTest.o \
 	xmlParser.o
 
 SOBJS=horus.@SO@
 
 
-all: $(SOBJS) biftest
+all: $(SOBJS) hcli
 
 # default rule
+%.o : $(srcdir)/%.cpp
+	$(CXX) -c $(CXXFLAGS) $< -o $@
+
+
 xmlParser.o : $(srcdir)/xmlParser/xmlParser.cpp
 	$(CXX) -c $(CXXFLAGS) $< -o $@
 
-%.o : $(srcdir)/%.cpp
-	$(CXX) -c $(CXXFLAGS) $< -o $@
 
 @DO_SECOND_LD@horus.@SO@: $(OBJS)
 @DO_SECOND_LD@	@SHLIB_CXX_LD@ -o horus.@SO@ $(OBJS) @EXTRA_LIBS_FOR_SWIDLLS@
 
-biftest: $(BIF_OBJS)
+
-	$(CXX) -o biftest $(BIF_OBJS)
+hcli: $(HCLI_OBJS)
+	$(CXX) -o hcli $(HCLI_OBJS)
 
 
 install: all
@@ -103,12 +124,12 @@ install: all
 
 
 clean:
-	rm -f *.o *~ $(OBJS) $(SOBJS) *.BAK biftest xmlParser/*.o
+	rm -f *.o *~ $(OBJS) $(SOBJS) *.BAK hcli xmlParser/*.o
 
 
 depend: $(HEADERS) $(CPP_SOURCES)
 	-@if test "$(GCC)" = yes; then\
-	  $(CC) -MM -MG $(CFLAGS) -I$(srcdir) -I$(srcdir)/../../../../include -I$(srcdir)/../../../../H $(CPP_SOURCES) >> Makefile;\
+	  $(CC) -std=c++0x -MM -MG $(CFLAGS) -I$(srcdir) -I$(srcdir)/../../../../include -I$(srcdir)/../../../../H $(CPP_SOURCES) >> Makefile;\
 	else\
           makedepend -f - -- $(CFLAGS) -I$(srcdir)/../../../../H -I$(srcdir)/../../../../include -- $(CPP_SOURCES) |\
           sed 's|.*/\([^:]*\):|\1:|' >> Makefile ;\

packages/CLPBN/clpbn/bp/SPSolver.cpp

@@ -0,0 +1,295 @@
+#include <cassert>
+#include <algorithm>
+#include <iostream>
+
+#include "SPSolver.h"
+#include "FactorGraph.h"
+#include "FgVarNode.h"
+#include "Factor.h"
+
+SPSolver* Link::klass = 0;
+
+
+SPSolver::SPSolver (const FactorGraph& fg) : Solver (&fg)
+{
+  fg_           = &fg;
+  accuracy_     = 0.0001;
+  maxIter_      = 10000;
+  //schedule_     = S_SEQ_FIXED;
+  //schedule_     = S_SEQ_RANDOM;
+  //schedule_     = S_SEQ_PARALLEL;
+  schedule_     = S_MAX_RESIDUAL;
+  Link::klass   = this;
+  FgVarSet vars = fg_->getFgVarNodes();
+  for (unsigned i = 0; i < vars.size(); i++) {
+    msgs_.push_back (new MessageBanket (vars[i]));
+  }
+}
+
+
+
+SPSolver::~SPSolver (void)
+{
+  for (unsigned i = 0; i < msgs_.size(); i++) {
+    delete msgs_[i];
+  }
+}
+
+
+
+void
+SPSolver::runSolver (void)
+{
+  nIter_ = 0;
+  vector<Factor*> factors = fg_->getFactors();
+  for (unsigned i = 0; i < factors.size(); i++) {
+    FgVarSet neighbors = factors[i]->getFgVarNodes();
+    for (unsigned j = 0; j < neighbors.size(); j++) {
+      updateOrder_.push_back (Link (factors[i], neighbors[j]));
+    }
+  }
+  
+  while (!converged() && nIter_ < maxIter_) {
+    if (DL >= 1) {
+      cout << endl;
+      cout << "****************************************" ;
+      cout << "****************************************" ;
+      cout << endl;
+      cout << " Iteration " << nIter_ + 1 << endl;
+      cout << "****************************************" ;
+      cout << "****************************************" ;
+      cout << endl;
+    }
+
+    switch (schedule_) {
+
+      case S_SEQ_RANDOM:
+        random_shuffle (updateOrder_.begin(), updateOrder_.end());
+        // no break
+
+      case S_SEQ_FIXED:
+        for (unsigned c = 0; c < updateOrder_.size(); c++) {
+          Link& link = updateOrder_[c];
+          calculateNextMessage (link.source, link.destination);
+          updateMessage (updateOrder_[c]);
+        }
+        break;
+
+      case S_PARALLEL:
+        for (unsigned c = 0; c < updateOrder_.size(); c++) {
+          Link link = updateOrder_[c];
+          calculateNextMessage (link.source, link.destination);
+        }
+        for (unsigned c = 0; c < updateOrder_.size(); c++) {
+          Link link = updateOrder_[c];
+          updateMessage (updateOrder_[c]);
+        }
+        break;
+
+      case S_MAX_RESIDUAL:
+        maxResidualSchedule();
+        break;
+    }
+
+    nIter_++;
+  }
+  cout << endl;
+  if (DL >= 1) {
+    if (nIter_ < maxIter_) {
+      cout << "Loopy Sum-Product converged in " ; 
+      cout << nIter_ << " iterations" << endl;
+    } else {
+      cout << "The maximum number of iterations was hit, terminating..." ;
+      cout << endl;
+    }
+  }
+}
+
+
+
+ParamSet
+SPSolver::getPosterioriOf (const Variable* var) const
+{
+  assert (var);
+  assert (var == fg_->getVariableById (var->getVarId()));
+  assert (var->getIndex() < msgs_.size());
+
+  ParamSet probs (var->getDomainSize(), 1);
+  if (var->hasEvidence()) {
+    for (unsigned i = 0; i < probs.size(); i++) {
+      if ((int)i != var->getEvidence()) {
+        probs[i] = 0;
+      }
+    }
+
+  } else {
+
+    MessageBanket* mb = msgs_[var->getIndex()]; 
+    const FgVarNode* varNode  = fg_->getFgVarNodes()[var->getIndex()];
+    vector<Factor*> neighbors = varNode->getFactors();
+    for (unsigned i = 0; i < neighbors.size(); i++) {
+      const Message& msg = mb->getMessage (neighbors[i]);
+      for (unsigned j = 0; j < msg.size(); j++) {
+        probs[j] *= msg[j];
+      }
+    }
+    Util::normalize (probs);
+  }
+
+  return probs;
+}
+
+
+
+bool
+SPSolver::converged (void)
+{
+  if (nIter_ == 0 || nIter_ == 1) {
+    return false;
+  }
+  bool converged = true;
+  for (unsigned i = 0; i < updateOrder_.size(); i++) {
+    double residual = getResidual (updateOrder_[i]);
+    if (DL >= 1) {
+      cout << updateOrder_[i].toString();
+      cout << " residual = " << residual << endl;
+    }
+    if (residual > accuracy_) {
+      converged = false;
+      if (DL == 0) {
+        break;
+      }
+    } 
+  }
+  return converged;
+}
+
+
+
+void
+SPSolver::maxResidualSchedule (void)
+{
+  if (nIter_ == 0) {
+    for (unsigned c = 0; c < updateOrder_.size(); c++) {
+      Link& l = updateOrder_[c];
+      calculateNextMessage (l.source, l.destination);
+      if (DL >= 1) {
+        cout << updateOrder_[c].toString() << " residual = " ;
+        cout << getResidual (updateOrder_[c]) << endl;
+      }
+    }
+    sort (updateOrder_.begin(), updateOrder_.end(), compareResidual);
+  } else {
+
+    for (unsigned c = 0; c < updateOrder_.size(); c++) {
+      Link& link = updateOrder_.front();
+      updateMessage (link); 
+      resetResidual (link);
+
+      // update the messages that depend on message source --> destination
+      vector<Factor*> fstLevelNeighbors = link.destination->getFactors();
+      for (unsigned i = 0; i < fstLevelNeighbors.size(); i++) {
+        if (fstLevelNeighbors[i] != link.source) {
+          FgVarSet sndLevelNeighbors;
+          sndLevelNeighbors = fstLevelNeighbors[i]->getFgVarNodes();
+          for (unsigned j = 0; j < sndLevelNeighbors.size(); j++) {
+            if (sndLevelNeighbors[j] != link.destination) {
+              calculateNextMessage (fstLevelNeighbors[i], sndLevelNeighbors[j]);
+            }
+          }
+        }
+      }  
+      sort (updateOrder_.begin(), updateOrder_.end(), compareResidual);
+    }
+  }
+}
+
+
+
+void
+SPSolver::updateMessage (const Link& link)
+{
+  updateMessage (link.source, link.destination);
+}
+
+
+
+void
+SPSolver::updateMessage (const Factor* src, const FgVarNode* dest)
+{
+  msgs_[dest->getIndex()]->updateMessage (src);
+/*  cout << src->getLabel() << " --> " << dest->getLabel() << endl; 
+  cout << "    m: " ;
+  Message msg = msgs_[dest->getIndex()]->getMessage (src);
+  for (unsigned i = 0; i < msg.size(); i++) {
+    if (i != 0) cout << ", " ;
+    cout << msg[i];
+  }
+  cout << endl;
+*/
+}
+
+
+
+void
+SPSolver::calculateNextMessage (const Link& link)
+{
+  calculateNextMessage (link.source, link.destination);
+}
+
+
+void
+SPSolver::calculateNextMessage (const Factor* src, const FgVarNode* dest)
+{
+  FgVarSet neighbors = src->getFgVarNodes(); 
+  // calculate the product of MessageBankets sended
+  // to factor `src', except from var `dest'
+  Factor result = *src;
+  for (unsigned i = 0; i < neighbors.size(); i++) {
+    if (neighbors[i] != dest) {
+      Message msg (neighbors[i]->getDomainSize(), 1);
+      calculateVarFactorMessage (neighbors[i], src, msg);
+      result *= Factor (neighbors[i], msg);
+    }
+  }
+  // marginalize all vars except `dest'
+  for (unsigned i = 0; i < neighbors.size(); i++) {
+    if (neighbors[i] != dest) {
+      result.marginalizeVariable (neighbors[i]);
+    }
+  }
+  msgs_[dest->getIndex()]->setNextMessage (src, result.getParameters());
+}
+
+
+
+void
+SPSolver::calculateVarFactorMessage (const FgVarNode* src,
+                                     const Factor* dest,
+                                     Message& placeholder) const
+{
+  assert (src->getDomainSize() == (int)placeholder.size());
+  if (src->hasEvidence()) {
+    for (unsigned i = 0; i < placeholder.size(); i++) {
+      if ((int)i != src->getEvidence()) {
+        placeholder[i] = 0.0;
+      } else {
+        placeholder[i] = 1.0;
+      }
+    }
+
+  } else {
+
+    MessageBanket* mb = msgs_[src->getIndex()];
+    vector<Factor*> neighbors = src->getFactors();
+    for (unsigned i = 0; i < neighbors.size(); i++) {
+      if (neighbors[i] != dest) {
+        const Message& fromFactor = mb->getMessage (neighbors[i]);
+        for (unsigned j = 0; j < fromFactor.size(); j++) {
+          placeholder[j] *= fromFactor[j];
+        }
+      }
+    }
+  }
+}
+

packages/CLPBN/clpbn/bp/SPSolver.h

@@ -0,0 +1,171 @@
+#ifndef BP_SPSOLVER_H
+#define BP_SPSOLVER_H
+
+#include <cmath>
+#include <map>
+#include <vector>
+#include <string>
+
+#include "Solver.h"
+#include "FgVarNode.h"
+#include "Factor.h"
+
+using namespace std;
+
+class FactorGraph;
+class SPSolver;
+
+struct Link
+{
+  Link (Factor* s, FgVarNode* d)
+  { 
+    source      = s;
+    destination = d;
+  }
+  string toString (void) const
+  {
+    stringstream ss;
+    ss << source->getLabel() << " --> " ;
+    ss << destination->getLabel();
+    return ss.str();
+  }
+  Factor*           source;
+  FgVarNode*        destination;
+  static SPSolver*  klass;
+};
+
+
+
+class MessageBanket
+{
+  public:
+    MessageBanket (const FgVarNode* var)
+    {
+      vector<Factor*> sources = var->getFactors();
+      for (unsigned i = 0; i < sources.size(); i++) {
+        indexMap_.insert (make_pair (sources[i], i));
+        currMsgs_.push_back (Message(var->getDomainSize(), 1));
+        nextMsgs_.push_back (Message(var->getDomainSize(), -10));
+        residuals_.push_back (0.0);
+      }
+    }
+
+    void updateMessage (const Factor* source)
+    {
+      unsigned idx = getIndex(source);
+      currMsgs_[idx] = nextMsgs_[idx];
+    }
+
+    void setNextMessage (const Factor* source, const Message& msg)
+    {
+      unsigned idx = getIndex(source);
+      nextMsgs_[idx] = msg;
+      residuals_[idx] = computeResidual (source);
+    }
+
+    const Message& getMessage (const Factor* source) const
+    {
+      return currMsgs_[getIndex(source)];
+    }
+
+    double getResidual (const Factor* source) const
+    {
+      return residuals_[getIndex(source)];
+    }
+    
+    void resetResidual (const Factor* source)
+    {
+      residuals_[getIndex(source)] = 0.0;
+    }
+
+  private:
+    double computeResidual (const Factor* source)
+    {
+      double change = 0.0;
+      unsigned idx  = getIndex (source);
+      const Message& currMessage = currMsgs_[idx];
+      const Message& nextMessage = nextMsgs_[idx];
+      for (unsigned i = 0; i < currMessage.size(); i++) {
+        change += abs (currMessage[i] - nextMessage[i]);
+      }
+      return change;
+    }
+ 
+    unsigned getIndex (const Factor* factor) const
+    {
+      assert (factor);
+      assert (indexMap_.find(factor) != indexMap_.end());
+      return indexMap_.find(factor)->second;
+    }
+ 
+    typedef map<const Factor*, unsigned> IndexMap;
+
+    IndexMap               indexMap_;
+    vector<Message>        currMsgs_;
+    vector<Message>        nextMsgs_;
+    vector<double>         residuals_;
+};
+
+
+
+class SPSolver : public Solver
+{
+  public:
+    SPSolver (const FactorGraph&);
+   ~SPSolver (void);
+
+    void            runSolver (void);
+    ParamSet        getPosterioriOf (const Variable* var) const;
+
+  private:
+    bool            converged (void);
+    void            maxResidualSchedule (void);
+    void            updateMessage (const Link&);
+    void            updateMessage (const Factor*, const FgVarNode*);
+    void            calculateNextMessage (const Link&);
+    void            calculateNextMessage (const Factor*, const FgVarNode*);
+    void            calculateVarFactorMessage (
+                         const FgVarNode*, const Factor*, Message&) const;
+    double          getResidual (const Link&) const;
+    void            resetResidual (const Link&) const;
+    friend bool     compareResidual (const Link&, const Link&);
+
+    const FactorGraph*          fg_;
+    vector<MessageBanket*>      msgs_;
+    Schedule                    schedule_;
+    int                         nIter_;
+    double                      accuracy_;
+    int                         maxIter_;
+    vector<Link>                updateOrder_;
+};
+
+
+
+inline double
+SPSolver::getResidual (const Link& link) const
+{
+  MessageBanket* mb = Link::klass->msgs_[link.destination->getIndex()];
+  return mb->getResidual (link.source);
+}
+
+
+
+inline void
+SPSolver::resetResidual (const Link& link) const
+{
+  MessageBanket* mb = Link::klass->msgs_[link.destination->getIndex()];
+  mb->resetResidual (link.source);
+}
+
+
+
+inline bool
+compareResidual (const Link& link1, const Link& link2)
+{
+  MessageBanket* mb1 = Link::klass->msgs_[link1.destination->getIndex()];
+  MessageBanket* mb2 = Link::klass->msgs_[link2.destination->getIndex()];
+  return mb1->getResidual(link1.source) > mb2->getResidual(link2.source);
+}
+
+#endif
+

packages/CLPBN/clpbn/bp/Shared.h

@@ -0,0 +1,203 @@
+#ifndef BP_SHARED_H
+#define BP_SHARED_H
+
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <cassert>
+#include <vector>
+#include <map>
+#include <unordered_map>
+
+// Macro to disallow the copy constructor and operator= functions
+#define DISALLOW_COPY_AND_ASSIGN(TypeName) \
+  TypeName(const TypeName&);               \
+  void operator=(const TypeName&)
+
+using namespace std;
+
+class Variable;
+class BayesNode;
+class FgVarNode;
+
+typedef double                             Param;
+typedef vector<Param>                      ParamSet;
+typedef vector<Param>                      Message;
+typedef vector<Variable*>                  VarSet;
+typedef vector<BayesNode*>                 NodeSet;
+typedef vector<FgVarNode*>                 FgVarSet;
+typedef vector<string>                     Domain;
+typedef vector<unsigned>                   DomainConf;
+typedef pair<unsigned, unsigned>           DomainConstr;
+typedef unordered_map<unsigned, unsigned>  IndexMap;
+
+
+//extern unsigned DL;
+static const unsigned DL = 0;
+
+// number of digits to show when printing a parameter
+static const unsigned PRECISION = 10;
+
+// shared by bp and sp solver
+enum Schedule
+{
+  S_SEQ_FIXED,
+  S_SEQ_RANDOM,
+  S_PARALLEL,
+  S_MAX_RESIDUAL
+};
+
+
+struct NetInfo
+{
+  NetInfo (unsigned c, double t)
+  { 
+    counting    = c;
+    solvingTime = t;
+  }
+  unsigned  counting;
+  double    solvingTime;
+};
+
+typedef map<unsigned, NetInfo> StatisticMap; 
+
+
+class Statistics
+{
+  public:
+
+    static void updateStats (unsigned size, double time)
+    {
+      StatisticMap::iterator it  = stats_.find(size);
+      if (it == stats_.end()) {
+        stats_.insert (make_pair (size, NetInfo (1, 0.0)));
+      } else {
+        it->second.counting ++;
+        it->second.solvingTime += time;
+      }
+    }
+
+    static unsigned getCounting (unsigned size)
+    {
+      StatisticMap::iterator it = stats_.find(size);
+      assert (it != stats_.end());
+      return it->second.counting;
+    }
+
+    static void updateIterations (unsigned nIters)
+    {
+      totalOfIterations += nIters;
+      if (nIters > maxIterations) { 
+        maxIterations = nIters;
+      }
+    }
+
+    static void writeStats (void)
+    {
+      ofstream out ("../../stats.txt");
+      if (!out.is_open()) {
+        cerr << "error: cannot open file to write at " ;
+        cerr << "Statistics:::updateStats()" << endl;
+        abort();
+      }
+      unsigned avgIterations = 0;
+      if (numSolvedLoopyNets > 0) {
+        avgIterations = totalOfIterations / numSolvedLoopyNets;
+      }
+      double totalSolvingTime = 0.0;
+      for (StatisticMap::iterator it = stats_.begin();
+          it != stats_.end(); it++) {
+        totalSolvingTime += it->second.solvingTime;
+      }
+      out << "created networks:              " << numCreatedNets     << endl;
+      out << "solver runs on polytrees:      " << numSolvedPolyTrees << endl;
+      out << "solver runs on loopy networks: " << numSolvedLoopyNets << endl;
+      out << "   unconverged:                " << numUnconvergedRuns << endl;
+      out << "   max iterations:             " << maxIterations      << endl;
+      out << "   average iterations:         " << avgIterations      << endl;
+      out << "total solving time             " << totalSolvingTime   << endl;
+      out << endl;
+      out << "Network Size\tCounting\tSolving Time\tAverage Time" << endl;
+      for (StatisticMap::iterator it = stats_.begin();
+          it != stats_.end(); it++) {
+        out << it->first;
+        out << "\t\t" << it->second.counting;
+        out << "\t\t" << it->second.solvingTime;
+        if (it->second.counting > 0) {
+          out << "\t\t" << it->second.solvingTime / it->second.counting;
+        } else {
+          out << "\t\t0.0" ;
+        }
+        out << endl;
+      }
+      out.close();
+    }
+
+    static unsigned numCreatedNets;
+    static unsigned numSolvedPolyTrees;
+    static unsigned numSolvedLoopyNets;
+    static unsigned numUnconvergedRuns;
+
+  private:
+    static StatisticMap stats_;
+    static unsigned maxIterations;
+    static unsigned totalOfIterations;
+
+};
+
+
+
+class Util
+{
+  public:
+    static void normalize (ParamSet& v)
+    {
+     double sum = 0.0;
+      for (unsigned i = 0; i < v.size(); i++) {
+        sum += v[i];
+      }
+      assert (sum != 0.0);
+      for (unsigned i = 0; i < v.size(); i++) {
+        v[i] /= sum;
+      }
+    }
+
+    static double getL1dist (const ParamSet& v1, const ParamSet& v2)
+    {
+      assert (v1.size() == v2.size());
+      double dist = 0.0;
+      for (unsigned i = 0; i < v1.size(); i++) {
+        dist += abs (v1[i] - v2[i]);
+      }
+      return dist;
+    }
+
+    static double getMaxNorm (const ParamSet& v1, const ParamSet& v2)
+    {
+      assert (v1.size() == v2.size());
+      double max = 0.0;
+      for (unsigned i = 0; i < v1.size(); i++) {
+        double diff = abs (v1[i] - v2[i]);
+        if (diff > max) {
+          max = diff;
+        }
+      }
+      return max;
+    }
+
+    static bool isInteger (const string& s)
+    {
+      stringstream ss1 (s);
+      stringstream ss2;
+      int integer;
+      ss1 >> integer;
+      ss2 << integer;
+      return (ss1.str() == ss2.str());
+    }
+};
+
+
+//unsigned Statistics::totalOfIterations  = 0;
+
+#endif
+

packages/CLPBN/clpbn/bp/Solver.h

@@ -0,0 +1,50 @@
+#ifndef BP_SOLVER_H
+#define BP_SOLVER_H
+
+#include <iomanip>
+
+#include "GraphicalModel.h"
+#include "Variable.h"
+
+using namespace std;
+
+class Solver
+{
+  public:
+    Solver (const GraphicalModel* gm)
+    {
+      gm_ = gm;
+    }
+    virtual void     runSolver (void) = 0;
+    virtual ParamSet getPosterioriOf  (const Variable*) const = 0;
+
+    void printPosterioriOf (const Variable* var) const
+    {
+      cout << endl;
+      cout << setw (20) << left << var->getLabel() << "posteriori" ;
+      cout << endl;
+      cout << "------------------------------" ;
+      cout << endl;
+      const Domain& domain = var->getDomain();
+      ParamSet results = getPosterioriOf (var);
+      for (int xi = 0; xi < var->getDomainSize(); xi++) {
+        cout << setw (20) << domain[xi];
+        cout << setprecision (PRECISION) << results[xi];
+        cout << endl;
+      }
+      cout << endl;
+    }
+
+    void printAllPosterioris (void) const
+    {
+      VarSet vars = gm_->getVariables();
+      for (unsigned i = 0; i < vars.size(); i++) {
+        printPosterioriOf (vars[i]);
+      }
+    }
+
+  private:
+     const GraphicalModel* gm_;
+};
+
+#endif

packages/CLPBN/clpbn/bp/Variable.h

@@ -0,0 +1,143 @@
+#ifndef BP_GENERIC_VARIABLE_H
+#define BP_GENERIC_VARIABLE_H
+
+#include <sstream>
+
+#include <algorithm>
+#include "Shared.h"
+
+using namespace std;
+
+class Variable
+{
+  public:
+
+    Variable (unsigned varId)
+    {
+      this->varId_      = varId;
+      this->dsize_      = 0;
+      this->evidence_   = -1;
+      this->label_      = 0;
+    }
+
+    Variable (unsigned varId, unsigned dsize, int evidence = -1)
+    {
+      assert (dsize != 0);
+      assert (evidence < (int)dsize);
+      this->varId_      = varId;
+      this->dsize_      = dsize;
+      this->evidence_   = evidence;
+      this->label_      = 0;
+    }
+
+    Variable (unsigned varId, const Domain& domain, int evidence = -1)
+    {
+      assert (!domain.empty());
+      assert (evidence < (int)domain.size());
+      this->varId_      = varId;
+      this->dsize_      = domain.size();
+      this->domain_     = domain;
+      this->evidence_   = evidence;
+      this->label_      = 0;
+    }
+
+   ~Variable (void)
+    {
+      delete label_;
+    }
+
+    unsigned       getVarId (void) const       { return varId_; }
+    unsigned       getIndex (void) const       { return index_; }
+    void           setIndex (unsigned idx)     { index_ = idx; }
+    int            getDomainSize (void) const  { return dsize_; }
+    bool           hasEvidence (void) const    { return evidence_ != -1; }
+    int            getEvidence (void) const    { return evidence_; }
+    bool           hasDomain (void)            { return !domain_.empty(); }
+    bool           hasLabel (void)             { return label_ != 0; }
+
+   bool isValidStateIndex (int index)
+   {
+     return index >= 0 && index < dsize_;
+   }
+
+    bool isValidState (const string& state)
+    {
+      return find (domain_.begin(), domain_.end(), state) != domain_.end();
+    }
+
+    Domain getDomain (void) const
+    {
+      assert (dsize_ != 0);
+      if (domain_.size() == 0) {
+        Domain d;
+        for (int i = 0; i < dsize_; i++) {
+          stringstream ss;
+          ss << "x" << i ;
+          d.push_back (ss.str());
+        }
+        return d;
+      } else {
+        return domain_;
+      }
+    }
+
+    void setDomainSize (unsigned dsize)
+    {
+      assert (dsize != 0);
+      dsize_ = dsize;
+    }
+
+    void setDomain (const Domain& domain)
+    {
+      assert (!domain.empty());
+      domain_ = domain; 
+      dsize_  = domain.size();
+    }
+
+    void setEvidence (int ev) 
+    { 
+      assert (ev < dsize_);
+      evidence_ = ev;
+    }
+
+    void setEvidence (const string& ev) 
+    { 
+      assert (isValidState (ev));
+      for (unsigned i = 0; i < domain_.size(); i++) {
+        if (domain_[i] == ev) {
+          evidence_ = i;
+        }
+      }
+    }
+
+    void setLabel (string label)
+    {
+      label_ = new string (label);
+    }
+
+    string getLabel (void) const
+    {
+      if (label_ == 0) {
+        stringstream ss;
+        ss << "v" << varId_;
+        return ss.str();
+      } else {
+        return *label_;
+      }
+    }
+
+  protected:
+    unsigned       varId_;
+    string*        label_;
+    unsigned       index_;
+    int            evidence_;
+
+  private:
+    DISALLOW_COPY_AND_ASSIGN (Variable);
+    Domain         domain_;
+    int            dsize_;
+
+};
+
+#endif // BP_GENERIC_VARIABLE_H
+

packages/CLPBN/clpbn/bp/callgrind.h

@@ -0,0 +1,147 @@
+
+/*
+   ----------------------------------------------------------------
+
+   Notice that the following BSD-style license applies to this one
+   file (callgrind.h) only.  The rest of Valgrind is licensed under the
+   terms of the GNU General Public License, version 2, unless
+   otherwise indicated.  See the COPYING file in the source
+   distribution for details.
+
+   ----------------------------------------------------------------
+
+   This file is part of callgrind, a valgrind tool for cache simulation
+   and call tree tracing.
+
+   Copyright (C) 2003-2010 Josef Weidendorfer.  All rights reserved.
+
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   1. Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+   2. The origin of this software must not be misrepresented; you must
+      not claim that you wrote the original software.  If you use this
+      software in a product, an acknowledgment in the product
+      documentation would be appreciated but is not required.
+
+   3. Altered source versions must be plainly marked as such, and must
+      not be misrepresented as being the original software.
+
+   4. The name of the author may not be used to endorse or promote
+      products derived from this software without specific prior written
+      permission.
+
+   THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS
+   OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+   WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+   ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
+   DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+   DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
+   GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+   INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+   WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+   ----------------------------------------------------------------
+
+   Notice that the above BSD-style license applies to this one file
+   (callgrind.h) only.  The entire rest of Valgrind is licensed under
+   the terms of the GNU General Public License, version 2.  See the
+   COPYING file in the source distribution for details.
+
+   ----------------------------------------------------------------
+*/
+
+#ifndef __CALLGRIND_H
+#define __CALLGRIND_H
+
+#include "valgrind.h"
+
+/* !! ABIWARNING !! ABIWARNING !! ABIWARNING !! ABIWARNING !!
+   This enum comprises an ABI exported by Valgrind to programs
+   which use client requests.  DO NOT CHANGE THE ORDER OF THESE
+   ENTRIES, NOR DELETE ANY -- add new ones at the end.
+
+   The identification ('C','T') for Callgrind has historical
+   reasons: it was called "Calltree" before. Besides, ('C','G') would
+   clash with cachegrind.
+ */
+
+typedef
+   enum {
+      VG_USERREQ__DUMP_STATS = VG_USERREQ_TOOL_BASE('C','T'),
+      VG_USERREQ__ZERO_STATS,
+      VG_USERREQ__TOGGLE_COLLECT,
+      VG_USERREQ__DUMP_STATS_AT,
+      VG_USERREQ__START_INSTRUMENTATION,
+      VG_USERREQ__STOP_INSTRUMENTATION
+   } Vg_CallgrindClientRequest;
+
+/* Dump current state of cost centers, and zero them afterwards */
+#define CALLGRIND_DUMP_STATS                                            \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__DUMP_STATS,                     \
+                            0, 0, 0, 0, 0);                             \
+   }
+
+/* Dump current state of cost centers, and zero them afterwards.
+   The argument is appended to a string stating the reason which triggered
+   the dump. This string is written as a description field into the
+   profile data dump. */
+#define CALLGRIND_DUMP_STATS_AT(pos_str)                                \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__DUMP_STATS_AT,                  \
+                            pos_str, 0, 0, 0, 0);                       \
+   }
+
+/* Zero cost centers */
+#define CALLGRIND_ZERO_STATS                                            \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__ZERO_STATS,                     \
+                            0, 0, 0, 0, 0);                             \
+   }
+
+/* Toggles collection state.
+   The collection state specifies whether the happening of events
+   should be noted or if they are to be ignored. Events are noted
+   by increment of counters in a cost center */
+#define CALLGRIND_TOGGLE_COLLECT                                        \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__TOGGLE_COLLECT,                 \
+                            0, 0, 0, 0, 0);                             \
+   }
+
+/* Start full callgrind instrumentation if not already switched on.
+   When cache simulation is done, it will flush the simulated cache;
+   this will lead to an artifical cache warmup phase afterwards with
+   cache misses which would not have happened in reality. */
+#define CALLGRIND_START_INSTRUMENTATION                                 \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__START_INSTRUMENTATION,          \
+                            0, 0, 0, 0, 0);                             \
+   }
+
+/* Stop full callgrind instrumentation if not already switched off.
+   This flushes Valgrinds translation cache, and does no additional
+   instrumentation afterwards, which effectivly will run at the same
+   speed as the "none" tool (ie. at minimal slowdown).
+   Use this to bypass Callgrind aggregation for uninteresting code parts.
+   To start Callgrind in this mode to ignore the setup phase, use
+   the option "--instr-atstart=no". */
+#define CALLGRIND_STOP_INSTRUMENTATION                                  \
+   {unsigned int _qzz_res;                                              \
+    VALGRIND_DO_CLIENT_REQUEST(_qzz_res, 0,                             \
+                            VG_USERREQ__STOP_INSTRUMENTATION,           \
+                            0, 0, 0, 0, 0);                             \
+   }
+
+#endif /* __CALLGRIND_H */

packages/CLPBN/clpbn/bp/examples/bayes-ball a.xml

@@ -0,0 +1,76 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Bayes-Ball: The Rational Pastime Network, Figure 4, a)</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>1</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>2</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>3</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>4</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>5</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>6</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>1</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>2</FOR>
+	<GIVEN>1</GIVEN>
+	<GIVEN>3</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>3</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>4</FOR>
+	<GIVEN>1</GIVEN>
+	<GIVEN>5</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>5</FOR>
+	<GIVEN>2</GIVEN>
+	<GIVEN>6</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>6</FOR>
+	<GIVEN>3</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/bayes-ball c.xml

@@ -0,0 +1,74 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Bayes-Ball: The Rational Pastime Network, Figure 4, c)</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>1</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>2</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>3</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>4</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>5</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>6</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>1</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>2</FOR>
+	<GIVEN>1</GIVEN>
+	<GIVEN>3</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>3</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>4</FOR>
+	<GIVEN>5</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>5</FOR>
+	<GIVEN>2</GIVEN>
+	<GIVEN>6</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>6</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/burglary-alarm.uai

@@ -0,0 +1,28 @@
+MARKOV
+5
+2 2 2 2 2
+5
+1 0
+1 1
+3 2 0 1
+2 3 2
+2 4 2
+
+2
+ .001 .009
+
+2
+ .002 .008 
+
+8
+ .95 .94 .29 .001
+ .05 .06 .71 .999
+
+4
+ .9 .05
+ .1 .95
+
+4
+ .7 .01 
+ .3 .99
+ 

packages/CLPBN/clpbn/bp/examples/burglary-alarm.xml

@@ -0,0 +1,81 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+    B     E
+     \   /
+      \ / 
+       A 
+      / \
+     /   \
+    J     M
+
+-->
+
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Simple Loop</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME>b1</OUTCOME>
+	<OUTCOME>b2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>E</NAME>
+	<OUTCOME>e1</OUTCOME>
+	<OUTCOME>e2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME>a1</OUTCOME>
+	<OUTCOME>a2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>J</NAME>
+	<OUTCOME>j1</OUTCOME>
+	<OUTCOME>j2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>M</NAME>
+	<OUTCOME>m1</OUTCOME>
+	<OUTCOME>m2</OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<TABLE> .001 .009 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>E</FOR>
+	<TABLE> .002 .008 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<GIVEN>B</GIVEN>
+	<GIVEN>E</GIVEN>
+	<TABLE> .95 .05 .94 .06 .29 .71 .001 .999 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>J</FOR>
+	<GIVEN>A</GIVEN>
+	<TABLE> .9 .1 .05 .95  </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>M</FOR>
+	<GIVEN>A</GIVEN>
+	<TABLE> .7 .3 .01 .99 </TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/burglary-alarm.yap

@@ -0,0 +1,54 @@
+
+:- use_module(library(clpbn)).
+
+:- set_clpbn_flag(solver, vel).
+
+%
+%    B     E
+%     \   /
+%      \ / 
+%       A 
+%      / \
+%     /   \
+%    J     M
+%
+
+
+b(B) :-
+    b_table(BDist),
+    { B = b with p([b1, b2], BDist) }.
+
+e(E) :-
+    e_table(EDist),
+    { E = e with p([e1, e2], EDist) }.
+
+a(A) :-
+    b(B),
+    e(E),
+    a_table(ADist),
+    { A = a with p([a1, a2], ADist, [B, E]) }.
+
+j(J):-
+    a(A),
+    j_table(JDist),
+    { J = j with p([j1, j2], JDist, [A]) }.
+
+m(M):-
+    a(A),
+    m_table(MDist),
+    { M = m with p([m1, m2], MDist, [A]) }.
+    
+
+b_table([0.001, 0.009]).
+
+e_table([0.002, 0.008]).
+
+a_table([0.95, 0.94, 0.29, 0.001,
+         0.05, 0.06, 0.71, 0.999]).
+
+j_table([0.9, 0.05,
+         0.1, 0.95]).
+
+m_table([0.7, 0.01,
+         0.3, 0.99]).
+

packages/CLPBN/clpbn/bp/examples/chain.xml

@@ -0,0 +1,58 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+  A
+  |
+  |
+  -
+  B
+  |
+  |
+  -
+  C
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Simple Chain</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME>a1</OUTCOME>
+	<OUTCOME>a2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME>b1</OUTCOME>
+	<OUTCOME>b2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>C</NAME>
+	<OUTCOME>c1</OUTCOME>
+	<OUTCOME>c2</OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<TABLE>0.3 0.7</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<GIVEN>A</GIVEN>
+	<TABLE>0.4 0.6 0.2 0.8</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>C</FOR>
+	<GIVEN>B</GIVEN>
+	<TABLE>0.9 0.1 0.25 0.75</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/convergence.xml

@@ -0,0 +1,51 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+  A     B
+   \   /
+    \ /
+     -
+     C
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Simple Convergence</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>C</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>C</FOR>
+	<GIVEN>A</GIVEN>
+	<GIVEN>B</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/divergence.xml

@@ -0,0 +1,51 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+     A
+    / \
+   /   \
+  -     -
+  B     C
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Simple Divergence</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>C</NAME>
+	<OUTCOME></OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<GIVEN>A</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>C</FOR>
+	<GIVEN>A</GIVEN>
+	<TABLE>1</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/john-mary-call.xml

@@ -0,0 +1,106 @@
+<?XML VERSION="1.0"?>
+
+
+<!--
+	Bayesian network in BIF (BayesNet Interchange Format)
+	Produced by JavaBayes (http://www.cs.cmu.edu/~javabayes/
+	Output created Fri Nov 14 13:14:15 GMT+00:00 1997
+-->
+
+
+
+<!-- DTD for the BIF format -->
+<!DOCTYPE BIF [
+	<!ELEMENT PROPERTY (#PCDATA)>
+	<!ELEMENT TYPE (#PCDATA)>
+	<!ELEMENT VALUE (#PCDATA)>
+	<!ELEMENT NAME (#PCDATA)>
+	<!ELEMENT NETWORK
+		( NAME, ( PROPERTY | VARIABLE | PROBABILITY )* )>
+	<!ELEMENT VARIABLE ( NAME, TYPE, ( VALUE |  PROPERTY )* ) >
+	<!ELEMENT PROBABILITY
+		( FOR | GIVEN | TABLE | ENTRY | DEFAULT | PROPERTY )* >
+	<!ELEMENT TABLE (#PCDATA)>
+	<!ELEMENT DEFAULT (TABLE)>
+	<!ELEMENT ENTRY ( VALUE* , TABLE )>
+]>
+
+
+<BIF>
+<NETWORK>
+<NAME>John-Mary-Call</NAME>
+
+<!-- Variables -->
+<VARIABLE>
+	<NAME>Burglary</NAME>
+	<TYPE>discrete</TYPE>
+	<VALUE>False</VALUE>
+	<VALUE>True</VALUE>
+	<PROPERTY>position = (145, 114)</PROPERTY>
+</VARIABLE>
+
+<VARIABLE>
+	<NAME>Earthquake</NAME>
+	<TYPE>discrete</TYPE>
+	<VALUE>False</VALUE>
+	<VALUE>True</VALUE>
+	<PROPERTY>position = (351, 110)</PROPERTY>
+</VARIABLE>
+
+<VARIABLE>
+	<NAME>Alarm</NAME>
+	<TYPE>discrete</TYPE>
+	<VALUE>False</VALUE>
+	<VALUE>True</VALUE>
+	<PROPERTY>position = (253, 224)</PROPERTY>
+</VARIABLE>
+
+<VARIABLE>
+	<NAME>JohnCalls</NAME>
+	<TYPE>discrete</TYPE>
+	<VALUE>False</VALUE>
+	<VALUE>True</VALUE>
+	<PROPERTY>position = (156, 343)</PROPERTY>
+</VARIABLE>
+
+<VARIABLE>
+	<NAME>MaryCalls</NAME>
+	<TYPE>discrete</TYPE>
+	<VALUE>False</VALUE>
+	<VALUE>True</VALUE>
+	<PROPERTY>position = (344, 341)</PROPERTY>
+</VARIABLE>
+
+
+<!-- Probability distributions -->
+<PROBABILITY>
+	<FOR>Burglary</FOR>
+	<TABLE>0.999 0.0010 </TABLE>
+</PROBABILITY>
+
+<PROBABILITY>
+	<FOR>Earthquake</FOR>
+	<TABLE>0.998 0.0020 </TABLE>
+</PROBABILITY>
+
+<PROBABILITY>
+	<FOR>Alarm</FOR>
+	<GIVEN>Burglary</GIVEN>
+	<GIVEN>Earthquake</GIVEN>
+	<TABLE>0.999 0.71 0.06 0.05 0.0010 0.29 0.94 0.95 </TABLE>
+</PROBABILITY>
+
+<PROBABILITY>
+	<FOR>JohnCalls</FOR>
+	<GIVEN>Alarm</GIVEN>
+	<TABLE>0.95 0.1 0.05 0.9 </TABLE>
+</PROBABILITY>
+
+<PROBABILITY>
+	<FOR>MaryCalls</FOR>
+	<GIVEN>Alarm</GIVEN>
+	<TABLE>0.99 0.3 0.01 0.7 </TABLE>
+</PROBABILITY>
+
+
+</BIF>

packages/CLPBN/clpbn/bp/examples/loop.xml

@@ -0,0 +1,81 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+        A     E
+       / \   /
+      /   \ /
+     B     C 
+      \   /
+       \ / 
+        D
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Loop</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME>a1</OUTCOME>
+	<OUTCOME>a2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME>b1</OUTCOME>
+	<OUTCOME>b2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>C</NAME>
+	<OUTCOME>c1</OUTCOME>
+	<OUTCOME>c2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>D</NAME>
+	<OUTCOME>d1</OUTCOME>
+	<OUTCOME>d2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>E</NAME>
+	<OUTCOME>e1</OUTCOME>
+	<OUTCOME>e2</OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<TABLE> .01 .09 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>B</FOR>
+    <GIVEN>A</GIVEN>
+	<TABLE> .03 .97 .6 .4 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>C</FOR>
+    <GIVEN>A</GIVEN>
+    <GIVEN>E</GIVEN>
+	<TABLE> .24 .76 .12 .88 .2 .4. 5. .6 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>D</FOR>
+	<GIVEN>B</GIVEN>
+	<GIVEN>C</GIVEN>
+	<TABLE> .2 .8 .7 .3 .45 .55 .22 .78 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>E</FOR>
+	<TABLE> .5 .6</TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/loop.yap

@@ -0,0 +1,53 @@
+
+:- use_module(library(clpbn)).
+
+:- set_clpbn_flag(solver, bp).
+
+%
+%        A     E
+%       / \   /
+%      /   \ /
+%     B     C 
+%      \   /
+%       \ / 
+%        D
+%
+
+a(A) :-
+    a_table(ADist),
+    { A = a with p([a1, a2], ADist) }.
+
+b(B) :-
+    a(A),
+    b_table(BDist),
+    { B = b with p([b1, b2], BDist, [A]) }.
+
+c(C) :-
+    a(A),
+    c_table(CDist),
+    { C = c with p([c1, c2], CDist, [A]) }.
+
+d(D) :-
+    b(B),
+    c(C),
+    d_table(DDist),
+    { D = d with p([d1, d2], DDist, [B, C]) }.
+
+e(E) :-
+    e_table(EDist),
+    { E = e with p([e1, e2], EDist) }.
+
+
+a_table([0.005, 0.995]).
+
+b_table([0.02, 0.97,
+         0.88, 0.03]).
+
+c_table([0.55, 0.94,
+         0.45, 0.06]).
+
+d_table([0.192, 0.98, 0.33, 0.013,
+         0.908, 0.02, 0.77, 0.987]).
+
+e_table([0.055, 0.945]).
+

packages/CLPBN/clpbn/bp/examples/neapolitan-A-F-joint.yap

@@ -0,0 +1,55 @@
+
+:- use_module(library(clpbn)).
+
+:- set_clpbn_flag(solver, bp).
+
+%
+%     B     F
+%      \   /
+%       \ / 
+%        A
+%
+
+b(B) :-
+    b_table(BDist),
+    { B = b with p([b1, b2], BDist) }.
+
+f(F) :-
+    f_table(FDist),
+    { F = f with p([f1, f2], FDist) }.
+
+a(A) :-
+    b(B),
+    f(F),
+    a_table(ADist),
+    { A = a with p([a1, a2], ADist, [B, F]) }.
+
+d(D) :-
+    a(A),
+    f(F),
+    d_table(DDist),
+    { D = d with p([d1, d2, d3, d4], DDist, [A, F]) }.
+ 
+
+b_table([0.005, 0.995]).
+
+f_table([0.03, 0.97]).
+
+a_table([0.992, 0.99, 0.2, 0.003,
+         0.008, 0.01, 0.8, 0.997]).
+
+d_table([1.0, 0.0, 0.0, 0.0,
+         0.0, 1.0, 0.0, 0.0,
+         0.0, 0.0, 1.0, 0.0,
+         0.0, 0.0, 0.0, 1.0]).
+
+%d_table([0.997, 0.001, 0.001, 0.001,
+%         0.001, 0.997, 0.001, 0.001,
+%         0.001, 0.001, 0.997, 0.001,
+%         0.001, 0.001, 0.001, 0.997]).
+
+%d_table([0.15, 0.1, 0.7, 0.5,
+%         0.25, 0.3, 0.2, 0.25,
+%         0.3, 0.15, 0.35, 0.2,
+%         0.3, 0.4, 0.2, 0.1]).
+

packages/CLPBN/clpbn/bp/examples/neapolitan.uai

@@ -0,0 +1,17 @@
+MARKOV
+3
+2 2 2
+3
+1 0
+1 1
+3 2 0 1
+
+2
+ 0.005 0.995
+
+2
+ 0.03 0.97
+
+8
+ 0.992 0.99 0.2 0.003
+ 0.008 0.01 0.8 0.997

packages/CLPBN/clpbn/bp/examples/neapolitan.xml

@@ -0,0 +1,53 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+
+     B     F
+      \   /
+       \ / 
+        A
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Neapolitan</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME>b1</OUTCOME>
+	<OUTCOME>b2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>F</NAME>
+	<OUTCOME>f1</OUTCOME>
+	<OUTCOME>f2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME>a1</OUTCOME>
+	<OUTCOME>a2</OUTCOME>
+</VARIABLE>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<TABLE> .005 .995 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>F</FOR>
+	<TABLE> .03 .97 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<GIVEN>B</GIVEN>
+	<GIVEN>F</GIVEN>
+	<TABLE> .992 .008 .99 .01 .2 .8 .003 .997 </TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/neapolitan.yap

@@ -0,0 +1,35 @@
+
+:- use_module(library(clpbn)).
+
+:- set_clpbn_flag(solver, bp).
+
+%
+%     B     F
+%      \   /
+%       \ / 
+%        A
+%
+
+
+b(B) :-
+    b_table(BDist),
+    { B = b with p([b1, b2], BDist) }.
+
+f(F) :-
+    f_table(FDist),
+    { F = f with p([f1, f2], FDist) }.
+
+a(A) :-
+    b(B),
+    f(F),
+    a_table(ADist),
+    { A = a with p([a1, a2], ADist, [B, F]) }.
+
+
+b_table([0.005, 0.995]).
+
+f_table([0.03, 0.97]).
+
+a_table([0.992, 0.99, 0.2, 0.003,
+         0.008, 0.01, 0.8, 0.997]).
+

packages/CLPBN/clpbn/bp/examples/several-parents-and-childs.xml

@@ -0,0 +1,128 @@
+<?xml version="1.0" encoding="US-ASCII"?>
+
+<!--
+ 
+ A   B   C
+  \  |  /
+   \ | /
+     D
+   / | \
+  /  |  \
+ E   F   G
+
+-->
+
+<BIF VERSION="0.3">
+<NETWORK>
+<NAME>Node with several parents and childs</NAME>
+
+<VARIABLE TYPE="nature">
+	<NAME>A</NAME>
+	<OUTCOME>a1</OUTCOME>
+	<OUTCOME>a2</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>B</NAME>
+	<OUTCOME>b1</OUTCOME>
+	<OUTCOME>b2</OUTCOME>
+	<OUTCOME>b3</OUTCOME>
+	<OUTCOME>b4</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>C</NAME>
+	<OUTCOME>c1</OUTCOME>
+	<OUTCOME>c2</OUTCOME>
+	<OUTCOME>c3</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>D</NAME>
+	<OUTCOME>d1</OUTCOME>
+	<OUTCOME>d2</OUTCOME>
+	<OUTCOME>d3</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>E</NAME>
+	<OUTCOME>e1</OUTCOME>
+	<OUTCOME>e2</OUTCOME>
+	<OUTCOME>e3</OUTCOME>
+	<OUTCOME>e4</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>F</NAME>
+	<OUTCOME>f1</OUTCOME>
+	<OUTCOME>f2</OUTCOME>
+	<OUTCOME>f3</OUTCOME>
+</VARIABLE>
+
+<VARIABLE TYPE="nature">
+	<NAME>G</NAME>
+	<OUTCOME>g1</OUTCOME>
+	<OUTCOME>g2</OUTCOME>
+</VARIABLE>
+
+
+<DEFINITION>
+	<FOR>A</FOR>
+	<TABLE> .1 .2 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>B</FOR>
+	<TABLE> .01 .02 .03 .04 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>C</FOR>
+	<TABLE> .11 .22 .33 </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>D</FOR>
+	<GIVEN>A</GIVEN>
+	<GIVEN>B</GIVEN>
+	<GIVEN>C</GIVEN>
+	<TABLE>
+      .522 .008 .99 .01 .2 .8 .003 .457 .423 .007 .92 .04 .5 .232 .033 .227 .112 .048 .91 .21 .24 .18 .005 .227
+      .212 .04 .59 .21 .6 .1 .023 .215 .913 .017 .96 .01 .55 .422 .013 .417 .272 .068 .61 .11 .26 .28 .205 .322 
+      .142 .028 .19 .11 .5 .67 .013 .437 .163 .067 .12 .06 .1 .262 .063 .167 .512 .028 .11 .41 .14 .68 .015 .92  
+    </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>E</FOR>
+	<GIVEN>D</GIVEN>
+	<TABLE> 
+      .111 .11 .1
+      .222 .22 .2
+      .333 .33 .3
+      .444 .44 .4
+    </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>F</FOR>
+	<GIVEN>D</GIVEN>
+	<TABLE> 
+      .112 .111 .110
+      .223 .222 .221
+      .334 .333 .332
+    </TABLE>
+</DEFINITION>
+
+<DEFINITION>
+	<FOR>G</FOR>
+	<GIVEN>D</GIVEN>
+	<TABLE>
+      .101 .102 .103
+      .201 .202 .203
+    </TABLE>
+</DEFINITION>
+
+</NETWORK>
+</BIF>
+

packages/CLPBN/clpbn/bp/examples/test.uai

@@ -0,0 +1,36 @@
+MARKOV
+5
+4 2 3 2 3
+7
+1 0
+1 1
+1 2
+1 3
+1 4
+2 0 1
+4 1 2 3 4
+
+4
+ 0.1 0.7 0.43 0.22
+
+2
+ 0.2 0.6
+
+3
+ 0.3 0.5 0.2
+
+2
+ 0.15 0.75
+
+3
+ 0.25 0.45 0.15
+
+8
+ 0.210 0.333 0.457 0.4
+ 0.811 0.000 0.189 0.89
+
+36
+ 0.1 0.15 0.2 0.25 0.3 0.45 0.5 0.55 0.65 0.7 0.75 0.9
+ 0.11 0.22 0.33 0.44 0.55 0.66 0.77 0.88 0.91 0.93 0.95 0.97
+ 0.42 0.22 0.33 0.44 0.15 0.36 0.27 0.28 0.21 0.13 0.25 0.17
+

packages/CLPBN/clpbn/graphviz.yap

@@ -22,7 +22,7 @@ output_var(Stream, V) :-
 	Parents = [_|_], !,
 	format(Stream, '	',[]),
 	output_parents(Stream, Parents),
-	format(' -> ',[]),
+	format(Stream,' -> ',[]),
 	output_key(Stream,Key),
 	nl(Stream).
 output_var(_, _).

packages/CLPBN/clpbn/table.yap

@@ -96,7 +96,7 @@ clpbn_table(F/N,M) :-
 	Key =.. L1,
 	atom_concat(F, '___tabled', NF),
 	L2 = [_|Args],
-	S1 =.. [NF|Args],
+	_S1 =.. [NF|Args],
 	L0 = [_|OArgs],
 	S2 =.. [NF|OArgs],
 	asserta(clpbn_table(S, M, S2)),

packages/CLPBN/clpbn/ve.yap

@@ -14,10 +14,10 @@
   
 *********************************/
 
-:- module(clpbn_vel, [vel/3,
+:- module(clpbn_ve, [ve/3,
-		check_if_vel_done/1,
+		check_if_ve_done/1,
-		init_vel_solver/4,
+		init_ve_solver/4,
-		run_vel_solver/3]).
+		run_ve_solver/3]).
 
 :- attribute size/1, all_diffs/1.
 
@@ -64,40 +64,41 @@
 	      [check_for_agg_vars/2]).
 
 
-check_if_vel_done(Var) :-
+check_if_ve_done(Var) :-
 	get_atts(Var, [size(_)]), !.
 
 %
 % implementation of the well known variable elimination algorithm
 %
-vel([[]],_,_) :- !.
+ve([[]],_,_) :- !.
-vel([LVs],Vs0,AllDiffs) :-
+ve([LVs],Vs0,AllDiffs) :-
-	init_vel_solver([LVs], Vs0, AllDiffs, State),
+	init_ve_solver([LVs], Vs0, AllDiffs, State),
 	% variable elimination proper
-	run_vel_solver([LVs], [LPs], State),
+	run_ve_solver([LVs], [LPs], State),
 	% bind Probs back to variables so that they can be output.
 	clpbn_bind_vals([LVs],[LPs],AllDiffs).
 
-init_vel_solver(Qs, Vs0, _, LVis) :-
+init_ve_solver(Qs, Vs0, _, LVis) :-
 	check_for_agg_vars(Vs0, Vs1),
 	% LVi will have a  list of CLPBN variables
 	% Tables0 will have the full data on each variable
 	init_influences(Vs1, G, RG),
-	init_vel_solver_for_questions(Qs, G, RG, _, LVis).
+	init_ve_solver_for_questions(Qs, G, RG, _, LVis).
 
-init_vel_solver_for_questions([], _, _, [], []).
+init_ve_solver_for_questions([], _, _, [], []).
-init_vel_solver_for_questions([Vs|MVs], G, RG, [NVs|MNVs0], [NVs|LVis]) :-
+init_ve_solver_for_questions([Vs|MVs], G, RG, [NVs|MNVs0], [NVs|LVis]) :-
 	influences(Vs, _, NVs0, G, RG),
 	sort(NVs0, NVs),
 %clpbn_gviz:clpbn2gviz(user_error, test, NVs, Vs),
-	init_vel_solver_for_questions(MVs, G, RG, MNVs0, LVis).
+	init_ve_solver_for_questions(MVs, G, RG, MNVs0, LVis).
 
 % use a findall to recover space without needing for GC
-run_vel_solver(LVs, LPs, LNVs) :-
+run_ve_solver(LVs, LPs, LNVs) :-
-	findall(Ps, solve_vel(LVs, LNVs, Ps), LPs).
+	findall(Ps, solve_ve(LVs, LNVs, Ps), LPs).
 
-solve_vel([LVs|_], [NVs0|_], Ps) :-
+solve_ve([LVs|_], [NVs0|_], Ps) :-
-%	length(NVs0, L), (L > 64 -> clpbn_gviz:clpbn2gviz(user_error,sort,NVs0,LVs) ; true ),
+%	length(NVs0, L), (L > 415 -> clpbn_gviz:clpbn2gviz(user_error,sort,NVs0,LVs) ; true ),
+%	length(NVs0, L), writeln(+LVs:L),
 	find_all_clpbn_vars(NVs0, NVs0, LV0, LVi, Tables0),
 	sort(LV0, LV),
 	% construct the graph
@@ -108,8 +109,8 @@ solve_vel([LVs|_], [NVs0|_], Ps) :-
 	% move from potentials back to probabilities
 	normalise_CPT(Dist,MPs),
 	list_from_CPT(MPs, Ps).
-solve_vel([_|MoreLVs], [_|MoreLVis], Ps) :-
+solve_ve([_|MoreLVs], [_|MoreLVis], Ps) :-
-	solve_vel(MoreLVs, MoreLVis, Ps).
+	solve_ve(MoreLVs, MoreLVis, Ps).
 
 exps([],[]).
 exps([L|LD],[O|LDE]) :-
@@ -133,7 +134,7 @@ find_all_clpbn_vars([V|Vs], NVs0, [Var|LV], ProcessedVars, [table(I,Table,Parent
 	% variables with evidence should not be processed.
 	(var(Ev) ->
 	    Var = var(V,I,Sz,Vals,Parents,Ev,_,_),
-	    vel_get_dist_size(V,Sz),
+	    ve_get_dist_size(V,Sz),
 	    ProcessedVars = [Var|ProcessedVars0]
 	;
 	    ProcessedVars = ProcessedVars0
@@ -191,7 +192,7 @@ compute_size([tab(_,Vs,_)|Tabs],Vs0,K) :-
 
 multiply_sizes([],K,K).
 multiply_sizes([V|Vs],K0,K) :-
-	vel_get_dist_size(V, Sz),
+	ve_get_dist_size(V, Sz),
 	KI is K0*Sz,
 	multiply_sizes(Vs,KI,K).
 
@@ -280,9 +281,9 @@ 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).
 
-vel_get_dist_size(V,Sz) :-
+ve_get_dist_size(V,Sz) :-
 	get_atts(V, [size(Sz)]), !.
-vel_get_dist_size(V,Sz) :-
+ve_get_dist_size(V,Sz) :-
 	clpbn:get_atts(V,dist(Id,_)), !,
 	get_dist_domain_size(Id,Sz),
 	put_atts(V, [size(Sz)]).

packages/CLPBN/learning/em.yap

@@ -11,6 +11,7 @@
 	      [clpbn_init_graph/1,
 	       clpbn_init_solver/5,
 	       clpbn_run_solver/4,
+	       clpbn_finalize_solver/1,
 	       clpbn_flag/2]).
 
 :- use_module(library('clpbn/dists'),
@@ -53,6 +54,7 @@
 em(Items, MaxError, MaxIts, Tables, Likelihood) :-
 	catch(init_em(Items, State),Error,handle_em(Error)),
 	em_loop(0, 0.0, State, MaxError, MaxIts, Likelihood, Tables),
+	clpbn_finalize_solver(State),
 	assert(em_found(Tables, Likelihood)),
 	fail.
 % get rid of new random variables the easy way :)
@@ -90,7 +92,7 @@ init_em(Items, state( AllDists, AllDistInstances, MargVars, SolverVars)) :-
 em_loop(Its, Likelihood0, State, MaxError, MaxIts, LikelihoodF, FTables) :-
 	estimate(State, LPs),
 	maximise(State, Tables, LPs, Likelihood),
-	writeln(Likelihood:Its:Likelihood0:Tables),
+%	writeln(Likelihood:Its:Likelihood0:Tables),
 	(
 	    (
 	     abs((Likelihood - Likelihood0)/Likelihood) < MaxError
@@ -205,7 +207,7 @@ compute_parameters([Id-Samples|Dists], [Id-NewTable|Tables],  MDistTable, Lik0,
 	empty_dist(Id, Table0),
 	add_samples(Samples, Table0, MDistTable),
 	soften_sample(Table0, SoftenedTable),
-	matrix:matrix_sum(Table0,TotM),
+%	matrix:matrix_sum(Table0,TotM),
 	normalise_counts(SoftenedTable, NewTable),
 	compute_likelihood(Table0, NewTable, DeltaLik),
 	dist_new_table(Id, NewTable),

packages/CLPBN/learning/learn_utils.yap

@@ -38,7 +38,7 @@ run_all(M:Gs) :-
 run_all([],_).
 run_all([G|Gs],M) :-
 %	(G = _:ge(ybr136w,t8,23,-1) -> nb_getval(clpbn_tables, Tab), writeln(Tab) ; true ),
-	( call(M:G) -> true ;  writeln(bad:M:G), start_low_level_trace, M:G ; halt ),
+	( call(M:G) -> true ;  throw(bad_call(M:G)) ),
 	run_all(Gs,M).
 
 clpbn_vars(Vs,BVars) :-

packages/chr

@@ -1 +1 @@
-Subproject commit bf6525f85cfcf3c08fff8cf91fb189fe71dc34fd
+Subproject commit b2eb894ce3e41925070215f800d6df3a356dc29d

packages/meld/examples/pagerank/l4.meld

@@ -1,7 +1,7 @@
 %edge(0,1).
-%edge(0,4).
+edge(0,4).
-%edge(1,4).
 %edge(1,2).
-edge(2,3).
+edge(1,4).
+%edge(2,3).
 edge(2,4).
 edge(3,4).

packages/meld/examples/pagerank/pagerank.meld

@@ -1,6 +1,5 @@
 
 type rank(node, int, float).
-type reachable(node, node).
 type calcRank(node, int, sum float).
 % type persistent numPages(node,  int).
 type persistent numPages(node, sum int).
@@ -8,21 +7,24 @@ type numLinks(node, sum int).
 type path(node, node).
 
 const damping = 0.85.
-const num_iterations = 4.
+const num_iterations = 100.
 % extern float to_float(int).
 % extern float float_abs(float).
 
 rank(A, 0, 1.0 / to_float(T)) :- numPages(A,T).
 rank(A, I, V) :-
+	numLinks(B,L),
 	numPages(A, Ps),
 	calcRank(A, I, T),
-	Before = I - 1,
+%	Before = I - 1,
-	rank(A, Before, VOld),
+%	rank(A, Before, VOld),
-	V = (damping + (1.0 - damping) * T)/to_float(Ps),
+	V = damping + (1.0 - damping) * T,
 	I =< num_iterations.
 
 % //float_abs((damping + (1.0 - damping) * T) - VOld) > 0.001.
 
+calcRank(A, I + 1, 0.0) :-
+	 rank(A, I, _).
 calcRank(A, I + 1, O / to_float(C)) :-
 	edge(B, A),
 	rank(B, I, O),