function testsubdivision(flag)
% testsubdivision( [flag] )
% Demonstrates how to use the subdivision and tjsr-package and tests if the subdivision package works correctly
% if flag is 1, then the functions in the subdivision-package are tested.
%
% Written by: tommsch, 2018

%#ok<*NOPRT>
%#ok<*ASGLU>
%#ok<*NASGU>

if(nargin>=1 && flag); selftestsubdivision; return; end;

vprintf('Either start the function with argument <1> or read the source code of this file.\n'); 
return; 
2; %#ok<UNRCH>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 1-D example of multiple subdivision scheme
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%make finite set of subdivision operators
fprintf('make finite set of subdivision operators\n'); 
S1=getS('a',1/4*[1;3;3;1],'M',-2);
S2=getS('1_rand'); %random scheme
S=[S1;S2]; %each subdivision operator is one line in the cell array. Due to strange matlab-behaviour, if the cell C shall be displayed, an error is thrown.

%plot blf with two different ordering
blf({[2 2],[1 1 2]},S); %ordering 2211211212112
blf(S); %random ordering

[T, Om]=transitionmatrix(S); %compute transition matrices and the set Omega

plotm(Om,'x') %plot the set Omega

V0=constructV(Om,0)                   %Construct space orthogonal to polynomial sequences
V0bar=constructVbar(Om,0);            %Construct space of differences of delta
if(size(V0,2)~=size(V0bar,2)); 
    fprintf('Set Omega is wrong.');   %Test ifV0==V0bar, if not the JSR is no characterization for regularity
end
TV0=restrictmatrix(T,V0); vdisp(TV0); %restrict transition matrices

[JSR, type]=tjsr(TV0);  %compute JSR of restricted matrices. If the algorithm does not terminate, set delta to some value less then one.
%[JSR, type]=tjsr(TV0,'maxsmpdepth',2,'delta',.999999); %compute the continuity regularity. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% s-D example of subdivision scheme
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S=getS('2_butterfly','sym');        %Butterfly-scheme
[~,~,xyzv]=blf(S,'iteration',3);   %plot the blf
tile(getS(S,'supp'));               %plot the support of the blf
[T, Om]=transitionmatrix(S);        %compute transition matrices and the set Omega
plotm(Om,'x');                      %plot the set Omega
V3=constructV(Om,3);                %Construct space of fourth differences
[~,TV3]=restrictmatrix(T,V3);           %restrict transition matrices

% The transition matrices of this scheme have a bunch of invariant subspaces. 
% Thus the invariant-polytope algorithm behaves badly.
% To encounter that problem, we attack the subspaces one by one
TV3=invariantsubspace(TV3)
TV3=[invariantsubspace(TV3{1}) invariantsubspace(TV3{2})]

%lets have a look at the subspaces sizes
cellfun(@(x) size(x{1},2),TV3)

[JSR1, type1]=tjsr(TV3{1}); %compute Hoelder regularity.
[JSR2, type2]=tjsr(TV3{2},'maxsmpdepth',1,'delta',.999999999); %compute Hoelder regularity.
[JSR3, type3]=tjsr(TV3{3},'maxsmpdepth',1,'delta',.999999999); %compute Hoelder regularity.
cellfun(@(x) max(abs([x{:}])),TV3(4:end)) %all the other spectral radii

end

