glmnet.m

function fit = glmnet(x, y, options)
% Speed up version by Zhe Zhu (10/28/2014)
%--------------------------------------------------------------------------
% glmnet.m: fit an elasticnet model path
%--------------------------------------------------------------------------
%
% DESCRIPTION:
%    Fit a regularization path for the elasticnet at a grid of values for
%    the regularization parameter lambda. Can deal with all shapes of data.
%    Fits linear, logistic and multinomial regression models.
%
% USAGE:
%    fit = glmnet(x, y)
%    fit = glmnet(x, y, family, options)
%
% EXTERNAL FUNCTIONS:
% options         = glmnetSet;                  provided with glmnet.m
%
% INPUT ARGUMENTS:
% x           Input matrix, of dimension nobs x nvars; each row is an
%             observation vector. Currently sparse format is NOT supported.
% y           Response variable. Quantitative for family =
%             'gaussian'. For family = 'binomial' should be either a vector
%             of two levels, or a two-column matrix of counts or
%             proportions. For family = 'multinomial', can be either a
%             vector of nc>=2 levels, or a matrix with nc columns of counts
%             or proportions.
% family      Reponse type. (See above). Default is 'gaussian'.
% options     A structure that may be set and altered by glmnetSet (type
%             help glmnetSet).
%
% OUTPUT ARGUMENTS:
% fit         A structure.
% fit.a0      Intercept sequence of length length(fit.lambda).
% fit.beta    For "elnet" and "lognet" models, a nvars x length(lambda)
%             matrix of coefficients. For "multnet", a list of nc such
%             matrices, one for each class.
% fit.lambda  The actual sequence of lambda values used.
% fit.dev     The fraction of (null) deviance explained (for "elnet", this
%             is the R-square).
% fit.nulldev Null deviance (per observation).
% fit.df      The number of nonzero coefficients for each value of lambda.
%             For "multnet", this is the number of variables with a nonzero
%             coefficient for any class.
% fit.dfmat   For "multnet" only. A matrix consisting of the number of
%             nonzero coefficients per class.
% fit.dim     Dimension of coefficient matrix (ices).
% fit.npasses Total passes over the data summed over all lambda values.
% fit.jerr    Error flag, for warnings and errors (largely for internal
%             debugging).
% fit.class   Type of regression - internal usage.
%
% DETAILS:
%    The sequence of models implied by lambda is fit by coordinate descent.
%    For family = 'gaussian' this is the lasso sequence if alpha = 1, else
%    it is the elasticnet sequence. For family = 'binomial' or family =
%    "multinomial", this is a lasso or elasticnet regularization path for
%    fitting the linear logistic or multinomial logistic regression paths.
%    Sometimes the sequence is truncated before options.nlambda values of
%    lambda have been used, because of instabilities in the logistic or
%    multinomial models near a saturated fit. glmnet(..., family =
%    'binomial') fits a traditional logistic regression model for the
%    log-odds. glmnet(..., family = 'multinomial') fits a symmetric
%    multinomial model, where each class is represented by a linear model
%    (on the log-scale). The penalties take care of redundancies. A
%    two-class "multinomial" model will produce the same fit as the
%    corresponding "binomial" model, except the pair of coefficient
%    matrices will be equal in magnitude and opposite in sign, and half the
%    "binomial" values. Note that the objective function for
%    "gaussian" is
%                1 / (2 * nobs) RSS + lambda * penalty
%    , and for the logistic models it is
%                1 / nobs - loglik + lambda * penalty
%
% LICENSE: GPL-2
%
% DATE: 14 Jul 2009
%
% AUTHORS:
%    Algorithm was designed by Jerome Friedman, Trevor Hastie and Rob Tibshirani
%    Fortran code was written by Jerome Friedman
%    R wrapper (from which the MATLAB wrapper was adapted) was written by Trevor Hasite
%    MATLAB wrapper was written and maintained by Hui Jiang, jiangh@stanford.edu
%    Department of Statistics, Stanford University, Stanford, California, USA.
%
% REFERENCES:
%    Friedman, J., Hastie, T. and Tibshirani, R. (2009)
%    Regularization Paths for Generalized Linear Models via Coordinate Descent.
%    Journal of Statistical Software, 33(1), 2010
%
% SEE ALSO:
%    glmnetSet, glmnetPrint, glmnetPlot, glmnetPredict and glmnetCoef methods.
%
% EXAMPLES:
%    x=randn(100,20);
%    y=randn(100,1);
%    g2=randsample(2,100,true);
%    g4=randsample(4,100,true);
%    fit1=glmnet(x,y);
%    glmnetPrint(fit1);
%    glmnetCoef(fit1,0.01) % extract coefficients at a single value of lambda
%    glmnetPredict(fit1,'response',x(1:10,:),[0.01,0.005]') % make predictions
%    fit2=glmnet(x,g2,'binomial');
%    fit3=glmnet(x,g4,'multinomial');
%
% DEVELOPMENT:
%    14 Jul 2009: Original version of glmnet.m written.
%    26 Jan 2010: Fixed a bug in the description of y, pointed out by
%                 Peter Rijnbeek from Erasmus University.
%    09 Mar 2010: Fixed a bug of printing "ka = 2", pointed out by
%                 Ramon Casanova from Wake Forest University.
%    25 Mar 2010: Fixed a bug when p > n in multinomial fitting, pointed
%                 out by Gerald Quon from University of Toronto
%    25 Jul 2010: Check for input matrix format and size
%    27 Sep 2010: Fixed a bug of undefined "df" in multinomial fitting,
%                 pointed by Jeff Howbert from Insilicos.
%    27 Jan 2013  Faster version made especially for CCDC (Zhe Zhu)

%  % Check input arguments
%  if nargin < 2
%      error('more input arguments needed.');
%  end
%
%  if nargin < 3
%      family = 'gaussian';
%  end
%
%  if nargin < 4
%      options = glmnetSet;
%  end

% fixed family
% family = 'gaussian';
% fixed options
% options = glmnetSet;

% Prepare parameters
% nlam = options.nlambda;

%  if (issparse(x) || issparse(y))
%      error('currently sparse matrix is NOT supported.');
%  end
%
%  if (~isa(x,'double') || ~isa(y,'double'))
%      error('only DOUBLE precision matrix is supported.');
%  end

[nobs,nvars] = size(x);

%  if (nobs <= 1)
%      error('at least two observations should be provided.');
%  end

weights = options.weights;
if isempty(weights)
    weights = ones(nobs,1);
end

% maxit = options.maxit;

%  if strcmp(family, 'binomial') || strcmp(family, 'multinomial')
%      [noo,nc] = size(y);
%      kopt = double(options.HessianExact);
%      if noo ~= nobs
%          error('x and y have different number of rows');
%      end
%      if nc == 1
%          classes = unique(y);
%          nc = length(classes);
%          indexes = eye(nc);
%          y = indexes(y,:);
%      end
%      if strcmp(family, 'binomial')
%          if nc > 2
%              error ('More than two classes; use multinomial family instead');
%          end
%          nc = 1; % for calling multinet
%      end
%      if ~isempty(weights)
%          % check if any are zero
%          o = weights > 0;
%          if ~all(o) %subset the data
%              y = y(o,:);
%              x = x(o,:);
%              weights = weights(o);
%              nobs = sum(o);
%          end
%          [my,ny] = size(y);
%          y = y .* repmat(weights,1,ny);
%      end
%      % Compute the null deviance
%      prior = sum(y,1);
%      sumw = sum(sum(y));
%      prior = prior / sumw;
%      nulldev = -2 * sum(sum(y .* (ones(nobs, 1) * log(prior)))) / sumw;
%  elseif strcmp(family, 'gaussian')
% Compute the null deviance
ybar = y' * weights/ sum(weights);
nulldev = (y' - ybar).^2 * weights / sum(weights);
%      if strcmp(options.type, 'covariance')
%          ka = 1;
%      elseif strcmp(options.type, 'naive')
%          ka = 2;
%      else
%          error('unrecognized type');
%      end
ka = 1;
%  else
%      error('unrecognized family');
%  end

% ne = options.dfmax;
%  if ne == 0
%      ne = nvars + 1;
%  end
ne = nvars + 1;

% nx = options.pmax;
%  if nx == 0
%      nx = min(ne * 1.2, nvars);
%  end
nx = min(ne * 1.2, nvars);

% exclude = options.exclude;
%  if ~isempty(exclude)
%      exclude = unique(exclude);
%  %      if ~all(exclude > 0 & exclude <= nvars)
%  %          error('Some excluded variables out of range');
%  %      end
%      jd = [length(exclude); exclude];
%  else
%      jd = 0;
%  end
jd = 0;

% vp = options.penalty_factor;
%  if isempty(vp)
%      vp = ones(nvars,1);
%  end
vp = ones(nvars,1);

isd = double(options.standardize);
thresh = options.thresh;
lambda = options.lambda;
% lambda_min = options.lambda_min;
% %  if lambda_min == 0
% if nobs < nvars
%     lambda_min = 5e-2;
% else
%     lambda_min = 1e-4;
% end
%  end
%  if isempty(lambda)
%      if (lambda_min >= 1)
%          error('lambda_min should be less than 1');
%      end
%      flmin = lambda_min;
%      ulam = 0;
%  else
flmin = 1.0;
%      if any(lambda < 0)
%          error ('lambdas should be non-negative');
%      end
ulam = -sort(-lambda);
nlam = length(lambda);
% end

parm = options.alpha;

%  if strcmp(family, 'gaussian')
[a0,ca,ia,nin,rsq,alm,nlp,jerr] = glmnetMex(parm,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,weights,ka);
%  else
%      [a0,ca,ia,nin,dev,alm,nlp,jerr] = glmnetMex(parm,x,y,jd,vp,ne,nx,nlam,flmin,ulam,thresh,isd,nc,maxit,kopt);
%  end

% Prepare output
lmu = length(alm);
ninmax = max(nin);
lam = alm;
if isempty(options.lambda)
    lam = fix_lam(lam); % first lambda is infinity; changed to entry point
end
% errmsg = err(jerr, maxit, nx);
%  if errmsg.n == 1
%      error(errmsg.msg);
%  elseif errmsg.n == -1
%      warning(errmsg.msg);
%  end

%  if strcmp(family, 'multinomial')
%      beta_list = {};
%      a0 = a0 - repmat(mean(a0), nc, 1);
%      dfmat=a0;
%      dd=[nvars, lmu];
%      if ninmax > 0
%          ca = reshape(ca, nx, nc, lmu);
%          ca = ca(1:ninmax,:,:);
%          ja = ia(1:ninmax);
%          [ja1,oja] = sort(ja);
%          df = any(abs(ca) > 0, 2);
%          df = sum(df, 1);
%          df = df(:);
%          for k=1:nc
%              ca1 = reshape(ca(:,k,:), ninmax, lmu);
%              cak = ca1(oja,:);
%              dfmat(k,:) = sum(sum(abs(cak) > 0));
%              beta = zeros(nvars, lmu);
%              beta(ja1,:) = cak;
%              beta_list{k} = beta;
%          end
%      else
%          for k = 1:nc
%              dfmat(k,:) = zeros(1,lmu);
%              beta_list{k} = zeros(nvars, lmu);
%          end
%          df = zeros(1,lmu);
%      end
%      fit.a0 = a0;
%      fit.beta = beta_list;
%      fit.dev = dev;
%      fit.nulldev = nulldev;
%      fit.dfmat = dfmat;
%      fit.df = df';
%      fit.lambda = lam;
%      fit.npasses = nlp;
%      fit.jerr = jerr;
%      fit.dim = dd;
%      fit.class = 'multnet';
%  else
dd=[nvars, lmu];
if ninmax > 0
    ca = ca(1:ninmax,:);
    df = sum(abs(ca) > 0, 1);
    ja = ia(1:ninmax);
    [ja1,oja] = sort(ja);
    beta = zeros(nvars, lmu);
    beta (ja1, :) = ca(oja,:);
else
    beta = zeros(nvars,lmu);
    df = zeros(1,lmu);
end

% if strcmp(family, 'binomial')
%     a0 = -a0;
%     fit.a0 = a0;
%     fit.beta = -beta; %sign flips make 2 arget class
%     fit.dev = dev;
%     fit.nulldev = nulldev;
%     fit.df = df';
%     fit.lambda = lam;
%     fit.npasses = nlp;
%     fit.jerr = jerr;
%     fit.dim = dd;
%     fit.class = 'lognet';
% else
fit.a0 = a0;
fit.beta = beta;
fit.dev = rsq;
fit.nulldev = nulldev;
fit.df = df';
fit.lambda = lam;
fit.npasses = nlp;
fit.jerr = jerr;
fit.dim = dd;
fit.class = 'elnet';
% end
%  end

%------------------------------------------------------------------
% End function glmnet
%------------------------------------------------------------------

function new_lam = fix_lam(lam)

new_lam = lam;
llam=log(lam);
new_lam(1)=exp(2*llam(2)-llam(3));

%------------------------------------------------------------------
% End private function fix_lam
%------------------------------------------------------------------

% function output = err(n,maxit,pmax)
% 
% if n==0
%     output.n=0;
%     output.msg='';
% elseif n>0 %fatal error
%     if n<7777
%         msg='Memory allocation error; contact package maintainer';
%     elseif n==7777
%         msg='All used predictors have zero variance';
%     elseif (8000<n) && (n<9000)
%         msg=sprintf('Null probability for class %d < 1.0e-5', n-8000);
%     elseif (9000<n) && (n<10000)
%         msg=sprintf('Null probability for class %d > 1.0 - 1.0e-5', n-9000);
%     elseif n==10000
%         msg='All penalty factors are <= 0';
%     end
%     output.n=1;
%     output.msg=['in glmnet fortran code - %s',msg];
% elseif n<0 %non fatal error
%     if n > -10000
%         msg=sprintf('Convergence for %dth lambda value not reached after maxit=%d iterations; solutions for larger lambdas returned', -n, maxit);
%     elseif n < -10000
%         msg=sprintf('Number of nonzero coefficients along the path exceeds pmax=%d at %dth lambda value; solutions for larger lambdas returned', pmax, -n-10000);
%     end
%     output.n=-1;
%     output.msg=['from glmnet fortran code - ',msg];
% end
% 
% %------------------------------------------------------------------
% % End private function err
% %------------------------------------------------------------------