wolffd@0: function bnet = mk_higher_order_dbn(intra, inter, node_sizes, varargin)
wolffd@0: % MK_DBN Make a Dynamic Bayesian Network.
wolffd@0: %
wolffd@0: % BNET = MK_DBN(INTRA, INTER, NODE_SIZES, ...) makes a DBN with arcs
wolffd@0: % from i in slice t to j in slice t iff intra(i,j) = 1, and 
wolffd@0: % from i in slice t to j in slice t+1 iff inter(i,j) = 1,
wolffd@0: % for i,j in {1, 2, ..., n}, where n = num. nodes per slice, and t >= 1.
wolffd@0: % node_sizes(i) is the number of values node i can take on.
wolffd@0: % The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
wolffd@0: % See also mk_bnet.
wolffd@0: %
wolffd@0: % Optional arguments [default in brackets]
wolffd@0: % 'discrete' - list of discrete nodes [1:n]
wolffd@0: % 'observed' - the list of nodes which will definitely be observed in every slice of every case [ [] ]
wolffd@0: % 'eclass1' - equiv class for slice 1 [1:n]
wolffd@0: % 'eclass2' - equiv class for slice 2 [tie nodes with equivalent parents to slice 1]
wolffd@0: %    equiv_class1(i) = j means node i in slice 1 gets its parameters from bnet.CPD{j},
wolffd@0: %    i.e., nodes i and j have tied parameters.
wolffd@0: % 'intra1' - topology of first slice, if different from others
wolffd@0: % 'names' - a cell array of strings to be associated with nodes 1:n [{}]
wolffd@0: %    This creates an associative array, so you write e.g.
wolffd@0: %     'evidence(bnet.names{'bar'}) = 42' instead of  'evidence(2} = 42' 
wolffd@0: %     assuming names = { 'foo', 'bar', ...}.
wolffd@0: %    
wolffd@0: % For backwards compatibility with BNT2, arguments can also be specified as follows
wolffd@0: %   bnet = mk_dbn(intra, inter, node_sizes, dnodes, eclass1, eclass2, intra1)
wolffd@0: %
wolffd@0: % After calling this function, you must specify the parameters (conditional probability
wolffd@0: % distributions) using bnet.CPD{i} = gaussian_CPD(...) or tabular_CPD(...) etc.
wolffd@0: 
wolffd@0: 
wolffd@0: n = length(intra);
wolffd@0: ss = n;
wolffd@0: bnet.nnodes_per_slice = ss;
wolffd@0: bnet.intra = intra;
wolffd@0: bnet.inter = inter;
wolffd@0: bnet.intra1 = intra;
wolffd@0: 
wolffd@0: % As this method is used to generate a higher order Markov Model
wolffd@0: % also connect from time slice t - i -> t with i > 1 has to be 
wolffd@0: % taken into account.
wolffd@0: 
wolffd@0: %inter should be a three dimensional array where inter(:,:,i)
wolffd@0: %describes the connections from time-slice t - i to t.  
wolffd@0: [rows,columns,order] = size(inter);
wolffd@0: assert(rows    == n);
wolffd@0: assert(columns == n);
wolffd@0: dag = zeros((order + 1)*n);
wolffd@0: 
wolffd@0: i = 0;
wolffd@0: while i <= order
wolffd@0:     j = i;
wolffd@0:     while j <= order
wolffd@0:         if j == i
wolffd@0:             dag(1 + i*n:(i+1)*n,1+i*n:(i+1)*n) = intra;
wolffd@0:         else
wolffd@0:             dag(1+i*n:(i+1)*n,1+j*n:(j+1)*n) = inter(:,:,j - i);
wolffd@0:         end
wolffd@0:         j = j + 1;
wolffd@0:     end;
wolffd@0:     i = i + 1;
wolffd@0: end;
wolffd@0: 
wolffd@0: bnet.dag = dag;
wolffd@0: bnet.names = {};
wolffd@0: 
wolffd@0: directed = 1;
wolffd@0: if ~acyclic(dag,directed)
wolffd@0:   error('graph must be acyclic')
wolffd@0: end
wolffd@0: 
wolffd@0: % Calculation of the equivalence classes
wolffd@0: bnet.eclass1 = 1:n;
wolffd@0: bnet.eclass = zeros(order + 1,ss);
wolffd@0: bnet.eclass(1,:) = 1:n;
wolffd@0: for i = 1:order
wolffd@0:     bnet.eclass(i+1,:) = bnet.eclass(i,:);
wolffd@0:     for j = 1:ss 
wolffd@0:         if(isequal(parents(dag,(i-1)*n+j)+ss,parents(dag,(i*n + j))))
wolffd@0: 	   %fprintf('%d has isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))
wolffd@0:         else
wolffd@0: 	   bnet.eclass(i + 1,j) = max(bnet.eclass(i+1,:))+1;
wolffd@0: 	   %fprintf('%d has non isomorphic parents, eclass %d \n',j,bnet.eclass(i,j))  
wolffd@0: 	end;
wolffd@0:     end;
wolffd@0: end;
wolffd@0: bnet.eclass1 = 1:n;
wolffd@0: 
wolffd@0: % To be compatible with whe rest of the code 
wolffd@0: bnet.eclass2 = bnet.eclass(2,:);
wolffd@0: 
wolffd@0: dnodes = 1:n;
wolffd@0: bnet.observed = [];
wolffd@0: 
wolffd@0: if nargin >= 4
wolffd@0:   args = varargin;
wolffd@0:   nargs = length(args);
wolffd@0:   if ~isstr(args{1})
wolffd@0:     if nargs >= 1 dnodes = args{1}; end
wolffd@0:     if nargs >= 2 bnet.eclass1 = args{2}; bnet.eclass(1,:) = args{2}; end
wolffd@0:     if nargs >= 3 bnet.eclass2 = args{3}; bnet.eclass(2,:) = args{2}; end
wolffd@0:     if nargs >= 4 bnet.intra1 = args{4}; end
wolffd@0:   else
wolffd@0:     for i=1:2:nargs
wolffd@0:       switch args{i},
wolffd@0:        case 'discrete', dnodes = args{i+1}; 
wolffd@0:        case 'observed', bnet.observed = args{i+1}; 
wolffd@0:        case 'eclass1',  bnet.eclass1 = args{i+1}; bnet.eclass(1,:) = args{i+1}; 
wolffd@0:        case 'eclass2',  bnet.eclass2 = args{i+1}; bnet.eclass(2,:) = args{i+1};
wolffd@0:        case 'eclass',   bnet.eclass = args{i+1};  
wolffd@0:        case 'intra1',  bnet.intra1 = args{i+1}; 
wolffd@0:        %case 'ar_hmm',  bnet.ar_hmm = args{i+1};  % should check topology
wolffd@0:        case 'names',  bnet.names = assocarray(args{i+1}, num2cell(1:n)); 
wolffd@0:        otherwise,  
wolffd@0: 	error(['invalid argument name ' args{i}]);       
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: bnet.observed = sort(bnet.observed); % for comparing sets
wolffd@0: ns = node_sizes;
wolffd@0: bnet.node_sizes_slice = ns(:)';
wolffd@0: bnet.node_sizes = repmat(ns(:),1,order + 1);
wolffd@0: 
wolffd@0: cnodes = mysetdiff(1:n, dnodes);
wolffd@0: bnet.dnodes_slice = dnodes;
wolffd@0: bnet.cnodes_slice = cnodes;
wolffd@0: bnet.dnodes = dnodes;
wolffd@0: bnet.cnodes = cnodes;
wolffd@0: % To adapt the function to higher order Markov models include dnodes for more 
wolffd@0: % time slices
wolffd@0: for i = 1:order
wolffd@0:     bnet.dnodes = [bnet.dnodes dnodes+i*n];
wolffd@0:     bnet.cnodes = [bnet.cnodes cnodes+i*n];
wolffd@0: end
wolffd@0: 
wolffd@0: % Generieren einer Matrix, deren i-te Spalte die Aequivalenzklassen
wolffd@0: % der i-ten Zeitscheibe enthaelt. 
wolffd@0: bnet.equiv_class = [bnet.eclass(1,:)]';
wolffd@0: for i = 2:(order + 1)
wolffd@0:     bnet.equiv_class = [bnet.equiv_class   bnet.eclass(i,:)'];
wolffd@0: end
wolffd@0: 
wolffd@0: bnet.CPD = cell(1,max(bnet.equiv_class(:)));
wolffd@0: 
wolffd@0: ss = n;
wolffd@0: onodes = bnet.observed;
wolffd@0: hnodes = mysetdiff(1:ss, onodes);
wolffd@0: bnet.hidden_bitv = zeros(1,(order + 1)*ss);
wolffd@0: for i = 0:order
wolffd@0:     bnet.hidden_bitv(hnodes +i*ss) = 1;
wolffd@0: end;
wolffd@0: 
wolffd@0: bnet.parents = cell(1, (order + 1)*ss);
wolffd@0: for i=1:(order + 1)*ss
wolffd@0:   bnet.parents{i} = parents(bnet.dag, i);
wolffd@0: end
wolffd@0: 
wolffd@0: bnet.auto_regressive = zeros(1,ss);
wolffd@0: % ar(i)=1 means (observed) node i depends on i in the  previous slice
wolffd@0: for o=bnet.observed(:)'
wolffd@0:   if any(bnet.parents{o+ss} <= ss)
wolffd@0:     bnet.auto_regressive(o) = 1;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: