#Given an undirected (symmetric) adjacency matrix adj,
#returns the matrix of all maximal cliques where each column is
#a maximal clique
function [cliques] = maxcliques(adj)
	numverts = size(adj, 1);
	assert(size(adj, 2) == numverts);
	adj |= speye(numverts);
#	assert(issymmetric(adj));
	cliques = sparse(numverts, 1);
	for n = 1:numverts
		edgesto = diag(adj(n,:)) * cliques;
		if nnz(edgesto) == 0
			cliques(:,end+1) = sparse(n,1,true,numverts,1);
			continue
		endif

#		MISSING FEATURE: &~ as a single operator
#		subsumed = ~any(cliques &~ edgesto, 1);

		subsumed = ~any(andnot(cliques, edgesto), 1);
		cliques(:,subsumed) = [];		

		edgesto(:,~any(edgesto,1)) = [];
		edgesto = unique_subsets(edgesto);
		num_new = size(edgesto,2);
		edgesto(n,:) = true;
		cliques(:,end+1:end+num_new) = edgesto;
	endfor
endfunction

# Given a logical matrix, removes every column which is a subset of an earlier column
# or a strict subset of a later column.
function [cols] = unique_subsets(cols)
	cols = remove_sub(cols,true);
	cols = remove_sub(cols,false);
endfunction

#If earlier_only, removes all columns which are a subset of an earlier column,
#otherwise removes all columns which are a subset of any other column
function [cols] = remove_sub(cols,earlier_only)
	numcols = size(cols,2);
	submat = get_submat(cols);

	if earlier_only
		# MISSING FEATURE:
		# Retain only terms above the diagonal of a sparse matrix
		# Workaround:
		[i,j] = find(submat);
		# BUG WORKAROUND: find returns 0x0 rather than 0x1 matrix
		if size(i,1) == 0
			i = zeros(0,1);
			j = zeros(0,1);
		endif
		pairs = [i,j];
		pairs(i >= j,:) = [];
		submat = sparse(pairs(:,1), pairs(:,2), true, size(cols, 1), size(cols, 2));
	else
		# MISSING FEATURE/BUG WORKAROUND:
		# Clear the diagonal of a sparse matrix
		# Note: m(logical(speye(size(m,1)))) = 0 fails because
		# logical indexing into sparse matrices is broken
		submat -= speye(numcols) .* submat;
	endif

	issub = any(submat, 1);
	cols(:,issub) = [];
endfunction

#Returns a matrix submat where submat(i,j) is true iff cols(:,j) is a subset of cols(:,i)
function [submat] = get_submat(cols)
	counts = sum(cols, 1);
	assert(counts > 0);

#	BUG WORKAROUND: bsxfun unsupported for @ge on sparse matrix with positive vector
#	submat = bsxfun(@ge, cols.' * cols, counts);
#	This workaround is very clever

	multmat = cols.' * cols;
	submat = multmat > spones(multmat) * diag(counts - 1);
endfunction

function [a] = andnot(a,b)
	a = spones(a);
	a -= a & b;
	a = logical(a);
endfunction