% Version 1.000
%
% Code provided by Ruslan Salakhutdinov and Geoff Hinton
%
% Permission is granted for anyone to copy, use, modify, or distribute this
% program and accompanying programs and documents for any purpose, provided
% this copyright notice is retained and prominently displayed, along with
% a note saying that the original programs are available from our
% web page.
% The programs and documents are distributed without any warranty, express or
% implied.  As the programs were written for research purposes only, they have
% not been tested to the degree that would be advisable in any important
% application.  All use of these programs is entirely at the user's own risk.

% Version 1.100
%
% Updated by Computational Cognitive Neuroscience Lab
% University of Padova
% ccnl.psy.unipd.it
%
% Implementation on Cluster using MPITB for Octave

if (rnk==0)
	partbatchdata = batchdata(:,:,1:numbatches/siz); % Prepare single processor's batch matrixes
    	% 1st MPI: distribute data size
    	MPI_Bcast(A, 0, MPI_COMM_WORLD);
    	% 2nd MPI: distribute parts of batchdata
	MPI_Scatter(batchdata, partbatchdata, 0, MPI_COMM_WORLD);   % This master distributes the matrix to each single processor
    	batchdata=partbatchdata;
else    
	A = [0, 0, 0];
    	% 1st MPI: get data size
	MPI_Bcast(A, 0, MPI_COMM_WORLD);
	A(1,3) = A(1,3)/siz; % Number of batches per processor.
    	% 2nd MPI: get a batch
	batchdata = zeros(A);
	MPI_Scatter([], batchdata, 0, MPI_COMM_WORLD); 	% 2nd: get a partial batch
end
    

% MASTER Initialize symmetric weights and biases & distributes them to SLAVES
[numcases numdims numbatches]=size(batchdata);
if(rnk==0)
	vishid = 0.1*randn(numdims, numhid);		
else
	vishid= zeros(numdims,numhid);
end
MPI_Bcast(vishid, 0, MPI_COMM_WORLD);

hidbiases  = zeros(1,numhid);
visbiases  = zeros(1,numdims);
poshidprobs = zeros(numcases,numhid);
neghidprobs = zeros(numcases,numhid);
posprods    = zeros(numdims,numhid);
negprods    = zeros(numdims,numhid);
vishidinc  = zeros(numdims,numhid);
hidbiasinc = zeros(1,numhid);
visbiasinc = zeros(1,numdims);
batchposhidprobs= zeros(numcases,numhid,numbatches);

for epoch = 1:maxepoch 
	for i = 1:numbatches
		%%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	  	data = batchdata(:,:,i);
	  	poshidprobs = 1./(1 + exp(-data*vishid - repmat(hidbiases,numcases,1)));
		if (epoch==maxepoch)    
	  		batchposhidprobs(:,:,i)=poshidprobs;
		end
	  	posprods    = data' * poshidprobs;
	  	poshidact   = sum(poshidprobs);
	  	posvisact = sum(data);
		%%%%%%%%% END OF POSITIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 	  	poshidstates = poshidprobs > rand(numcases,numhid);
		%%%%%%%%% START NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	  	negdata = 1./(1 + exp(-poshidstates*vishid' - repmat(visbiases,numcases,1)));
	  	neghidprobs = 1./(1 + exp(-negdata*vishid - repmat(hidbiases,numcases,1)));    
	  	negprods  = negdata'*neghidprobs;
	  	neghidact = sum(neghidprobs);
	  	negvisact = sum(negdata); 
		%%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		if epoch>5
			momentum=finalmomentum;
		else
		  	momentum=initialmomentum;
        	end

		X1 = (posprods-negprods)/numcases;
	    	X2 = (posvisact-negvisact)/numcases;
	    	X3 = (poshidact-neghidact)/numcases;
                
        	if (rnk==0)
			Y1 = ones(size(X1));
		    	Y2 = ones(size(X2));
		    	Y3 = ones(size(X3));
			% master receive the slave's weights update 
		    	MPI_Reduce(X1, Y1, MPI_SUM, 0, MPI_COMM_WORLD);	 
		    	MPI_Reduce(X2, Y2, MPI_SUM, 0, MPI_COMM_WORLD);	
		    	MPI_Reduce(X3, Y3, MPI_SUM, 0, MPI_COMM_WORLD);
            		Y1 /= siz;		
            		Y2 /= siz;
            		Y3 /= siz;
			%%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
			vishidinc = momentum*vishidinc + epsilonw*( Y1 - weightcost*vishid);
			visbiasinc = momentum*visbiasinc + epsilonvb*Y2;
			hidbiasinc = momentum*hidbiasinc + epsilonhb*Y3;

		    	vishid = vishid + vishidinc;
		    	visbiases = visbiases + visbiasinc;
		    	hidbiases = hidbiases + hidbiasinc;
        	else
            		% slaves send the weights update
            		MPI_Reduce(X1, [], MPI_SUM, 0, MPI_COMM_WORLD);	 
            		MPI_Reduce(X2, [], MPI_SUM, 0, MPI_COMM_WORLD);	
            		MPI_Reduce(X3, [], MPI_SUM, 0, MPI_COMM_WORLD);           
        	end
		% distribute new weights and bias
		MPI_Bcast(vishid, 0, MPI_COMM_WORLD);
		MPI_Bcast(visbiases, 0, MPI_COMM_WORLD);
		MPI_Bcast(hidbiases, 0, MPI_COMM_WORLD);
	end % End of Batch
end % End of Epoch

% Get together all outputs: (1) make room for data from all cores (2) Gather the data (3) rename this var
if (rnk==0)
	batchposhidprobstot=zeros(size(batchposhidprobs,1),size(batchposhidprobs,2),size(batchposhidprobs,3)*siz);
end
MPI_Gather(batchposhidprobs, batchposhidprobstot, 0, MPI_COMM_WORLD);
if (rnk==0) 
  	batchposhidprobs=batchposhidprobstot;    % Output from the positive phase
end
