function [ GenLog ] = gen_cdbn ( CDBN , InitData, numframes, fr)
% Version 2.0 
%
% Modified as a function by Omid Alemi
%
% Code provided by Graham Taylor, Geoff Hinton and Sam Roweis 
%
% For more information, see:
%     http://www.cs.toronto.edu/~gwtaylor/publications/nips2006mhmublv
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.
%
% This program uses the n-level CRBM to generate data
% More efficient version than the original
%
% The program assumes that the following variables are set externally:
% numframes    -- number of frames to generate
% fr           -- a starting frame from initdata (for initialization)

toplayer = CDBN.NumLayers;

Atopflat = reshape(CDBN.layer(toplayer).model.A,CDBN.layer(toplayer-1).numhid,CDBN.layer(toplayer).order*CDBN.layer(toplayer-1).numhid);
Btopflat = reshape(CDBN.layer(toplayer).model.B,CDBN.layer(toplayer).numhid,CDBN.layer(toplayer).order*CDBN.layer(toplayer-1).numhid);

numGibbs = 30; %number of alternating Gibbs iterations

initdata = InitData.batchdata;

%We have saved some initialization data in "initdata"
%How many frames of it do we need to make a prediction for the first h2
%frame? we need n1 + n2 + ... frames of data

max_clamped = 0;
for ll = 1:size(CDBN.layer,2)
    max_clamped = max_clamped + CDBN.layer(ll).order;
end

%use this data to get the posteriors at layer top - 1
numcases = CDBN.layer(2).order; %number of hidden units to generate
%numdims = CDBN.layer(1).numdims;

%initialize visible data
visible = zeros(numframes,CDBN.layer(1).numdims);
visible(1:max_clamped,:) = initdata(fr:fr+max_clamped-1,:);

data = zeros(numcases,CDBN.layer(1).numdims,CDBN.layer(1).order+1);
dataindex = max_clamped-CDBN.layer(1).order+1:max_clamped;

a = visible(dataindex,:);
b = data(:,:,1);
data(:,:,1) = visible(dataindex,:); %store current data
%store delayed data
for hh=1:CDBN.layer(1).order
  data(:,:,hh+1) =visible(dataindex-hh,:);
end

%Calculate contributions from directed visible-to-hidden connections
bjstar = zeros(CDBN.layer(1).numhid,numcases);
for hh = 1:CDBN.layer(1).order
    bjstar = bjstar + CDBN.layer(1).model.B(:,:,hh)*data(:,:,hh+1)';
end

%Calculate "posterior" probability -- hidden state being on
%Note that it isn't a true posterior
eta = CDBN.layer(1).model.w*(data(:,:,1)./CDBN.layer(1).gsd)' + ...   %bottom-up connections
    repmat(CDBN.layer(1).model.bj, 1, numcases) + ...    %static biases on unit
    bjstar;                           %dynamic biases

hposteriors = 1./(1 + exp(-eta));     %logistic

%initialize hidden layer 1
%first CDBN.layer(1).order frames are just padded
hidden{1} = ones(numframes,CDBN.layer(1).numhid);
hidden{1}(CDBN.layer(1).order+1:CDBN.layer(1).order+CDBN.layer(2).order,:) = hposteriors';

%%%%%%%%%%%%%%

numcases = CDBN.layer(3).order; %number of hidden units to generate

