function pot2 = marginalize_cgpot(pot1, onto)


dom1 = union(pot1.cdom, pot1.ddom);
sumover = setdiff(dom1, onto);
csumover = intersect(sumover, pot1.cdom);
dsumover = intersect(sumover, pot1.ddom);
keep = onto;
dkeep = intersect(keep, pot1.ddom);
ckeep = intersect(keep, pot1.cdom);
%ns = sparse(1, max(dom1));
ns = zeros(1, max(dom1));
ns(pot1.ddom) = pot1.dsizes;
ns(pot1.cdom) = pot1.csizes;

% sum(ns(csumover))==0 is like isempty(csumover) but handles clamped nodes.
% Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)

% If not marginalizing over any hidden nodes, there is nothing to do except change the domain
if sum(ns(csumover))==0 & prod(ns(dsumover))==1
  pot2 = mk_cgpot(dkeep, ckeep, ns, pot1.subtype);
  pot2.pot = pot1.pot;
  return;
end

% Marginalize the cts parts.
% If we are in canonical form, we stay that way, since moment form might not exist.
% Besides, we would like to minimize the number of conversions.
if sum(ns(csumover)) > 0
  if pot1.subtype == 'm'
    for i=1:pot1.dsize
      [pot1.pot{i}.logp, pot1.pot{i}.mu, pot1.pot{i}.Sigma] = ...
	  marginalize_gaussian(pot1.pot{i}.logp, pot1.pot{i}.mu, pot1.pot{i}.Sigma, ...
			       ckeep, csumover, pot1.cdom, ns(pot1.cdom));
    end
  else
    for i=1:pot1.dsize
      [pot1.pot{i}.g, pot1.pot{i}.h, pot1.pot{i}.K] = ...
	  marginalize_canonical(pot1.pot{i}.g, pot1.pot{i}.h, pot1.pot{i}.K, ...
				ckeep, csumover, pot1.cdom, ns);
    end
  end
end

% Not marginalizing over any discrete nodes, so we are done.
if prod(ns(dsumover))==1
  pot2 = mk_cgpot(dkeep, ckeep, ns, pot1.subtype);
  pot2.pot = pot1.pot;
  return;
end

% To marginalize the discrete parts, we must be in moment form.
pot1 = cg_can_to_mom(pot1);

I = prod(ns(dkeep));
J = prod(ns(dsumover));
C = sum(ns(ckeep));
pot2 = mk_cgpot(dkeep, ckeep, ns, 'm');

% Reshape pot1 into the form mu1(:,j,i), where i is in dkeep, j is in dsumover
T1 = zeros(I,J);
mu1 = zeros(C,J,I);
Sigma1 = zeros(C,C,J,I);
sum_map = find_equiv_posns(dsumover, pot1.ddom);
keep_map = find_equiv_posns(dkeep, pot1.ddom);
iv = zeros(1, length(pot1.ddom));
for i=1:I
  keep_iv = ind2subv(ns(dkeep), i);
  iv(keep_map) = keep_iv;
  for j=1:J
    sum_iv = ind2subv(ns(dsumover), j);
    iv(sum_map) = sum_iv;
    k = subv2ind(ns(pot1.ddom), iv);
    T1(i,j) = exp(pot1.pot{k}.logp);
    if C > 0 % so mu1 and Sigma1 are non-empty
      mu1(:,j,i) = pot1.pot{k}.mu;
      Sigma1(:,:,j,i) = pot1.pot{k}.Sigma;
    end
  end
end

% Collapse the mixture of Gaussians
if I==1
  %coef = T1;
  coef = normalise(T1); % coef must be convex combination
else
  coef = mk_stochastic(T1);
end
T2 = sum(T1,2);
T2 = T2 + (T2==0)*eps;
%if C > 0, disp('collapsing onto '); disp(onto); end
for i=1:I
  if C > 0
    [pot2.pot{i}.mu, pot2.pot{i}.Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
  end
  pot2.pot{i}.logp = log(T2(i));
end






