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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 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,'iterations',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.','cpr',[.1 .4 .1]);
    vprintf('If there is some error or un-commented output, the test fails.\n','cpr',[.1 .4 .1]);

    S1=getS('a',[1 3 3 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/4*[1 3 3 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_rand');
    S=getS([S1;S2]);
    
    figure; blf(S,'iterations',4)
    vprintf('A figure should show a 1d-blf.\n\n');
    
    S2=getS('2_butterfly');
    figure; blf(S2,'iterations',3,'diff',1)
    vprintf('A figure should show two 2d-blf.\n\n','cpr',[.1 .4 .1]);
    
    T=transitionmatrix(S);
    
    %from here on the functions are called in lexicographic ordering
    
    [a,amin]=characteristic([1 0; 2 0; 1 -1]','amin',[-1; -1])
    if(~isequal(a,[0 0;0 0;1 1;0 1]) || ~isequal(amin,[-1;-1])); vprintf('characteristic() 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=constructdigits( [2 0; 1 2] );
    if(~isequal(d,[0 0 1 1;0 1 1 2])); vprintf('constructdigits() is broken.\n\n','cpr','err'); end;
    
    mu=constructmu(2,2);
    if(~isequal(mu,[2 1 0;0 1 2])); vprintf('constructmu() is broken.\n\n','cpr','err'); end;
    
    Om=constructOmega('1_1133','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
    V=constructV([0 1 2 3],1);
    
    Vbar=constructVbar([0 1 2 3]);
    
    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;
    
    %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;    
    
    mask2symbol([1 2 3; 4 5 6]);
    
    Sn=normalizeS(S,'equal');
    
    

    
end


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