function [ Mret, B ] = invariantsubspace(varargin)
% [ Mret, B ] = invariantsubspace( M, ['type'], [options] )
% Searches for invariant subspaces of matrices M.
%
% Input:
%   M       Cell array of matrices
%   'type'  A string of the following
%               'none'      Nothing happens,  Mret={M}, B=eye(dim)
%               'auto'      (default)
%               'perm'      Tries to find a permutation such that M is in block-diagonal form, using permTriangul from the JSR-toolbox by Jungers
%               'basis'     Tries to find a basis such that M is in block-diagonal form, using jointTriangul from the JSR-toolbox by Jungers
%               'trans'     If M are transition matrices, tries to find subspaces of differences U, as described in <'Smoothness of anisotropic wavelets, frames and subdivision schemes', M. Charina and V.Yu. Protasov> 
%                           First tries numerically, then symbolically, then with vpa.
%   If type is not given, then the algorithm tries in the following order: 'perm', 'basis', 'trans' (only numerically)
%
% Options:
%   'verbose',val           Verbose level, Default=0
%
% Output:
%   Mret    cell array of the blocks in the block diagonal. 
%   B       The change of basis. I.e. B'*M{i}*B has the block diagonal form
%
% E.g.: [M,B]=invariantsubspace({[1 1 1; 1 0 1; 0 1 0], [1 -1 0; 0 2 0; 1 1 2]},'verbose',2)
%
% See also: jointTriangul, permTriangul, constructU, restrictmatrix
%
% Written by: tommsch, 2018, and Jungers (JSR-toolbox)

%Rename: 'trans' to 'diff'


[verbose,varargin]=parsem('verbose',varargin,0);
M=varargin{1};
dim=size(M{1},2);
if(length(varargin)==1); type='auto'; else; type=varargin{2}; end;


info=identifymatrix(M);
if(isequal(type,'auto') && info.nonneg); 
    vprintf('Since the matrices are nonnegative, we only search permutations which triangulises. \n','imp',[2,verbose]);
    type='perm'; 
end;

EYE=eye(dim);
vprintf('Search for joint triangulizations: ','imp',[2,verbose])
triang=false;
switch type
    case 'none'
        B=EYE;
        Mret={M};    
        
    case 'perm'
        searchperm();
        
    case 'basis'
        searchbasis();
        
    case 'trans'
        searchtrans(1); %flag says: symbolically
        
    otherwise;
        searchperm();
        if(~triang); searchbasis(); end
        if(~triang); searchtrans(0); end; %flag says: not symbolically
end
vprintf('\n','imp',[1,verbose])


    function searchperm;
        vprintf('Search for permutation: ','imp',[1,verbose])
        [triang,Mret,B]=permTriangul(M);
        if(~triang);
            vprintf('no permutation found. ','imp',[1,verbose])
            B=EYE;
            Mret={M};
            triang=false;
        else
            vprintf('permutation found. ','imp',[1,verbose])
            B=EYE(B,:);
            triang=true;
        end;
    end

    function searchbasis;
        vprintf('Search for basis: ','imp',[1,verbose])
        [~,Mret,B]=jointTriangul(M);
        if(~isequal(B,eye(dim)));
            vprintf('basis found. ','imp',[1,verbose])
            triang=true;
        else
            vprintf('no basis found. ','imp',[1,verbose])
            triang=false;
        end      
        
    end

    function searchtrans(symflag);
        try
            vprintf('Search for difference subspace ','imp',[1,verbose]);
            if(symflag); vprintf('(symolically): ','imp',[1,verbose]); else; vprintf('(numerically): ','imp',[1,verbose]); end;
            U=constructU(M,1,'verbose',symflag);
            val=identifymatrix(M); 
            if(val.real); U=real(U); end;
            [MU,~,MR,NN,BASIS]=restrictmatrix(M,U,'verbose',verbose+symflag-2);
            if(symflag)
                if(max(max([NN{:}]))>10e-12); 
                    U=constructU(M,1,'sym','verbose',verbose+symflag-2);
                    [MU,~,MR,NN,BASIS]=restrictmatrix(M,U,'verbose',verbose+symflag-2);
                end;
                if(max(max([NN{:}]))>10e-12); 
                    U=constructU(M,1,'vpa',200);
                    [MU,~,MR,NN,BASIS]=restrictmatrix(M,U,'verbose',verbose+symflag-2);
                end;
            end
            if(max(max([NN{:}]))<10e-12)
                vprintf('difference subspace found. \n','imp',[1,verbose]);
                Mret{1}=MU; Mret{2}=MR; B=BASIS;
                triang=true;
            else
                vprintf('no difference subspace found. \n','imp',[1,verbose]);
                B=EYE;
                Mret={M};
                triang=false;
            end
        catch
            vprintf('could not search for difference subspace. \n','imp',[1,verbose]);
            B=EYE;
            Mret={M};
            triang=false;
        end
    end
        
end
    
function dummy; end %#ok<DEFNU> %Generates an error, if the 'end' of a function is missing.   
