function [ flag, C, T ] = checktile( varargin)
% [ flag, C, T ] = checktile( M ,D, [options])
% [ flag, C, T ] = checktile( S, [options] )
% Computes if a pair dilation-matrix/digit-set generates a tile (i.e. an attractor with area 1).
%
% Input:
%   M               dilation matrix
%   D               digit set for M
%
% Options:
%   'verbose',val   Verbose level
%   'legacy',val    Use old algorithm to compute the contact-vectors T, with set of integers of size [-val val]^dim.
%                   If val==0, a value for val will be determined automatically.
%
% Output:
%   flag            true:tile, false:no tile
%   C               Contact matrix
%   T               Contact vectors
%
% Info: Lattice = \ZZ^dim
%
% E.g.: M1=[1 2; -2 -2];  D1=[0 1; 0 -1]; 
%       M2=[-1 1; -1 -2]; D2=[0 0 0; -2 -1 0]; 
%       M12=M2*M1; D12=setplus(M2*D1,D2);
%       C1=checktile(M1,D1);
%       C2=checktile(M2,D1);
%       C12=checktile(M12,D12);
%
% References:
% K.H. Groechenig and A. Haas, Self-Similar Lattice Tilings
%
% See also: tilearea
%
% Written by: tommsch, 2018

%Parse input
[verbose,varargin]=parsem('verbose',varargin,1);
[legacy,varargin]=parsem('legacy',varargin,[]);

if(isS(varargin{1})); 
    M=varargin{1}{2};
    D=varargin{1}{3};
    varargin(1)=[];
else
    M=varargin{1};
    D=varargin{2};
    varargin(1:2)=[];
end

if(~isempty(varargin)); vprintf('Unkown options: \n %v',varargin,'cpr','err','imp',[0 verbose]); end;

%Make contact vectors
if(isempty(legacy)); 
    T=contactvectors( M, D ); 
else
    T=contactvectors_old( M ,D, legacy, verbose ); 
end

C=contactmatrix(M, D, T);

SPEC=max(abs(eig(C)));
if SPEC<abs(det(M));
    vprintf('(M,D) generates a tile.\n','imp',[1 verbose],'cpr',[0.1 .5 .0]);
    flag=1;
else
    vprintf('(M,D) generates NO tile.\n','imp',[1 verbose],'cpr','err');
    flag=0;
end

end

function [ C ] = contactmatrix( M, D, T );
nT=size(T,2);
C=zeros(nT,nT);
for k=1:nT
    for l=1:nT
        val=intersect( setplus(M*T(:,k),D).', setplus(T(:,l),D).', 'rows').';
        C(k,l)=size(val,2);
    end
end
end

function [ T ] = contactvectors( M ,D )
%Computes the matrices $\mathcal{T}_n$ and $\mathcal{T}$ as described in 
%Algortihm from KH Groechnig and A Haas, Self-Similar Lattice Tilings, 
%Thm 2.2 and Lemma 4.5
dim=size(M,1);

%Make T0 set
T0=unique([eye(dim) -eye(dim)].','rows').'; %start-cector for lattice = \ZZ^2, T0=\mathcal{T}_0  %I use unique to sort the vector

Tn{1}=T0;
detM=abs(det(M));

while(true)
    %vprintf('%v\n',Tn{end});
    val=M\setplus(Tn{end},D,-D);
    idx=allm(abs(val-round(val))<1/(2*detM),1);
    Tn{end+1}=unique([round(val(:,idx)) Tn{end}].','rows').'; %#ok<AGROW> %In the paper it is only written, the Tn stabilize. They are not necessarily contained in each other!
    if isequal(Tn{end},Tn{end-1});
        T=unique([Tn{:}]','rows')';
        break;
    end;
end
idx=ismember(T.',zerosm(1,dim),'rows');
T(:,idx)=[]; %remove zero vector %See Definition of \mathcal{T}_n^* before Lemma 4.5
    

end


function [ T ] = contactvectors_old( M ,D, Zval, verbose ) 
%Computes the matrices $\mathcal{T}_n$ and $\mathcal{T}$ as described in 
%Algortihm from KH Groechnig and A Haas, Self-Similar Lattice Tilings, 
%Thm 2.2 and Lemma 4.5
dim=size(M,1);
INDEX=[];

%Make Z lattice without origin
if(Zval==0); 
    Zval=(maxm(D(:))+1)*(max(abs(M(:)))+1)*size(D,1); 
    vprintf('Zval=%d\n',Zval,'imp',[1 verbose]);
end;
    
ZZ=mixvector(-Zval: Zval,dim); %represents the integers
zeroindex = ismember(ZZ',zeros(1,dim),'rows');
ZZ=ZZ(:,~zeroindex);%remove (0,0)

%Make T0 set
T0=unique([eye(dim) -eye(dim)],'rows'); %Start Vector for Lattice = Z^2, T0=\mathcal{T}_0  %I use unique to sort the vector

Tn{1}=T0;

while(true)
    for i=1:size(ZZ,2);
        INTERSECT=intersect(setplus(M*ZZ(:,i),D)',setplus(Tn{end},D)','rows');
        if(~isempty(INTERSECT));
            INDEX=[INDEX i]; %#ok<AGROW>
        end
    end
    Tn{end+1}=unique(ZZ(:,INDEX)','rows')'; %#ok<AGROW> %I use unique to sort the vector
    if isequal(Tn{end},Tn{end-1});
        T=unique([Tn{:}]','rows')';
        break;
    end;
end
    

end

function dummy; end %#ok<DEFNU> %Generates an error, if the 'end' of a function is missing. 

