function [opsi,singletons] = factorizeClique(psi)
  lpsi = log(psi./min(reshape(psi,prod(size(psi)),1)));
spsi = size(lpsi);
N = prod(spsi);
n = sum(spsi);
f = zeros(N+n+1,1);
f(1:(N+n))=1;
lb = -250*ones(N+n+1,1);
lb(1:(N+n))=0;
b = [reshape(lpsi,prod(spsi),1)];
A = zeros(N,N+n+1);
A(1:N,1:N)=eye(N,N);
for i=1:N
	q=ind2subT(spsi,i);
for j=1:length(spsi)
	A(i,N+sum(spsi(1:(j-1)))+q(j))=1;
  end
end
  A(:,end)=1;
ub = 250*ones(N+n,1);
X = linprog(f,[],[],A,b,lb,ub);
opsi = exp(reshape(X(1:N),spsi));
ind = N+1;
for i=1:length(spsi)
	singletons{i} = exp(X(ind:((ind+spsi(i))-1)));
ind=ind+spsi(i);
end
siz = sz;
nout = length(siz);
siz = double(siz);
