function [w,b, invH]= logist_phi(phi,Y,alpha, inind,w0)
% logistic regression on class-dependent features
% phi is a cell array of length C of features (C is the # of classes)
% each element is of size nxk_i (i=1..C)
% Y is in the range of 1..C
% let kall=sum_i k_i
% alpha, inind, and w0 are all of length kall
% OUTPUT:
% w is of length kall
% b is of length C
% invH is of size kallxkall
  
% check inputs and fill in defaults
  C=length(phi); assert(max(Y)==C);
  ns=cellfun('size',phi,1); assert(all(diff(ns)==0)); n=ns(1); assert(n==length(Y));
  ks=cellfun('size',phi,2);
  k2=[0; cumsum(ks(:))]; kall=k2(end);
  if (~exist('alpha','var'))
    alpha=zeros(kall,1);
  end
  if (exist('alpha','var') & length(alpha)==1 & kall>1)
    alpha=alpha*ones(kall,1);
  end
  if (~exist('inind','var'))
    inind=true(size(alpha));
  end
  if (~exist('w0','var')),
    w0=zeros(size(alpha));
  end
  assert(islogical(inind));
  inind=inind(:);
  for i=1:C
    k_in{i}=inind(k2(i)+1:k2(i+1));
  end



  % augment constant to input and prune input according to inind
  for i=1:C
    phi{i}=[phi{i}(:,k_in{i}) ones(n,1)];
  end
  
  % prepare for sub indices
  ks=zeros(C,1);
  for i=1:C
    k_in{i}=[k_in{i}; 1];
    ks(i)=sum(k_in{i});
  end
  b_ind=cumsum(ks);
  k2=[0;b_ind];
  kall=b_ind(end);
  w_ind=setdiff([1:kall],b_ind);

  % set optimization constants
  opts=optimset('MaxIter',100);
  opts=optimset(opts,'Display','iter');
  alpha_const=2*1E-5;
  
  % initialize w and alpha
  wb0=zeros(kall,1);
  wb0(w_ind)=w0(inind); 
  alpha_wb=zeros(kall,1);
  alpha_wb(w_ind)=alpha(inind);
  alpha_wb(b_ind)=alpha_const;

  % now do the newton steps
  [wb,H]=newton(@logist_phi_fgH, wb0, opts, phi,Y,alpha_wb);

% return values
  invH=inv(H(w_ind,w_ind));
  % update the w
  w=w0;
  w(inind)=wb(w_ind);
  b=wb(b_ind);