function Binder = BinderFinder(Phi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Binder Finder
%%%%%% Code to find the embedded simplices (ETFs for their spans with one
%%%%%% more vector than their rank) in a given ETF.
%%%%%%
%%%%%% Input: Phi, the synthesis matrix of an equiangular tight frame.
%%%%%%
%%%%%% Output: Binder, the incidence matrix (rows = simplices, cols = frame
%%%%%% vectors) of the simplices embedded in the given ETF
%%%%%%
%%%%%% Citation: "Equiangular tight frames that contain regular simplices,"
%%%%%% Matthew Fickus, John Jasper, Emily J. King, Dustin G. Mixon, 2017
%%%%%%
%%%%%% Code created: July 2016
%%%%%% Last updated: November 22, 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[d n] = size(Phi); 
k = sqrt(d*(n-1)/(n-d))+1;

Gram = Phi'*Phi;
Gram = n*(k-1)*Gram/trace(Gram); % Normalization so that triple products equal 1

% Initializing incidence matrix of triples
TripleCounter = 0;
for ii = 1:n-2
    for jj = ii+1:n-1
        for kk = jj+1:n
            if abs(Gram(ii,jj)*Gram(jj,kk)*Gram(kk,ii)+1) < 10^(-2)
                TripleCounter = TripleCounter + 1;
            end
        end
    end
end

% Filling incidence matrix of triples
Triple = zeros(TripleCounter,n);
TripleCounter = 0;
for ii = 1:n-2
    for jj = ii+1:n-1
        for kk = jj+1:n
            if abs(Gram(ii,jj)*Gram(jj,kk)*Gram(kk,ii)+1) < 10^(-2)
                TripleCounter = TripleCounter + 1;
                Triple(TripleCounter,:) = zeros(1,n);
                Triple(TripleCounter,[ii jj kk]) = 1;
            end
        end
    end
end

% Incrementally finding j-tuples for j>3
jTuple = Triple;
for jTupleSize = 3:k-1
    tic
    for ii = 1:size(jTuple,1)
        Indices = find(jTuple(ii,:));
        Indicator = sum(double(jTuple(:,Indices)*(ones(jTupleSize)-eye(jTupleSize)) == jTupleSize-1),2)'*jTuple;
        Indicator(Indices) = 0;
        NewIndices = find(Indicator == jTupleSize);
        A = kron(ones(length(NewIndices),1),jTuple(ii,:));
        A(:,NewIndices) = eye(length(NewIndices));
        if ii == 1
            NewjTuple = A;
        else
            NewjTuple(size(NewjTuple,1)+(1:size(A,1)),:) = A;
        end
        jTuple(ii,:) = zeros(1,size(jTuple,2));
    end
    jTuple = NewjTuple;
    toc
end

Binder = jTuple;

% Visualization of (incidence matrix of) binder and Gram matrix of binder
figure('Name','Binder'), imagesc(Binder)
figure('Name','Gram Matrix of Binder'), imagesc(Binder'*Binder), colorbar

end
