function [inter CPTnode]= learn_struct_dbn_reveal_f(seqs, ns, max_fan_in)
%function inter = learn_struct_dbn_reveal_f(seqs, ns, max_fan_in, penalty)
% LEARN_STRUCT_DBN_REVEAL Learn inter-slice adjacency matrix given fully observable discrete time series
% inter = learn_struct_dbn_reveal(seqs, node_sizes, max_fan_in, penalty)
%
% seqs{l}{i,t} = value of node i in slice t of time-series l.
%   If you have a single time series in an N*T array D, use
%      seqs = { num2cell(D) }.
%   If you have L time series, each of length T, in an N*T*L array D, use
%      seqs= cell(1,L); for l=1:L, seqs{l} = num2cell(D(:,:,l)); end
%   or, in vectorized form,
%      seqs = squeeze(num2cell(num2cell(D),[1 2]));
% Currently the data is assumed to be discrete (1,2,...)
%
% node_sizes(i) is the number of possible values for node i
% max_fan_in is the largest number of parents we allow per node (default: N)
%
% inter(i,j) = 1 iff node in slice t connects to node j in slice t+1
%
% The parent set for each node in slice 2 is computed by a K2-like search
% strategy.
% Since all the nodes are observed, we do not need to use an inference engine.
% And since we are only learning the inter-slice matrix, we do not need to check for cycles.
%
% This algorithm is described in
% - "REVEAL: A general reverse engineering algorithm for inference of genetic network
%      architectures", Liang et al. PSB 1998
% - "Extended dependency analysis of large systems",
%       Roger Conant, Intl. J. General Systems, 1988, vol 14, pp 97-141
% - "Learning the structure of DBNs", Friedman, Murphy and Russell, UAI 1998.

n = length(ns);

if nargin < 3, max_fan_in = n; end

inter = zeros(n,n);

if ~iscell(seqs)
    data{1} = seqs;
end

nseq = length(seqs);
nslices = 0;
data = cell(1, nseq);
for l=1:nseq
    nslices = nslices + size(seqs{l}, 2);
    data{l} = cell2num(seqs{l})'; % each row is a case
end
ndata = nslices - nseq; % subtract off the initial slice of each sequence


% We concatenate the sequences as in the following example.
% Let there be 2 sequences of lengths 4 and 5, with n nodes per slice,
% and let i be the target node.
% Then we construct following matrix D
%
% s{1}{1,1} ... s{1}{1,3}     s{2}{1,1} ... s{2}{1,4}
% ....
% s{1}{n,1} ... s{1}{n,3}     s{2}{n,1} ... s{2}{n,4}
% s{1}{i,2} ... s{1}{i,4}     s{2}{i,2} ... s{2}{i,5}
%
% D(1:n, i) is the i'th input and D(n+1, i) is the i'th output.
%
% We concatenate each sequence separately to avoid treating the transition
% from the end of one sequence to the beginning of another as a "normal" transition.

%progression bar
global bar;
set(bar,'maximum',n);
set(bar,'StringPainted',1);
set(bar,'String','Please wait...');

CPTnode=cell(n,1);

for i=1:n

    set(bar,'value',i);
    calin=java.util.Date();

    D = [];
    for l=1:nseq
        T = size(seqs{l}, 2);
        A = cell2num(seqs{l}(:, 1:T-1));
        B = cell2num(seqs{l}(i, 2:T));
        C = [A;B];
        D = [D C];
    end

    target = n+1;

    %1 parent
    bic_score = zeros(1, n);
    CPTn=cell(n,1);
    ns2 = [ns ns(i)];
    for j=1:n
        dom = [j target];
        counts = compute_counts(D(dom, :), ns2(dom));
        CPT = mk_stochastic(counts);
        bic_score(j) = bic_score_family(counts, CPT, ndata);
        CPTn{j}=CPT;
    end

    jmax = argmax(bic_score);
    curr_score=bic_score(jmax);
    par=jmax;
    CPTnode{i}=CPTn{jmax};   

    %K2-like search: try adding 1 parent at a time
    done=0;

    while ~done && (length(par)<max_fan_in)
  
        %possible parents
        pps=mysetdiff(1:n,par);
        nps=length(pps);

        if nps>0

            score=zeros(nps,1);
            CPTn=cell(nps,1);
            %try adding 1 par at a time
            for k=1:nps
                ps=[par pps(k)];
                dom = [ps target];
                counts = compute_counts(D(dom, :), ns2(dom));
                CPT = mk_stochastic(counts);
                score(k) = bic_score_family(counts, CPT, ndata);
                CPTn{k}=CPT;

            end

            kmax=argmax(score);
            newscore=score(kmax);

            if newscore>curr_score
                par=[par pps(kmax)];
                curr_score=newscore;
                CPTnode{i}=CPTn{kmax};
            else
                done=1;

            end


        else
            done=1;

        end
   
    end

    inter(par,i)=1;
    
    %calculate approximate time to the end of algorithm
    calfine=java.util.Date();
    time(i)=calfine.getTime()-calin.getTime(); %time for variable i in milliseconds
    wait=round((mean(time))*(n-i)/1000); %approximate wait time in seconds
    if (wait<=120)
        set(bar,'String',strcat('Please wait...',num2str(wait),' seconds left'));
    else
        set(bar,'String',strcat('Please wait...',num2str(ceil(wait/60)),' minutes left'));
    end


end


