function selidx = pls_multiclass_v2( Data_Tr, Labels_Tr, K, D)
% concatenate K groups of features
%
% input: 
%       Data_Tr         p-by-NTr, original representation of training data
%       Data_Te         p-by-NTe
%       Labels_Tr       1-by-NTr, labels of training data, taking values 1, ..., K
%       K               number of class
%       D               intended reduced dimension per-class
%
% output:
%       Feat_Tr         Dall-by-NTr, Representation in reduced dimension,
%                           Dall is the total number of nonrepetitive
%                           selected indeces
%       Feat_Te         Dall-by-NTe
%       meta.D          dimension per-class
%
% OLS (Parital Least Square) supervised dimension reduction for K-classes
%
% train the selected subspace for each class in one-vs-all mode, and
% concatenate the features
%
% Xiuyuan CHENG 03/10/2014
%
%%%%%

%addpath ./olscode/

[p, NTr] = size( Data_Tr);


% OLS per class 

meta.select_index = cell(numel(K),1);

Feat_Tr = zeros( numel(K)*D, NTr);

Q_cell_tr=cell(numel(K),1);

select_index2=cell(numel(K),1);

selidx = zeros(p,1);

for icls = 1: numel(K)
    fprintf('PLS processing class %d:', icls);
    
    %%% linear regression: y = X * beta
    
    allidx = 1:length(Labels_Tr);
    curidx = allidx(Labels_Tr == K(icls));
    xtr = Data_Tr(:,curidx).'; % training data, N-by-p
    
    [~,maxidx] = max(sum(xtr.^2,1));
    y = xtr(:, maxidx);
    xtr(:,maxidx) = 0;
    
%     y = zeros( NTr, 1); y( Labels_Tr == K(icls) ) = numel(K); %training label, N-by-1 before it was K
%     y=y-mean(y);
    %%% firstly select the bias term/the constant feature, equivalently remove the mean per-coordinate
    %xtr1 = bsxfun( @minus, xtr, mean( xtr, 1));
    xtr1 = xtr;
    
    %%% OLS on training data to learn the selected D indeces
    atom_ind = ols( y, xtr1, D-1);
    if numel( atom_ind) < D-1
        warning(sprintf('selected %d dimensions only.', numel( atom_ind) ));
        %pause();
    end
    
    atom_ind = [maxidx atom_ind];
    
    select_index2{icls}=atom_ind;
    %meta.select_index{icls} = atom_ind;

    % The new representation is the Q_tr so that Q_tr*R_tr = (X_tr|_selected)
    % thus, Q_te = (X_te|_selected) * inv(R_tr)
    [Q,R] = qr( xtr(:, atom_ind), 0); %if N > D, then Q is N-by-D, and R is D-by-D

   % alpha=Q' * y;
    Q_tr = Q;
    Q_cell_tr{icls}=Q_tr';
    %Feat_Tr( D*(icls-1)+1: D*icls, :) = Q_tr';
    %Feat_Te( D*(icls-1)+1: D*icls, :) = Q_te.';

end


k1=1;
k2=1;

% meta.select_index=zeros(numel(K),D);
% for icls=1:numel(K)
%     meta.select_index(icls,1:numel(select_index2{icls}))=select_index2{icls};
%      Feat_Tr(k1:(k1-1+size(Q_cell_tr{icls},1)), :) = Q_cell_tr{icls};
%     k1=k1+size(Q_cell_tr{icls},1);
%     k2=k2+size(Q_cell_tr{icls},1);
% end

for icls=1:numel(K)
    selidx(select_index2{icls}) = 1;
end
