wolffd@0: function engine = belprop_gdl_inf_engine(gdl, varargin) 
wolffd@0: % BELPROP_GDL_INF_ENGINE Make a belief propagation inference engine for a GDL graph
wolffd@0: % engine = belprop_gdl_inf_engine(gdl_graph, ...)
wolffd@0: %
wolffd@0: % If the GDL graph is a tree, this will give exact results.
wolffd@0: %
wolffd@0: % The following optional arguments can be specified in the form of name/value pairs:
wolffd@0: % [default in brackets]
wolffd@0: % e.g., engine = belprop_inf_engine(gdl, 'tol', 1e-2, 'max_iter', 10)
wolffd@0: %
wolffd@0: % protocol - 'tree' means send messages up then down the tree,
wolffd@0: %            'parallel' means use synchronous updates ['parallel']
wolffd@0: % max_iter - max. num. iterations [ 2*num_nodes ]
wolffd@0: % momentum - weight assigned to old message in convex combination (useful for damping oscillations) [0]
wolffd@0: % tol - tolerance used to assess convergence [1e-3]
wolffd@0: % maximize - 1 means use max-product, 0 means use sum-product [0]
wolffd@0: 
wolffd@0: 
wolffd@0: engine = init_fields;
wolffd@0: engine = class(engine, 'belprop_gdl_inf_engine');
wolffd@0: 
wolffd@0: % set default params
wolffd@0: N = length(gdl.G);
wolffd@0: engine.protocol = 'parallel';
wolffd@0: engine.max_iter = 2*N;
wolffd@0: engine.momentum = 0;
wolffd@0: engine.tol = 1e-3;
wolffd@0: engine.maximize = 0;
wolffd@0: 
wolffd@0: engine = set_params(engine, varargin);
wolffd@0: 
wolffd@0: engine.gdl = gdl;
wolffd@0: 
wolffd@0: if strcmp(engine.protocol, 'tree')
wolffd@0:   % Make a rooted tree, so there is a fixed message passing order.
wolffd@0:   root = N;
wolffd@0:   [engine.tree, engine.preorder, engine.postorder, height, cyclic] = mk_rooted_tree(gdl.G, root);
wolffd@0:   assert(~cyclic);
wolffd@0: end
wolffd@0: 
wolffd@0: % store results computed by enter_evidence here
wolffd@0: ndoms = length(gdl.doms);
wolffd@0: nvars = length(gdl.vars);
wolffd@0: engine.marginal_domains = cell(1, ndoms);
wolffd@0: 
wolffd@0: % to compute the marginal on each variable, we need to know which domain to marginalize
wolffd@0: % and we want to choose the lightest. We compute the weight once we have seen the evidence.
wolffd@0: engine.dom_weight = [];
wolffd@0: engine.evidence = [];
wolffd@0: 
wolffd@0: 
wolffd@0: %%%%%%%%%
wolffd@0: 
wolffd@0: function engine = init_fields()
wolffd@0: 
wolffd@0: engine.protocol = [];
wolffd@0: engine.gdl = [];
wolffd@0: engine.max_iter = [];
wolffd@0: engine.momentum = [];
wolffd@0: engine.tol = [];
wolffd@0: engine.maximize = [];
wolffd@0: engine.marginal_domains = [];
wolffd@0: engine.evidence = [];
wolffd@0: engine.tree = [];
wolffd@0: engine.preorder = [];
wolffd@0: engine.postorder = [];
wolffd@0: engine.dom_weight = [];