function [rbm,output] = 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)
% output= hidden activation probabilities for every input (after training)
% input = dimension
% nh    = number of hidden units
e=0.01;

nv=size(input);   output = zeros(nh);
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;  
    %Sample hidden units given input, then reconstruct visible and hidden.
    v = input;
    whos
    h = 1./(1 + exp(-(w *v + hb)));
    hs= h > rand(nh);
    vr= 1./(1 + exp(-(w'*hs+ vb)));
    hr= 1./(1 + exp(-(w *vr+ hb)));

    %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 ));
    output(:,:,batch)=h;  %
  %fprintf(1, ' error %.4f  \n', err/(nv*batsz*nbat));
  plot( epoch, err/nv, '.'); drawnow;
end
rbm = {w,hb,vb}; %alternatively, rbm = struct('w',w,'hb',hb,'vb',vb);