function [a_bst,P_bst,covest,eta_bst,etas,ind] = penscm(X,varargin)
% PENSCM Computes different penalized (shrinkage) sample covariance matrix 
% (PSCM) estimators given an n x p data matrix X (rows are observations).
% The method computes the solution (minimizer) of the penalized negative 
% Gaussian likelihood function that uses an additive orthogonally
% invariant penalty. Current choises are  Riemannian distance or the 
% scaled Riemannian distance. 
%
% [a_bst,P_bst,covest,eta_bst,ind] = PENSCM(X,varargin)
%
% PENSCM can be called with numerous optional arguments. Optional
% arguments are given in parameter pairs, so that first argument is
% the name of the parameter and the next argument is the value for
% that parameter. Optional parameter pairs can be given in any order.
% 
% Positional parameters:
%
%      X                A numeric matrix (dimension,nxp), where rows are
%                       observations
%
% Optional Parameters:
%
% Parameter name        Values and description
%==========================================================================
% --Basic parameter is the choise of the estimator to use
%
% 'penalty'             (string) which orthogonally invariant penalty is to 
%                       be used. Current choises are: 
%                        's-riemann'  to use scaled Riemannian penalty
%                        'riemann'    to use Riemannian penalty                         
%==========================================================================
% -- Optional parameter values
% 
% 'eta'                 (nonnegative vector) regularization parameter
%                       values. If not given, then default set of parameter 
%                       values are calculated by PENSCM
% 'CV   '               (positive integer) when CV is positive interger K, 
%                       then K-fold crossvalidation i used. 
% 'verbose'             (string) Either 'on' (default) or 'off'. When 'on' 
%                       then the algorithm prints the progress of the 
%                       function in text format.
% 'c'                   (positive scalar) The covariance matrix estimator is 
%                       shrinked  towards c*I, where c>0. This option is
%                       used only when using non-scaled penalty such as
%                       Riemannian distance. The default value is c=1. When
%                       'c' is 'mean' the method uses the mean of the
%                       eigenvalues of the sample covariance matrix
%                       (=Trace(S)/p). 
% 'preconly'            logical (default:false). If true, then compute only 
%                       the precesion matrix.
% 
% -- Return parameters:
%
%  a_bst, P_bst         (p x 1 vector) the log of the eigenvalues and 
%                       orthogonal  p xp matrix  of the eigenvectors
%                       corresponding to the solution usin the penalty 
%                       parameter value eta_bst found by cross-validation
%  P_bst                orthogonal p x p matrix of eigenvectors of the
%                       solution that uses the 
%  covest               struct that contains the actual covariance matrix  
%                       and precesion matrix estimates with fields
%   'Sigma_hat'         covariance matrix estimator = P*diag(exp(a_bst))*P' 
%   'Theta_hat'         precision matrix estimator = P*diag(exp(-a_bst))*P' 
%  eta_bst              best penalty parameter found by cross-validation 
%  etas                 array of the used penalty parameters
%  ind                  index of the the best penalty on the grid
%
% See also PENSCM_SRIE0 PENSCM_RIE0
%
%% Details and references:
%         
%  The algorithm uses the NR algorithm and the CV approach described in 
% 
%   Tyler, David E., and Mengxi Yi. "Shrinking the Sample Covariance Matrix 
%   using Convex Penalties on the Matrix-Log Transformation." 
%   ArXiv preprint arXiv:1903.08281 (2019).
% 
% toolbox: SCATTER
% authors: Copyright by Esa Ollila and Michael Muma, 2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-- Check requirements on the input data matrix X
argin = inputParser;
valX = @(X) assert(isreal(X) && (numel(size (X))==2),['''X'' must be a ' ...
    'real-valued matrix having at least two columns (variables)']);
addRequired(argin,'X',valX);

if nargin == 0
  error ('You must supply the data matrix  as input argument.');
end

if any (isnan (X))
  error ('Input data contains NaN''s.');
end

if ~isa (X, 'double')
  X = double(X);
end

[n,p] = size(X); % data matrix dimensions

%-- optional input parsing rules

valCV = @(x) assert(x > 0 && isreal(x) && (x-round(x)==0),['''CV'' must be a' ...
    ' positive integer']);
