Daniel@0: function bnet = mk_bnet(dag, node_sizes, varargin)
Daniel@0: % MK_BNET Make a Bayesian network.
Daniel@0: %
Daniel@0: % BNET = MK_BNET(DAG, NODE_SIZES, ...) makes a graphical model with an arc from i to j iff DAG(i,j) = 1.
Daniel@0: % Thus DAG is the adjacency matrix for a directed acyclic graph.
Daniel@0: % The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
Daniel@0: %
Daniel@0: % node_sizes(i) is the number of values node i can take on,
Daniel@0: %   or the length of node i if i is a continuous-valued vector.
Daniel@0: % node_sizes(i) = 1 if i is a utility node.
Daniel@0: % 
Daniel@0: % Below are the names of optional arguments [and their default value in brackets].
Daniel@0: % Pass as 'PropertyName1', PropertyValue1, 'PropertyName2', PropertyValue2, ...
Daniel@0: % 
Daniel@0: % discrete - the list of nodes which are discrete random variables [1:N]
Daniel@0: % equiv_class - equiv_class(i)=j  means node i gets its params from CPD{j} [1:N]
Daniel@0: % observed - the list of nodes which will definitely be observed in every case [ [] ]
Daniel@0: % 'names' - a cell array of strings to be associated with nodes 1:n [{}]
Daniel@0: %    This creates an associative array, so you write e.g.
Daniel@0: %     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42' 
Daniel@0: %     assuming names = { 'foo', 'bar', ...}.
Daniel@0: %
Daniel@0: % e.g., bnet = mk_bnet(dag, ns, 'discrete', [1 3])
Daniel@0: %
Daniel@0: % For backwards compatibility with BNT2, you can also specify the parameters in the following order
Daniel@0: %   bnet = mk_bnet(dag, node_sizes, discrete_nodes, equiv_class)
Daniel@0: 
Daniel@0: n = length(dag);
Daniel@0: 
Daniel@0: % default values for parameters
Daniel@0: bnet.equiv_class = 1:n;
Daniel@0: bnet.dnodes = 1:n; % discrete 
Daniel@0: bnet.observed = [];
Daniel@0: bnet.names = {};
Daniel@0: 
Daniel@0: if nargin >= 3
Daniel@0:   args = varargin;
Daniel@0:   nargs = length(args);
Daniel@0:   if ~isstr(args{1})
Daniel@0:     if nargs >= 1, bnet.dnodes = args{1}; end
Daniel@0:     if nargs >= 2, bnet.equiv_class = args{2}; end
Daniel@0:   else    
Daniel@0:     for i=1:2:nargs
Daniel@0:       switch args{i},
Daniel@0:        case 'equiv_class', bnet.equiv_class = args{i+1}; 
Daniel@0:        case 'discrete',    bnet.dnodes = args{i+1}; 
Daniel@0:        case 'observed',    bnet.observed = args{i+1}; 
Daniel@0:        case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n)); 
Daniel@0:        otherwise,  
Daniel@0: 	error(['invalid argument name ' args{i}]);       
Daniel@0:       end
Daniel@0:     end
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: bnet.observed = sort(bnet.observed); % for comparing sets
Daniel@0: bnet.hidden = mysetdiff(1:n, bnet.observed(:)');
Daniel@0: bnet.hidden_bitv = zeros(1,n);
Daniel@0: bnet.hidden_bitv(bnet.hidden) = 1;
Daniel@0: bnet.dag = dag;
Daniel@0: bnet.node_sizes = node_sizes(:)';
Daniel@0: 
Daniel@0: bnet.cnodes = mysetdiff(1:n, bnet.dnodes);
Daniel@0: % too many functions refer to cnodes to rename it to cts_nodes - 
Daniel@0: % We hope it won't be confused with chance nodes!
Daniel@0: 
Daniel@0: bnet.parents = cell(1,n);
Daniel@0: for i=1:n
Daniel@0:   bnet.parents{i} = parents(dag, i);
Daniel@0: end
Daniel@0: 
Daniel@0: E = max(bnet.equiv_class);
Daniel@0: mem = cell(1,E);
Daniel@0: for i=1:n
Daniel@0:   e = bnet.equiv_class(i);
Daniel@0:   mem{e} = [mem{e} i];
Daniel@0: end
Daniel@0: bnet.members_of_equiv_class = mem;
Daniel@0: 
Daniel@0: bnet.CPD = cell(1, E);
Daniel@0: 
Daniel@0: bnet.rep_of_eclass = zeros(1,E);
Daniel@0: for e=1:E
Daniel@0:   mems = bnet.members_of_equiv_class{e};
Daniel@0:   bnet.rep_of_eclass(e) = mems(1);
Daniel@0: end
Daniel@0: 
Daniel@0: directed = 1;
Daniel@0: if ~acyclic(dag,directed)
Daniel@0:   error('graph must be acyclic')
Daniel@0: end
Daniel@0: 
Daniel@0: bnet.order = topological_sort(bnet.dag);