diff --git a/packages/CLPBN/Makefile.in b/packages/CLPBN/Makefile.in
index 2f908d017..b49f60cfc 100644
--- a/packages/CLPBN/Makefile.in
+++ b/packages/CLPBN/Makefile.in
@@ -24,7 +24,8 @@ INSTALL_DATA=@INSTALL_DATA@
 INSTALL_PROGRAM=@INSTALL_PROGRAM@
 srcdir=@srcdir@
 
-CLPBN_TOP= $(srcdir)/clpbn.yap
+CLPBN_TOP= $(srcdir)/clpbn.yap \
+	$(srcdir)/pfl.yap
 
 CLPBN_SRCDIR = $(srcdir)/clpbn
 
@@ -41,9 +42,11 @@ CLPBN_PROGRAMS= \
 	$(CLPBN_SRCDIR)/display.yap \
 	$(CLPBN_SRCDIR)/dists.yap \
 	$(CLPBN_SRCDIR)/evidence.yap \
+	$(CLPBN_SRCDIR)/fove.yap \
 	$(CLPBN_SRCDIR)/gibbs.yap \
 	$(CLPBN_SRCDIR)/graphs.yap \
 	$(CLPBN_SRCDIR)/graphviz.yap \
+	$(CLPBN_SRCDIR)/ground_factors.yap \
 	$(CLPBN_SRCDIR)/hmm.yap \
 	$(CLPBN_SRCDIR)/jt.yap \
 	$(CLPBN_SRCDIR)/matrix_cpt_utils.yap \
@@ -66,6 +69,7 @@ CLPBN_SCHOOL_EXAMPLES= \
 	$(CLPBN_EXDIR)/School/README \
 	$(CLPBN_EXDIR)/School/evidence_128.yap \
 	$(CLPBN_EXDIR)/School/schema.yap \
+	$(CLPBN_EXDIR)/School/parschema.yap \
 	$(CLPBN_EXDIR)/School/school_128.yap \
 	$(CLPBN_EXDIR)/School/school_32.yap \
 	$(CLPBN_EXDIR)/School/school_64.yap \
diff --git a/packages/CLPBN/clpbn.yap b/packages/CLPBN/clpbn.yap
index 8dbcebc92..ddbc66003 100644
--- a/packages/CLPBN/clpbn.yap
+++ b/packages/CLPBN/clpbn.yap
@@ -46,6 +46,14 @@
 	       finalize_bp_solver/1
 	      ]).
 
