% DPS Correntropy based, hierarchical density preserving data split
% 
%   R	      = DPS(A,LEVELS,CLASSWISE)
%   [R H]   = DPS(A,LEVELS,CLASSWISE)
%
% INPUT
%   A			    Input dataset
%   LEVELS		Number of split levels, default: 3
%   CLASSWISE	Use (1, default) or ignore (0) label information
%
% OUTPUT
%   R			Index array with rotation set with 2^LEVELS folds
%   H			Hierarchy of splits
%
% DESCRIPTION
% Density Preserving Sampling (DPS) divides the input dataset into a given
% number of folds (2^LEVELS) by maximizing the correntropy between the folds
% and can be used as an alternative for cross-validation. The procedure is
% deterministic, so unlike cross-validation it does not need to be repeated.
%
% REFERENCE
% M. Budka, B. Gabrys, Correntropy-based density-preserving data sampling
% as an alternative to standard cross-validation, IJCNN2010, 1-8
% http://www.budka.co.uk/
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, CROSSVAL
  
function [R H] = dps(A,levels,classwise)
	
 	if (nargin<3) || isempty(classwise), classwise = 1; end
 	if (nargin<2) || isempty(levels), levels = 3; end
	
	islabtype(A,'crisp');

	H = zeros(levels,size(A,1));

	idxs = cell(1,levels+1);
	idxs{1} = {1:size(A,1)};
	for i = 1:levels
		for j = 1:2^(i-1)
			t = helper(A(idxs{i}{j},:),classwise);
			idxs{i+1}{2*j-1} = idxs{i}{j}(t{1});
			idxs{i+1}{2*j} = idxs{i}{j}(t{2});
		end
		for j = 1:length(idxs{i+1})
			H(i,idxs{i+1}{j}) = j;
		end
	end
	
	R = H(end,:);

end

function idxs = helper(A,classwise)

	% if the classes aro too small to be divided further, switch to classless mode
	if classwise && all(classsizes(A)>2), c = getsize(A,3);
	else classwise = 0;	c = 1;
	end
	 
	siz = zeros(2,1);
	idx = cell(1,c);

	for i = 1:c
 		if classwise, [B BI] = seldat(A,i);					% select i-th class
		else B = A; BI = (1:size(A,1));
		end
		
		m = length(BI);
		mask = true(1,m);									% mask is used for counting remaining objects
		D = +distm(B) + diag(inf(m,1));						% working distance matrix
		Dorg = D;											% original distance matrix
		idx{i} = nan(2,ceil(m/2));
		
		for j = 1:floor(m/2)
			[mD,I] = min(D,[],1);							% \
			[mmD,J] = min(mD);								%   find two closest objects
			I = I(J(1)); J = J(1);							% /

			mask(I) = 0; mask(J) = 0;						% mark them as used
			
			% split the objects to maximally increase coverage of both subsets
			if (mean(Dorg(I,idx{i}(1,1:j-1))) + mean(Dorg(J,idx{i}(2,1:j-1))) < ...
				mean(Dorg(I,idx{i}(2,1:j-1))) + mean(Dorg(J,idx{i}(1,1:j-1))))
				idx{i}(1,j) = J;
				idx{i}(2,j) = I;
			else  
				idx{i}(1,j) = I;
				idx{i}(2,j) = J;
			end
			
			% remove used objects from the distance matrix
			D(I,:) = inf; D(:,I) = inf;
			D(J,:) = inf; D(:,J) = inf;
		end
		if isempty(j), j = 0; end							% in case the loop is not entered at all
		
		% odd number of points in class
		if sum(mask)>0
			I = find(mask);
			if siz(1)<siz(2)
				idx{i}(1,end) = I;
			elseif siz(1)>siz(2)
				idx{i}(2,end) = I;
			else
				if (mean(Dorg(I,idx{i}(1,1:j))) < mean(Dorg(I,idx{i}(2,1:j))))
					idx{i}(2,j+1) = I;
				else
					idx{i}(1,j+1) = I;
				end
			end
		end
		
		% convert indexes from class-specific to dataset-specific
		idx{i}(1,~isnan(idx{i}(1,:))) = BI(idx{i}(1,~isnan(idx{i}(1,:))));
		idx{i}(2,~isnan(idx{i}(2,:))) = BI(idx{i}(2,~isnan(idx{i}(2,:))));

		% update fold sizes
		siz(1) = siz(1) + sum(~isnan(idx{i}(1,:)));
		siz(2) = siz(2) + sum(~isnan(idx{i}(2,:)));

	end

	idx = cell2mat(idx);
	
	idxs = cell(1,2);
	idxs{1} = idx(1,~isnan(idx(1,:)));
	idxs{2} = idx(2,~isnan(idx(2,:)));
	
end