function CPT = mk_CPT_from_mux2(rgate, cgate, input, arity)
% MK_CPT_FROM_MUX2 Make a CPT representation of a 2D multiplexer (can be quite large!)
% CPT = mk_CPT_from_mux2(rgate, cgate, input, arity)
% 
% e.g., for a 2x2 multiplexer, CPT(M(1,1), M(1,2), M(2,1), M(2,2), r, c, i) = p
% where p = 1 if M(r,c) = i, and p = 0 otherwise

[nrows ncols] = size(input);
n = nrows*ncols;
an = arity^n;
if 1
  CPT = zeros(an, nrows, ncols, arity);
  U = ind2subv(arity*ones(1,n), 1:an);
  for r=1:nrows
    for c=1:ncols
      g = (r-1)*nrows + c;
      tmp = sparse(1:an, U(:,g), ones(an,1), an, arity);
      % tmp(i, U(i,g)) = 1, where tmp is A^n * A in size
      CPT(:,r,c,:) = full(tmp);
    end
  end
  CPT = reshape(CPT, [arity*ones(1, n) nrows ncols arity]);
  % We may need to re-order the dimensions to be consistent with the node numbering.
  % e.g., rgate = 1, cgate = 2, input = [4 5; 6 7],
  % Then we need to put the gate's dimension (assumed to be n+1=3 above) first.
  if rgate < input | cgate < input
    perm = [n+1 n+2 1:n n+3];
    CPT = permute(CPT, perm);
  end
else
  ns = zeros(1, n+3);
  ns(rgate) = nrows;
  ns(cgate) = ncols;
  ns(end) = arity;
  ns(input(:)) = arity;
  CPT = myzeros(ns);
  for r=1:nrows
    for c=1:ncols
      for ui=1:arity^n
	u = ind2subv(ns(input(:)), ui);
	u = reshape(u, [nrows ncols]);
	v = zeros(1,n+2);
	v(rgate)=r;
	v(cgate)=c;
	v(input)=u(:);
	j = subv2ind(ns, [v u(r,c)]);
	CPT(j)=1.0;
      end 
    end
  end
end
