function rbm = trainrbm(input, nh, epochs)
%Given a series of inputs, trains an RBM with nh hidden units.
%Plotted error is the average across visible units and inputs.
% rbm{1}= weights (w), rbm{2}= hidden biases (hb), rbm{3}= visible biases (vb)
% input = [dimesion, batch_size, num_batches]
% nh    = number of hidden units
e=0.01; % learning rate

[nv,npat]=size(input);   output = zeros(nh,npat);
vb =zeros(nv,1);  hb =zeros(nh,1);  w =.1*randn(nh,nv);
dvb=zeros(nv,1);  dhb=zeros(nh,1);  dw=   zeros(nh,nv);  %d for delta

figure; hold on; xlabel epoch; ylabel 'avg error'; xlim([0 epochs]);
for epoch=1:epochs;
  err=0;  
  for pat=1:npat
    %Sample hidden units given input, then reconstruct visible and hidden.
    v = input(:,pat);
    h = 1./(1 + exp(-(w *v + hb)));
    hs= h > rand(nh,1);
    vr= 1./(1 + exp(-(w'*hs+ vb)));
    hr= 1./(1 + exp(-(w *vr+ hb)))-0.0001*sum(h);

    %Contrastive Divergence rule: dw ~ h*v - hr*vr
    dw  =  e*( (h*v'-hr*vr') );  w = w +dw;
    dvb =  e*(sum(v,2)-sum(vr,2));   vb= vb+dvb;
    dhb =  e*(sum(h,2)-sum(hr,2));   hb= hb+dhb;
    
    err = err   + sum(sum( (v-vr).^2 ));
  end
  %fprintf(1, ' error %.4f  \n', err/(nv*batsz*nbat));
  plot( epoch, err/(nv*npat), '.'); drawnow;
end
rbm = {w,hb,vb}; %alternatively, rbm = struct('w',w,'hb',hb,'vb',vb);