function [ GenLog ] = gen_crbm_relu ( CRBM , 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 a single level CRBM (or the first level of a multi-
% level CRBM) to generate data

% The program assumes that the following variables are set externally:
% numframes    -- number of frames to generate
% fr           -- a starting frame from initdata (for initialization)

numGibbs = 30; %number of alternating Gibbs iterations 

% rootindx = [1:6];
% lowidx = [rootindx 7:22];
% upidx = [23:52];
   

initdata = InitData.batchdata;

numdims = size(initdata,2);

%initialize visible layer
visible = zeros(numframes,CRBM.numdims);
visible(1:CRBM.order,:) = initdata(fr:fr+CRBM.order-1,:);

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

%initialize hidden layer
hidden{1} = ones(numframes,CRBM.numhid);

for tt=CRBM.order+1:numframes
  
  %initialize using the last frame + noise
  visible(tt,:) = visible(tt-1,:) + 0.01*randn(1,numdims);
  
  %Dynamic biases aren't re-calculated during Alternating Gibbs
  %First, add contributions from autoregressive connections 
  bistar = zeros(numdims,1);
  for hh=1:CRBM.order
    %should modify to data * A'
    bistar = bistar +  CRBM.model.A(:,:,hh)*visible(tt-hh,:)' ;
  end

  %Next, add contributions to hidden units from previous time steps
  bjstar = zeros(CRBM.numhid,1);
  for hh = 1:CRBM.order
    bjstar = bjstar + CRBM.model.B(:,:,hh)*visible(tt-hh,:)';
  end

  %Gibbs sampling
  for gg = 1:numGibbs
    %Calculate posterior probability -- hidden state being on (estimate)
    %add in bias
    bottomup =  CRBM.model.w*(visible(tt,:)./CRBM.gsd)';
    eta = bottomup + ...                   %bottom-up connections
      CRBM.model.bj + ...                            %static biases on unit
      bjstar;                              %dynamic biases
    
    hposteriors =relu(-eta);% 1./(1 + exp(-eta));      %logistic
    
    hidden{1}(tt,:) = double(hposteriors' > rand(1,CRBM.numhid));
    
    %Downward pass; visibles are Gaussian units
    %So find the mean of the Gaussian    
    topdown = CRBM.gsd.*(hidden{1}(tt,:)*CRBM.model.w);
    
    %Mean-field approx
    visible(tt,:) = topdown + ...            %top down connections
      CRBM.model.bi' + ...                             %static biases
      bistar';                               %dynamic biases     
  
  end

  %If we are done Gibbs sampling, then do a mean-field sample
  %(otherwise very noisy)  
  topdown = CRBM.gsd.*(hposteriors'*CRBM.model.w);                 

  visible(tt,:) = topdown + ...              %top down connections
      CRBM.model.bi' + ...                             %static biases
      bistar';                               %dynamic biases


end

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

end

  

