function [W,bvis,bhid] = RBMtrain(pat,nhid,nepoch,e,momentum)
%train a restricted Boltzmann machine on a set of training patterns
%pat: training patterns; pattern x input unit
%nhid: number of hidden units in the network
%nepoch: number of times to run through the training set
%e: learning rate
%momentum: carryover of gradients between epochs, for momentum learning
%returns final network weights (W=weights, bvis=visible biases, bhid=hidden biases) 

%set up network
nvis = size(pat,2); %number of visible units
W = zeros(nvis,nhid); %initialize weight matrix; visible x hidden
bvis = zeros(nvis,1); %initialize biases for visibile units
bhid = zeros(nhid,1); %initialize biases for hidden units

%initialize variables for the learning routine
npat = size(pat,1); %number of training patterns
DeltaW = zeros(nvis,nhid); %for storing weight updates
Deltabv = zeros(nvis,1); %for storing visible bias updates
Deltabh = zeros(nhid,1); %for storing hidden bias updates

%some variables for tracking convergence
corrs = zeros(nepoch,1); %correlation of weight updates on each epoch with previous updates
rmse = zeros(nepoch,1); %reconstruction error by epoch

for epoch = 1:nepoch %loop through epochs
    for i = 1:npat %loop through patterns
        vpos = pat(i,:)'; %input pattern
        hpos = rand(nhid,1) < 1./(1+exp(-W'*vpos-bhid)); %sample values for hidden units (positive phase)
        vneg = 1./(1+exp(-W*hpos-bvis)); %new values for visible units (negative phase); using expected values rather than binary samples, for efficiency (following Hinton, 2002, sec 8)
        hneg = rand(nhid,1) < 1./(1+exp(-W'*vneg-bhid)); %sample values for hidden units (negative phase)
        rmse(epoch) = rmse(epoch) + mean((vpos-vneg).^2); %add mean squared reconstruction error for this pattern (vpos vs vneg)
        DeltaW = DeltaW + vpos*hpos' - vneg*hneg'; %increment the updates to the weights
        Deltabv = Deltabv + vpos - vneg; %increment the updates to the visible biases
        Deltabh = Deltabh + hpos - hneg; %increment the updates to the hidden biases
    end
    rmse(epoch) = sqrt(rmse(epoch)/npat); %convert to root mean squared error (averaged across all patterns and all hidden units)
    W = W + e*DeltaW; %update weights once at end of epoch
    bvis = bvis + e*Deltabv; %update visible biases
    bhid = bhid + e*Deltabh; %update hidden biases
    disp(['completed epoch ' num2str(epoch)])
    if epoch>1
        corrs(epoch) = corr(DeltaW(:)-DeltaWprev(:)*momentum,DeltaWprev(:)); %correlation of new and old weight updates, across weights; indicates whether learning is progressing in a consistent direction
        disp(['rmse = ' num2str(rmse(epoch)) ', update correlation = ' num2str(corrs(epoch))]) %progress report
    end
    DeltaWprev = DeltaW; %save current updates to use for correlation measure on next epoch
    DeltaW = DeltaW * momentum; %depreciate updates for next epoch (momentum learning)
    Deltabv = Deltabv * momentum; %depreciate updates for next epoch (momentum learning)
    Deltabh = Deltabh * momentum; %depreciate updates for next epoch (momentum learning)
    figure(1),plot([corrs(2:epoch) rmse(2:epoch)]),drawnow %graphical tracking of convergence; at asymptote, rmse should be constant and corrs should be zero
end