function result = glmnetPredict(object,s)

%--------------------------------------------------------------------------
% glmnetPredict.m: make predictions from a "glmnet" object.
%--------------------------------------------------------------------------
%
% DESCRIPTION:
%    Similar to other predict methods, this functions predicts fitted
%    values, logits, coefficients and more from a fitted "glmnet" object.
%
% USAGE: 
%    glmnetPredict(object)
%    glmnetPredict(object, type)
%    glmnetPredict(object, type, newx)
%    glmnetPredict(object, type, newx, s)
%
% INPUT ARGUMENTS:
% fit         Fitted "glmnet" model object.
% type        Type of prediction required. Type "link" gives the linear
%             predictors for "binomial" or "multinomial" models; for
%             "gaussian" models it gives the fitted values. Type "response"
%             gives the fitted probabilities for "binomial" or
%             "multinomial"; for "gaussian" type "response" is equivalent
%             to type "link". Type "coefficients" computes the coefficients
%             at the requested values for s. Note that for "binomial"
%             models, results are returned only for the class corresponding
%             to the second level of the factor response. Type "class"
%             applies only to "binomial" or "multinomial" models, and
%             produces the class label corresponding to the maximum
%             probability. Type "nonzero" returns a list of the indices of
%             the nonzero coefficients for each value of s.
% newx        Matrix of new values for x at which predictions are to be 
%             made. Must be a matrix; This argument is not used for 
%             type=c("coefficients","nonzero")
% s           Value(s) of the penalty parameter lambda at which predictions
%             are required. Default is the entire sequence used to create
%             the model.
%
% DETAILS:
%    The shape of the objects returned are different for "multinomial"
%    objects. glmnetCoef(fit, ...) is equivalent to glmnetPredict(fit, "coefficients", ...)
%
% 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:
%    glmnet, glmnetSet, glmnetPrint, glmnetPlot 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);
%    glmnetPredict(fit1,'link',x(1:5,:),[0.01,0.005]') % make predictions
%    glmnetPredict(fit1,'coefficients')
%    fit2=glmnet(x,g2,'binomial');
%    glmnetPredict(fit2, 'response', x(2:5,:))
%    glmnetPredict(fit2, 'nonzero')
%    fit3=glmnet(x,g4,'multinomial');
%    glmnetPredict(fit3, 'response', x(1:3,:), 0.01)
%
% DEVELOPMENT: 
%    14 Jul 2009: Original version of glmnet.m written.
%    20 Oct 2009: Fixed a bug in bionomial response, pointed out by Ramon
%                 Casanov from Wake Forest University.
%    26 Jan 2010: Fixed a bug in multinomial link and class, pointed out by
%                 Peter Rijnbeek from Erasmus University.
%    23 Jun 2010: Fixed a bug in multinomial with s, pointed out by 
%                 Robert Jacobsen from Aalborg University.

% if nargin < 2
%     type = 'link';
% end
% 
% if nargin < 3
%     newx = [];
% end
% 
% if nargin < 4
%     s = object.lambda;
% end

% if strcmp(object.class, 'elnet')
a0=transpose(object.a0);
nbeta=[a0; object.beta];
% if nargin == 4
lambda=object.lambda;
lamlist=lambda_interp(lambda,s);
nbeta=nbeta(:,lamlist.left).*repmat(lamlist.frac',size(nbeta,1),1) +nbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(nbeta,1),1));
% end
% if strcmp(type, 'coefficients')
result = nbeta;
% elseif strcmp(type, 'link')
%     result = [ones(size(newx,1),1), newx] * nbeta;
% elseif strcmp(type, 'response')
%     result = [ones(size(newx,1),1), newx] * nbeta;
% elseif strcmp(type, 'nonzero')
%     result = nonzeroCoef(nbeta(2:size(nbeta,1),:), true);
% else
%     error('Unrecognized type');
% end
% elseif strcmp(object.class, 'lognet')
%     
%     a0=transpose(object.a0);
%     nbeta=[object.a0; object.beta];
%     if nargin == 4    
%         lambda=object.lambda;
%         lamlist=lambda_interp(lambda,s);
%         nbeta=nbeta(:,lamlist.left).*repmat(lamlist.frac',size(nbeta,1),1) +nbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(nbeta,1),1));
%     end
%     %%%   remember that although the fortran lognet makes predictions
%     %%%   for the first class, we make predictions for the second class
%     %%%   to avoid confusion with 0/1 responses.
%     %%%   glmnet flipped the signs of the coefficients 
%     if strcmp(type,'coefficients')
%         result = nbeta;
%     elseif strcmp(type,'nonzero')
%         result = nonzeroCoef(nbeta(2:size(nbeta,1),:), true);
%     else
%         nfit = [ones(size(newx,1),1), newx] * nbeta;
%         
%         if strcmp(type,'response')
%             pp=exp(-nfit);
%             result = 1./(1+pp);
%         elseif strcmp(type,'link')
%             result = nfit;
%         elseif strcmp(type,'class')
%             result = (nfit > 0) * 2 + (nfit <= 0) * 1;
%         else
%             error('Unrecognized type');
%         end
%     end
% elseif strcmp(object.class, 'multnet')
%     a0=object.a0;
%     nbeta=object.beta;
%     nclass=size(a0,1);
%     nlambda=length(s);
%     if nargin == 4
%         lambda=object.lambda;
%         lamlist=lambda_interp(lambda,s);
%         for i=1:nclass
%             kbeta=[a0(i,:); nbeta{i}];
% %            kbeta=kbeta(:,lamlist.left)*lamlist.frac +kbeta(:,lamlist.right)*(1-lamlist.frac);
%             kbeta=kbeta(:,lamlist.left).*repmat(lamlist.frac',size(kbeta,1),1)+kbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(kbeta,1),1));
%             nbeta{i}=kbeta;
%         end
%     else
%         for i=1:nclass
%             nbeta{i} = [a0(i,:);nbeta{i}];
%         end
%     end
%     if strcmp(type, 'coefficients')
%         result = nbeta;
%     elseif strcmp(type, 'nonzero')
%         for i=1:nclass
%             result{i}=nonzeroCoef(nbeta{i}(2:size(nbeta{i},1),:),true);
%         end
%     else    
%         npred=size(newx,1);
%         dp = zeros(nclass,nlambda,npred);
%         for i=1:nclass
%             fitk = [ones(size(newx,1),1), newx] * nbeta{i};
%             dp(i,:,:)=dp(i,:,:)+reshape(transpose(fitk),1,nlambda,npred);
%         end
%         if strcmp(type, 'response')
%             pp=exp(dp);
%             psum=sum(pp,1);
%             result = permute(pp./repmat(psum,nclass,1),[3,1,2]);
%         elseif strcmp(type, 'link')
%             result=permute(dp,[3,1,2]);
%         elseif strcmp(type, 'class')
%             dp=permute(dp,[3,1,2]);
%             result = [];
%             for i=1:size(dp,3)
%                 result = [result, softmax(dp(:,:,i))];
%             end
%         else
%             error('Unrecognized type');
%         end
%     end
% else
%     error('Unrecognized class');
% end

%-------------------------------------------------------------
% End private function glmnetPredict
%-------------------------------------------------------------

function result = lambda_interp(lambda,s)
% lambda is the index sequence that is produced by the model
% s is the new vector at which evaluations are required.
% the value is a vector of left and right indices, and a vector of fractions.
% the new values are interpolated bewteen the two using the fraction
% Note: lambda decreases. you take:
% sfrac*left+(1-sfrac*right)
  
  if length(lambda)==1 % degenerate case of only one lambda
    nums=length(s);
    left=ones(nums,1);
    right=left;
    sfrac=ones(nums,1);
else
      s(s > max(lambda)) = max(lambda);
      s(s < min(lambda)) = min(lambda);
      k=length(lambda);
      sfrac =(lambda(1)-s)/(lambda(1) - lambda(k));
      lambda = (lambda(1) - lambda)/(lambda(1) - lambda(k));
      coord = interp1(lambda, 1:length(lambda), sfrac);
      left = floor(coord);
      right = ceil(coord);
      sfrac=(sfrac-lambda(right))./(lambda(left) - lambda(right));
      sfrac(left==right)=1;
  end
  result.left = left;
  result.right = right;
  result.frac = sfrac;

%-------------------------------------------------------------
% End private function lambda_interp
%-------------------------------------------------------------
%   
% function result = softmax(x, gap) 
% if nargin < 2
%     gap = false;
% end
% d = size(x);
% maxdist = x(:, 1);
% pclass = repmat(1, d(1), 1);
% for i =2:d(2)
%     l = x(:, i) > maxdist;
%     pclass(l) = i;
%     maxdist(l) = x(l, i);
% end
% if gap
%     x = abs(maxdist - x);
%     x(1:d(1), pclass) = x * repmat(1, d(2));
%     gaps = pmin(x);
% end
% if gap
%     result = {pclass, gaps};
% else
%     result = pclass;
% end

%-------------------------------------------------------------
% End private function softmax
%-------------------------------------------------------------
  