function [ GenLog ] = gen_cdbnUpLow ( JointsCRBM, CRBMBody, CRBMBodyTop , 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)

n1 = JointsCRBM{1}.order; % Currently Up and Low have the same order 
n2 = CRBMBody.order;
toplayer = 3;

Atopflat = reshape(CRBMBodyTop.model.A,CRBMBody.numhid,CRBMBodyTop.order*CRBMBody.numhid);
Btopflat = reshape(CRBMBodyTop.model.B,CRBMBodyTop.numhid,CRBMBodyTop.order*CRBMBody.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 = JointsCRBM{1}.order + CRBMBody.order + CRBMBodyTop.order;


numcases = max_clamped - JointsCRBM{1}.order;
numdims = size(initdata,2);

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


data = zeros(numcases,numdims,n1+1);
dataindex = n1+1:max_clamped;

data(:,:,1) = visible(dataindex,:); %store current data
%store delayed data
for hh=1:n1
  data(:,:,hh+1) =visible(dataindex-hh,:);
end


%%%%% joint specific
lasthididx = 0;

for jj = 1:length(JointsCRBM)
    
    JointsCRBM{jj}.hididx = lasthididx+1:lasthididx+1+JointsCRBM{jj}.numhid-1;
    lasthididx = lasthididx+1+JointsCRBM{jj}.numhid-1;
    
    %Calculate contributions from directed visible-to-hidden connections
    bjstar = zeros(JointsCRBM{jj}.numhid,numcases);
    for hh = 1:n1
        bjstar = bjstar + JointsCRBM{jj}.model.B(:,:,hh)*data(:,JointsCRBM{jj}.jointidx,hh+1)';
    end
    
    %Calculate "posterior" probability -- hidden state being on
    %Note that it isn't a true posterior
    eta =  JointsCRBM{jj}.model.w*(data(:,JointsCRBM{jj}.jointidx,1)./JointsCRBM{jj}.gsd)' + ...   %bottom-up connections
        repmat(JointsCRBM{jj}.model.bj, 1, numcases) + ...    %static biases on unit
        bjstar;                           %dynamic biases
    
    hposteriors = 1./(1 + exp(-eta));     %logistic
    
    %initialize hidden layer 1
    %first n1 frames are just padded
    hidden1{jj} = ones(numframes,JointsCRBM{jj}.numhid);
    hidden1{jj}(n1+1:max_clamped,:) = hposteriors';
    
end

%%%%%
hidden1 = cell2mat(hidden1);
hidden{1} = hidden1;


CRBMBody.jointidx = [];
CRBMBodyTop.jointidx = [];
CRBMBody.hididx = [];
CRBMBodyTop.hididx = [];

CDBN.layer(1) = JointsCRBM{1};
CDBN.layer(2) = CRBMBody;
CDBN.layer(3) = CRBMBodyTop;

layerN = 2;
hidden{2} = gen_init(CDBN, layerN, numframes, hidden{layerN-1}, max_clamped);


%initialize second layer (first n1+n2 frames padded)
hidden{3} = ones(numframes,CDBN.layer(3).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{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
    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 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
   for jj = 1:length(JointsCRBM)
        
        %Add contributions from autoregressive connections
        bistar = zeros(JointsCRBM{jj}.numdims,1);
        for hh=1:n1
            bistar = bistar +  JointsCRBM{jj}.model.A(:,:,hh)*visible(tt-hh,JointsCRBM{jj}.jointidx)' ;
        end
        
        %Mean-field approx; visible units are Gaussian
        %(filtering distribution is the data we just generated)
        topdown = JointsCRBM{jj}.gsd.*(hidden1(tt,JointsCRBM{jj}.hididx)*JointsCRBM{jj}.model.w);
        visible(tt,JointsCRBM{jj}.jointidx) = topdown + ...             %top down connections
            JointsCRBM{jj}.model.bi' + ...                            %static biases
            bistar';                              %dynamic biases
        
    end
end

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