% test double preconditioning
% by XXL (Peter Balazs)
% contact: xxl@kfs.oeaw.ac.at
% double preconditioning algorithm at lines 180-200
% for documentation see 
% [xxl] Peter Balazs, 'Irregular And Regular Gabor Multipliers With 
%       Application To Psychoacutsic Masking'
%
% Known Problems: Some Strings are repeated because of the usage of sprintf
% and %s. This is necessary for Layout of this help, but seems to make some
% problems.


iter0 = 10; % number of different lattice parameters
iter1 = 1; % number of different tests (= different windows) for every 
            %set of parameter
t0 = cputime; % for measuring calculation time

if exist('circessbet') == 1
    % incorporate past experiment in this one
    if show == 1
        disp('Using old data to get bigger basis.');
    end;
    n0 = n0+iter0*iter1;
else
	if exist('show') ~= 1
        show = 0; % print output? Or be quite (==0)
        % use it outside to set it
    end;

	maxn  = 1000; % maximal signal length
    
    % initiate variables
	diagsmal = 0;
	noframe = 0;
	diagesssmal = 0;
	diagnonconv = 0;
	diagmat = 0;
	diagmaterr = 0;
	diagessbet = 0;
	rrd = maxn+1;
    rrde = maxn+1;
    rrnd = maxn+1;
	nodouble = 0;
    diagsupp = 0;
    doubleessbet = 0;
    essnodouble = 0;
	
	circsmal = 0;
	circesssmal = 0;
	circnonconv = 0;
	circmat = 0;
	circmaterr = 0;
	noiter = 0;
	circessbet = 0;
	rrc = maxn+1;
    rrce = maxn+1;
    n0 = iter0*iter1;
    
    methstr = ...
        'First choosing b, then a and supp. So condition on a and supp.';
    winstr = 'Hamming'; % change with window !!!!

end;
    
for jj = 1:iter0
     n = 1; %reset
     nAlph = 1;
     % make sure that we don't choose a prime number (or one)
     while isprime(n) | n == 1 
         n = ceil(maxn*rand(1,1));
         % choose a random n in [1,maxn]         
     end;
  
     alph = alldiv(n); % divisors or n
     nAlph = length(alph);
   
     % choose a support length (now done later)
     % supp = alph(ceil(nAlph*rand(1,1)));
     % supp = n;
     % nAlph = length(find(alph <= supp));
     % a <= nAlpg % else no frame is possible.
     % choose arbitrary lattice parameter 

     b = alph(ceil(nAlph*rand(1,1)));
	 M = n/b;
	
	 % choose b such that a*b <= n, otherwise no frame is possible
	
     naaa = length(find(alph <= M));
	 a = alph(ceil(naaa*rand(1,1)));
	 red = n/(a*b); %redundancy
   	 N = n/a;
     
     % choose a support length
     % bigger than M such that we don't have an diagonal frame operator
     % alsu = alph(find(alph >= M));
     % bigger than a such that we don't have a frame
