function [ent,re] = wentropy(x,varargin)
%Wavelet entropy
% Private MATLAB code generation function

%#codegen
%   Copyright 2022 The MathWorks, Inc.
% Validate input
validFormat = (iscell(x) || isnumeric(x)) && ~isstruct(x);
coder.internal.assert(validFormat,'Wavelet:entropy:MustBeCellOrNumeric');
IsCellInput = iscell(x);
coder.internal.errorIf(IsCellInput && isempty(x),'Wavelet:entropy:EmptyCell');
if ~IsCellInput
    validateattributes(x,{'double','single'},{'real','finite','nonempty','2d'},'wentropy','x');
else
    for ii = 1:length(x)
        validateattributes(x{ii},{'double','single'},{'real','finite','nonempty','vector'},'wentropy','x');
    end
end
coder.internal.prefer_const(varargin{:});
parms = struct( ...
    'Wavelet', uint32(0),...
    'Level',uint32(0),...
    'Entropy', uint32(0),...
    'Transform',uint32(0), ...
    'Distribution',uint32(0), ...
    'Exponent',uint32(0),...
    'EnergyThreshold',uint32(0),...
    'Scaled',uint32(0));
% Select parsing options.
poptions = struct( ...
    'CaseSensitivity',false, ...
    'PartialMatching','unique', ...
    'StructExpand',false, ...
    'IgnoreNulls',true, ...
    'SupportOverrides',false);
% Parse the inputs.
pstruct = coder.internal.parseParameterInputs(parms,poptions,varargin{:});
% Retrieve parameter values.
tmpEntropy = coder.internal.getParameterValue(pstruct.Entropy,"Shannon",varargin{:});
tmpTransform = coder.const(coder.internal.getParameterValue(pstruct.Transform,'',varargin{:}));
% Validate the entropy
pEntropy = validatestring(tmpEntropy,{'Shannon','Renyi','Tsallis'},'wentropy','Entropy');
% Normalization
tmpDistribution = coder.const(coder.internal.getParameterValue(pstruct.Distribution,"scale",varargin{:}));
pNormalization = coder.const(validatestring(tmpDistribution,{'global','scale'},'wentropy','Distribution'));
% wavelet
tmpWavelet = coder.const(coder.internal.getParameterValue(pstruct.Wavelet,coder.const(''),varargin{:}));
% Initial exponent must be real, finite, and scalar.
inExponent = coder.internal.getParameterValue(pstruct.Exponent,zeros(0,1),varargin{:});
tmpEnergyThreshold = coder.internal.getParameterValue(pstruct.EnergyThreshold,1e-8,varargin{:});
pScaled = coder.internal.getParameterValue(pstruct.Scaled,true,varargin{:});
validateattributes(pScaled,{'double','single','logical'}, {'scalar',...
    'real','finite'},'wentropy','Scaled');
% We will cast to match data type
validateattributes(tmpEnergyThreshold,{'double','single'}, {'scalar',...
    'nonnegative'},'wentropy','EnergyThreshold');
tmpLevel = coder.internal.getParameterValue(pstruct.Level,[],varargin{:});
if ~isempty(tmpLevel)
    validateattributes(tmpLevel,{'double','single'},{'scalar','positive'},'wentropy','Level');
end
if IsCellInput
    dataType = coder.const(underlyingType(x{1}));    
else
    dataType = coder.const(underlyingType(x));    
end
coder.internal.errorIf(~isempty(inExponent) && startsWith(pEntropy,'S'),...
    'Wavelet:entropy:expNoShannon');
if ~isempty(inExponent)
    validateattributes(inExponent,{'double','single'},{'scalar','real','finite'},'wentropy','Exponent');
    tmpExponent = cast(inExponent(1,1),dataType);
else
    tmpExponent = Inf(1,1,dataType);
end
if ~isfinite(tmpExponent) && (startsWith(pEntropy,'R') || startsWith(pEntropy,'T'))
    defaultExponent = cast(2,dataType);
elseif ~isfinite(tmpExponent) && startsWith(pEntropy,'S')
    defaultExponent = ones(1,1,dataType);
else
    defaultExponent = tmpExponent;
end
% Make exponent and level in the data type of the input
% Handle limiting case when exponent goes to 1.
if ~isempty(defaultExponent) && ~startsWith(pEntropy,'S')
    [entmethod,expvalue] = validateExponent(defaultExponent,pEntropy,dataType);
else
    entmethod = 'Shannon';
    expvalue = cast(1,dataType);
end
% Is the input time series data
tsData = (coder.internal.isConst(size(x,1)) || coder.internal.isConst(size(x,2))) ...
    && (size(x,1) == 1 || size(x,2) == 1) && ~IsCellInput;
coder.internal.errorIf(tsData && numel(x) < 4,'Wavelet:entropy:FourSamples');
coder.internal.errorIf(~tsData && ~IsCellInput && size(x,1) < 2,'Wavelet:entropy:TwoRows');
thresh = cast(tmpEnergyThreshold,dataType);
coder.internal.errorIf(~tsData && (~isempty(tmpTransform) || ~isempty(tmpWavelet)...
    || ~isempty(tmpLevel)), 'Wavelet:entropy:InvalidWithTransform');