%initialize visible data
hidden{1} = zeros(numframes,CDBN.layer(1).numhid);
% hidden{1}(1:max_clamped,:) = hidden{(fr:fr+max_clamped-1,:);

data = zeros(numcases,CDBN.layer(2).numdims,CDBN.layer(2).order+1);
dataindex = max_clamped-CDBN.layer(2).order+1:max_clamped;


data(:,:,1) = hidden{1}(dataindex,:); %store current data
%store delayed data
for hh=1:CDBN.layer(2).order
  data(:,:,hh+1) =hidden{1}(dataindex-hh,:);
end

%Calculate contributions from directed visible-to-hidden connections
bjstar = zeros(CDBN.layer(2).numhid,numcases);
for hh = 1:CDBN.layer(2).order
    bjstar = bjstar + CDBN.layer(2).model.B(:,:,hh)*data(:,:,hh+1)';
end

%Calculate "posterior" probability -- hidden state being on
%Note that it isn't a true posterior
eta = CDBN.layer(2).model.w*(data(:,:,1)./CDBN.layer(2).gsd)' + ...   %bottom-up connections
    repmat(CDBN.layer(2).model.bj, 1, numcases) + ...    %static biases on unit
    bjstar;                           %dynamic biases

hposteriors2 = 1./(1 + exp(-eta));     %logistic


%%%%%%%%%


%initialize hidden layer 1
%first CDBN.layer(1).order frames are just padded
hidden{2} = ones(numframes,CDBN.layer(2).numhid);
hidden{2}(CDBN.layer(1).order+1:CDBN.layer(2).order+CDBN.layer(3).order,:) = hposteriors2';



hidden{toplayer} = ones(numframes,CDBN.layer(toplayer).numhid);

%initialize next hidden layers (first n1+n2 + ... frames padded)
% for ll=2:toplayer
%     hidden{ll} = ones(numframes,CDBN.layer(ll).numhid);
% end


%keep the recent past in vector form
%for input to directed links
%it's slightly faster to pre-compute this vector and update it (by
%shifting) at each time frame, rather than to fully recompute each time
%frame
past = reshape(hidden{toplayer-1}(max_clamped:-1:max_clamped+1-CDBN.layer(toplayer).order,:)',CDBN.layer(toplayer-1).numhid*CDBN.layer(toplayer).order,1);

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%First generate a hidden sequence (top layer)
%Then go down through the first CRBM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Generating hidden states\n');
for tt=max_clamped+1:numframes  
  %initialize using the last frame
  %noise is not added for binary units
 
  
%   for ll=1:toplayer-1
%      hidden{ll}(tt,:) = hidden{ll}(tt-1,:);
%   end

  hidden{toplayer-1}(tt,:) = hidden{toplayer-1}(tt-1,:);

  %Dynamic biases aren't re-calculated during Alternating Gibbs
  bistar = Atopflat*past;
  bjstar = Btopflat*past;

  %Gibbs sampling
  for gg = 1:numGibbs
      
      
      %Calculate posterior probability -- hidden state being on (estimate)
      %add in bias
      bottomup =  CDBN.layer(toplayer).model.w*hidden{toplayer-1}(tt,:)';
      eta = bottomup + ...                   %bottom-up connections
          CDBN.layer(toplayer).model.bj + ...                            %static biases on unit
          bjstar;                              %dynamic biases
      
      hposteriors = 1./(1 + exp(-eta));      %logistic
      
      hidden{toplayer}(tt,:) = double(hposteriors' > rand(1,CDBN.layer(toplayer).numhid)); %Activating hiddens
      
      
      %Downward pass; visibles are binary logistic units
      topdown = hidden{toplayer}(tt,:)*CDBN.layer(toplayer).model.w;
      
      eta = topdown + ...                      %top down connections
          CDBN.layer(toplayer).model.bi' + ...                             %static biases
          bistar';                               %dynamic biases
      
      hidden{toplayer-1}(tt,:) = 1./(1 + exp(-eta));      %logistic
      
      
  end

  %If we are done Gibbs sampling, then do a mean-field sample
%   
%   topdown = hposteriors'*CDBN.layer(toplayer).model.w;  %Very noisy if we don't do mean-field here
% 
%   eta = topdown + ...                        %top down connections
%       CDBN.layer(toplayer).model.bi' + ...                             %static biases
%       bistar';                               %dynamic biases
%   
%   hidden{toplayer-1}(tt,:) = 1./(1 + exp(-eta));

  %update the past
  past(CDBN.layer(toplayer-1).numhid+1:end) = past(1:end-CDBN.layer(toplayer-1).numhid); %shift older history down
  past(1:CDBN.layer(toplayer-1).numhid) = hidden{toplayer-1}(tt,:); %place most recent frame at top
  
  
  if mod(tt,10)==0
    fprintf('Finished frame %d\n',tt);
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Now that we've decided on the "filtering distribution", generate visible
%data through CRBM 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Generating intermediate hidden units and visible data\n');
for tt=max_clamped+1:numframes
    
    % intermediate hidden layers
    for ll=toplayer-1:2
        %Add contributions from autoregressive connections
        bistar = zeros(CDBN.layer(ll).numdims,1);
        for hh=1:CDBN.layer(ll).order
            bistar = bistar +  CDBN.layer(ll).model.A(:,:,hh)*hidden{ll-1}(tt-hh,:)' ;
        end
        
        
        %Mean-field approx; visible units are Gaussian
        %(filtering distribution is the data we just generated)
        topdown = CDBN.layer(ll).gsd.*(hidden{ll}(tt,:)*CDBN.layer(ll).model.w);
        hidden{ll-1}(tt,:) = topdown + ...             %top down connections
            CDBN.layer(ll).model.bi' + ...                            %static biases
            bistar';                              %dynamic biases
        
        eta = topdown + ...                      %top down connections
          CDBN.layer(ll).model.bi' + ...                             %static biases
          bistar';                               %dynamic biases
      
      hidden{ll-1}(tt,:) = 1./(1 + exp(-eta));      %logistic
      
    end
    
    
    % visible layer
     %Add contributions from autoregressive connections
    bistar = zeros(CDBN.layer(1).numdims,1);
    for hh=1:CDBN.layer(1).order
        bistar = bistar + CDBN.layer(1).model.A(:,:,hh)*visible(tt-hh,:)' ;
    end

    %Mean-field approx; visible units are Gaussian
    %(filtering distribution is the data we just generated)
    topdown = CDBN.layer(1).gsd.*(hidden{1}(tt,:)*CDBN.layer(1).model.w);
    visible(tt,:) = topdown + ...             %top down connections
        CDBN.layer(1).model.bi' + ...                            %static biases
        bistar'; 
    
end

GenLog.visible = visible;
GenLog.hidden = hidden;

end
