wolffd@0: function engine = jtree_inf_engine(bnet, varargin)
wolffd@0: % JTREE_INF_ENGINE Junction tree inference engine
wolffd@0: % engine = jtree_inf_engine(bnet, ...)
wolffd@0: %
wolffd@0: % The following optional arguments can be specified in the form of name/value pairs:
wolffd@0: % [default value in brackets]
wolffd@0: %
wolffd@0: % clusters  - a cell array of sets of nodes we want to ensure are in the same clique (in addition to families) [ {} ]
wolffd@0: % root      - the root of the junction tree will be a clique that contains this set of nodes [N]
wolffd@0: % stages    - stages{t} is a set of nodes we want to eliminate before stages{t+1}, ... [ {1:N} ]
wolffd@0: %
wolffd@0: % e.g., engine = jtree_inf_engine(bnet, 'maximize', 1);
wolffd@0: %
wolffd@0: % For more details on the junction tree algorithm, see
wolffd@0: % - "Probabilistic networks and expert systems", Cowell, Dawid, Lauritzen and Spiegelhalter, Springer, 1999
wolffd@0: % - "Inference in Belief Networks: A procedural guide", C. Huang and A. Darwiche, 
wolffd@0: %      Intl. J. Approximate Reasoning, 15(3):225-263, 1996.
wolffd@0: 
wolffd@0: 
wolffd@0: % set default params
wolffd@0: N = length(bnet.dag);
wolffd@0: clusters = {};
wolffd@0: root = N;
wolffd@0: stages = { 1:N };
wolffd@0: maximize = 0;
wolffd@0: 
wolffd@0: if nargin >= 2
wolffd@0:   args = varargin;
wolffd@0:   nargs = length(args);
wolffd@0:   if ~isstr(args{1})
wolffd@0:     error('the interface to jtree has changed; now, onodes is not allowed and all optional params must be passed by name')
wolffd@0:   end
wolffd@0:   for i=1:2:nargs
wolffd@0:     switch args{i},
wolffd@0:      case 'clusters', clusters = args{i+1}; 
wolffd@0:      case 'root',     root = args{i+1}; 
wolffd@0:      case 'stages',   stages = args{i+1}; 
wolffd@0:      case 'maximize', maximize = args{i+1};
wolffd@0:      otherwise,  
wolffd@0:       error(['invalid argument name ' args{i}]);       
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: engine = init_fields;
wolffd@0: engine = class(engine, 'jtree_inf_engine', inf_engine(bnet));
wolffd@0: 
wolffd@0: engine.maximize = maximize;
wolffd@0: 
wolffd@0: onodes = bnet.observed;
wolffd@0: 
wolffd@0: %[engine.jtree, dummy, engine.cliques, B, w, elim_order, moral_edges, fill_in_edges, strong] = ...
wolffd@0: %    dag_to_jtree(bnet, onodes, stages, clusters);
wolffd@0: 
wolffd@0: porder = determine_elim_constraints(bnet, onodes);
wolffd@0: strong = ~isempty(porder);
wolffd@0: ns = bnet.node_sizes(:);
wolffd@0: ns(onodes) = 1; % observed nodes have only 1 possible value
wolffd@0: [engine.jtree, root2, engine.cliques, B, w] = ...
wolffd@0:     graph_to_jtree(moralize(bnet.dag), ns, porder, stages, clusters);
wolffd@0: 
wolffd@0: 
wolffd@0: engine.cliques_bitv = B;
wolffd@0: engine.clique_weight = w;
wolffd@0: C = length(engine.cliques);
wolffd@0: engine.clpot = cell(1,C);
wolffd@0: 
wolffd@0: % Compute the separators between connected cliques.
wolffd@0: [is,js] = find(engine.jtree > 0);
wolffd@0: engine.separator = cell(C,C);
wolffd@0: for k=1:length(is)
wolffd@0:   i = is(k); j = js(k);
wolffd@0:   engine.separator{i,j} = find(B(i,:) & B(j,:)); % intersect(cliques{i}, cliques{j});
wolffd@0: end
wolffd@0: 
wolffd@0: % A node can be a member of many cliques, but is assigned to exactly one, to avoid
wolffd@0: % double-counting its CPD. We assign node i to clique c if c is the "lightest" clique that
wolffd@0: % contains i's family, so it can accomodate its CPD.
wolffd@0: 
wolffd@0: engine.clq_ass_to_node = zeros(1, N);
wolffd@0: for i=1:N
wolffd@0:   %c = clq_containing_nodes(engine, family(bnet.dag, i));
wolffd@0:   clqs_containing_family = find(all(B(:,family(bnet.dag, i)), 2)); % all selected columns must be 1
wolffd@0:   c = clqs_containing_family(argmin(w(clqs_containing_family)));  
wolffd@0:   engine.clq_ass_to_node(i) = c; 
wolffd@0: end
wolffd@0: 
wolffd@0: % Make the jtree rooted, so there is a fixed message passing order.
wolffd@0: if strong
wolffd@0:   % the last clique is guaranteed to be a strong root
wolffd@0:   % engine.root_clq = length(engine.cliques);
wolffd@0:   
wolffd@0:   % --- 4/17/2010, by Wei Sun (George Mason University):
wolffd@0:   % It has been proved that the last clique is not necessary to be the  
wolffd@0:   % strong root, instead, a clique called interface clique, that contains
wolffd@0:   % all discrete parents and at least one continuous node from a connected
wolffd@0:   % continuous component in a CLG, is guaranteed to be a strong root.
wolffd@0:   engine.root_clq = findroot(bnet, engine.cliques) ;
wolffd@0: else
wolffd@0:   % jtree_dbn_inf_engine requires the root to contain the interface.
wolffd@0:   % This may conflict with the strong root requirement! *********** BUG *************
wolffd@0:   engine.root_clq = clq_containing_nodes(engine, root);
wolffd@0:   if engine.root_clq <= 0
wolffd@0:     error(['no clique contains ' num2str(root)]);
wolffd@0:   end
wolffd@0: end  
wolffd@0: 
wolffd@0: [engine.jtree, engine.preorder, engine.postorder] = mk_rooted_tree(engine.jtree, engine.root_clq);
wolffd@0: 
wolffd@0: % collect 
wolffd@0: engine.postorder_parents = cell(1,length(engine.postorder));
wolffd@0: for n=engine.postorder(:)'
wolffd@0:   engine.postorder_parents{n} = parents(engine.jtree, n);
wolffd@0: end
wolffd@0: % distribute
wolffd@0: engine.preorder_children = cell(1,length(engine.preorder));
wolffd@0: for n=engine.preorder(:)'
wolffd@0:   engine.preorder_children{n} = children(engine.jtree, n);
wolffd@0: end
wolffd@0: 
wolffd@0:   
wolffd@0: 
wolffd@0: %%%%%%%%
wolffd@0: 
wolffd@0: function engine = init_fields()
wolffd@0: 
wolffd@0: engine.jtree = [];
wolffd@0: engine.cliques = [];
wolffd@0: engine.separator = [];
wolffd@0: engine.cliques_bitv = [];
wolffd@0: engine.clique_weight = [];
wolffd@0: engine.clpot = [];
wolffd@0: engine.clq_ass_to_node = [];
wolffd@0: engine.root_clq = [];
wolffd@0: engine.preorder = [];
wolffd@0: engine.postorder = [];
wolffd@0: engine.preorder_children = [];
wolffd@0: engine.postorder_parents = [];
wolffd@0: engine.maximize = [];
wolffd@0: engine.evidence = [];
wolffd@0: