function [Aeq,Beq] = constrain_marg(m,n,Pmarg)
mn = m^n;

Aeq = zeros(0,mn);
Beq = zeros(0,1);
rind = 0;
for i=1:n
  % enforce constraint: marg(X_i) = Pmarg(:,i)
  for y = 1:m
    rind = rind + 1;
    Beq(rind,1) = Pmarg(y,i);
    for xi = 1:m^n
      x = i2c(xi,m,n);
      if x(i) == y
        Aeq(rind,xi) = 1;
      end
    end
  end
end

return

function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return