function net = fastTrainRBM(data,para)
%% parameter specification
% epsilonw      = 0.001;   % Learning rate for weights 
% epsilonvb     = 0.001;   % Learning rate for biases of visible units 
% epsilonhb     = 0.001;   % Learning rate for biases of hidden units 

epsilonw      = 0.01;   % Learning rate for weights 
epsilonvb     = 0.01;   % Learning rate for biases of visible units 
epsilonhb     = 0.01;   % Learning rate for biases of hidden units 

weightcost  = 0.0002;   
initialmomentum  = 0.5;
finalmomentum    = 0.9;

nIn = para.visibleSize; nHid = para.hiddenSize;
maxEpoch = para.maxepoch;


Winc = zeros(nIn, nHid); hidbinc = zeros(1,nHid); visbinc = zeros(1,nIn);
W = 0.001*randn(nIn, nHid); 
hidbias = zeros(1,nHid); visbias = zeros(1,nIn);
% hidbias = -1+0.01*randn(1,nHid); visbias = 0.01*randn(1,nIn);

%% training in action

% batchid = genBatch(data, para.nbatch);
batchid = genBatchFixNum(size(data,1),para.nbatch);
nBatch = size(batchid,2);
batchSize = batchid(2,1) - batchid(1,1) + 1

for epoch = 1: maxEpoch
   errsum = 0;
   seq = randperm(size(data,1));      
   trdata = data(seq,:);
   
   
   for bid = 1: nBatch
       numcases = batchid(2,bid) - batchid(1,bid) + 1;
       V0 = trdata(batchid(1,bid):batchid(2,bid),:);
       
       %% One step of Gibbs sampling (CD-1)     
       
       h0prob = sigmoid(V0*W+repmat(hidbias,numcases,1));              
       h0sample = h0prob > rand(numcases,nHid);
       
       v1prob = h0sample*W' + repmat(visbias,numcases,1);
       h1prob = sigmoid(v1prob*W+repmat(hidbias,numcases,1));
       
       posprods = V0' * h0prob;
       negprods = v1prob' * h1prob;
       
       poshidact   = sum(h0prob);
       posvisact = sum(V0);
       
       neghidact = sum(h1prob);
       negvisact = sum(v1prob);
             
       
       %% Update weights and biases
       errsum = errsum + sum(sum((v1prob - V0).^2));
       if epoch>5
         momentum=finalmomentum;
       else
         momentum=initialmomentum;
       end
       

        Winc = momentum*Winc + ...
                    epsilonw*( (posprods-negprods)/numcases - weightcost*W);
        visbinc = momentum*visbinc + (epsilonvb/numcases)*(posvisact-negvisact);
        hidbinc = momentum*hidbinc + (epsilonhb/numcases)*(poshidact-neghidact);

        W = W + Winc;
        visbias = visbias + visbinc;
        hidbias = hidbias + hidbinc;       
       
              
   end
   
   disp(['Reconstruction error at epoch ', int2str(epoch), ' is ', num2str(errsum)]);
   
end

net.W = W'; net.W1 = net.W; 
net.visbias = visbias; net.hidbias = hidbias;
net.b = hidbias'; net.b1 = net.b;
