% Calculation of gradient and objective for multinomial MMMF.
%
% function [obj,grad,lossobj,regobj] = m3fmultinomial(v,Y,lambda,varargin)
% v - vector of parameters [n*p+m*p,1]
% Y - target matrix [n,m]
% lambda - regularization parameter [scalar]
% obj - value of objective at v [scalar]
% grad - gradient at v [n*p+m*p,1]
% lossobj - loss component of objective [scalar]
% regobj - regularization component of objective [scalar]
% 
% Written by Jason Rennie, February 2006
% Last modified: Wed Mar 29 21:49:06 2006

function [obj,grad,lossobj,regobj] = m3fmultinomial(v,Y,lambda,varargin)
  fn = mfilename;
  if nargin < 3
    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)./m;
  p = length(v)./(n+m);
  if p ~= floor(p) | p < 1
    error('dimensions of v and Y don''t match');
  end
  U = reshape(v(1:n*p),n,p); % [n,p]
  %U = ones(n,p);
  %V = reshape(v(1:m*p),m,p); % [m,p]
  V = reshape(v(n*p+1:n*p+m*p),m,p); % [m,p]
  %V = ones(m,p);

  X = U*V';
  eX = exp(X); % [n,m]
  Yp = sum(Y,2)*ones(1,p); % [n,p]
  Ym = sum(Y,2)*ones(1,m); % [n,m]
  eXp = sum(eX,2)*ones(1,p); % [n,p]
  eXm = sum(eX,2)*ones(1,m); % [n,m]
  dU = lambda.*U + (eX*V)./eXp.*Yp - Y*V; % [n,p]
  dV = lambda.*V + (eX.*Ym./eXm)'*U - Y'*U; % [m,p]
  regobj = lambda./2.*(sum(U(:).^2) + sum(V(:).^2)); % [scalar]
  lossobj = log(sum(eX,2))'*sum(Y,2) - sum(sum(Y.*X)); % [scalar]
  obj = regobj + lossobj;
  grad = [dU(:); dV(:)];
  %grad = dV(:);
  if verbose
    fprintf(1,'lambda=%.2e obj=%.2e grad''*grad=%.2e time=%.1f\n',lambda,obj,grad'*grad,etime(clock,t0));
  end