wolffd@0: function seq = sample_dbn(bnet, varargin)
wolffd@0: % SAMPLE_DBN Generate a random sequence from a DBN.
wolffd@0: % seq = sample_dbn(bnet, ...)
wolffd@0: %
wolffd@0: % seq{i,t} contains the values of the i'th node in the t'th slice.
wolffd@0: %
wolffd@0: % Optional arguments:
wolffd@0: %
wolffd@0: % length - length of sequence  to be generated (can also just use sample_dbn(bnet,T))
wolffd@0: % stop_test - name of a function which is used to decide when to stop;
wolffd@0: %   This will be called as feval(stop_test, seq(:,t))
wolffd@0: %   i.e., stop_test is passed a cell array containing all the nodes in the current slice.   
wolffd@0: % evidence - initial evidence; if evidence{i,t} is non-empty, this node won't be sampled.
wolffd@0: 
wolffd@0: args = varargin;
wolffd@0: nargs = length(args);
wolffd@0: 
wolffd@0: if (nargs == 1) & ~isstr(args{1})
wolffd@0:   % Old syntax: sample_dbn(bnet, T)
wolffd@0:   T = args{1};
wolffd@0: else
wolffd@0:   % get length
wolffd@0:   T = 1;
wolffd@0:   for i=1:2:nargs
wolffd@0:     switch args{i},
wolffd@0:      case 'length',      T = args{i+1}; 
wolffd@0:      case 'evidence',    T = size(args{i+1}, 2);
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: ss = length(bnet.intra);
wolffd@0: % set default arguments
wolffd@0: seq = cell(ss, T);
wolffd@0: stop_test = [];
wolffd@0: for i=1:2:nargs
wolffd@0:   switch args{i},
wolffd@0:    case 'evidence',    seq = args{i+1}; % initialise observed nodes
wolffd@0:    case 'stop_test',  stop_test = args{i+1};
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: t = 1;
wolffd@0: for i=1:ss
wolffd@0:   if ~isempty(stop_test) | isempty(seq{i,t})
wolffd@0:     ps = parents(bnet.dag, i);
wolffd@0:     e = bnet.equiv_class(i,1);
wolffd@0:     pvals = seq(ps);
wolffd@0:     seq{i,t} = sample_node(bnet.CPD{e}, pvals);
wolffd@0:     %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
wolffd@0:   end
wolffd@0: end
wolffd@0: t = 2;
wolffd@0: done = 0;
wolffd@0: while ~done
wolffd@0:   for i=1:ss
wolffd@0:     if ~isempty(stop_test) | isempty(seq{i,t})
wolffd@0:       ps = parents(bnet.dag, i+ss) + (t-2)*ss;
wolffd@0:       e = bnet.equiv_class(i,2);
wolffd@0:       pvals = seq(ps);
wolffd@0:       seq{i,t} = sample_node(bnet.CPD{e}, pvals);
wolffd@0:       %fprintf('sample i=%d,t=%d,val=%d,ps\n', i, t, seq(i,t)); pvals(:)'
wolffd@0:     end
wolffd@0:   end
wolffd@0:   if ~isempty(stop_test)
wolffd@0:     done = feval(stop_test, seq(:,t));
wolffd@0:   else
wolffd@0:     if t==T
wolffd@0:       done = 1;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   t = t + 1;
wolffd@0: end