function valTab = testDLplot(inn, makePlot, verbose)
% testDLplot      Show some dictionary properties and results on testdata,
% Three plots of development during learning may be plotted (if available
% and if 'makePlot' is true. An array of some properties is returned.
% Also atoms may be plottet for a 64x256 dictionary.
%
% (This file is in test-mode, and quite messy)
% For dictionary properties rather use: dictprop.m
% but if SRC is needed this is the file to use (as dictprop does not
% include this special sparse representation capability measure)
% 
% valTab = testDLplot(res, makePlot, verbose)
% testDLplot(res);    % dictionary and other fields in struct res
% testDLplot(fn);     % struct res should be stored in this file
%----------------------------------------------------------------------
% res      : output from testDL_*.m
% fn       : a filename, ex: fn = 'testDL_RLS_May111630.mat';
% makePlot : if 0 (false) plots are not made
% verbose  : if 0 (false) nothing is displayed
% valTab   : array with elements: propD.A, propD.B, propD.mumse, propD.mu, 
%            propD.muavg, propD.mugap, , , SRC, achieved PSNR, 
%----------------------------------------------------------------------
% ex:
% res.D = get64x256();
% testDLplot(res,0);   % only field 'D', and no plot

%----------------------------------------------------------------------
% Karl Skretting. 
% University of Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  20.01.2011  KS: test started
% Ver. 1.1  25.01.2011  KS: minor changes
% Ver. 1.2  10.02.2011  KS: may give input as a filename
% Ver. 1.3  18.04.2011  KS: updated output a little bit
% Ver. 1.4  02.09.2011  KS: search catalogs for mat-file
%----------------------------------------------------------------------

mfile = 'testDLplot';
% It may be convenient to store the learned dictionaries in some
% subcatalogs, thus to be able to find them possible locations should be
% listed in the cell array below.
catalogs = {pwd,[pwd,'/d3'],[pwd,'/d1']};

valTab = zeros(1,12);
if nargin < 2
    makePlot = ispc();
end
if nargin < 3
    verbose = true;
end

if ischar(inn)
    for j=1:numel(catalogs)
        filename = [catalogs{j},'/',inn];
        if exist(filename,'file'); break; end;
    end
    if exist(filename, 'file')
        load(filename); 
    else
        disp([mfile,': does not find ',inn]); 
        return;
    end
elseif isstruct(inn)
    res = inn;
end
if ~exist('res', 'var') 
    disp([mfile,': variable res was not given as expected.']); 
    return;
end
if ~isfield(res, 'D')
    disp([mfile,': input argument has no field for the dictionary, no res.D.']); 
    return
end        

if verbose
    disp(' ');
    if isfield(res, 'inputParam')
        disp([mfile,' plot results generated by ',res.inputParam.mfile,...
            ' ',res.inputParam.vecSel,' ',res.startTime]);
    else
        disp([mfile,' plot results for dictionary in res.D ']);
        res.inputParam.dataSet = 'A';
    end
    if isfield(res, 'resultFile')
        disp([' and stored in ',res.resultFile,' ',res.endTime]);
    end
    if isfield(res, 'vectorsProcessed') && isfield(res, 'totalTime')
        hour = floor(res.totalTime/3600);
        minu = floor(res.totalTime/60) - hour*60;
        sec = floor(res.totalTime+0.5) - hour*3600 - minu*60;
        disp([' ',sprintf('%6g',res.vectorsProcessed),' vectors processed during ',...
            sprintf('%02i.%02i.%02i',hour,minu,sec),' (hh.mm.ss) of learning.']);
        fprintf(' which is %6.1f vectors per second ', res.vectorsProcessed/res.totalTime);
        if isfield(res,'computer')
            fprintf('on a %s computer.', res.computer);
        end
        if isfield(res.inputParam,'useMOD')
            fprintf(' (useMOD)\n');
        else
            fprintf('\n');
        end
    end
end
%
if ~isfield(res, 'inputParam') 
    if size(res.D,1) == 64
        res.inputParam.dataSet = 'A';
    end
end
if strcmpi(res.inputParam.dataSet,'lena')
    res.inputParam.dataSet = 'A';
end

% dictionary properties
[N,K] = size(res.D);
propD = dictprop(res.D);
if verbose
    fprintf(' Dictionary is %ix%i and frame bounds are A: %6.4f  B: %6.2f and betamin: %5.2f  betaavg: %5.2f \n', ...
        N, K, propD.A, propD.B, propD.betamin, propD.betaavg );
    fprintf(' betagap >= %5.2f  and betamse: %5.2f  mumse: %6.4f  \n', ...
        propD.betagap, propD.betamse, propD.mumse );
end
valTab = zeros(1,12);   % values that may be in output table
valTab(1:6) = [ propD.A, propD.B, propD.mumse, propD.mu, propD.muavg, propD.mugap ];

%%  test the final dictionary in a standard way
testDataFile = ['testDLdata',res.inputParam.dataSet,'.mat'];  
if res.inputParam.dataSet == 'A'
    tPSNR = 35;  % results in actual PSNR in range 36.00 to 36.30
elseif res.inputParam.dataSet == 'B'
    tPSNR = 28;  % results in actual PSNR in range 29.10 to 29.20
elseif res.inputParam.dataSet == 'C'
    tPSNR = 28;  % results in actual PSNR in range 29.40 to 29.60
end    
errLim = 255*sqrt(N)*10^(-tPSNR/20);  % ORMP exit when  ||r||_2 < errLim 
sparseLim = min(floor(N/2),100);
Xtest = load(testDataFile);
if isfield(res.inputParam, 'transform') && strcmpi(res.inputParam.transform, 'dct')
    Xtest = ks_dct(Xtest.X);
else
    Xtest = Xtest.X;
end
avg_x0 = mean(sum(Xtest ~= 0));
avg_x1 = mean(sum(abs(Xtest)));
avg_x2 = mean(sqrt(sum(Xtest.*Xtest)));
avg_w0 = 0; avg_w1 = 0; avg_w2 = 0;
avg_r0 = 0; avg_r1 = 0; avg_r2 = 0; avg_rr = 0;

if exist('mexOMP.m','file')
    param.paramOMP = struct('eps',errLim.^2, 'L',sparseLim  );
    if verbose
        fprintf(' Test dictionary on X in %s using mexOMP with sparseLim (L) = %i and epsilon = %6.2f \n', ...
            testDataFile, sparseLim, errLim.^2 );
    end
    W = mexOMP(Xtest, res.D, param.paramOMP);   % use mexOMP in test also for mexLasso in learning
    R = full(Xtest - res.D*W);
    avg_w0 = full(mean(sum(W ~= 0))); 
    avg_w1 = full(mean(sum(abs(W)))); 
    avg_w2 = full(mean(sqrt(sum(W.*W)))); 
    avg_r0 = full(mean(sum(R ~= 0)));
    avg_r1 = full(mean(sum(abs(R))));
    avg_r2 = full(mean(sqrt(sum(R.*R)))); 
    avg_rr = full(mean(sum(R.*R))); 
    clear W R
    %
else
    if verbose
        fprintf(' Test dictionary on X in %s using javaORMP with sparseLim = %i and errLim = %6.2f \n', ...
            testDataFile, sparseLim, errLim );
    end
    java_access();
    jD = mpv2.SimpleMatrix(res.D);
    jDD = mpv2.SymmetricMatrix(size(res.D, 2), size(res.D, 2));
    jDD.eqInnerProductMatrix(jD);
    jMP = mpv2.MatchingPursuit(jD, jDD);
    jMP.setNormalized();
    %
    for j = 1:size(Xtest,2);    % for each element
        x = Xtest(:,j);
        if sqrt(x'*x) < errLim
            w = zeros(K,1);
        else
            w = jMP.vsORMP(x, sparseLim, errLim/sqrt(x'*x));
        end
        r = x - res.D*w;
        rr = r'*r;
        avg_w0 = avg_w0 + sum(w ~= 0);
        avg_w1 = avg_w1 + sum(abs(w));
        avg_w2 = avg_w2 + sqrt(w'*w);
        avg_r0 = avg_r0 + sum(r ~= 0);
        avg_r1 = avg_r1 + sum(abs(r));
        avg_r2 = avg_r2 + sqrt(rr);
        avg_rr = avg_rr + rr;
    end
    avg_w0 = avg_w0/size(Xtest,2);
    avg_w1 = avg_w1/size(Xtest,2);
    avg_w2 = avg_w2/size(Xtest,2);
    avg_r0 = avg_r0/size(Xtest,2);
    avg_r1 = avg_r1/size(Xtest,2);
    avg_r2 = avg_r2/size(Xtest,2);
    avg_rr = avg_rr/size(Xtest,2);
end
if verbose
    fprintf(' Mean values for testdata      0-norm: %6.2f  1-norm: %7.2f  2-norm: %7.2f \n', ...
        avg_x0, avg_x1, avg_x2 );
    fprintf(' Mean values for coefficients  0-norm: %6.2f  1-norm: %7.2f  2-norm: %7.2f   (SRC: %5.2f)\n', ...
        avg_w0, avg_w1, avg_w2, avg_w0);
    fprintf(' Mean values for errors        0-norm: %6.2f  1-norm: %7.2f  2-norm: %7.2f  (PSNR: %5.2f)\n', ...
        avg_r0, avg_r1, avg_r2, 10*log10((255^2)/(avg_rr/N)) );
end
valTab(9) = avg_w0;
valTab(10) = 10*log10((255^2)/(avg_rr/N));


%% make figures if ispc

if makePlot && isfield(res,'fbA')
    figure(1); clf;
    x = (0:(numel(res.fbA)-1));
    subplot(221);
    plot(x,res.fbA);
    if isfield(res,'resultFile')
        title(['  ',res.resultFile,' :'],'Interpreter','none');
        ylabel('lower A');
    else
        title('Lower frame bound A');
    end
    textFinalValue('%6.4f', res.fbA(end));
    subplot(222);
    plot(x,res.fbB);
    if isfield(res,'resultFile')
        title('Frame bounds');
        ylabel('upper B');
    else
        title('Upper frame bound B');
    end
    textFinalValue('%6.2f', res.fbB(end));
    subplot(223);
    if isfield(res, 'mu')
        plot(x,res.mu,'r-');
        textFinalValue('%5.2f', res.mu(end));
    elseif isfield(res, 'betamin')
        plot(x,cos(res.betamin*pi/180),'r-');
        textFinalValue('%5.2f', cos(res.betamin(end)*pi/180));
    end
    title('\mu = cos \beta_m_i_n');
    xlabel('Main iteration number');
    % ylabel('Angle in degrees.');
    subplot(224);
    % several names have been used for this property, mse is latest
    if isfield(res, 'mumse')
        plot(x,res.mumse,'r-');
        textFinalValue('%5.2f', res.mumse(end));
    elseif isfield(res, 'mueig')
        plot(x,res.mueig,'r-');
        textFinalValue('%5.2f', res.mueig(end));
    elseif isfield(res, 'betamse')
        plot(x,cos(res.betamse*pi/180),'b-');
        textFinalValue('%5.2f', cos(res.betamse(end)*pi/180));
    elseif isfield(res, 'betaeig')
        plot(x,cos(res.betaeig*pi/180),'b-');
        textFinalValue('%5.2f', cos(res.betaeig(end)*pi/180));
    end
    title('\mu_m_s_e = cos \beta_m_s_e ');
    xlabel(['each with ',int2str(res.inputParam.dataLength),' vectors.']);
    % title('betamin (red) and betaavg (blue) for dictioary during training');
    
    if isfield(res, 'delta0') && isfield(res, 'deltai')
        valTab(7) = res.delta0(end);
        valTab(8) = res.deltai(end);
        figure(2); clf;
        if (numel(res.delta0) < 40)  % 'circle' plot
            xy = zeros(numel(res.delta0)+1, 2);
            xy(2,1) = res.delta0(1);
            th = 0;  % radians
            for i=2:numel(res.delta0)
                r1 = res.delta0(i-1);
                r2 = res.delta0(i);
                d = res.deltai(i);
                th = th + acos( (r1*r1+r2*r2-d*d)/(2*r1*r2) );
                xy(i+1,:) = r2*[cos(th), sin(th)];
            end
            plot(xy(:,1), xy(:,2), 'b-o');
            axis equal;
            if isfield(res,'resultFile')
                title([res.resultFile,' : dictionary difference.'],'Interpreter','none')
            else
                title('Change (Frobenius norm) in D_i relative to D_i_-_1 and D_0');
            end
            V = axis();
            h = text(0, (V(4)-V(3))/25, 'D_0');
            set(h, 'HorizontalAlignment', 'center');
            td = ['D',subscriptInFig(int2str(numel(res.delta0)))];
            h = text(xy(end,1)+(V(2)-V(1))/40, xy(end,2), td);
            set(h, 'VerticalAlignment', 'middle');
            t = ['norm( ',td,' - D_0 ) = ',sprintf('%6.2f',res.delta0(end))];
            xlabel(t);
        else
            i1 = 1:floor(0.5*numel(res.delta0));
            i2 = i1(end):numel(res.delta0);
            b1 = ( (abs(max(res.delta0(i1))-max(res.delta0(i2))) + ...
                    abs(min(res.delta0(i1))-min(res.delta0(i2)))) < ...
                   0.5*(max(res.delta0)-min(res.delta0)) );
            b2 = ( (abs(max(res.deltai(i1))-max(res.deltai(i2))) + ...
                    abs(min(res.deltai(i1))-min(res.deltai(i2)))) < ...
                   0.5*(max(res.deltai)-min(res.deltai)) );
            % b1 = (numel(i1) > 5); b2 = b1;   % always true ?   
            if b1  % one plot for each value
                subplot(211)
                plot(1:numel(res.delta0), res.delta0);
                title('Distance between D_i and D_0');
                ylabel('Angle in degrees.');
                t = ['|D',subscriptInFig(int2str(numel(res.delta0))),...
                    sprintf(' - D_0| = %5.2f',res.delta0(end))];
                textFinalValue(t, res.delta0(end));
            else                % divide x-axis in two plots
                subplot(221)
                plot(i1, abs(res.delta0(i1)));
                title('Distance between D_i and D_0');
                ylabel('Angle in degrees.');
                V = axis; V(2)=i1(end); axis(V);
                subplot(222)
                plot(i2, res.delta0(i2));
                t = ['|D',subscriptInFig(int2str(numel(res.delta0))),...
                    sprintf(' - D_0| = %5.2f',res.delta0(end))];
                title(t);
                V = axis; V(1)=i2(1); axis(V);
            end
            if b2 %  one plot for each value
                subplot(212)
                plot(1:numel(res.deltai), res.deltai);
                title('Distance between D_i and D_i_-_1');
                t = ['|D',subscriptInFig(int2str(numel(res.deltai))),...
                    ' - D',subscriptInFig(int2str(numel(res.deltai)-1)),...
                     sprintf('| = %6.4f',res.deltai(end))];
                textFinalValue(t, res.deltai(end));
                if isfield(res,'resultFile')
                    xlabel(res.resultFile,'Interpreter','none')
                end
            else
                subplot(223)
                plot(i1, abs(res.deltai(i1)));
                V = axis; V(2)=i1(end); axis(V);
                title('Distance between D_i and D_i_-_1');
                ylabel('Angle in degrees.');
                if isfield(res,'resultFile')
                    xlabel(res.resultFile,'Interpreter','none')
                end
                subplot(224)
                plot(i2, res.deltai(i2));
                V = axis; V(1)=i2(1); axis(V);
                t = ['|D',subscriptInFig(int2str(numel(res.deltai))),...
                     ' - D',subscriptInFig(int2str(numel(res.deltai)-1)),...
                     sprintf('| = %6.4f',res.deltai(end))];
                title(t);
            end
        end
    elseif isfield(res, 'deltai')  % plot only deltai
        valTab(8) = res.deltai(end);
        figure(2); clf;
        plot(1:numel(res.deltai), abs(res.deltai));
        V = axis; V(1)=1; V(2)=numel(res.deltai); axis(V);
        title('Distance between D_i and D_i_-_1');
        ylabel('Angle in degrees.');
        if isfield(res,'resultFile')
            xlabel(res.resultFile,'Interpreter','none')
        end
    end
    
    if isfield(res, 'avg_w0') || isfield(res, 'avg_r2') || isfield(res, 'avg_rr')
        figure(3); clf;
        if isfield(res, 'avg_w0') && (isfield(res, 'avg_r2') || isfield(res, 'avg_rr'))
            subplot(211)
        end
        if isfield(res, 'avg_w0')
            valTab(11) = res.avg_w0(end);
            x = 1:numel(res.avg_w0);
            plot(x, res.avg_w0);
            if isfield(res,'resultFile')
                title([res.resultFile,' : mean sparseness.'],'Interpreter','none')
            else
                title('Average of ||w||_0 for each main iteration');
            end
            textFinalValue('Final ||w||_0 : %5.2f', res.avg_w0(end));
        end
        if isfield(res, 'avg_w0') && (isfield(res, 'avg_r2') || isfield(res, 'avg_rr'))
            subplot(212);
        end
        if isfield(res, 'avg_rr')
            x = 1:numel(res.avg_rr);
            y = 10*log10( size(res.D,1)*(res.inputParam.dataPeak^2) ./ res.avg_rr ); % PSNR
            plot(x, y);
            title(sprintf('PSNR for each main iteration, target = %5.2f.',res.inputParam.targetPSNR));
            textFinalValue('Final PSNR : %5.2f', y(end));
            valTab(12) = y(end);
        elseif isfield(res, 'avg_r2')
            x = 1:numel(res.avg_r2);
            plot(x, res.avg_r2);
            title('Average of ||r||_2 for each main iteration');
            textFinalValue('Final ||r||_2 : %5.2f', res.avg_r2(end));
        end
        if res.inputParam.useTestData
            xlabel('The test data set is used.');
        else
            xlabel('Training data is used.');
        end
    end
else   % not Plot, only set valTab
    if isfield(res, 'delta0') && isfield(res, 'deltai')
        valTab(7) = res.delta0(end);
        valTab(8) = res.deltai(end);
    end
    if isfield(res, 'avg_w0')
        valTab(11) = res.avg_w0(end);
    end
    if isfield(res, 'avg_rr')
        y = 10*log10( size(res.D,1)*(res.inputParam.dataPeak^2) ./ res.avg_rr ); % PSNR
        valTab(12) = y(end);
    end
end

if verbose 
    if isfield(res,'resultFile')
        t1  = res.resultFile((end-12):(end-4));
    else
        t1 = 'given D';
    end
    fprintf('\n File          Frame bounds               mu-prop               dictdiff(...)     TestData        Learning \n');
    fprintf('                A        B      mse     mu      avg     gap    to D0   to prev   SRC     PSNR    w0      PSNR \n');
    fprintf(' %9s & %6.4f & %6.2f & %5.2f & %5.2f & %5.2f & %5.2f & %5.2f & %6.4f & %5.2f & %5.2f & %5.2f & %5.2f  \n\n', ...
             t1, valTab );
    %
    % may display atoms in figure 4
    if makePlot && (N == 64) && (K == 256) 
        D = res.D;
        dv = zeros(1,K);
        dh = zeros(1,K);
        dd = zeros(1,K);
        for i=1:K
            atom = D(:,i);
            ma = max(abs(atom));
            atom = reshape(0.5+0.499*atom/ma, 8, 8);
            for i1=2:8; for i2=1:8; dv(i) = dv(i) + abs(atom(i1,i2)-atom(i1-1,i2)); end; end;
            for i1=1:8; for i2=2:8; dh(i) = dh(i) + abs(atom(i1,i2)-atom(i1,i2-1)); end; end;
            for i1=2:8; for i2=2:8; dd(i) = dd(i) + abs(atom(i1,i2)-atom(i1-1,i2-1)) + abs(atom(i1,i2-1)-atom(i1-1,i2)); end; end;
        end
        [t,I] = sort(dd+dv+dh);   % sort by 'total' differences.
        m1 = 12; m2 = ceil(K/m1);
        % sort more
        [t,Iv] = sort(dv);   % sort 'vertically
        for i1=1:m1
            i3 = ((i1-1)*m2+1):min(K,i1*m2);
            [t,Ih] = sort(dh(Iv(i3)));
            I(i3) = Iv(i3(Ih));
        end
        D = D(:,I);
        
        A = ones(9*m1+1, 9*m2+1);
        i = 1;
        for i1=1:m1
            for i2=1:m2
                atom = D(:,i);
                ma = max(abs(atom));
                atom = reshape(0.5+0.499*atom/ma, 8, 8);
                A( (9*i1-7):(9*i1), (9*i2-7):(9*i2) ) = atom;
                i = i+1;
                if (i>K); break; end;
            end
        end
        
        figure(4); clf;
        imagesc(A);
        colormap(gray);
        axis equal;
        axis off;
        title( sprintf('Atoms in dictionary %s. Properties: A=%6.2f, B=%6.2f, mumse=%6.4f, muavg=%6.4f, mu=%6.4f, SRC=%5.2f', ...
            t1, valTab([1:3,5,4,9])), 'Interpreter','none' );
    end
    
end

return

function t2 = subscriptInFig(t1)
t2 = blanks(2*numel(t1));
for i=1:numel(t1)
    t2(2*i-1) = '_';
    t2(2*i) = t1(i);
end
return;

function textFinalValue(t1, v)
V = axis();
if (v > (0.5*V(3)+0.5*V(4)))
    y = v - 0.12*(V(4)-V(3));
else
    y = v + 0.1*(V(4)-V(3));
end
text( 0.05*V(1)+0.95*V(2), y, sprintf(t1, v), 'HorizontalAlignment', 'right');
return

% Alternatively errLim and sparseLim may be found from input parameters
% function [errLim, sparseLim] = limitsFromInput(res)
% if isfield(res.inputParam, 'sparseNum'); 
%     sparseLim = res.inputParam.sparseNum;
% elseif isfield(res.inputParam, 'paramOMP'); 
%     sparseLim = res.inputParam.paramOMP.L;
% else
%     sparseLim = floor(size(res.D, 1)/2);
% end
% if isfield(res.inputParam, 'absErrorLimit'); 
%     errLim = res.inputParam.absErrorLimit;
% elseif isfield(res.inputParam, 'paramOMP'); 
%     errLim = sqrt( res.inputParam.paramOMP.eps );
% elseif isfield(res.inputParam, 'targetMSE');
%     errLim = sqrt( res.inputParam.targetMSE );
% else
%     errLim = 50;
% end
%    