%     alsu = alph(find(alph > M));

     nsss = nAlph-naaa;
     if nsss == 0
         % supp <= M
         % see tschurtschenthaler or qiu
         % in this case it must be diagonal
      
         % We could not avoid it
         % either cound it
         if 1 == 1
            diagmat = diagmat+1;
            diagsupp = diagsupp+1;
            diagsmal = diagsmal+1;
            if show == 1
                disp('The operator is diagonal.');
            end;
        else
            n0 = n0-1;
        end
         % or just continue
            continue;
   else
        supp = alph(naaa+ceil(nsss*rand(1,1)));        
   end;
                               
    if show == 1
        disp(sprintf(...
            'n = %g, a = %g, b = %g, red = %g; supp = %g',n,a,b,red,supp));
    end;
        
    for ii = 1:iter1 
        % use this inner loop only if the window can change for the 
        % same n,a,b and supp for example with random windows.
        g = zeros(1,n);
        
        % (1) random window (white noise)
        w = rand(1,supp);
        % g(1:(floor(supp/2)+1))=w(ceil(supp/2):supp);
        % g(n-ceil(supp/2)+2:n)=w(1:(ceil(supp/2)-1));

        % (2) hanning window
        % w = hanning(supp+2);
        % g(1:(floor(supp/2)+1))=w(ceil(supp/2)+1:supp+1);
        % g(n-ceil(supp/2)+2:n)=w(2:ceil(supp/2));
        % winstr='Hanning'
        % change above !!!!

        % (2b) hamming window
        % winstr = 'Hamming';
        % change above
        % w = hamming(supp);
        g(1:(floor(supp/2)+1))=w(ceil(supp/2):supp);      
        g(n-ceil(supp/2)+2:n)=w(1:ceil(supp/2)-1);
        %        g(M+2) = 1     
        % in this case use iter1 = 1
        % blackman(n)
        % kaiser(n,6)
      
        % (3) full Gaussian
        % w = gaussnk(supp);
        % g(1:(floor(supp/2)))=w(1:floor(supp/2));
        % g(n-ceil(supp/2)+1:n)=w(floor(supp/2)+1:supp);
        % in this case use iter1 = 1
        
        %% (4) cut-off Gaussian
        % supp2 = ceil(maxn*rand(1,1));
        % w = gaussnk(supp2)
        % g(1:(floor(supp/2)+1))=w(supp2-ceil(supp/2):supp2); (??)
        % g(n-ceil(supp/2)+2:n)=w(1:(ceil(supp/2)-1));
        % in this case use iter1 = 1
        %% !!! Use Kaiser-Bessel with beta < 1 instead ?
        
              
        B = blockxxl(g,a,b); % b x a block matrix
        
        % mde = min(abs(B(1,:)));
        if min(abs(B(1,:))) == 0 
            % disp(sprintf('Minimal diagonal entry = %g',mde));
            if show == 1
                disp(sprintf('No frame'));
            end;
            noframe = noframe+1;
            continue;
        end;
                          
        % Id=eye(n);
        Id=zeros(b,a);
        Id(1,:)=ones(1,a);

        % Single Preconditioning: diagonal
        P1=zeros(b,a);
        P1(1,:)=1./(B(1,:));   
        % Block-Matrix with inverse diagonal elements
        D=blockm(B,P1,n);      
        % Preconditioning by multiplication of blocks
        diagnorm = walnorm(D-Id); 
        % 'distance' of preconditioned matrix to the identiy

        % Single Preconditioning: circulant
        v=mean(B.');      % mean of rows of B
        w=ifft(oneover(fft(v))); % deconvolution 
        P1=w.'*ones(1,a);
        A=blockm(B,P1,n);         
        circnorm = walnorm(A-Id);
        
        % Double Preconditioning
        v=mean(D.'); 
        w=ifft(oneover(fft(v)));
        P1=w.'*ones(1,a);
        C=blockm(D,P1,n); 
        doublenorm = walnorm(C-Id);

        % full n x n matrix
        % S = blo2walxxl(B,n);
        % rS = rank(S);
        rS = n;
        
        if rS ~= n
            % n elements are linear dependent, only create a space 
            % with dimension = rank, so cannot be a frame.
            noframe = noframe+1;
            if show == 1
                disp(sprintf('No frame'));
            end;
            continue;
        else
            if show == 1
                disp(sprintf(...
                 'Rank of S: %g, Diagonal Dominance Norm: %g, %s %g; %s : %g', ...
                 '  Circulant Norm:', '  Double Preconditioning', ...
                 rS, diagnorm, circnorm, doublenorm));
            end;
        end;
        if doublenorm >= 0.99
            if diagnorm >= 0.99 & circnorm >=0.99
                % to prevent calculation erros
                noiter = noiter+1;
                if show == 1
                    disp(...
                        'The system is a frame, but no iterative scheme converges.');
                end;
            else
                nodouble = nodouble+1;
                if show == 1
                   disp('-----------------------');
                   disp('!!! Bad Moon Rising !!!');
                   disp('-----------------------');
                end;
                     aand = a;
                     bbnd = b;
                     ggnd = g;
                     rrnd = min(rrnd,red);
                 %    return;
                 % end if this bad case occurs
                 if diagnorm <= 0.9 | circnorm <= 0.9
                     essnodouble = essnodouble+1;
                 end;
            end;
        else
            if show == 1
               disp('Double Preconditioning is convergent!!');
            end;
            if diagnorm < doublenorm & diagnorm < circnorm
                diagsmal = diagsmal+1;
                if diagnorm <= 10*eps
                    diagmat = diagmat+1;
                    diagmaterr = max(diagmaterr,doublenorm);
                else
                    relerr = diagnorm/doublenorm;
                    %if relerr > 0.01
                    if relerr < 0.1
                         diagessbet = diagessbet+1;
                        % if show == 1
                             disp('Diagnorm is essentially smaller');
                             %end;
                         aaed = a;
                         bbed = b;
                         gged = g;
                         rred = min(rrd,red);
                     end;
                     if relerr > diagesssmal
                        diagesssmal = relerr;
		%                     if show == 1
		%                         disp('Diagnorm is smaller');
		%                     end;
		%                     aad = a;
		%                     bbd = b;
		%                     ggd = g;
		%                     rrd = min(rrd,red):
                    end;  
                end;
            elseif circnorm < doublenorm
                circsmal = circsmal+1;
                if circnorm >= 1
                    circnonconv = circnonconv+1;
                elseif circnorm <= 10*eps
                    circmat = circmat+1;
                    circmaterr = max(circmaterr,doublenorm);
                else
                  relerr = circnorm/doublenorm;
                  %  if relerr > 0.01
                  if relerr < 0.1
                        circessbet = circessbet+1;
                         %if show == 1
                             disp('Circnorm is essentially smaller');
                             %end;
                         aace = a;
                         bbce = b;
                         ggce = g;
                         rrce = min(rrd,red);
                    end;
                    if relerr > circesssmal
                        circesssmal = relerr;
		%                      aac = a
		%                      bbc = b
		%                      ggc = g
		%                      rrc = min(rrc,red)
                    end;
                end;
            elseif 10*doublenorm < diagnorm & 10*doublenorm < circnorm
                        doubleessbet = doubleessbet+1;
                         % if show == 1
                             disp(sprintf(...
                   'Doublenorm is essentially smaller, a = %g, b=%g',a,b));
                             % end;
                         adde = a;
                         bdde = b;
                         gdde = g;
                         rdde = min(rrd,red);
            end;
		end;
    end;
end;
t1 = cputime;
disp(sprintf('time = %g s',t1-t0));
disp(methstr);
disp(sprintf('We have tested %g times using a %s window',n0,winstr));
disp(sprintf(...
    'In %g cases (%g percent) the Gabor system was no frame.',noframe,...
    noframe/n0*100));
disp(sprintf(...
    'In %g cases (%g percent) no norm was smaller than 1.',noiter,...
    noiter/n0*100));
disp(sprintf(...
    '%g times the diagonal norm was smaller than the double precond norm',...
    diagsmal));
disp(sprintf('1.) In %g cases the Frame Matrix is already diagonal', ...
    diagmat));
disp(sprintf(...
    '    and the difference is only due to precision errors (%g).',...
    diagmaterr));
disp(sprintf('    In %g cases this is due to the support property.',...
    diagsupp));
disp(sprintf(...
 '2.) the rest (%g cases): here the maximal difference was (relative) %g.',...
  diagsmal-diagmat,diagesssmal));
disp(sprintf('Only in %g cases the difference was bigger than 1 percent.',...
    diagessbet));
% if rrde ~= maxn+1
%     disp(sprintf(...
%    'The minimal redundancy for which the diagonal norm 
%     was essential smaller was %g.',rrde));
% end;
% disp('---------------------------------------------------------------');
disp(sprintf(...
    '%g times the circulant norm was smaller than the double precond norm',...
    circsmal));
disp(sprintf(...
    '1.) In %g cases the Frame Matrix is already circulant',circmat));
disp(sprintf(...
    '    and the difference is only due to precision errors (%g).',...
    circmaterr));
disp(sprintf(...
    '2.) the rest (%g cases): here the maximal difference was (relative) %g.',...
    circsmal-circmat,circesssmal));
disp(sprintf(...
    'Only in %g cases the difference was bigger than 1 percent.',...
    circessbet));
% if rrc ~= maxn+1
%   disp(sprintf(...
%     'The minimal redundancy for the circulant case was %g, for which 
%       the error was 
%        essential in the circulant case was %g.',rrc,rrce));
%    if rrce ~= maxn+1
%         disp(sprintf('The minimal redundancy for which the error was essential 
%        in the circulant case was %g.',rrce));  
%    end;
% end;
if nodouble ~= 0
    disp(sprintf('*ATTENTION* There was %g case, %s',...
     'when a single preconditiner would have converged, but the double didnt',...
     nodouble));
    disp(sprintf('But only in %g cases, the smaller norm was < 0.9.',...
        essnodouble));
end;
disp(sprintf('Heureka! In %g cases the double preconditioning norm %s',...
    'was essentially smaller (factor:10!).',doubleessbet)); 