addParameter(argin,'CV',5, valCV);

valVerbose = @(x) assert(any(strcmpi(x,{'on','off'})),['''verbose'' must be ' ...
    'string equal to ''on'' or ''off''']);
addParameter(argin,'verbose','on',valVerbose);

valEtas = @(x) assert( numel(x>0)==length(x) && isreal(x) && isvector(X), ...
    ['''etas'' must be a real-valued vector of positive reals']);
addParameter(argin,'etas',[], valEtas);

addParameter(argin,'penalty','riemann', ... 
    @(x) logical(sum(strcmp(x,{'riemann','s-riemann'}))));

addParameter(argin,'c',1);

addParameter(argin,'preconly',false, @(x) assert(islogical(x),['input ' ...
    ' ''prec_only'' must be logical (false/true)']));

% must be riemann or s-riemann

%---  parse inputs
parse(argin,X,varargin{:});
print_info = strcmpi(argin.Results.verbose,'on');
penalty = argin.Results.penalty;
CV = argin.Results.CV;
etas = argin.Results.etas;
c = argin.Results.c;


if isempty(etas)
    %% Set-up the grid of eta values
    betas = fliplr([0.01 [1:6 8]/10]); % small beta --> large eta (for PSCM)
    etas = (1-betas)./betas; % eta grid for PSCM
%    betas = fliplr([0.01 (1:8)/10]); % small beta --> large eta
%    etas = (1-betas)./betas;
    % NOTE: If eta is very small and n < p, then the NR algorithm may have
    % convergence problems
end

if (isa(c,'char') && ~strcmpi(c,'mean'))
   error(message('stats:RSCM:Invalid_c_type'));
elseif  isreal(c) 
    if c < 0
        error(message('stats:RSCM:negative_c'));
    end
end

%-- print information about data
if print_info
  fprintf('Number of variables           : %d\n', p);
  fprintf('Number of samples             : %d\n', n);
  fprintf('Computing the estimator with %d-fold cross validation\n',CV);
  fprintf('Using %s penalty\n',penalty);   
  if strcmpi(penalty,'riemann') 
      if strcmpi(c,'mean')
        fprintf('Shrinkage towards the mean of the eigenvalues\n');
      else
        fprintf('Shrinkage towards %.1f times the identity matrix\n',c);
      end
  end
end
   

%-- Now compute the optimal shrinkage parameter value and the estimator
% using cross validation
      
CVO = cvpartition(n,'KFold',CV);
err = zeros(length(etas),CV);
      
for i = 1:CVO.NumTestSets
       trIdx = CVO.training(i);
       teIdx = CVO.test(i);
       if strcmpi(penalty,'s-riemann')
           [~,~,stats] = penscm_srie0(X(trIdx,:),X(teIdx,:),etas,false,false); 
       else              
           [~,~,stats] = penscm_rie0(X(trIdx,:),X(teIdx,:),c,etas,false,false);
       end
       err(:,i) = stats.lik_eta;
       if print_info
           fprintf('  fold %d done\n',i); 
       end
 end
 etas = stats.etas;     
 [~,ind] = min(mean(err,2));
 eta_bst = etas(ind);
 
 %-- Then compute the final estimator based on the eta value chosen by 
 % cross validation
 if strcmpi(penalty,'s-riemann')
    [a_bst,P_bst,~,covest] = penscm_srie0(X,[],eta_bst,false,true,true);
 else              
    [a_bst,P_bst,~,covest] = penscm_rie0(X,[],c,eta_bst,false,true,true);
 end
     