+:- use_module('clpbn/fove',
+	      [fove/3,
+               check_if_fove_done/1,
+	       init_fove_solver/4,
+	       run_fove_solver/3,
+	       finalize_fove_solver/1
+	      ]).
+
 
 :- use_module('clpbn/jt',
 	      [jt/3,
@@ -65,14 +73,6 @@
 	       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,
@@ -96,7 +96,7 @@
 :- use_module('clpbn/evidence',
 	      [
 	       store_evidence/1,
-	       add_evidence/2,
+	       add_stored_evidence/2,
 	       incorporate_evidence/2,
 	       check_stored_evidence/2,
 	       put_evidence/2
@@ -110,7 +110,11 @@
 :- use_module('clpbn/graphviz',
 	[clpbn2gviz/4]).
 
-:- dynamic solver/1,output/1,use/1,suppress_attribute_display/1, parameter_softening/1, em_solver/1.
+:- use_module(clpbn/ground_factors,
+	[generate_bn/2]).
+
+
+:- dynamic solver/1,output/1,use/1,suppress_attribute_display/1, parameter_softening/1, em_solver/1, use_parfactors/1.
 
 solver(ve).
 em_solver(ve).
@@ -122,6 +126,7 @@ em_solver(ve).
 output(no).
 suppress_attribute_display(false).
 parameter_softening(m_estimate(10)).
+use_parfactors(off).
 
 clpbn_flag(Flag,Option) :-
 	clpbn_flag(Flag, Option, Option).
@@ -153,7 +158,9 @@ clpbn_flag(suppress_attribute_display,Before,After) :-
 clpbn_flag(parameter_softening,Before,After) :-
 	retract(parameter_softening(Before)),
 	assert(parameter_softening(After)).
-
+clpbn_flag(use_factors,Before,After) :-
+	retract(use_parfactors(Before)),
+	assert(use_parfactors(After)).
 
 {_} :-
 	solver(none), !.
@@ -204,8 +211,8 @@ add_evidence(V,Key,Distinfo,NV) :-
 	store_var(NV),
 	clpbn:put_atts(NV,evidence(Pos)).
 add_evidence(V,K,_,V) :-
-	store_var(V),
-	add_evidence(K,V).
+	add_stored_evidence(K,V),
+	store_var(V).
 
 clpbn_marginalise(V, Dist) :-
 	attributes:all_attvars(AVars),
@@ -216,8 +223,9 @@ clpbn_marginalise(V, Dist) :-
 % called by top-level
 % or by call_residue/2
 %
-project_attributes(GVars, AVars) :-
+project_attributes(GVars, AVars0) :-
 	suppress_attribute_display(false),
+	generate_vars(GVars, AVars0, AVars),
 	AVars = [_|_],
 	solver(Solver),
 	( GVars = [_|_] ; Solver = graphs), !,
@@ -235,10 +243,21 @@ project_attributes(GVars, AVars) :-
 	).
 project_attributes(_, _).
 
+generate_vars(GVars, _, NewAVars) :-
+	use_parfactors(on), !,
+	generate_bn(GVars, NewAVars).
+generate_vars(_GVars, AVars, AVars).
+
 clpbn_vars(AVars, DiffVars, AllVars) :-
 	sort_vars_by_key(AVars,SortedAVars,DiffVars),
 	incorporate_evidence(SortedAVars, AllVars).
 
+get_clpbn_vars([V|GVars],[V|CLPBNGVars]) :-
+	get_atts(V, [key(_)]), !,
+	get_clpbn_vars(GVars,CLPBNGVars).
+get_clpbn_vars([_|GVars],CLPBNGVars) :-
+	get_clpbn_vars(GVars,CLPBNGVars).
+
 get_clpbn_vars([],[]).
 get_clpbn_vars([V|GVars],[V|CLPBNGVars]) :-
 	get_atts(V, [key(_)]), !,
@@ -276,6 +295,8 @@ write_out(gibbs, GVars, AVars, DiffVars) :-
 	gibbs(GVars, AVars, DiffVars).
 write_out(bnt, GVars, AVars, DiffVars) :-
 	do_bnt(GVars, AVars, DiffVars).
+write_out(fove, GVars, AVars, DiffVars) :-
+	fove(GVars, AVars, DiffVars).
 
 get_bnode(Var, Goal) :-
 	get_atts(Var, [key(Key),dist(Dist,Parents)]),
diff --git a/packages/CLPBN/clpbn/bp.yap b/packages/CLPBN/clpbn/bp.yap
index 2c9ab9bc9..ea02fbc99 100644
--- a/packages/CLPBN/clpbn/bp.yap
+++ b/packages/CLPBN/clpbn/bp.yap
@@ -7,7 +7,6 @@
 
 :- module(clpbn_bp,
           [bp/3,
-           check_if_bp_done/1,
            set_solver_parameter/2,
            use_log_space/0,
            init_bp_solver/4,
@@ -72,7 +71,7 @@ bp([QueryVars], AllVars, Output) :-
 	clpbn_bind_vals([QueryVars], LPs, Output).
 
 
-init_bp_solver(_, AllVars0, _, bp(BayesNet, DistIds, AllParFactors)) :-
+init_bp_solver(_, AllVars0, _, bp(BayesNet, DistIds, _AllParFactors)) :-
 	check_for_agg_vars(AllVars0, AllVars),
 	%inc_network_counting,
 %writeln_clpbn_vars(AllVars),
@@ -100,126 +99,6 @@ parents_to_keys(Var.Parents, Key.Keys) :-
 	clpbn:get_atts(Var, [key(Key)]),
 	parents_to_keys(Parents, Keys).
 
-generate_parfactors(AllVars, ParFactors) :-
-	generate_factors(AllVars, Factors),
-%writeln(Factors),
-	% sort factors by distribution
-%	sort(Factors, DistFactors),
-%writeln(DistFactors),
-	group(DistFactors, ParFactors).
-%writeln(ParFactors).
-
-generate_factors(Var.AllVars, f(Dist,[Var|Parents]).AllFactors) :-
-	clpbn:get_atts(Var, [dist(Dist,Parents)]),
-	generate_factors(AllVars, AllFactors).
-generate_factors([], []).
-
-group([], []).
-group(f(Dist,Vs).DistFactors, phi(Dist,NConstraints,Domain).ParFactors) :-
-	number(Dist),
-	grab_similar_factors(Dist, Vs, f(Dist,Vs).DistFactors, RemainingDistFactors, Constraints),
-	simplify_constraints(Constraints, NConstraints, Domain), !,
-	group(RemainingDistFactors, ParFactors).
-
-group(f(Dist,Vs).DistFactors, phi(Dist,Constraints,[]).ParFactors) :-
-	grab_similar_factors(Dist, Vs, f(Dist,Vs).DistFactors, RemainingDistFactors, Constraints),
-	group(RemainingDistFactors, ParFactors).
-
-simplify_constraints([[1=El]|Constraints], [NEl], [in(1,NEl,Domain)]) :-
-	functor(El,Name,1), !,
-	functor(NEl,Name,1),
-	constraints_to_domain(1,[[1=El]|Constraints],Domain0),
-	sort(Domain0, Domain).
-simplify_constraints(Constraints, NewConstraints, Ds) :-
-        Constraints = [Constraint|_],
-	generate_domains(Constraint, Constraints, Ds), !,
-	normalize_constraints(Ds, Constraints, NewConstraints).
-simplify_constraints(Constraints, Constraints, []).
-
-normalize_constraints(Ds, Constraints, [T|GeneralizedConstraints]) :-
-	unique(Ds, I, T, RemDs), !,
-	remove_i(Constraints, I, ConstraintsI),
-	normalize_constraints(RemDs, ConstraintsI, GeneralizedConstraints).
-normalize_constraints(Ds, Constraints, [(S1,S2)|GeneralizedConstraints]) :-
-	equal(Ds, I, J, RemDs, S1, S2),
-	arg(1,S1,V),
-	arg(1,S2,V),
-%writeln(start:Ds:I:J),
-	remove_eqs(Constraints, I, J, ConstraintsI), !,
-	normalize_constraints(RemDs, ConstraintsI, GeneralizedConstraints).
-normalize_constraints(_Ds, Constraints, []) :-
-	Constraints = [[_]|_], !.
-normalize_constraints(_, Constraints, Constraints).
-
-unique([in(I,T,[_])|Ds], I, T, Ds).	
-unique([D|Ds], I, T, D.NewDs) :-
-	unique(Ds, I, T, NewDs).	
-
-equal([in(I,S1,Vals)|Ds], I, J, Ds, S1, S2) :-
-	equal2(Ds, Vals, J, S2), !.
-equal([D|Ds], I, J, D.NewDs, S1, S2) :-
-	equal(Ds, I, J, NewDs, S1, S2).	
-
-equal2([in(J,S2,Vals)|Ds], Vals, J, S2).
-equal2([D|Ds], Vals, J, S2) :-
-	equal2(Ds, Vals, J, S2).	
-
-remove_i([], _I, []).
-remove_i(C.Constraints, I, NewC.ConstraintsI) :-
-	remove_ic(C,I,NewC),
-	remove_i(Constraints, I, ConstraintsI).
-
-remove_ic([I=_|C], I, C) :- !.
-remove_ic(El.C, I, El.NewC) :-
-	remove_ic(C, I, NewC).
-
-remove_eqs([], _I, _J, []).
-remove_eqs(C.Constraints, I, J, NewC.ConstraintsI) :-
-	remove_eqs2(C, I, J, NewC),
-	remove_eqs(Constraints, I, J, ConstraintsI).
-
-remove_eqs2([I=V|C], I, J, C) :- !,
-	arg(1,V,A),
-	check_match(C, J, A).
-remove_eqs2(El.C, I, J, El.NewC) :-
-	remove_eqs2(C, I, J, NewC).
-
-check_match([J=V1|C], J, V) :- !,
-	arg(1,V1,V).
-check_match(El.C, J, V) :-
-	check_match(C, J, V).
-
-
-generate_domains([], _Constraints, []).
-generate_domains([I=El|Constraint], Constraints, in(I,NEl,Domain).Ds) :-
-	functor(El,Name,1), !,
-	functor(NEl,Name,1),
-	constraints_to_domain(I,Constraints,Domain0),
-	sort(Domain0, Domain),
-	generate_domains(Constraint, Constraints, Ds).	
-
-
-constraints_to_domain(_,[],[]).
-constraints_to_domain(I,[Constraint|Constraints],El.Domain) :-
-	add_constraint_to_domain(I, Constraint, El),
-	constraints_to_domain(I,Constraints,Domain).
-
-add_constraint_to_domain(I, [I=El|_], A) :- !,
-	arg(1, El, A).
-add_constraint_to_domain(I, _.Constraint, El) :-
-	add_constraint_to_domain(I, Constraint, El).
-
-
-grab_similar_factors(Dist, Vs, f(Dist,DVs).DistFactors, RemainingDistFactors, Constraint.Constraints) :-
-	grab_similar_factor(DVs, 1, Constraint), !,
-	grab_similar_factors(Dist, Vs, DistFactors, RemainingDistFactors, Constraints).
-grab_similar_factors(_Dist, _Vs, DistFactors, DistFactors, []).
-
-grab_similar_factor([], _Arg, []).
-grab_similar_factor(V.VDVs, Arg, (Arg=Key).Constraint) :-
-	clpbn:get_atts(V,key(Key)),
-	Arg1 is Arg+1,
-	grab_similar_factor(VDVs, Arg1, Constraint).
 
 
 process_ids([], _, []).
diff --git a/packages/CLPBN/clpbn/dists.yap b/packages/CLPBN/clpbn/dists.yap
index 362837947..ebd2454a6 100644
--- a/packages/CLPBN/clpbn/dists.yap
+++ b/packages/CLPBN/clpbn/dists.yap
@@ -179,6 +179,10 @@ add_dist(Domain, Type, CPT, Parents, Key, Id) :-
 
 record_parent_sizes([], Id, [], DSizes) :-
 	recordz(clpbn_dist_psizes,db(Id, DSizes),_).
+record_parent_sizes([P|Parents], Id, [Size|Sizes], DSizes) :-
+	integer(P), !,
+	Size = P,
+	record_parent_sizes(Parents, Id, Sizes, DSizes).
 record_parent_sizes([P|Parents], Id, [Size|Sizes], DSizes) :-
 	clpbn:get_atts(P,dist(Dist,_)), !,
 	get_dist_domain_size(Dist, Size),
diff --git a/packages/CLPBN/clpbn/evidence.yap b/packages/CLPBN/clpbn/evidence.yap
index 68fadc7e4..62b221858 100644
--- a/packages/CLPBN/clpbn/evidence.yap
+++ b/packages/CLPBN/clpbn/evidence.yap
@@ -8,7 +8,7 @@
 	   store_evidence/1,
 	   incorporate_evidence/2,
 	   check_stored_evidence/2,
-	   add_evidence/2,
+	   add_stored_evidence/2,
 	   put_evidence/2
 	  ]).
 
@@ -19,7 +19,7 @@
     ]).
 
 :- use_module(library('clpbn/dists'), [
-       get_evidence_position/3
+	get_dist/4
     ]).
 
 :- use_module(library(rbtrees), [
@@ -30,7 +30,7 @@
 
 :- meta_predicate store_evidence(:).
 
-:- dynamic node/4, edge/2, evidence/2.
+:- dynamic node/3, edge/2, evidence/2.
 
 %
 % new evidence storage algorithm. The idea is that instead of 
@@ -46,28 +46,42 @@ store_evidence(G) :-
 	compute_evidence(G, PreviousSolver).
 
 compute_evidence(G, PreviousSolver) :-
-	catch(call_residue(G, Vars), Ball, evidence_error(Ball,PreviousSolver)), !,
-	store_graph(Vars),
+	catch(get_clpbn_vars(G, Vars), Ball, evidence_error(Ball,PreviousSolver)), !,
+	store_graph(Vars), !,
 	set_clpbn_flag(solver, PreviousSolver).
 compute_evidence(_,PreviousSolver) :-
 	set_clpbn_flag(solver, PreviousSolver).
 
+get_clpbn_vars(G, Vars) :-
+%	attributes:all_attvars(Vars0),	
+	once(G),
+	attributes:all_attvars(Vars).	
 
 evidence_error(Ball,PreviousSolver) :-
 	set_clpbn_flag(solver,PreviousSolver),
 	throw(Ball).
 
 store_graph([]).
-store_graph([_-node(K,Dom,CPT,TVs,Ev)|Vars]) :-
-	\+ node(K,_,_,_), !,
-	assert(node(K,Dom,CPT,TVs)),
-	( nonvar(Ev) -> assert(evidence(K,Ev)) ; true),
+store_graph([V|Vars]) :-
+	clpbn:get_atts(V,[key(K),dist(Id,Vs)]),        
+	\+ node(K, Id, _), !,
+	translate_vars(Vs,TVs),
+	assert(node(K,Id,TVs)),
+	( clpbn:get_atts(V,[evidence(Ev)]) -> assert(evidence(K,Ev)) ; true),
 	add_links(TVs,K),
 	store_graph(Vars).
 store_graph([_|Vars]) :-
 	store_graph(Vars).
 
+translate_vars([],[]).
+translate_vars([V|Vs],[K|Ks]) :-
+	clpbn:get_atts(V, [key(K)]),
+	translate_vars(Vs,Ks).
+
 add_links([],_).
+add_links([K0|TVs],K) :-
+	edge(K,K0), !,
+	add_links(TVs,K).
 add_links([K0|TVs],K) :-
 	assert(edge(K,K0)),
 	add_links(TVs,K).
@@ -82,7 +96,7 @@ incorporate_evidence(Vs,AllVs) :-
 create_open_list([], L, L, C, C).
 create_open_list([V|Vs], [K-V|OL], FL, C0, CF) :-
 	clpbn:get_atts(V,[key(K)]),
-	add_evidence(K, V),
+	add_stored_evidence(K, V),
 	rb_insert(C0, K, V, CI),
 	create_open_list(Vs, OL, FL, CI, CF).
 
@@ -91,12 +105,26 @@ do_variables([K-V|Vs], Vf, C0) :-
 	check_for_evidence(K, V, Vf, Vff, C0, Ci),
 	do_variables(Vs, Vff, Ci).
 
-create_new_variable(K, V, Vf0, Vff, C0, Cf) :-
-	node(K,Dom, CPT, TVs),
-	{ V = K with p(Dom, CPT, NTVs) },
-	add_evidence(K, V),
-	add_variables(TVs, NTVs, Vf0, Vff, C0, Cf).
+extract_vars([], []).
+extract_vars([_-V|Cache], [V|AllVs]) :-
+	extract_vars(Cache, AllVs).
 
+%make sure that we are consistent
+check_stored_evidence(K, Ev) :-
+	evidence(K, Ev0), !,
+	Ev0 = Ev.
+check_stored_evidence(_, _).
+
+add_stored_evidence(K, V) :-
+	evidence(K, Ev), !,
+	put_evidence(Ev, V).
+add_stored_evidence(_, _).
+
+check_for_evidence(_, V, Vf, Vf, C, C) :-
+	clpbn:get_atts(V, [evidence(_)]), !.
+check_for_evidence(K, _, Vf0, Vff, C0, Ci) :-
+	findall(Rt,edge(Rt,K),Rts),
+	add_variables(Rts, _, Vf0, Vff, C0, Ci).
 
 add_variables([], [], Vf, Vf, C, C).
 add_variables([K|TVs], [V|NTVs], Vf0, Vff, C0, Cf) :-
@@ -107,30 +135,22 @@ add_variables([K|TVs], [V|NTVs], [K-V|Vf0], Vff, C0, Cf) :-
 	create_new_variable(K, V, Vf0, Vf1, C1, C2),
 	add_variables(TVs, NTVs, Vf1, Vff, C2, Cf).
 
+create_new_variable(K, V, Vf0, Vff, C0, Cf) :-
+	node(K, Id, TVs),
+writeln(add:K:Id),
+	get_dist(Id,_,Dom,CPT), !,
+	{ V = K with p(Dom, CPT, NTVs) },
+	add_stored_evidence(K, V),
+	add_variables(TVs, NTVs, Vf0, Vff, C0, Cf).
+create_new_variable(K, V, Vf0, Vff, C0, Cf) :-
+	node(K, Id, TVs),
+	Id =.. [Na,Dom],
+	Dist =.. [Na,Dom,NTVs],
+	{ V = K with Dist },
+writeln(done),
+	add_stored_evidence(K, V),
+	add_variables(TVs, NTVs, Vf0, Vff, C0, Cf).
 
-extract_vars([], []).
-extract_vars([_-V|Cache], [V|AllVs]) :-
-	extract_vars(Cache, AllVs).
-
-%make sure that we are 
-check_stored_evidence(K, Ev) :-
-	evidence(K, Ev0), !, Ev0 = Ev.
-check_stored_evidence(_, _).
-
-add_evidence(K, V) :-
-	evidence(K, Ev), !,
-	store_evidence(V, Ev),
-	clpbn:put_atts(V, [evidence(Ev)]).
-add_evidence(_, _).
-
-check_for_evidence(_, V, Vf, Vf, C, C) :-
-	clpbn:get_atts(V, [evidence(_)]), !.
-check_for_evidence(K, _, Vf0, Vff, C0, Ci) :-
-	findall(Rt,edge(Rt,K),Rts),
-	add_variables(Rts, _, Vf0, Vff, C0, Ci).
-
-put_evidence(K, V) :-
-	clpbn:get_atts(V, [dist(Id,_)]),
-	get_evidence_position(K, Id, Ev),
+put_evidence(Ev, V) :-
 	clpbn:put_atts(V, [evidence(Ev)]).
 
diff --git a/packages/CLPBN/clpbn/fove.yap b/packages/CLPBN/clpbn/fove.yap
new file mode 100644
index 000000000..6bc8c94b7
--- /dev/null
+++ b/packages/CLPBN/clpbn/fove.yap
@@ -0,0 +1,94 @@
+:- module(clpbn_fove,
+          [fove/3,
+           set_solver_parameter/2,
+           init_fove_solver/4,
+           run_fove_solver/3,
+           finalize_fove_solver/1
+          ]).
+
+:- use_module(library(pfl), [
+          factor/5,
+	  skolem/2]).
+
+%
+% support fove method
+%
+
+fove([[]],_,_) :- !.
+fove([QueryVars], AllVars, Output) :-
+	init_fove_solver(_, AllVars, _, GraphicalNet),
+	run_fove_solver([QueryVars], LPs, GraphicalNet),
+	finalize_fove_solver(GraphicalNet),
+	clpbn_bind_vals([QueryVars], LPs, Output).
+
+%
+% set up network, add evidence, and query all marginals at the same time?
+%
+init_fove_solver(_, AllAttVars, _, fove(ParNet, EvidenceVariables)) :-
+	all_factors(Factors),
+	all_domains(Domains),
+	evidence_variables(AllAttVars, EvidenceVariables),
+writeln(ev:EvidenceVariables),
+	% c-code, just receives the par factors
+	init_fove(Factors, Domains, ParNet).	
+
+evidence_variables([], []).
+evidence_variables(V.AllAttVars, [K:E|EvidenceVariables]) :-
+	clpbn:get_atts(V,[key(K),evidence(E)]), !,
+	evidence_variables(AllAttVars, EvidenceVariables).
+evidence_variables(_V.AllAttVars, EvidenceVariables) :-
+	evidence_variables(AllAttVars, EvidenceVariables).
+
+all_domains(Domains) :-
+	findall(X:Y, skolem(X,Y), Domains).
+
+
+
+:- table all_factors/1.
+
+%
+% enumerate all par-factors and enumerate their domain as tuples.
+%
+% output is list of pf(
+%        ID: an unique number
+%        Ks: a list of keys, also known as the pf formula [a(X),b(Y),c(X,Y)]
+%        Vs: the list of free variables [X,Y]
+%        Phi: the table following usual CLP(BN) convention
+%        Tuples: tuples with all ground bindings for variables in Vs, of the form [fv(x,y)]
+%
+all_factors(Factors) :-
+	findall(F, is_factor(F), Factors).
+
+is_factor(pf(Id, Ks, Vs, Phi, Tuples)) :-
+	factor(Id, Ks, Vs, Table, Constraints),
+	Table \= avg,
+	gen_table(Table, Phi),
+	all_tuples(Constraints, Vs, Tuples).
+
+gen_table(Table, Phi) :-
+	( is_list(Table)
+          -> 
+	  Phi = Table
+          ;
+	  call(user:Table, Phi)
+        ).
+
+all_tuples(Constraints, Tuple, Tuples) :-
+	setof(Tuple, Constraints^run(Constraints), Tuples).
+
+run([]).
+run(Goal.Constraints) :-
+	user:Goal,
+	run(Constraints).
+
+%
+% ask probability of a single variable
+%
+run_fove_solver(QueryVars, LPs, fove(ParFactors, EvidenceVariables)) :-
+	fove(QueryVars, EvidenceVariables, ParFactors, LPs).	
+
+
+
+	
+
+
diff --git a/packages/CLPBN/clpbn/ground_factors.yap b/packages/CLPBN/clpbn/ground_factors.yap
new file mode 100644
index 000000000..966b74681
--- /dev/null
+++ b/packages/CLPBN/clpbn/ground_factors.yap
@@ -0,0 +1,232 @@
+
+%parfactor(
+%	[ability(P),grade(C,S), satisfaction(C,S,P)],
+%        \phi = [....], 
+%       [P,C,S],
+%	[P \in [p1,p2,p4], C \in [c1,c3], S \in [s2,s3]]).
+%	[S \= s2])
+
+
+:- module(clpbn_ground_factors, [
+          generate_bn/2,
+          ground_parfactors/1]).
+
+:- use_module(library(bhash), [
+          b_hash_new/1,
+          b_hash_lookup/3,
+	  b_hash_insert/4]).
+
+:- use_module(library(pfl), [
+          factor/5,
+	  skolem/2]).
+
+:- use_module(library(clpbn/dists), [
+          dist/4]).
+
+
+%
+% generate a CLP(BN) network that can be run in CLP(BN).
+%
+generate_bn(QueryVars, AllAttVars) :-
+	attributes:all_attvars(AVars),
+	b_hash_new(H0),
+	check_for_evidence(AVars, EVars),
+	run_through_factors(QueryVars, H0, H1, AllAttVars, IVars),
+	run_through_factors(EVars, H1, _HF, IVars, []).
+
+check_for_evidence(V.AVars, V.EVars) :-
+	clpbn:get_atts(V,[evidence(_E)]), !,
+	check_for_evidence(AVars, EVars).
+check_for_evidence(_V.AVars, EVars) :-
+	check_for_evidence(AVars, EVars).
+check_for_evidence([], []).
+
+run_through_factors([], H, H) --> [].
+run_through_factors(V.Vars, H0, HF) -->
+	{ clpbn:get_atts(V,[key(K)]), 
+	  b_hash_lookup(K,V,H0)
+        }, !,
+	run_through_factors(Vars, H0, HF).
+run_through_factors(V.Vars, H0, HF) -->
+	% it is a new clpbn variable 
+	[V],
+	{
+          % should already have a key
+          clpbn:get_atts(V,[key(K)]),
+	  % insert it into a table of seen variables
+	  b_hash_insert(H0,K,V,HI),
+	  construct_clpbn_node(K, V, HI, MoreVars, Vars)
+        },
+	run_through_factors(MoreVars, HI, HF).
+
+% aggregates are special.
+construct_clpbn_node(K, V, HI) -->
+	% and get the corresponding factor
+	{ factor(Id, [K|Ks], _, avg, Constraints) }, !,
+	{
+          skolem(K, Domain),
+	  dist(avg(Domain, Parents), DistId, K, Parents),
+          clpbn:put_atts(V,[dist(DistId,Parents)]), 
+	  % we also need to set the constraints
+	  % this should set all the keys to rights
+	  run(Constraints)
+        },
+	% now let's look up parents and set up the graph
+	run_bayesian_factor(Ks, HI, Parents, []).
+construct_clpbn_node(K, V, HI) -->
+	{
+	  % and get the corresponding factor
+	  factor(Id, [K|Ks], _, _Phi, Constraints),
+	  factor_to_dist(Id, DistId),
+	  % and the dist constraint
+          clpbn:put_atts(V,[dist(DistId,Parents)]), 
+	  % we also need to set the constraints
+	  % this should set all the keys to rights
+	  run(Constraints)
+        },
+	% now let's look up parents and set up the graph
+	run_bayesian_factor(Ks, HI, Parents, []).
+
+factor_to_dist(Id, NewId) :-
+	factor(Id, [K|Ks], _, Phi, _Constraints),
+	skolem(K, Domain),
+	( is_list(Phi)
+          -> 
+	  CPT = Phi
+          ;
+	  call(user:Phi, CPT)
+        ),
+	keys_to_sizes(Ks, Szs),
+	dist(p(Domain, CPT, Szs), NewId, K, Szs).
+	
+keys_to_sizes([], []).
+keys_to_sizes(K.Ks, Sz.Szs) :-
+	skolem(K, Domain),
+	length(Domain, Sz),
+	keys_to_sizes(Ks, Szs).
+
+run([]).
+run(Goal.Constraints) :-
+	user:Goal, !,
+	run(Constraints).
+
+run_bayesian_factor([], _H, Vs, Vs) --> [].
+run_bayesian_factor(K.Ks, H, Vs, Vs0) -->
+	run_var(K, H, Vs, Vs1),
+	run_bayesian_factor(Ks, H, Vs1, Vs0).
+
+%
+% this function returns a list of *new* variables
+%
+% collection of random variables
+run_var(avg(Els), H, Vs, Vs0) --> !,
+	run_vars(Els, H, Vs, Vs0).
+% existing random variable
+run_var(K, H, V.Vs, Vs) -->
+	{  b_hash_lookup(K,V,H) }, !.
+% new random variable
+run_var(K, _H, V.Vs, Vs) -->
+	[V],
+	{
+          clpbn:put_atts(V,[key(K)])
+        }.
+
+run_vars([], _H, Vs, Vs) --> [].
+run_vars(K.Els, H, Vs, Vs0) -->
+	run_var(K, H, Vs, VsI),
+	run_vars(Els, H, VsI, Vs0).
+
+ground_parfactors(ParFactors) :-
+	findall(Factor, factor(Factor), SourceFactors),
+	run_all_parfactors(SourceFactors, ParFactors).
+
+factor(Factor) :- 
+	user:parfactor(Factor).
+factor(Factor) :- 
+	user:bayes(Factor).
+
+run_all_parfactors([], []).
+run_all_parfactors(Source.SourceFactors, Factor.ParFactors) :-
+	run_parfactors(Source, Factor),
+	run_all_parfactors(SourceFactors, ParFactors).
+
+run_parfactors((Formula ; Phi ; ConstraintGenerator), parfactor(Formula, Phi, FV, Domain, NewConstraints)) :-
+	term_variables(Formula, FreeVars),
+	FV =.. fv(FreeVars),
+	evaluate_constraints(FV, ConstraintGenerator, NewConstraints, Domain).
+
+evaluate_constraints(FreeVars, Constraint.ConstraintGenerators, NC, Domain) :-
+	functor(FreeVars, fv, NOf),
+	setof(FreeVars, user:Constraint, Bindings),
+	run_free_vars(0, NOf, FreeVars, Bindings, Domain, Domain0),
+	get_list_of_conditions(Domain, 0, N, Conditions),
+	add_extra_constraints(N, Conditions, Bindings, NC, NC0),
+	evaluate_constraints(FreeVars, ConstraintGenerators, NC0, Domain0).
+evaluate_constraints(_FreeVars, [], []).
+
+run_free_vars(N, N, _FreeVars, _Bindings) --> !.
+run_free_vars(I0, N, FreeVars, Bindings) -->
+	{ I is I0+1,
+	 arg(I, FreeVars, V),
+	 Bindings = B._,
+	 arg(I, B, C), ground(C)
+        }, !,
+	{ setof(C, check_val(Bindings, I, C), Dom) },
+	[domain(I,V,Dom)],
+	run_free_vars(I, N, FreeVars, Bindings).
+run_free_vars(I0, N, FreeVars, Bindings) -->
+	I is I0+1,
+	run_free_vars(I, N, FreeVars, Bindings).
+
+add_extra_constraints(0, [], _Bindings) --> !.
+add_extra_constraints(1, _Conditions, _Bindings) --> !.
+add_extra_constraints(N, Conditions, Bindings) -->
+	{ extract_unique(Conditions, NewConditions) }, !,
+	{ N1 is N-1 },
+	add_extra_constraints(N1, NewConditions, Bindings).
+add_extra_constraints(N, [dom(I1,V1,Dom1),dom(I2,V2,Dom2)|Conditions], Bindings) -->
+	{ length(Dom1, Sz), length(Dom2, Sz) }, !,
+	{ N1 is N-2 },
+	{ generate_map(Bindings, I1, I2, Mapping) },
+	[map(V1,V2,Mapping)],
+	add_extra_constraints(N1, dom(I1,V1,Dom1).Conditions, Bindings).
+add_extra_constraints(_N, Conditions, Bindings) -->
+	[or(Vs,Or)],
+	{ gather_vs(Conditions, Vs, Indices),
+	  generate(Bindings, Indices, Or) }.
+
+% domain is a singleton constant
+extract_unique(domain(_,_,[_]).Conditions, Conditions) :- !.
+extract_unique(_.Conditions, NewConditions) :-
+	extract_unique(Conditions, NewConditions).
+
+get_list_of_conditions([], N, N, []).
+get_list_of_conditions(Dom._, N, N, _Conditions) :- 
+	var(Dom), !.
+get_list_of_conditions(Dom.Domain, I0, N, Dom.Conditions) :-
+	I is I0+1,
+	get_list_of_conditions(Domain, I, N, Conditions).
+
+check_val(B._Bindings, I, C) :-
+	arg(I, B, C).
+check_val(_.Bindings, I, C) :-
+	check_val(Bindings, I, C).
+	
+generate_map(B.Bindings, I1, I2, [[A1|A2]|Mapping]) :-
+	arg(I1, B, A1),
+	arg(I2, B, A2),
+	generate_map(Bindings, I1, I2, Mapping).
+
+gather_vs([], [], []).
+gather_vs(domain(I,V,_).Conditions, V.Vs, I.Indices) :-
+	gather_vs(Conditions, Vs, Indices).
+
+generate([], _, []).
+generate(B.Bindings, Indices, O.Or) :-
+	generate_el(B, Indices, O),
+	generate(Bindings, Indices, Or).
+
+generate_el(_B, [], []).
+generate_el(B, I.Indices, A.O) :-
+	arg(I, B, A),
+	generate_el(B, Indices, O).
diff --git a/packages/CLPBN/clpbn/matrix_cpt_utils.yap b/packages/CLPBN/clpbn/matrix_cpt_utils.yap
index 053ec4595..26d217416 100644
--- a/packages/CLPBN/clpbn/matrix_cpt_utils.yap
+++ b/packages/CLPBN/clpbn/matrix_cpt_utils.yap
@@ -36,6 +36,7 @@
 	       matrix_dims/2,
 	       matrix_sum/2,
 	       matrix_sum_logs_out/3,
+	       matrix_sum_out/3,
 	       matrix_sum_logs_out_several/3,
 	       matrix_op_to_all/4,
 	       matrix_to_exps2/1,
@@ -61,7 +62,9 @@ project_from_CPT(V,tab(Table,Deps,_),tab(NewTable,NDeps,NSzs)) :-
 	matrix_dims(NewTable, NSzs).
 project_from_CPT(V,tab(Table,Deps,_),tab(NewTable,NDeps,NSzs)) :-
 	vnth(Deps, 0, V, N, NDeps),
+%	matrix_to_exps2(Table),
 	matrix_sum_logs_out(Table, N, NewTable),
+%	matrix_to_logs(NewTable),
 	matrix_dims(NewTable, NSzs).
 
 evidence(V, Pos) :-
diff --git a/packages/CLPBN/examples/School/evidence_128.yap b/packages/CLPBN/examples/School/evidence_128.yap
index d3d72799f..c7402b510 100644
--- a/packages/CLPBN/examples/School/evidence_128.yap
+++ b/packages/CLPBN/examples/School/evidence_128.yap
@@ -1,5 +1,5 @@
 
-:- [school_128].
+:- ensure_loaded(school_128).
 
 professor_popularity(p0,h) :- {}.
 professor_popularity(p3,h) :- {}.
diff --git a/packages/CLPBN/examples/School/parschema.yap b/packages/CLPBN/examples/School/parschema.yap
new file mode 100644
index 000000000..8e7834a0d
--- /dev/null
+++ b/packages/CLPBN/examples/School/parschema.yap
@@ -0,0 +1,66 @@
+
+:- use_module(library(pfl)).
+
+/* base file for school database. Supposed to be called from school_*.yap */
+
+%
+% bayes is a parfactor for a bayesian network,
+% first argument is target of other arguments pop(K) <- abi(K)
+% second argument is the name of a predicate to call for \phi (CPT)
+% last argument is a list of goals defining the constraints over the elements
+% of the 
+%
+
+%
+% these states that skolem variables abi(K) are in a parametric factor with
+% with \phi defined by abi_table(X) and whose domain and constraints
+% is obtained from professor/1.
+%
+bayes abi(K)::[h,m,l] ; abi_table ; [professor(K)].
+
+bayes pop(K)::[h,m,l], abi(K) ; pop_table ; [professor(K)].
+
+bayes grade(C,S)::[a,b,c,d], int(S), diff(C) ; grade_table ; [registration(_,C,S)].
+
+bayes sat(C,S,P)::[h,m,l], abi(P), grade(C,S) ; sat_table ; [reg_sat(C,S,P)].
+
+bayes rat(C) :: [h,m,l], avg(Sats) ; avg ; [course_rating(C, Sats)].
+
+bayes diff(C) :: [h,m,l] ; diff_table ; [course(C,_)].
+
+bayes int(S) :: [h,m,l] ; int_table ; [student(S)].
+
+bayes rank(S) :: [a,b,c,d], avg(Grades) ; avg ; [student_ranking(S,Grades)].
+
+
+grade(Key, Grade) :-
+	registration(Key, CKey, SKey),
+	grade(CKey, SKey, Grade).
+
+reg_sat(CKey, SKey, PKey) :-
+	registration(_Key, CKey, SKey),
+	course(CKey, PKey).
+
+course_rating(CKey, Sats) :-
+	course(CKey,  _),
+	setof(sat(CKey,SKey,PKey),
+	   reg_sat(CKey, SKey, PKey),
+          Sats).
+
+student_ranking(SKey, Grades) :-
+	student(SKey),
+	setof(grade(CKey,SKey), RKey^registration(RKey,CKey,SKey), Grades).
+
+:- ensure_loaded(tables).
+
+%
+% evidence
+%
+abi(p0, h).
+
+pop(p1, m).
+pop(p2, h).
+
+% Query
+% ?- abi(p0, X).
+
diff --git a/packages/CLPBN/examples/School/school_128.yap b/packages/CLPBN/examples/School/school_128.yap
index 260be5302..034df01d9 100644
--- a/packages/CLPBN/examples/School/school_128.yap
+++ b/packages/CLPBN/examples/School/school_128.yap
@@ -18428,5 +18428,5 @@ registration(r13919,c221,s4095).
 registration(r13920,c39,s4095).
 
 
-% :- [evidence_128].
+ :- [evidence_128].
 
diff --git a/packages/CLPBN/examples/School/school_32.yap b/packages/CLPBN/examples/School/school_32.yap
index a948a2e64..1c436668c 100644
--- a/packages/CLPBN/examples/School/school_32.yap
+++ b/packages/CLPBN/examples/School/school_32.yap
@@ -15,9 +15,7 @@ total_students(256).
 
 :- yap_flag(write_strings,on).
 
-:- use_module(library(clpbn)).
-
-:- [-schema].
+:- [-parschema].
 
 professor(p0).
 professor(p1).
diff --git a/packages/CLPBN/examples/School/tables.yap b/packages/CLPBN/examples/School/tables.yap
index 21aff50e6..1bf6c14fa 100644
--- a/packages/CLPBN/examples/School/tables.yap
+++ b/packages/CLPBN/examples/School/tables.yap
@@ -1,27 +1,31 @@
 
-int_table(_, [0.5,
-	      0.4,
-	      0.1],[h, m, l]).
+int_table(_,T ,[h, m, l]) :- int_table(T).
+
+int_table([0.5,
+	0.4,
+	0.1]).
+
+	   /* h h   h m  h l  m h  m m  m l  l h  l m l l */
+grade_table([
+             0.2,  0.7, 0.85, 0.1, 0.2, 0.5, 0.01, 0.05,0.1 ,
+	     0.6, 0.25, 0.12, 0.3, 0.6,0.35,0.04, 0.15, 0.4 ,
+             0.15,0.04, 0.02, 0.4,0.15,0.12, 0.5, 0.6, 0.4,
+	     0.05,0.01, 0.01, 0.2,0.05,0.03, 0.45, 0.2, 0.1 ]).
 
 grade_table(I, D, 
 	   /* h h   h m  h l  m h  m m  m l  l h  l m l l */
-	    p([a,b,c,d],
-		[ 0.2,  0.7, 0.85, 0.1, 0.2, 0.5, 0.01, 0.05,0.1 ,
-	          0.6, 0.25, 0.12, 0.3, 0.6,0.35,0.04, 0.15, 0.4 ,
-                  0.15,0.04, 0.02, 0.4,0.15,0.12, 0.5, 0.6, 0.4,
-	          0.05,0.01, 0.01, 0.2,0.05,0.03, 0.45, 0.2, 0.1 ], [I,D])).
+	    p([a,b,c,d], T, [I,D])) :- grade_table(T).
     
-	
+sat_table(
+       /*   h a  h b  h c   h d  m a   m b  m c  m d   l a   l b  l c  l d */
+/*h*/     [0.98, 0.9, 0.8 , 0.6, 0.9,  0.4, 0.2, 0.01, 0.5,  0.2, 0.01, 0.01,
+/*m*/	   0.01, 0.09,0.15, 0.3, 0.05, 0.4, 0.3, 0.04, 0.35, 0.3, 0.09, 0.01 ,
+/*l*/	   0.01, 0.01,0.05, 0.1, 0.05, 0.2, 0.5, 0.95, 0.15, 0.5, 0.9,  0.98]).	
 /*
   A: professor's ability;
   B: student's grade (for course registration).
 */
-satisfaction_table(A, G,
-       /* h a  h b  h c  h d  m a  m b  m c  m d  l a  l b  l c  l d */
-	p([h,m,l],
-/*h*/     [0.98, 0.9,0.8 , 0.6, 0.9, 0.4, 0.2, 0.01, 0.5, 0.2, 0.01, 0.01,
-/*m*/	   0.01, 0.09,0.15, 0.3,0.05, 0.4, 0.3, 0.04,0.35, 0.3, 0.09, 0.01 ,
-/*l*/	   0.01, 0.01,0.05, 0.1,0.05, 0.2, 0.5, 0.95,0.15, 0.5, 0.9, 0.98], [A,G])).
+satisfaction_table(A, G, p([h,m,l], T, [A,G])) :- sat_table(T).
 
 
 % The idea is quite simple:
@@ -35,11 +39,18 @@ rating_prob_table([0.9,0.05,0.01,
 		   0.09,0.9,0.09,
 		   0.01,0.05,0.9]).
 
-abi_table( _, [0.50, 0.40,  0.10]).
+abi_table( [0.50, 0.40,  0.10]).
+
+abi_table( _, T) :- abi_table(T).
 
 
-pop_table(_, [0.9, 0.2, 0.01,
-              0.09, 0.6, 0.09,
-	      0.01, 0.2, 0.9]).
+pop_table( [0.9, 0.2, 0.01,
+           0.09, 0.6, 0.09,
+           0.01, 0.2, 0.9]).
+
+pop_table(_, T) :- pop_table(T). 
+
+diff_table([0.25, 0.50, 0.25]).
+
+dif_table(_, T) :- diff_table(T). 
 
-dif_table( _, [0.25, 0.50, 0.25]).
diff --git a/packages/CLPBN/pfl.yap b/packages/CLPBN/pfl.yap
new file mode 100644
index 000000000..43949abbe
--- /dev/null
+++ b/packages/CLPBN/pfl.yap
@@ -0,0 +1,92 @@
+%
+% This module defines PFL, the prolog factor language.
+%
+%
+
+:- module(clpbn_parfactors, [
+          factor/5,
+	  skolem/2,
+	  op(550,yfx,::),
+	  op(1150,fx,bayes),
+	  op(1150,fx,parfactor)]).
+
+:- use_module(library(lists),
+	[nth0/3,
+	 append/3]).
+
+:- dynamic factor/5, skolem_in/2, skolem/2, preprocess/3, evidence/2, id/1.
+
+:- reexport(library(clpbn),
+	[clpbn_flag/3,
+	 set_clpbn_flag/2]).
+
+:- set_clpbn_flag(use_factors,on).
+
+user:term_expansion( bayes((Formula ; Phi ; Constraints)), clpbn_parfactors:factor(Id,FList,FV,Phi,Constraints)) :-
+	!,
+	term_variables(Formula, FreeVars),
+	FV =.. [fv|FreeVars],
+	new_id(Id),
+	process_args(Formula, Id, 0, _, FList, []).
+user:term_expansion( Goal, [] ) :-
+	preprocess(Goal, Sk,Var), !,
+	(ground(Goal) -> true ; throw(error('non ground evidence',Goal))),
+%	prolog_load_context(module, M),
+	assert(evidence(Sk,Var)).
+
+id(0).
+
+new_id(Id) :-
+	retract(id(Id0)),
+	Id is Id0+1,
+	assert(id(Id)).
+
+process_args((Arg1,Arg2), Id, I0, I ) --> !,
+	process_args(Arg1, Id, I0, I1),
+	process_args(Arg2, Id, I1, I).
+process_args(Arg1, Id, I0, I ) -->
+	{ I is I0+1 },
+	process_arg(Arg1, Id, I).
+
+process_arg(Sk::D, Id, _I) -->
+	!,
+	{
+         new_skolem(Sk,D),
+	 assert(skolem_in(Sk, Id))
+        },
+	[Sk].
+process_arg(Sk, Id, _I) -->
+	!,
+	{
+	 assert(skolem_in(Sk, Id))
+        },
+	[Sk].
+
+new_skolem(Sk,D) :-
+	copy_term(Sk, Sk1),
+	skolem(Sk1, D1),
+	Sk1 =@= Sk,
+	!,
+	D1 = D.
+new_skolem(Sk,D) :-
+	interface_predicate(Sk),
+	assert(skolem(Sk, D)).
+
+interface_predicate(Sk) :-
+	Sk =.. SKAs,
+	append(SKAs, [Var], ESKAs),
+	ESk =.. ESKAs,
+	assert(preprocess(ESk, Sk, Var)),
+	assert_static((user:ESk :-
+	     var(Var) -> insert_atts(Var,Sk) ; add_evidence(Sk,Var) )
+	).
+
+insert_atts(Var,Sk) :-
+	clpbn:put_atts(Var,[key(Sk)]).
+
+add_evidence(Sk,Var) :-
+	skolem(Sk,D),
+	once(nth0(E,D,Var)),
+	clpbn:put_atts(_V,[key(Sk),evidence(E)]).
+
+