% Configuration
% -------------
nv = 784;   % number of visible units
nh = 400;   % number of hidden units
mb = 100;   % size of the minibatch
it = 100;   % number of iterations (WARNING: set to 10000 or higher to produce the final results)

% Initialize data
% ---------------
data = load('mnist.mat');       % Load the dataset
X = single(data.x)/255.0 > 0.5; % Preprocess the dataset and convert it into a matrix of binary values
W = normrnd(0,0.001,[nv,nh]);   % Weight matrix
a = zeros(1,nv);                % Bias vector for visible units
b = zeros(1,nh);                % Bias vector for hidden units

% Algorithm
% ---------
for i=1:it,
	
	% Take a minibatch from the dataset
	P = ceil(unifrnd(0,10000,[mb,1]));
	v = X(P(:),:);
	
	% Perform three steps of Gibbs sampling
	h = % TO COMPLETE
	V = % TO COMPLETE
	H = % TO COMPLETE
	
	% Compute the directions of parameter updates
	dW = % TO COMPLETE
	da = % TO COMPLETE
	db = % TO COMPLETE
	
	% Update the parameters
	W = W + 0.001 * dW;
	b = b + 0.001 * db;
	a = a + 0.001 * da;
	
	% Print iterations and some information on the parameters of the RBM and the reconstruction error
	sprintf('i=%5d   std(W)=%.3f   mean(a)=%.3f   mean(b)=%.3f   e=%.3f\n',i,std(W(:)),mean(a(:)),mean(b(:)),mean(mean(abs(v-V))))
end

% Display the weights and the reconstructed digits
% ------------------------------------------------
wdisp = reshape(permute(reshape(W,[28,28,20,20]),[2,4,1,3]),[28*20,28*20]); % Create a human-readable view of matrix W
wdisp = wdisp - min(wdisp(:));
wdisp = wdisp / max(wdisp(:));
imwrite(wdisp,'w.png')

xdisp = reshape(permute(reshape(V,[10,10,28,28]),[4,2,3,1]),[28*10,28*10]);
xdisp = xdisp - min(xdisp(:));
xdisp = xdisp / max(xdisp(:));
imwrite(xdisp,'x.png')

