function [ GenLog ] = gen_cdbnUpLow ( CRBMUp, CRBMLow, CRBMBody , 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 = CRBMUp.order; % Currently Up and Low have the same order 
n2 = CRBMBody.order;

numhiduplow = CRBMUp.numhid + CRBMLow.numhid;

rootindx = InitData.rootindx;
lowidx = [InitData.rootindx InitData.lowidx];
upidx = [InitData.rootindx InitData.upidx];


lowhididx = 1:CRBMLow.numhid;
uphididx = lowhididx(end)+1:lowhididx(end)+CRBMUp.numhid;


A2flat = reshape(CRBMBody.model.A,numhiduplow,n2*numhiduplow);
B2flat = reshape(CRBMBody.model.B,CRBMBody.numhid,n2*numhiduplow);

initdata = InitData.batchdata;


numGibbs = 30; %number of alternating Gibbs iterations

%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 = n1+n2;

%use this data to get the posteriors at layer 1
numcases = n2; %number of hidden units to generate
numdims = size(initdata,2);

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

%visible(1:max_clamped,lowidx) = initdata(fr:fr+max_clamped-1,lowidx);
%fr = fr + 15; % tease the network
%visible(1:max_clamped,InitData.upidx) = initdata(fr:fr+max_clamped-1,InitData.upidx);

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




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

%Calculate "posterior" probability -- hidden state being on
%Note that it isn't a true posterior
eta =  CRBMUp.model.w*(data(:,upidx,1)./CRBMUp.gsd)' + ...   %bottom-up connections
    repmat(CRBMUp.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
hidden1up = ones(numframes,CRBMUp.numhid);
hidden1up(n1+1:n1+n2,:) = hposteriors';


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

%Calculate "posterior" probability -- hidden state being on
%Note that it isn't a true posterior
eta =  CRBMLow.model.w*(data(:,lowidx,1)./CRBMLow.gsd)' + ...   %bottom-up connections
    repmat(CRBMLow.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
hidden1low = ones(numframes,CRBMLow.numhid);
hidden1low(n1+1:n1+n2,:) = hposteriors';

%%%%%

hidden1 = [hidden1low hidden1up];



%initialize second layer (first n1+n2 frames padded)
hidden2 = ones(numframes,CRBMBody.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(hidden1(max_clamped:-1:max_clamped+1-n2,:)',numhiduplow*n2,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
  hidden1(tt,:) = hidden1(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 =  CRBMBody.model.w*hidden1(tt,:)';
    eta = bottomup + ...                   %bottom-up connections
      CRBMBody.model.bj + ...                            %static biases on unit
      bjstar;                              %dynamic biases
    
    hposteriors = 1./(1 + exp(-eta));      %logistic
    
    hidden2(tt,:) = double(hposteriors' > rand(1,CRBMBody.numhid)); %Activating hiddens
    
    %Downward pass; visibles are binary logistic units     
    topdown = hidden2(tt,:)*CRBMBody.model.w;
        
    eta = topdown + ...                      %top down connections
      CRBMBody.model.bi' + ...                             %static biases
      bistar';                               %dynamic biases   
    
    hidden1(tt,:) = 1./(1 + exp(-eta));      %logistic 
  end

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

  eta = topdown + ...                        %top down connections
      CRBMBody.model.bi' + ...                             %static biases
      bistar';                               %dynamic biases
  hidden1(tt,:) = 1./(1 + exp(-eta));

  %update the past
  past(numhiduplow+1:end) = past(1:end-numhiduplow); %shift older history down
  past(1:numhiduplow) = hidden1(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
    
    % Lower Body
    %Add contributions from autoregressive connections
    bistar = zeros(CRBMUp.numdims,1);
    for hh=1:n1
        bistar = bistar +  CRBMUp.model.A(:,:,hh)*visible(tt-hh,upidx)' ;
    end
    
    %Mean-field approx; visible units are Gaussian
    %(filtering distribution is the data we just generated)
    topdown = CRBMUp.gsd.*(hidden1(tt,uphididx)*CRBMUp.model.w);
    visible(tt,upidx) = topdown + ...             %top down connections
        CRBMUp.model.bi' + ...                            %static biases
        bistar';                              %dynamic biases
    
    % Upper Body
    %Add contributions from autoregressive connections
    bistar = zeros(CRBMLow.numdims,1);
    for hh=1:n1
        bistar = bistar +  CRBMLow.model.A(:,:,hh)*visible(tt-hh,lowidx)' ;
    end
    
    %Mean-field approx; visible units are Gaussian
    %(filtering distribution is the data we just generated)
    topdown = CRBMLow.gsd.*(hidden1(tt,lowhididx)*CRBMLow.model.w);
    visible(tt,lowidx) = topdown + ...             %top down connections
        CRBMLow.model.bi' + ...                            %static biases
        bistar';
    
    
end


GenLog.visible = visible;
GenLog.hidden{1} = hidden1;
GenLog.hidden{2} = hidden2;
end