if isempty(tmpTransform) && tsData
    pTransform = coder.const('modwt');
elseif ~isempty(tmpTransform) && ~IsCellInput
    pTransform = coder.const(validatestring(tmpTransform,{'dwt','dwpt','modwt','modwpt'},'wentropy','Transform'));    
else
    pTransform = coder.const(tmpTransform);
end
%Common orientation
if tsData && isrow(x)
    xcol = reshape(x,numel(x),1);
else
    xcol = x;
end
% length of data, or cell array, or number of rows
N = size(xcol,1);
pLevel = setEntropyLevel(N,tmpLevel,pTransform,dataType);

if isempty(tmpWavelet) && tsData
    switch pTransform
        case {'modwpt','dwpt'}
            pWavelet = coder.const('fk18');
        otherwise
            pWavelet = coder.const('sym4');
    end
else
    pWavelet = tmpWavelet;
end
if tsData
    wt = getEntropyTransform(xcol,N,pTransform,pWavelet,pLevel);
else
    wt = xcol;
end
IsCellTransform = iscell(wt);
if IsCellTransform
    wtcol = transposeCell(wt,dataType);    
end
ProbCell = {};
Prob = cast([],dataType);
Nest = 0;
if IsCellTransform && startsWith(pNormalization,'s')
    [ProbCell,Nest] = getProbCellScale(wtcol,thresh,dataType);
elseif IsCellTransform && startsWith(pNormalization,'g')
    [Prob,Nest] = getProbCellGlobal(wtcol,thresh,dataType);
elseif ~IsCellTransform
    [Prob,Nest] = getProbabilityNumeric(wt,pNormalization,thresh);
end
if IsCellTransform && startsWith(pNormalization,'s')
    ent = getEntropyCell(ProbCell,entmethod,pScaled,expvalue,Nest,dataType);
elseif IsCellTransform && startsWith(pNormalization,'g')
    ent = getEntropyNumeric(Prob,entmethod,pScaled,expvalue,Nest,dataType);
else
    ent = getEntropyNumeric(Prob,entmethod,pScaled,expvalue,Nest,dataType);
end
if IsCellTransform && startsWith(pNormalization,'s')
    re = ProbCell;
else
    re = Prob;
end
%-------------------------------------------------------------------------
function se = getEntropyNumeric(Pij,entmethod,scaledTF,expvalue,N,datatype)

if startsWith(entmethod,'S')
    se = shannonEntropyNumeric(Pij,scaledTF,N,datatype);
elseif startsWith(entmethod,'R')
    se = renyiEntropyNumeric(Pij,expvalue,scaledTF,N,datatype);
else
    se = tsallisEntropyNumeric(Pij,expvalue,scaledTF,N,datatype);
end
%--------------------------------------------------------------------------
function se = getEntropyCell(Pij,entmethod,scaledTF,expvalue,N,datatype)

if startsWith(entmethod,'S')
    se = shannonEntropyCell(Pij,scaledTF,N,datatype);
elseif startsWith(entmethod,'R')
    se = renyiEntropyCell(Pij,expvalue,scaledTF,N,datatype);
else
    se = tsallisEntropyCell(Pij,expvalue,scaledTF,N,datatype);
end
%--------------------------------------------------------------------------
function se = shannonEntropyNumeric(Pij,scaledTF,N,datatype)
if size(Pij,2) == 1
    se = squeeze(-sum(Pij.*log(Pij),1,'omitnan'));
    se = se';
    se(se<eps(datatype)) = zeros(1,1,datatype);
else
    se = squeeze(-sum(Pij.*log(Pij),2,'omitnan'));
    se(se < eps(datatype)) = zeros(1,1,datatype);
end
se(~isfinite(se)) = zeros(1,1,datatype);
if scaledTF
    se = se./log(N);
    se(N==1) = zeros(1,1,datatype);
end
if isrow(se)
    se = reshape(se,[],1);
end
%--------------------------------------------------------------------------
function re = renyiEntropyNumeric(Pij,alpha,scaledTF,N,datatype)
if size(Pij,2) == 1
    re = 1/(1-alpha)*log(sum(Pij.^alpha,1,'omitnan'));
    re = reshape(re,1,[]);
else
    re = 1/(1-alpha)*log(sum(Pij.^alpha,2,'omitnan'));
end

re(~isfinite(re)) = zeros(1,1,datatype);
if scaledTF
    re = re./log(N);
    re(N==1) = zeros(1,1,'like',re);
end
if isrow(re)
    re = reshape(re,[],1);
end
%--------------------------------------------------------------------------
function te = tsallisEntropyNumeric(Pij,q,scaledTF,N,datatype)
qMinus1 = q-1;
if size(Pij,2) == 1
    te = 1/qMinus1*sum(Pij-Pij.^q,1,'omitnan');
    te = reshape(te,1,[]);
