function [B, mu,stats] = CRDA_coefmat(X,y,covtype,mu)
% Computes the CRDA coefficinent matrix B = Sigma^-1 * M 
% 
% [B, mu] = CRDA_coefmat(X,y)
% [B, mu] = CRDA_coefmat(X,y,covtype,mu)
%  
% Input
%  X        data matrix p x N (p = # variables, N = # observations)
%  y        class matrix N x 1 (each element in [1,..,G], where G ist the
%           # of classes).
%  covtype  string either equal to 'ell2' (default), 'ell1' or 'riemann'
%  mu       matrix of p x G, of sample mean vectors   
%
% See also CRDACV and CRDA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p = size(X,1);
n  = histcounts(y);	% class labels counts
G = length(n);

if nargin < 3
    covtype = 'ell2';
end

if nargin < 4
    % Default: calculate the class means (mu = []);
    mu = [];
end

%% Center the data
if isempty(mu)
    % need to calculate the class means and center the data 
    
    % in case of Ell2-RSCM and Rie-PSCM, centering is w.r.t. to sample mean
    %if strcmpi(covtype,{'ell2','riemann'})
    I = eye(max(y));    
    Y = I(y, :);	% nxG indicator matrix => Yij = 1, i belong to group j
    mu = X*Y ./ n; 
    %else
    % in case of Ell1-RSCM centering is w.r.t. to spatial median

    
    for ii=1:G
        X(:,y==ii) = X(:,y==ii) - repmat(mu(:,ii),1,n(ii)); 
    end         
end

%% Calculate the pooled covariance matrix estimate 
stats=[];
switch covtype
    case 'ell2'
        [~,iCOV,stats] = regscm(X.','approach','ell2','inverse',true, ... 
            'centered',true,'verbose','off'); % pooled ELL2-RSCM   
    case 'ell1'
        [~,iCOV,stats] = regscm(X.','approach','ell1','inverse',true, ... 
            'centered',true,'verbose','off'); % pooled ELL1-RSCM   
    case 'riemann'
        [~,~,covest] = penscm(X','penalty','riemann','CV',5,'c','mean',...
                'verbose','off','preconly',true);   
        iCOV = covest.Theta_hat;
    otherwise
        error('sorry: wrong covariance matrix type');
end

%% Calculate the coefficient matrix
B = iCOV*mu;  




