function [F, alpha] = varElim(factors, query, elimOrder)
% varElim Variable elimination algorithm
% function [F, alpha] = varElim(factors, query, order)

% Infer number and sizes of variables from the factors
NF = length(factors);
scope = [];
for f=1:NF
  scope = union(scope, factors{f}.scope);
end
N = length(scope);
if ~isequal(scope, 1:N)
  error('variables should be ordered 1,2,...') 
end
sizes = ones(1,N);
for f=1:NF
  sizes(match(factors{f}.scope, scope)) = factors{f}.sizes;
end

if nargin < 3, elimOrder = 1:N; end
% make sure query nodes is last to be elminated
% order = [query setdiff(elimOrder, query)]; % setdiff sorts the order!
elimOrder(match(query, elimOrder)) = [];
order = [query elimOrder];

% Initialize the buckets with the product of the factors assigned to them
% (could be done later if we want to be lazy)
B = cell(1,N+1); 
% N+1 is a dummy bucket to hold the normalization constant
for b=1:N+1
  B{b} = factorMake([], 1, 1);
end
for f=1:NF
  b = bucket_num(factors{f}.scope, order);
  if isempty(B{b})
    B{b} = factorMake(factors{f}.scope, factors{f}.sizes, factors{f}.T);
  else
    B{b} = factorMultiply(B{b}, factors{f});
  end
end
B{N+1} = factorMake([], 1, 1); % dummy bucket to hold the normalization constant

% Do the marginalization
sum_over = setdiff(scope, query);
sum_over = sum_over(length(sum_over):-1:1); % reverse

for i=sum_over(:)'
  % summing over variable i which occurs in bucket bsrc
  bsrc = bucket_num(i, order);
  newFactor = factorMarginalize(B{bsrc}, i); % tau(i)
  B{bsrc} = [];
  % put result in bdest
  bdest = bucket_num(newFactor.scope, order);
  %fprintf('summing over bucket %d (var %d) into %d\n',  bsrc, i, bdest);
  B{bdest} = factorMultiply(B{bdest}, newFactor);
  %B{bdest}.T
end

% Combine all the remaining buckets into one
result = factorMake(query, sizes(query));
for b=1:N+1
  if ~isempty(B{b})
    result = factorMultiply(result, B{b});
  end
end
[F, alpha] = factorNormalize(result);

%%%%%%

function b = bucket_num(domain, order)
if isempty(domain)
  b = length(order)+1; % use dummy bucket at end
else
  b = max(match(domain, order));
end
