function [ A, oo ] = testnormdecay(varargin)
% testnormdecay( A, delta, [options] )
% Constructs sequence of matrices with maximal spectral radius below a given threshold.
% More specific: Constructs a sequence (i_n)_n s.t. \ro(A_{i_k}\cdots A_{i_1})^{1/k}\leq\delta \forall k\in\NN
%
% Input:
%   A               Cell array, matrices to use. If empty matrices are random.
%
% Options:
%   'J',val         Number of matrices to generate (if M is not given), default=2
%   'rho',val       Spectral radius of matrices to generate (if M is not given), default=.9
%   'dim',val       Dimension (if M is not given), default=2
%   'RNG',val       Random number generator seed
%   'delta',val     Threshold as described above. If empty, delta=1.
%   'l',val         Maximal length of sequence to generate. default=30
%   'verbose',val   Verbose level
%   'start',val     start matrix
%
% Output:
%   A               Matrices used
%   oo              Sequence generated
%   nice pictures 
%
% E.g.: testnormdecay('delta',.9999,'rho',.9,'l',200,'RNG',82159266)
%
% Written by: tommsch, 2018

[l,varargin]=parsem('l',varargin,30);
[J,varargin]=parsem('J',varargin,2);
[RNG,varargin]=parsem('RNG',varargin,0); if(RNG); rng(RNG); else; RNG=randi(100000000); rng(RNG); end;
[rhoval,varargin]=parsem('rho',varargin,0.9);
[dim,varargin]=parsem('dim',varargin,3);
[delta,varargin]=parsem('delta',varargin,1);
[verbose,varargin]=parsem('verbose',varargin,1);

if(isequal(length(varargin),0) || isequal(length(varargin),1) && isempty(varargin{1}));
    if(isequal(length(varargin),1)); varargin(1)=[]; end;
    A=cell(1,J); for j=1:J; A{j}=randn(dim); A{j}=A{j}/rho(A{j})*rhoval; end;
    vprintf('Matrices: \n%v\n',A,'imp',[1 verbose]);
else
    A=varargin{1};
    varargin(1)=[];
    dim=size(A{1},1);
    if(any(cellfun(@rho, A)>delta)); 
        vprintf('Spectral radius of matrices bigger then delta.\n','cpr','err','imp',[1 verbose]); 
    end;
end

[X,varargin]=parsem('start',varargin,eye(dim));

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


NORM=zeros(1,0);
RHO=zeros(1,0);
oo=zeros(1,0);
%TJSR=tjsr(A,'delta',.9999,'maxtime',120,'verbose',1);
%REFERENCE=zeros(1,0); %the decay in the 1-d case

i=0;
while(true)
    i=i+1;
    
    Xnew=cellfun(@(x) X*x,A,'UniformOutput',0);
    rhonew=cellfun(@rho, Xnew).^(1/i);
    rhonew(rhonew>delta)=-1;
    [val,idx]=max(rhonew);
    if(i>l); break; end;
    if(val<0); break; end;
    
    X=X*A{idx};
    NORM(i)=norm(X);
    RHO(i)=rho(X);
    oo(i)=idx;
	%REFERENCE(i)=TJSR^i; 
end
%REFERENCE(1)=[];

if(~isempty(oo));
    oo=findperiod(oo);

    clf;
    %semilogy(RHO,'r.-'); hold on;
    %semilogy(RHO.^(1./(1:(i-1))),'r.-'); hold on;
    val=RHO.^(1./(1:(i-1)));
    semilogy([liminf(val);limsup(val)].','r.-'); hold on;
    %semilogy(NORM,'k.-');
    %semilogy(NORM.^(1./(1:(i-1))),'k.-');
    val=NORM.^(1./(1:(i-1)));
    semilogy([liminf(val);limsup(val)].','k.-');
    %semilogy(REFERENCE,'g.-'); 
    vprintf('Product: { %v , %v }\n', oo{1}, oo{2} ,'imp',[1 verbose]);
    vprintf('Productlength: %i\n', i-1 ,'imp',[1 verbose]);
    %legend('rho','rho-normalized','norm','norm-normalized');
end

end


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