function alpha = dirpolya_mackay(file,prec,maxiter)
% alpha = dirpolya_mackay(file,[prec,maxiter])
% calculates a MacKay approximation to Dirichlet-Multinomial/Polya
% distribution.
% file    : file of h * w sparse feature format.
% maxiter : number of maximum iteration.
% prec    : precision of convergence. (default 0.0001)
% $Id: dirpolya_mackay.m,v 1.2 2005/07/04 15:27:53 dmochiha Exp $
fprintf(1,'loading data..\n');
N = fmatrix(file)';
F = full(sum(N,1));
D = size(N,1);
L = size(N,2);
T = 100;
alpha = F / L;
thresh = 1.0e-4;
% optional parameters
if nargin > 1
  thresh = prec;
  if nargin > 2
    T = maxiter;
  end
end
% statistics cache
for i = 1:L
  fprintf(1,'initializing cache %d/%d..\r',i,L);
  v_cache(i) = func_v(i,N);
  g_cache(i) = func_g(i,N);
  h_cache(i) = func_h(i,N);
end
fprintf(1,'\n');
% main
for t = 1:T
  % fprintf(1,'iteration %d..\n',t);
  alpha_prev = alpha;
  for i = 1:L
    fprintf(1,'iteration %d : lexicon %d/%d..\r',t,i,L);
    v = v_cache(i);
    g = g_cache(i);
    h = h_cache(i);
    k = func_k(alpha,F);
    alpha(i) = 2 * v / (k - g + sqrt((k - g)^2 + 4 * h * v));
  end
  if converged(alpha_prev,alpha,thresh)
    fprintf(1,'\nconverged.\n');
    break;
  else
    fprintf(1,'\t\t\t\t\tdifference = %g%%\r', ...
	    diff_vec(alpha_prev,alpha) * 100);
  end
end

%% functions

function n = func_v(i,N)
n = size(find(N(:,i)),1);

function z = func_g(i,N)
z = 0;
M = full(max(N(:,i)));
for f = 2:M
  n = size(find(N(:,i) >= f),1);
  z = z + n / (f - 1);
end

function z = func_h(i,N)
z = 0;
M = full(max(N(:,i)));
for f = 2:M
  n = size(find(N(:,i) >= f),1);
  z = z + n / ((f - 1)^2);
end

function z = func_k(alpha,F)
z = 0;
L = size(F,2);
alpha0 = sum(alpha);
log_alpha0 = log(alpha0);
for i = 1:L
  f = F(i);
  z = z + log(f + alpha0) - log_alpha0 ...
        + 0.5 * (f / (alpha0 * (f + alpha0)));
end

% $Id: dirpolya_mackay.m,v 1.2 2005/07/04 15:27:53 dmochiha Exp $
