% Calculation of gradient and objective for All-Threshold Logistic
% MMMF.
%
% function [obj,grad,lossobj,regobj] = m3flogistic(v,Y,lambda,l,varargin)
% v - vector of parameters [n*p+m*p+n*(l-1),1]
% Y - rating matrix (labels) [n,m]
% lambda - regularization parameter [scalar]
% l - # of unique rating values (1..l)
% obj - value of objective at v [scalar]
% grad - gradient at v [n*p+m*p+n*(l-1),1]
% lossobj - loss component of objective [scalar]
% regobj - regularization component of objective [scalar]
% 
% Written by Jason Rennie, January 2005
% Last modified: Fri Dec 15 17:29:00 2006

function [obj,grad,lossobj,regobj] = m3flogistic(v,Y,lambda,l,varargin)
  fn = mfilename;
  if nargin < 4
    error('insufficient parameters')
  end
  % Parameters that can be set via varargin
  verbose = 1;
  % Process varargin
  paramgt;
  
  t0 = clock;
  [n,m] = size(Y);
  p = (length(v)-n.*(l-1))./(n+m);
  if p ~= floor(p) | p < 1
    error('dimensions of v and Y don''t match l');
  end
  U = reshape(v(1:n*p),n,p);
  V = reshape(v(n*p+1:n*p+m*p),m,p);
  theta = reshape(v(n*p+m*p+1:n*p+m*p+n*(l-1)),n,l-1);

  X = U*V';
  Ygt0 = Y>0;
  BX = X.*Ygt0;
  %reallyzero = (X==0).*Ygt0;
  %fprintf(1,'%s: sum(sum(reallyzero))=%d\n',fn,full(sum(sum(reallyzero))));
  %end
  clear X;
  dU = lambda.*U; % [n,p]
  dV = lambda.*V; % [m,p]
  dtheta = zeros(n,l-1); % [n,l-1]
  regobj = lambda.*(sum(U(:).^2)+sum(V(:).^2))./2; % [scalar]
  lossobj = 0;
  for k=1:l-1
    S = Ygt0-2.*(Y>k);
    BZ = S.*(theta(:,k)*ones(1,m))-S.*BX; % [n,m]
    issparse(BZ);
    tmp = S.*hprime(BZ); % [n,m]
    dU = dU - tmp*V; % [n,p]
    dV = dV - tmp'*U; % [m,p]
    dtheta(:,k) = tmp*ones(m,1);
    lossobj = lossobj + sum(sum(h(BZ)));
  end
  obj = regobj + lossobj;
  grad = [dU(:); dV(:); dtheta(:)];
  if verbose
    fprintf(1,'%s: lambda=%.2e obj=%.2e grad''*grad=%.2e time=%.1f\n',fn,lambda,obj,grad'*grad,etime(clock,t0));
  end

function [ret] = h(z)
  [m,n] = size(z);
  [i,j,v] = find(z);
  hv = log(1+exp(-v));
  ret = sparse(i,j,hv,m,n);

% ret = (z>0).*(z<1).*(z-1) - (z<=0);
function [ret] = hprime(z)
  [m,n] = size(z);
  [i,j,v] = find(z);
  hprimev = -(exp(-v)./(1+exp(-v)));
  ret = sparse(i,j,hprimev,m,n);

% ChangeLog
% 12/8/06 - remove reallyzero (print warning if it would be necessary)
% 3/23/05 - made calcultions take better advantage of sparseness
% 3/18/05 - fixed bug in objective (wasn't squaring fro norms)
% 3/1/05 - added objective calculation
% 2/23/05 - fixed bug in hprime()