% function x = logistic(a, y, w, ridge, param)
%
% Logistic regression.  Design matrix A, targets Y, optional instance
% weights W, optional ridge term RIDGE, optional parameters object PARAM.
%
% W is a vector with length equal to the number of training examples; RIDGE
% can be either a vector with length equal to the number of regressors, or
% a scalar (the latter being synonymous to a vector with all entries the
% same).
%
% PARAM has fields PARAM.MAXITER (an iteration limit), PARAM.VERBOSE
% (whether to print diagnostic information), PARAM.EPSILON (used to test
% convergence), and PARAM.MAXPRINT (how many regression coefficients to
% print if VERBOSE==1).
%
% Model is 
%
%   E(Y) = 1 ./ (1+exp(-A*X))
%
% Outputs are regression coefficients X.
%
% Copyright 2007 Geoffrey J. Gordon
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or (at 
% your option) any later version. 
%
% This program is distributed in the hope that it will be useful, but 
% WITHOUT ANY WARRANTY; without even the implied warranty of 
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function x = logistic(a, y, w, ridge, param)

% process parameters
[n, m] = size(a);

if ((nargin < 3) || (isempty(w)))
  w = ones(n, 1);
end

if ((nargin < 4) || (isempty(ridge)))
  ridge = 1e-5;
end

if (nargin < 5)
  param = [];
end

if (length(ridge) == 1)
    ridgemat = speye(m) * ridge;
elseif (length(ridge(:)) == m)
    ridgemat = spdiags(ridge(:), 0, m, m);
else
    error('ridge weight vector should be length 1 or %d', m);
end

if (~isfield(param, 'maxiter'))
  param.maxiter = 200;
end

if (~isfield(param, 'verbose'))
  param.verbose = 0;
end

if (~isfield(param, 'epsilon'))
  param.epsilon = 1e-10;
end

if (~isfield(param, 'maxprint'))
  param.maxprint = 5;
end

% do the regression
% rng('default');
x = zeros(m,1);
oldexpy = -ones(size(y));
for iter = 1:param.maxiter

  adjy = a * x;
  expy = 1 ./ (1 + exp(-adjy));
  deriv = expy .* (1-expy); 
  wadjy = w .* (deriv .* adjy + (y-expy));
  weights = spdiags(deriv .* w, 0, n, n);

  x = inv(a' * weights * a + ridgemat) * a' * wadjy;

  if (param.verbose)
    len = min(param.maxprint, length(x));
    fprintf('%3d: [',iter);
    fprintf(' %g', x(1:len));
    if (len < length(x))
      fprintf(' ... ');
    end
    fprintf(' ]\n');
  end

  if (sum(abs(expy-oldexpy)) < n*param.epsilon)
    if (param.verbose)
      fprintf('Converged.\n');
    end
    return;
  end
  
  oldexpy = expy;

end

warning('logistic:notconverged', 'Failed to converge');
