function holes = findholeofsizeD(EorA,D,timeout,silent)
%FINDHOLEofsizeD Finds holes (chordless cycles) in graph of length D.
% holes = findholeofsizeD(E or A, D, timeout)
%     E = |E(G)|-by-2 matrix enumerating the edges of graph G
%     A = |V(G)|-by-|V(G)| binary adjacency matrix for graph G
%     D = vector containing hole lengths to search for
%         e.g. D=[5 7 9 11 13] to find some obvious odd holes
%     timeout = number of seconds to search for holes before giving up
%
% THIS ALGORITHM HAS EXPONENTIAL RUNNING TIME, SO IT WILL RARELY BE
% FEASIBLE TO RUN TO COMPLETION.  ITS OUTPUT IS A LIST OF HOLES FOUND --
% RUN FOR SEVERAL MINUTES AND IT'S LIKELY THAT THE GRAPH CONTAINS NO HOLES
% OF LENGTH D.

% (by Delbert Dueck, delbert@cs.columbia.edu, 2009)

if size(EorA,2)==2, E=EorA; end;
if size(E,1)==size(E,2), A=EorA; end;
if nargin<3, timeout=60; end; % keep looking for a minute
% if nargin<4, silent=0; else silent=strcmp(silent,'silent'); end;

if exist('E','var'), V=max(max(E));
    for v=1:V, N{v}=[E(find(E(:,1)==v),2)',E(find(E(:,2)==v),1)']; end; % enumerate neighbors
elseif exist('A','var'), V=size(E,1);
    for v=1:V, N{v}=find(A(v,:)); end; % enumerate neighbors
end;
start = clock;

holes = cell(0);
findhole_(N,D,[]);

    function findhole_(N,D,v_)
        if etime(clock,start)>timeout, return; end;
        if isempty(v_), vs=1:numel(N); else vs=N{v_(end)}; end;
        for v=vs,
            if numel(v_)>=2, if v==v_(end-1), continue; end; end; % don't backtrack
            if ~isempty(intersect(N{v},v_(2:end-1))), continue; end; % chord encountered
            if numel(v_)>=2, % at least three deep
                if any(N{v}==v_(1)), % hole encountered
                    if any(D==numel(v_)+1), holes{end+1,1}=[v_,v,v_(1)]; % relevant hole
                        if ispc, disp([v_,v,v_(1)]); end; % fprintf(['%d-hole: %d',repmat('-%d',[1 numel(v_)]),'-%d\n',
                    end;
                    continue;
                end;
            end;
            if numel(v_)<=max(D), findhole_(N,D,[v_,v]); end;
        end;

    end

end


% i=find(E(:,1)>E(:,2)); E(i,[2 1])=E(i,[1 2]); E_=setdiff(nchoosek(1:max(max(E)),2),E,'rows'); % misc code to find graph complement
