function [ GenLog ] = gen_cdbn0 ( CDBN , clamped_h, 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 2-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)

A2flat = reshape(CDBN.layer(2).model.A,CDBN.layer(1).numhid,CDBN.layer(2).order*CDBN.layer(1).numhid);
B2flat = reshape(CDBN.layer(2).model.B,CDBN.layer(2).numhid,CDBN.layer(2).order*CDBN.layer(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:2%size(CDBN.layer,2)
    max_clamped = max_clamped + CDBN.layer(ll).order;
end

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

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

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

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';





%initialize second layer (first n1+n2 frames padded)
hidden{2} = ones(numframes,CDBN.layer(2).numhid);

%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{1}(max_clamped:-1:max_clamped+1-CDBN.layer(2).order,:)',CDBN.layer(1).numhid*CDBN.layer(2).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
  hidden{1}(tt,:) = hidden{1}(tt-1,:);
  
  %Dynamic biases aren't re-calculated during Alternating Gibbs
  bistar = A2flat*past;
  bjstar = B2flat*past;

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

  %If we are done Gibbs sampling, then do a mean-field sample
  
  topdown = hposteriors'*CDBN.layer(2).model.w;  %Very noisy if we don't do mean-field here

  eta = topdown + ...                        %top down connections
      CDBN.layer(2).model.bi' + ...                             %static biases
      bistar';                               %dynamic biases
  
  hidden{1}(tt,:) = 1./(1 + exp(-eta));

  hidden{1}(:,:) = 0;
    hidden{1}(:,clamped_h) = 1;
  
  %update the past
  past(CDBN.layer(1).numhid+1:end) = past(1:end-CDBN.layer(1).numhid); %shift older history down
  past(1:CDBN.layer(1).numhid) = hidden{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 visible data\n');
for tt=max_clamped+1:numframes
    %Add contributions from autoregressive connections
    bistar = zeros(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';                              %dynamic biases
end

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

end