function [lambda] = lambdaSequence(Xo,y,ds,type)
% LAMBDASEQUENCE Function used to compute lamda values such that the maximum
% lambda value produces all 0s for the Beta coefficients.
% [LAMBDA] = LAMBDASEQUENCE(X0,Y,DS,TYPE) return a vector LAMBDA containing
% lamda values generated for predictors X0, response variable Y, DS
% containing the data summary for the given TYPE ('classification' or
% 'regression') problem.

%    Copyright 2020-2021 The MathWorks, Inc.

% remove rows with NaNs in response
nanY = isnan(y);
if all(nanY)
     error(message('stats:featlearn:genfeatures:NoGoodYData'))
else
   Xo(nanY,:)   = [];
   y(nanY)     = [];
end

X = classreg.learning.internal.expandCategorical(Xo,ds.CategoricalPredictors,ds.OrdinalPredictors,ds.VariableRange);
[lambdaMax] = computeLambdaMax(X, y, type);
lambda = computeLambdaSequence(lambdaMax, 100, 1e-4);

end

function [lambdaMax] = computeLambdaMax(X, Y, type)
% lambdaMax is the penalty term (lambda) beyond which coefficients
% are guaranteed to be all zero: there is no benefit to calculating
% penalized fits with lambda > lambdaMax.

[N,~] = size(X);
if strcmp(type, 'classification')
    constantTerm = ones(length(Y),1);
    coeffs = glmfit(constantTerm,Y,'binomial','constant','off', ...
        'link','identity', 'weights',[]);
    predictedMu = glmval(coeffs,constantTerm,'identity','constant','off');
    muX = mean(X,1);
    X0 = X - muX;
    dotp = abs(X0' * (Y - predictedMu));
else % for regression
    muX = mean(X,1);
    X0 = X - muX;
    Y0 = Y - mean(Y);
    dotp = abs(X0' * Y0);
end

lambdaMax = max(dotp) / N;

end

function lambda = computeLambdaSequence(lambdaMax, nLambda, lambdaRatio)
% Fill in the log-spaced sequence of lambda values.        
% Fill in a number "nLambda" of smaller values, on a log scale.
lambdaMin = lambdaMax * lambdaRatio;
loghi = log(lambdaMax);
loglo = log(lambdaMin);
lambda = exp(linspace(loghi,loglo,nLambda));
lambda(end) = lambdaMin;
lambda = fliplr(lambda);

end