matthiasm@8: function inter = learn_struct_dbn_reveal(seqs, ns, max_fan_in, penalty)
matthiasm@8: % LEARN_STRUCT_DBN_REVEAL Learn inter-slice adjacency matrix given fully observable discrete time series
matthiasm@8: % inter = learn_struct_dbn_reveal(seqs, node_sizes, max_fan_in, penalty)
matthiasm@8: % 
matthiasm@8: % seqs{l}{i,t} = value of node i in slice t of time-series l.
matthiasm@8: %   If you have a single time series in an N*T array D, use
matthiasm@8: %      seqs = { num2cell(D) }.
matthiasm@8: %   If you have L time series, each of length T, in an N*T*L array D, use
matthiasm@8: %      seqs= cell(1,L); for l=1:L, seqs{l} = num2cell(D(:,:,l)); end
matthiasm@8: %   or, in vectorized form,
matthiasm@8: %      seqs = squeeze(num2cell(num2cell(D),[1 2]));
matthiasm@8: % Currently the data is assumed to be discrete (1,2,...)
matthiasm@8: %
matthiasm@8: % node_sizes(i) is the number of possible values for node i
matthiasm@8: % max_fan_in is the largest number of parents we allow per node (default: N)
matthiasm@8: % penalty is weight given to the complexity penalty (default: 0.5)
matthiasm@8: %  A penalty of 0.5 gives the BIC score.
matthiasm@8: %  A penalty of 0 gives the ML score.
matthiasm@8: %  Maximizing likelihood is equivalent to maximizing mutual information between parents and child.
matthiasm@8: %
matthiasm@8: % inter(i,j) = 1 iff node in slice t connects to node j in slice t+1
matthiasm@8: %
matthiasm@8: % The parent set for each node in slice 2 is computed by evaluating all subsets of nodes in slice 1,
matthiasm@8: % and picking the largest scoring one. This takes O(n^k) time per node, where n is the num. nodes
matthiasm@8: % per slice, and k <= n is the max fan in.
matthiasm@8: % Since all the nodes are observed, we do not need to use an inference engine.
matthiasm@8: % And since we are only learning the inter-slice matrix, we do not need to check for cycles.
matthiasm@8: %
matthiasm@8: % This algorithm is described in
matthiasm@8: % - "REVEAL: A general reverse engineering algorithm for inference of genetic network
matthiasm@8: %      architectures", Liang et al. PSB 1998
matthiasm@8: % - "Extended dependency analysis of large systems",
matthiasm@8: %       Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141
matthiasm@8: % - "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998.
matthiasm@8: 
matthiasm@8: n = length(ns);
matthiasm@8: 
matthiasm@8: if nargin < 3, max_fan_in = n; end
matthiasm@8: if nargin < 4, penalty = 0.5; end
matthiasm@8: 
matthiasm@8: inter = zeros(n,n);
matthiasm@8: 
matthiasm@8: if ~iscell(seqs)
matthiasm@8:   data{1} = seqs;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: nseq = length(seqs);
matthiasm@8: nslices = 0;
matthiasm@8: data = cell(1, nseq);
matthiasm@8: for l=1:nseq
matthiasm@8:   nslices = nslices + size(seqs{l}, 2);
matthiasm@8:   data{l} = cell2num(seqs{l})'; % each row is a case
matthiasm@8: end
matthiasm@8: ndata = nslices - nseq; % subtract off the initial slice of each sequence
matthiasm@8: 
matthiasm@8: % We concatenate the sequences as in the following example.
matthiasm@8: % Let there be 2 sequences of lengths 4 and 5, with n nodes per slice,
matthiasm@8: % and let i be the target node.
matthiasm@8: % Then we construct following matrix D 
matthiasm@8: %
matthiasm@8: % s{1}{1,1} ... s{1}{1,3}     s{2}{1,1} ... s{2}{1,4}
matthiasm@8: % ....
matthiasm@8: % s{1}{n,1} ... s{1}{n,3}     s{2}{n,1} ... s{2}{n,4}
matthiasm@8: % s{1}{i,2} ... s{1}{i,4}     s{2}{i,2} ... s{2}{i,5}
matthiasm@8: %
matthiasm@8: % D(1:n, i) is the i'th input and D(n+1, i) is the i'th output.
matthiasm@8: % 
matthiasm@8: % We concatenate each sequence separately to avoid treating the transition
matthiasm@8: % from the end of one sequence to the beginning of another as a "normal" transition.
matthiasm@8: 
matthiasm@8: 
matthiasm@8: for i=1:n
matthiasm@8:   D = [];
matthiasm@8:   for l=1:nseq
matthiasm@8:     T = size(seqs{l}, 2);
matthiasm@8:     A = cell2num(seqs{l}(:, 1:T-1));
matthiasm@8:     B = cell2num(seqs{l}(i, 2:T));
matthiasm@8:     C = [A;B];
matthiasm@8:     D = [D C];
matthiasm@8:   end
matthiasm@8:   SS = subsets(1:n, max_fan_in, 1); % skip the empty set 
matthiasm@8:   nSS = length(SS);
matthiasm@8:   bic_score = zeros(1, nSS);
matthiasm@8:   ll_score = zeros(1, nSS);
matthiasm@8:   target = n+1;
matthiasm@8:   ns2 = [ns ns(i)];
matthiasm@8:   for h=1:nSS
matthiasm@8:     ps = SS{h};
matthiasm@8:     dom = [ps target];
matthiasm@8:     counts = compute_counts(D(dom, :), ns2(dom));
matthiasm@8:     CPT = mk_stochastic(counts);
matthiasm@8:     [bic_score(h), ll_score(h)] = bic_score_family(counts, CPT, ndata);
matthiasm@8:   end
matthiasm@8:   if penalty == 0
matthiasm@8:     h = argmax(ll_score);
matthiasm@8:   else
matthiasm@8:     h = argmax(bic_score);
matthiasm@8:   end
matthiasm@8:   ps = SS{h};
matthiasm@8:   inter(ps, i) = 1;
matthiasm@8: end