function alpha = dirpolya(file,prec,maxiter)
% alpha = dirpolya(file,[prec,maxiter])
% calculates an exact solution to Dirichlet-Multinomial/Polya distribution.
% Thomas Minka version.
% file    : file of sparse feature format,
% prec    : desired precision,
% maxiter : maximum number of iteration.
% $Id: dirpolya.m,v 1.3 2005/07/04 15:27:43 dmochiha Exp $
N = fmatrix(file)';
F = full(sum(N,1));
I = size(N,1);
K = size(N,2);
T = 100;
alpha = F / I;
thresh = 1.0e-4;
% optional arguments
if nargin > 1
  thresh = prec;
  if nargin > 2
    T = maxiter;
  end
end
% main
for t = 1:T
  alpha_old = alpha;
  z1 = sum(psi(F + sum(alpha)) - psi(sum(alpha)));
  for k = 1:K
    fprintf(1,'iteration [%d] : lexicon %4d/%d..\r', t, k, K);
    z2 = 0;
    [i,j,c] = find(N(:,k));
    for n = c'
      z2 = z2 + psi(n + alpha(k)) - psi(alpha(k));
    end
    if z1 > 0
      alpha(k) = alpha(k) * z2 / z1;
    end
  end
  if converged(alpha_old,alpha,thresh)
    fprintf(1,'\nconverged.\n');
    break;
  else
    fprintf(1,'\t\t\t\t\tdifference = %g%%\r', ...
	    diff_vec(alpha_old,alpha) * 100);
  end
end
% $Id: dirpolya.m,v 1.3 2005/07/04 15:27:43 dmochiha Exp $