function selftestsubdivision
    vprintf('The subdivision-package will be tested. All output on the screen should be followed by an explanation.\n','cpr',[.1 .4 .1]);
    vprintf('If there is some error or un-commented output, the test fails.\n\n','cpr',[.1 .4 .1]);

    S1=getS('a',[1 2 1],'M',2);
    vprintf('Error that the mask does not sum up should be here.\n','cpr',[.1 .4 .1]);
    vprintf('Error that the mask is transposed should be here.\n\n','cpr',[.1 .4 .1]);
    S1=getS('a',1/2*[1 2 1]','M',2);
    
    S2=getS('2_rand');
    S=getS([S1;S2]);
    vprintf('Error that some dimension or size is wrong.\n\n','cpr',[.1 .4 .1]);
    S2=getS('1_1133');
    S=getS([S1;S2]);
    
    T=transitionmatrix(S);
    
    %from here on the functions are called in lexicographic ordering
    
    figure; blf(S,'iteration',4);
    vprintf('FIGURE: A figure should show a 1d-blf.\n\n','cpr',[.4 .2 0]);
    figure; blf(getS('2_butterfly'),'iteration',3,'diff',1)
    vprintf('FIGURE: A figure should show two 2d-blf.\n\n','cpr',[.4 .2 0]);
    
    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,'verbose',0);
    C12=checktile(M12,D12,'verbose',0);
    if(C1~=true || C2~=true || C12~=false); vprintf('checktile() is broken.\n\n','cpr','err'); end;
    
    cc=round(3*compresscoordinates([10 20 ; 40 50],[2 1; -1 1]));
    if(~isequal(cc,[0 1 -1 0;0 1 2 3;30 120 60 150])); vprintf('compresscoordinates() is broken.\n\n','cpr','err'); end;
    
    d=constructdigit( [2 0; 1 2] );
    if(~isequal(d,[0 0 1 1;0 1 1 2])); vprintf('constructdigit() is broken.\n\n','cpr','err'); end;
    
    Om=constructOmega(S,'Omega',[0],'stable');
    
    oo=constructordering([1 2],[3 ],[10 20],[30]);
    if(~isequal(oo,{[1 2;10 20], [3;30]})); vprintf('constructordering() is broken.\n\n','cpr','err'); end;
    
    %constructPhi %untested
    
    %constructU   %untested
    V0=constructV(Om,0);
    V1=constructV(Om,1);
    
    Vbar=constructVbar(Om);
    
    daubechiesmask( 5 );
    
    %differencescheme %untested
    
    dim=dimVVbar([1 1 0 0;1 0 1 1;0 1 0 0],'01','verbose',0);
    if(~isequal(dim,[5 3; 3 0])); vprintf('constructVbar()/constructV()/dimVVbar() is broken.\n\n','cpr','err'); end;
    
    [of]=findperiod([1 2 1 2 1 2; 4 5 5 5 5 5]);
    if(~isequal(of{1},[1;4]) || ~isequal(of{2},[2 1; 5 5])); vprintf('findperiod() is broken.\n\n','cpr','err'); end;
    
    %getS %already tested at the beginning
    
    if(isordering(oo)~=true || isordering([2 3])~=false); vprintf('isordering() is broken.\n\n','cpr','err'); end;
    
    if(isS(S1)~=true || isS([2 3]~=false)); vprintf('isS() is broken.\n\n','cpr','err'); end;    
    
    if(isT(T)~=true || isS([2 3]~=false)); vprintf('isT() is broken.\n\n','cpr','err'); end;    
    
    Sn=normalizeS(S,'equal');
    if(~isequal(2*Sn{1,1}.c,[1 2 1]') || ~isequal(2*Sn{2,1}.c,[1 1 1 1]')); vprintf('normalizeS() is broken.\n\n','cpr','err'); end;    
    
    num=ordering2num({[2 1 ],[]},getS({[],2}));
    if(~isequal(num,1/2)); vprintf('ordering2num() is broken.\n\n','cpr','err'); end;
    
    figure;
    Snum=getS({[],[1 2; -1 2]});
    num=ordering2num(num2ordering(Snum,[.24;.3],'l',50,'verbose',0,'plot',1),Snum);
    if(abs(num-[0.24;0.3])>10^-12); vprintf('ordering2num()/num2ordering() is broken.\n\n','cpr','err'); end;
    vprintf('FIGURE: A figure should show a tile and a connected lines in it.\n\n','cpr',[.4 .2 0]);
    
    vec=ordering2vector({ [10 20 30] , 1:2 },10);
    if(~isequal(vec,[10 20 30 1 2 1 2 1 2 1])); vprintf('ordering2vector() is broken.\n\n','cpr','err'); end;
    
    peter(rand(2,30),4);
    
    restrictmatrix(T,V0);
    restrictmatrix(T,V1);
    vprintf('An error should tell that two matrices are not invariant.\n\n','cpr',[.1 .4 .1]);
    
    %testsubdivision %this file
    
    figure;
    tile(getS('2_rand'))
    vprintf('FIGURE: A figure should show a fractal.\n\n','cpr',[.4 .2 0]);
    figure;
    tile([getS('2_rand');getS('2_rand')],'supertile',6);
    vprintf('FIGURE: A figure should show a Rorschach test.\n\n','cpr',[.4 .2 0]);
    
    tilearea(tile(S,'plot',0));
    vprintf('Some text output should state something.\n\n','cpr',[.1 .4 .1]);
    
    %transitionmatrix %tested above
    
    [o1, o2]=vector2ordering([1 2 1 2 1 2; 4 5 5 5 5 5]);
    if(~isequal(o1,[1;4]) || ~isequal(o2,[2 1; 5 5])); vprintf('vector2ordering()/findperiod() is broken.\n\n','cpr','err'); end;
end


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