%% X_r is datamatrix reduced to a 2D space from p-dims.
%  Thus, X_r is N x 2 matrix, y is response variable.

function plot_decision_surface (X_r, y)
	figure();
	hold on
	
	% find plot range
	R = [min(X_r); max(X_r)];
	% find step size
	s = diff(R) /500;
	
	X_range = R(1, 1):s(1):R(2, 1);
	Y_range = R(1, 2):s(1):R(2, 2);
	[xx, yy] = meshgrid(X_range, Y_range);
	
	Grid = [xx(:), yy(:)];
	
	% Train the model 
	model = fitcdiscr(X_r, y, 'DiscrimType','pseudoquadratic');
	labels = unique(y);
	predictions = predict(model, Grid);
	
	% Class colors for labels 1, 2, and 3
	fillcolor={'magenta', 'cyan', 'yellow'};
	
	% Color meshrid predictions - P stores index of 
	% points to be colored with fillcolor{j}
	for j = 1:length(labels)
	    P{j} = ismember(predictions, labels(j));
	    scatter(xx(P{j}), yy(P{j}), 'o' , ...
	    	    'MarkerEdgeColor', fillcolor{j}, ...
	            'MarkerFaceColor', fillcolor{j});
	end
	
	% Place the actual data points with outline edge in black
	for j = 1:length(labels)
		Q{j} = ismember(y, labels(j));
		scatter(X_r(Q{j},1), X_r(Q{j},2), 'o' , ...
		       'MarkerEdgeColor', 'black', ...
	           'MarkerFaceColor', fillcolor{j});
	end
end 
