% Calculation of gradient and objective for All-Threshold SHC
% MMMF.  The V matrix is fixed (and passed in as a separate
% argument).  This function calculates derivatives for U and theta.
%
% function [obj,grad,lossobj,regobj] = m3fshc(v,Y,V,lambda,l,varargin)
% v - vector of parameters [n*p+n*(l-1),1]
% Y - rating matrix (labels) [n,m]
% V - the "movie" feature matrix [m,p]
% 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: Sat Dec 16 13:58:01 2006

function [obj,grad,lossobj,regobj] = m3fshcFixedV(v,Y,V,lambda,l,varargin)
  fn = mfilename;
  if nargin < 5
    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;
  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);
  theta = reshape(v(n*p+1:n*p+n*(l-1)),n,l-1);

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

% ret = (z>0).*(z<1).*((1-z).^2)/2.0 + (z<=0).*(0.5-z);
function [ret] = h(z,reallyzero)
  zin01 = (z>0)-(z>=1);
  zle0 = z<0 + reallyzero;
  ret = zin01./2 - zin01.*z + zin01.*z.^2./2 + zle0./2 - zle0.*z;

% ret = (z>0).*(z<1).*(z-1) - (z<=0);
function [ret] = hprime(z,reallyzero)
  zin01 = (z>0)-(z>=1);
  zle0 = z<0 + reallyzero;
  ret = zin01.*z - zin01 - zle0;

% ChangeLog
% 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()