% Copyright (c) 2016 by Miguel A. Carreira-Perpinan
% for use in CSE176 Introduction to Machine Learning at UC Merced

% LDA for MNIST data
% This version adds small values to the diagonal of SW (to make it full rank)
% and then does LDA.

% ------------------------------- Load data -------------------------------
load MNIST.mat; [N,D] = size(train_fea); DD = sqrt(D);
train_fea = double(train_fea)/255; train_gnd = double(train_gnd);

% ---------------------------------- LDA ----------------------------------

% Subset of points to compute LDA on
I = find(train_gnd<5);			% 1s to 4s
%I=[find(train_gnd==4);find(train_gnd==9)]; % 4s & 9s (differ in top stroke)
%I=[find(train_gnd==7);find(train_gnd==9)]; % 7s & 9s (differ in 1-2 strokes)
I = I(randperm(length(I),1000));
%I = randperm(N,1000);	% rnd subset of all digit classes

% Instances X is of DxN, labels Y is of 1xN, C = classes, K = number of classes
X = train_fea(I,:)'; Y = train_gnd(I)'; C = unique(Y); K = length(C);

% Mean for each class and total within-class scatter matrix SW
D = size(X,1); M = zeros(D,K); SW = zeros(D,D); Nc = zeros(1,K);
for k=1:K
  J = find(Y==C(k)); Nc(k) = length(J);	% points in class C(k)
  M(:,k) = mean(X(:,J),2); SW = SW + Nc(k)*cov(X(:,J)',1);
end
% Total mean and between-class scatter matrix SB
m = mean(M,2); SB = zeros(D,D);
for k=1:K SB = SB + Nc(k)*(M(:,k)-m)*(M(:,k)-m)'; end

% SW should be full rank (= D) and rank(SB) <= K
SW = SW + 1e-10*mean(trace(SW))*eye(D,D);
rank(SB), rank(SW)

% Eigenvectors and LDA
[U,d] = eig(SB,SW); [d,J] = sort(diag(d),1,'descend'); U = U(:,J); W = U;

% Project training set (we obtain all projections 1..D at once)
Z = W'*X;

% ---------------------------------- Plots ----------------------------------

col = {'rgbkcmy','o+x.+sdv^<>ph'};	% colors and markers for points

% Plot eigenvalues and explained variance
figure(1); clf;
subplot(2,1,1); plot(d,'bo-'); axis([0 D+1 0 1.05*max(d)]);
xlabel('LDA component'); ylabel('\lambda'); title('eigenvalues');
subplot(2,1,2); plot(cumsum(d)/sum(d),'ro-'); axis([0 D+1 0 1.05]);
xlabel('LDA component'); ylabel('var'); title('explained variance');
set(gcf,'Name','eigenvalues of inv(SB)*SW and explained variance');

% Plot class mean as images
figure(2); clf; colormap(gray(256));
m1 = bsxfun(@plus,M,m);
k2 = ceil(sqrt(K)); k1 = ceil(K/k2);
for k=1:K
  subplot(k1,k2,k); imagesc(reshape(m1(:,k),DD,DD),[0 1]);
  set(gca,'XTick',[],'YTick',[]); axis image;
  title(['mean of class ' num2str(C(k))]);
end
set(gcf,'Name','mean of each class');

% Plot projections in 2D
c1 = 1; c2 = 2;				% LDA components to project on
figure(3); clf; hold on;
for k=1:K
  [k1,k2] = ind2sub([length(col{1}) length(col{2})],k);
  J = find(Y==C(k)); plot(Z(c1,J),Z(c2,J),[col{1}(k1) col{2}(k2)]);
end
hold off; daspect([1 1 1]); box on;
xlabel(['LDA' num2str(c1)]); ylabel(['LDA' num2str(c2)]);
legend(num2str(C'));
AX = axis; P = get(gca,'Position');
set(gcf,'Name','projection of the dataset on first 2 LDAs');

% Plot projections in 2D as digit images
c1 = 1; c2 = 2;				% LDA components to project on
w = 0.05;				% width of image relative to [0,1]
ax = P(3)/(AX(2)-AX(1)); bx = P(1) - ax*AX(1);	% map coordinate -> figure
ay = P(4)/(AX(4)-AX(3)); by = P(2) - ay*AX(3);
figure(30); clf; axis(AX); daspect([1 1 1]); box on;
xlabel(['LDA' num2str(c1)]); ylabel(['LDA' num2str(c2)]);
hold on;
for i=randperm(length(Y),50)
  axes('position',[ax*Z(c1,i)+bx-w/2 ay*Z(c2,i)+by-w/2 w w]);
  Xi = reshape(X(:,i)+m,DD,DD); colormap(gray(256));
  h = imagesc(1-Xi,[0 1]); alpha(h,1*(Xi>0)); axis off; daspect([1 1 1]);
end
hold off;
set(gcf,'Name','projection of a data subset on first 2 LDAs (as images)');

% Plot projections in 3D
c1 = 1; c2 = 2; c3 = 3;			% LDA components to project on
figure(4); clf; hold on;
for k=1:K
  [k1,k2] = ind2sub([length(col{1}) length(col{2})],k);
  J = find(Y==C(k)); plot3(Z(c1,J),Z(c2,J),Z(c3,J),[col{1}(k1) col{2}(k2)]);
end
hold off; daspect([1 1 1]); box on;
xlabel(['LDA' num2str(c1)]); ylabel(['LDA' num2str(c2)]);
zlabel(['LDA' num2str(c3)]); legend(num2str(C'));
set(gcf,'Name','projection of a data subset on first 3 components');
% Fly around to obtain different views
fr = 40;
el = linspace(-90,90,fr); el = [el el([end-1:-1:2])];
az = 30+linspace(0,360,length(el)+1); az = az([1:end-1]);
for k=1:length(az)
  view(az(k),el(k)); title(['az=' num2str(az(k),3) '; el=' num2str(el(k),3)]);
  pause(0.2);
end

% Plot top K-1 eigenvectors as images ("Fisherdigits")
L = K-1; tmp = U(:,1:L); c1 = min(tmp(:)); c2 = max(tmp(:));
figure(5); clf; colormap(parula(256));
k2 = ceil(sqrt(L+1)); k1 = ceil((L+1)/k2);
for i=1:L
  subplot(k1,k2,i); imagesc(reshape(U(:,i),DD,DD));
  set(gca,'XTick',[],'YTick',[]); axis image; title(['L = ' num2str(i)]);
end
subplot(k1,k2,L+1); imagesc(linspace(c2,c1,256)',[c1 c2]); daspect([1 50 1]);
set(gca,'XTick',[],'YTick',[1 128 256],'YTickLabel',num2str([c2;0;c1],2));
title('colorbar');
%colorbar;
set(gcf,'Name','top K-1 eigenvectors as images ("Fisherdigits")');

