wolffd@0: function mlr_plot(X, Y, W, D)
wolffd@0: 
wolffd@0:     [d, n] = size(X);
wolffd@0:     %%%
wolffd@0:     % Color-coding
wolffd@0:     CODES = {   'b.', 'r.', 'g.', 'c.', 'k.', 'm.', 'y.', ...
wolffd@0:                 'b+', 'r+', 'g+', 'c+', 'k+', 'm+', 'y+', ...
wolffd@0:                 'b^', 'r^', 'g^', 'c^', 'k^', 'm^', 'y^', ...
wolffd@0:                 'bx', 'rx', 'gx', 'cx', 'kx', 'mx', 'yx', ...
wolffd@0:                 'bo', 'ro', 'go', 'co', 'ko', 'mo', 'yo'};
wolffd@0: 
wolffd@0:     %%%
wolffd@0:     % First, PCA-plot of X
wolffd@0: 
wolffd@0:     figure;
wolffd@0: 
wolffd@0:     z = sum(X,3);
wolffd@0:     subplot(3,2,[1 3]), pcaplot(z, eye(d), Y, CODES), title('Native');
wolffd@0:     subplot(3,2,[2 4]), pcaplot(X, W, Y, CODES), title('Learned');
wolffd@0: 
wolffd@0:     [vecs, vals] = eig(z * z');
wolffd@0:     subplot(3,2,5), bar(sort(real(diag(vals)), 'descend')), title('X*X'' spectrum'), axis tight;
wolffd@0: 
wolffd@0:     if size(X,3) == 1
wolffd@0:         [vecs, vals] = eig(W);
wolffd@0:         vals = real(diag(vals));
wolffd@0:     else
wolffd@0:         if size(W,3) == 1
wolffd@0:             vals = real(W(:));
wolffd@0:         else
wolffd@0:             vals = [];
wolffd@0:             for i = 1:size(W,3)
wolffd@0:                 [vecs, vals2] = eig(W(:,:,i));
wolffd@0:                 vals = [vals ; real(diag(vals2))];
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0:     subplot(3,2,6), bar(sort(vals, 'descend')), title('W spectrum'), axis tight;
wolffd@0: 
wolffd@0:     if nargin < 4
wolffd@0:         return;
wolffd@0:     end
wolffd@0: 
wolffd@0:     %%%
wolffd@0:     % Now show some diagnostics
wolffd@0: 
wolffd@0:     figure;
wolffd@0:     subplot(2,1,1), semilogy(D.f), title('Objective');
wolffd@0:     subplot(2,1,2), barh([D.time_SO, D.time_solver D.time_total]), ...
wolffd@0:                     title('Time% in SO/solver/total');
wolffd@0: 
wolffd@0: 
wolffd@0: function pcaplot(X, W, Y, CODES)
wolffd@0:     if size(X,3) == 1
wolffd@0:         A = X' * W * X;
wolffd@0:     else
wolffd@0:         A = 0;
wolffd@0:         if size(W,3) == 1
wolffd@0:             for i = 1:size(X,3)
wolffd@0:                 A = A + X(:,:,i)' * bsxfun(@times, W(:,i), X(:,:,i));
wolffd@0:             end
wolffd@0:         else
wolffd@0:             for i = 1:size(X,3)
wolffd@0:                 A = A + X(:,:,i)' * W(:,:,i) * X(:,:,i);
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0:     [v,d] = eigs(A, 3);
wolffd@0:     X2 = d.^0.5 * v';
wolffd@0:     hold on;
wolffd@0:     for y = 1:max(Y)
wolffd@0:         z = Y == y;
wolffd@0:         scatter3(X2(1, z), X2(2, z), X2(3,z), CODES{y});
wolffd@0:     end
wolffd@0: 
wolffd@0:     axis equal;