else
    te = 1/qMinus1*sum(Pij-Pij.^q,2,'omitnan');
end
normFactor = qMinus1./(1-N.^-qMinus1);
if scaledTF
    te = normFactor.*te;
    te(N==1) = zeros(1,1,datatype);
end
if isrow(te)
    te = reshape(te,[],1);
end
%--------------------------------------------------------------------------
function [entmethod,expvalue] = validateExponent(exponent,entropy,type)
entmethod = 'Shannon';
expvalue = ones(1,1,type);
tmpexp = cast(exponent,type);
switch entropy
    case 'Renyi'
        validateattributes(exponent,{'numeric'},...
            {'scalar','real','nonnegative','nonempty','finite'},'WENTROPY','Exponent');
        if abs(expvalue - tmpexp) > eps(type)
            expvalue = tmpexp;
            entmethod = 'Renyi';
        end
    case 'Tsallis'
        validateattributes(exponent,{'numeric'},...
            {'scalar','real', 'nonempty','finite', '>=',cast(-1/2,type)},'WENTROPY','Exponent');
        if abs(expvalue - tmpexp) > eps(type)
            expvalue = tmpexp;
            entmethod = 'Tsallis';
        end
end
%-------------------------------------------------------------------------
function [Pij,N] = getProbabilityNumeric(cfs,pNormalization,energyThresh)
if startsWith(pNormalization,'s')
    Energy = sum(abs(cfs).*abs(cfs),2);
    zeroEnergy = Energy < energyThresh;
    cfs(zeroEnergy,:) = zeros(1,1,'like',energyThresh);
    Energy(zeroEnergy) = ones(1,1,'like',energyThresh);
    Pij = (abs(cfs).*abs(cfs))./Energy;
    N = size(Pij,2);
else
    TimeEnergy = sum(abs(cfs).*abs(cfs),2);
    globalEnergy = sum(TimeEnergy,1);
    if globalEnergy < energyThresh
        cfs = zeros(size(cfs),'like',energyThresh);
        globalEnergy = ones(1,1,'like',energyThresh);
    end
    Pij = sum(abs(cfs).*abs(cfs)./globalEnergy,2);
    if isrow(Pij)
        Pij = Pij(:);
    end
    N = size(Pij,1);
end
%-------------------------------------------------------------------------
function tf = mustBeValidCell(x)
N = length(x);
tf = true;
for ii = 1:N
    coder.internal.errorIf(~isvector(x{ii}) || numel(x{ii}) < 4,'Wavelet:entropy:FourSamples');
    if any(~isfinite(x{ii}))
        tf = false;
    end
    if isempty(x{ii})
        tf = false;
    end
    if any(imag(x{ii}))
        tf = false;
    end
    if ~strcmp(underlyingType(x{ii}),'double') && ~strcmp(underlyingType(x{ii}),'single')
        tf = false;
    end
end
%-------------------------------------------------------------------------
function tf = mustBeValidNumeric(x)
tf = true;
if any(~isfinite(x))
    tf  = false;
end
if isempty(x)
    tf = false;
end
if any(imag(x))
    tf = false;
end
dataType = underlyingType(x);
if ~strcmp(dataType,'double') && ~strcmp(dataType,'single')
    tf = false;
end
%-------------------------------------------------------------------------
function wt = transposeCell(x,dataType)
N = length(x);
if isrow(x)
    wt = repmat({zeros(0,1,dataType)},1,N);
else
    wt = repmat({zeros(0,1,dataType)},N,1);
end
for ii = 1:N
    if isrow(x{ii})
        wt{ii} = x{ii}';
    else
        wt{ii} = x{ii};
    end
end
%-------------------------------------------------------------------------
function [Pij,N] = getProbCellScale(wt,thresh,dataType)
Nt = length(wt);
if isrow(wt)
    Pij = repmat({zeros(0,1,dataType)},1,Nt);
else
    Pij = repmat({zeros(0,1,dataType)},Nt,1);
end
energy = ones(Nt,1);
N = ones(Nt,1);
for ii = 1:Nt
    energy(ii) = norm(wt{ii},2)*norm(wt{ii},2);
    if energy(ii) < thresh
        wt{ii} = zeros(1,1,dataType);
        energy(ii) = ones(1,1,dataType);
    end
    Pij{ii} = abs(wt{ii}.*wt{ii})./energy(ii);
    N(ii) = numel(Pij{ii});
end
%-------------------------------------------------------------------------
function [Pij,N] = getProbCellGlobal(wt,thresh,dataType)
N = length(wt);
Pij = zeros(N,1,dataType);
absSqSum = zeros(N,1,dataType);
globalEnergy = zeros(1,1,dataType);
for ii = 1:N
    absSqSum(ii) = sum(wt{ii}.*wt{ii});
    globalEnergy = globalEnergy+absSqSum(ii);
end
if globalEnergy > thresh
    Pij = absSqSum./globalEnergy;    
end










