function Ds = texdictlearn(varargin)
% texdictlearn     Learn a dictionary to represent a texture by a minibatch RLS-DLA variant. 
%                  i.e. make a model of the texture. The textures are available
% through the m-file teximage.m (see that one). Dictionary can also be stored 
% in a catalog given in variable 'cat_D' which must exist, or in a given
% named file ('dictfile' option and 'saveDs' option equal 3)
% 
%   Use of function:
%   ----------------
%   Ds = TEXDICTLEARN(varargin)
%     Ds is the dictionary as a struct, Ds.D is the actual dictionary
%     Arguments may be given as pairs: 'option',value, ...
%     (or options may be given in one (or more) struct or cell-array)
% 
%   Options for training texture, i.e. how to make the training vectors: 
%   --------------------------------------------------------
%   im, imno, or imfile   input image as used in teximage.m 
%   imrange   may be [-1,1] (default), [0,1], [-128,127], [0,255]
%       the image is scaled to be within this range, and imrange is then
%       set to the actual range of the input image (may be narrower)
%   vecMet    giving the neighborhood for each pixel as im makevector.m 'a'-'j'
%
%   Options for sparse approximation: 
%   --------------------------------------------------------
%   samet    the sparse approximation method to use. It can be the method
%       argument in sparseapprox.m. Default 'javaormp' 
%   saopt    the additional options to use in sparseapprox.m as a struct. 
%       Default is: struct('tnz',3, 'verbose',0); 
%   K   the number of atoms in the dictionary, default 100
%
%   Options for minibatch variant of RLS-DLA: 
%   --------------------------------------------------------
%   mb, minibatch   an array with number of batches and number of vectors in each batch.
%       Default: repmat([1000,1],5,1).*[2,10; 2,40; 3,100; 3,200; 2,500];
%       Total number of training vectors is: sum( minibatch(:,1).*minibatch(:,2) )
%   lam0, lambda0   The initial value of the forgetting factor lambda, default 0.996.
%   lam1, lambda1   when lambda should be increased to 1, given as a number relative
%       to the scheduled number of vectors to process, default 0.9.
%   outlim   limit for outlayers, limit for approximation error relative to
%       mean of the last errors. Default: inf (which is no limit)
%   dlim     limit for range of dictionary atoms. Only relevant when
%       the dictionary is not normalized, ||d_k||_2 == 1. Using dlim = [a,b] 
%       mean that  a <= D(n,k) <= b. Default dlim = [];   i.e. no limit
%   D    initial dictionary, default: random training vectors
%   A    initial A matrix, default: identity matrix
%
%   Options for mat-files to store dictionary in (or retrieve from): 
%   --------------------------------------------------------
%   dictin, optionfile   name of mat-file where options are stored.
%       This can be a dictionary file.
%   saveDs   0: the dictionary (Ds) is not saved in a mat-file. 
%       1: Ds is saved in a mat-file, name given as Dnnxxxx, 
%       2: Ds is saved in a mat-file, name given as Dnnxxxx-ddhhmm, 
%          nn is imno, xxxx is vecMet plus two chars given by samet, ddhhmm is date-time
%          if mat-file exists it will not be overwritten and (the old already
%          stored) Ds will be returned (no learning done).
%       3: save Ds in file given by dictfile, overwrite if dictfile exists 
%   dictfile  file to store Ds in if saveDs == 3
%
%   Options for logging and displaying results: 
%   --------------------------------------------------------
%   ptc, PropertiesToCheck   a cell array with names of dictionary properties to chech
%       at regular intervals during learning given by checkrate. None: {}.
%       Default: {{'mtvp','r2avg','traceA','fbA','fbB','betamse','d2min'}}
%       May also use: 'nofr', 'rcondA', 'nofo'
%       Note {{ and }} used to define a cell array as a field in a struct!
%       It is good to have 'mtvp' (or 'tvp') as first property (used as x-axis in plots)
%   checkrate   how often properties are reported and dictionary checked, default 1000
%   v, verbose  verbose level given as 0 (only warnings) 1 (normal) 2 (much)
%
%   Examples:
%   ---------
%   Ds = texdictlearn('imno',16);     % use default values
%   Ds = texdictlearn('imno',16, 'saveDs',1);     % and store Ds in a mat-file
%   Ds = texdictlearn('imno',16,'samet','cafixed','v',2);
%   opt = struct('verbose',2,'checkrate',1000,'samet','javaormp','dlim',2,'outlim',2);
%   % processing only 0.2 mill training vectors
%   opt.minibatch = repmat([100,1],5,1).*[1,20; 2,40; 3,100; 3,200; 2,500]

%--------------------------------------------------------------------------
% Copyright (c) 2012.  Karl Skretting.  All rights reserved.
% University of Stavanger, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  06.01.2012  Made function (based on maketm.m)
% Ver. 1.1  12.01.2012  Added several options
% Ver. 1.2  16.02.2012  Added exception handling in try block (line 590 ...)
% Ver. 1.3  25.06.2012  Add samet 'convexsets' and 'greedy'
%           31.07.2012  and corrected some errors introduced in ver. 1.3
% Ver. 1.4  08.08.2012  Add samet 'afformp', no normalization of D!
% Ver. 1.5  14.01.2013  Some minor (mainly help text) updates
% Ver. 2.0  10.11.2014  Added options to use mapping gMap, see makeGmap.m
%                       remove convex approximation methods: cafixed, convexsets, greedy
% Ver. 2.1  14.12.2016  Removed options to use mapping gMap, see version 2.0
%--------------------------------------------------------------------------

mfile = 'TexDictLearn';

%% defaults, initial values
tstart = tic;
text = [mfile,' started ',datestr(now)];  

cat_train = my_matlab_path('textures');
cat_D = my_matlab_path('.d\jan13');   % default catalog to store results in
if isunix
    cat_sep = '/';
else
    cat_sep = '\';
end

imno = 0;
imrange = [-1,1];
vecMet = 'd';
samet = 'javaORMP';
K = 100;
fac2K = 1000;   % only used to set minibatch below
minibatch = repmat([fac2K,1],5,1).*[2,10; 2,40; 3,100; 3,200; 2,500];  
lam0 = 0.996;
lam1 = 0.9;
tvp = 0;    % number of training vectors processed
outnum = 0; 
outlim = inf;
dlim = [];
D = [];
A = [];
saveDs = 0;
ptc = {'mtvp','r2avg','traceA','fbA','fbB','betamse','d2min'};
checkrate = 1000;
verbose = 1;
downscaling_factor = 0.95;   % scale down to be (well) inside leagal region, given by dlim

%% get the options
nofOptions = nargin;
optionNumber = 1;
fieldNumber = 1;
while (optionNumber <= nofOptions)
    if isstruct(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        sNames = fieldnames(sOptions);
        opName = sNames{fieldNumber};
        opVal = sOptions.(opName);
        % next option is next field or next (pair of) arguments
        fieldNumber = fieldNumber + 1;  % next field
        if (fieldNumber > numel(sNames)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    elseif iscell(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        opName = sOptions{fieldNumber};
        opVal = sOptions{fieldNumber+1};
        % next option is next pair in cell or next (pair of) arguments
        fieldNumber = fieldNumber + 2;  % next pair in cell
        if (fieldNumber > numel(sOptions)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    else
        opName = varargin{optionNumber};
        opVal = varargin{optionNumber+1};
        optionNumber = optionNumber + 2;  % next pair of options
    end
    % interpret opName and opVal
    if ( strcmpi(opName,'im') || strcmpi(opName,'image') || ...
         strcmpi(opName,'imfile') || strcmpi(opName,'imagefile') || ...
         strcmpi(opName,'imname') || strcmpi(opName,'imagename') || ...
         strcmpi(opName,'imno') || strcmpi(opName,'imnumber') || ...
         strcmpi(opName,'imageno') || strcmpi(opName,'imagenumber') )
        if ischar(opVal) || isnumeric(opVal)
            im = opVal;
        else
            error([mfile,': illegal type (class) of value for option ',opName]);
        end
    end
    if ( strcmpi(opName,'imrange') )
        if isnumeric(opVal) && (numel(imrange)==2)
            imrange = opVal;
            % perhaps the following tests should be removed? 
            if ((-128 <= imrange(1)) && (imrange(1) < -65)); imrange(1) = -128; end;
            if ((-1 <= imrange(1)) && (imrange(1) < -0.6)); imrange(1) = -1; end;
            if (0 <= imrange(1)); imrange(1) = 0; end;
            if ((130 < imrange(2)) && (imrange(2) <= 255)); imrange(2) = 255; end;
            if ((0.6 < imrange(2)) && (imrange(2) <= 1)); imrange(2) = 1; end;
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'vecMet') || strcmpi(opName,'letter')   % 
        if ischar(opVal)
            vecMet = lower(opVal(1));
            if (sum(vecMet == 'abcdefghij') == 0); vecMet = 'd'; end;
        elseif isnumeric(opVal) || islogical(opVal)
            vecMet = opVal;
        else
            error([mfile,': not char for option ',opName]);
        end
    end 
    %
    if ( strcmpi(opName,'samet') )
        if ischar(opVal)
            samet = opVal;
        else
            error([mfile,': not character value (string) for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'saopt') )
        if isstruct(opVal)
            saopt = opVal;
        else
            error([mfile,': not struct for option ',opName]);
        end
    end 
    if strcmpi(opName,'K')    
        if isnumeric(opVal)
            K = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    %
    if ( strcmpi(opName,'minibatch') || strcmpi(opName,'mb') )
        if isnumeric(opVal)
            minibatch = opVal;
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'lam0') || strcmpi(opName,'lambda0') )
        if isnumeric(opVal)
            lam0 = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'lam1') || strcmpi(opName,'lambda1') )
        if isnumeric(opVal)
            lam1 = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'r2avg') || strcmpi(opName,'r2average') )
        if isnumeric(opVal)
            r2avg = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'outlim') || strcmpi(opName,'outlimit') )
        if isnumeric(opVal)
            outlim = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'outnum') || strcmpi(opName,'outnumber') )
        if isnumeric(opVal)
            outnum = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'dlim') || strcmpi(opName,'dlimit') )
        if isnumeric(opVal)
            if (numel(opVal) == 1)
                dlim = [-opVal(1),opVal(1)];
            elseif (numel(opVal) > 1)
                dlim = [min(opVal(:)),max(opVal(:))];
            end
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'D') 
        if isnumeric(opVal)
            D = opVal;
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'A') 
        if isnumeric(opVal)
            A = opVal;
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    %
    if strcmpi(opName,'saveDs')
        if (islogical(opVal) && opVal); saveDs = 1; end;
        if isnumeric(opVal); saveDs = opVal(1); end;
    end
    if strcmpi(opName,'dictin') || strcmpi(opName,'optionfile')
        if ischar(opVal)
            dictin = opVal;
            if ((numel(dictin) < 5) || ~strcmpi(dictin((end-3):end),'.mat'))
                dictin = [dictin,'.mat']; %#ok<AGROW>
            end
            if (~exist(dictin, 'file') &&  (numel(strfind(dictin,cat_sep)) == 0))
                dictin = [cat_D,cat_sep,dictin];   %#ok<AGROW>
            end
            if exist(dictin, 'file')
                disp([mfile,': use all fields in ',dictin,' as options.']);
                % but some few should be ignored
                mfileHERE = mfile;
                tstartHERE = tstart;
                textHERE = text;
                clear text
                load(dictin);
                mfile = mfileHERE;  
                tstart = tstartHERE;
                if exist('text','var')
                    text = char(text, textHERE);
                else
                    text = textHERE;
                end
                dictfile = dictin;
                % the more careful code would be like
                %  Ds = load(dictin);  
                %  if isfield(Ds,'imno'); imno = Ds.imno; end;
                %  if isfield(Ds,'imfile'); imfile = Ds.imfile; end;
                %  ....
            else
                error([mfile,': can not find ',dictin]);
            end
        else
            error([mfile,': not char for option ',opName]);
        end
    end
    if strcmpi(opName,'dictfile') || strcmpi(opName,'dictout')
        if ischar(opVal)
            dictfile = opVal;
        else
            error([mfile,': not char for option ',opName]);
        end
    end
    %
    if strcmpi(opName,'tvp') 
        if isnumeric(opVal)
            tvp = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'text') 
        if ischar(opVal)
            text = char(opVal, text);
        else
            error([mfile,': not char for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'PropertiesToCheck') || strcmpi(opName,'ptc') )
        if iscell(opVal)
            ptc = opVal;
        else
            error([mfile,': not cell for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'PropertiesTable') || strcmpi(opName,'ptab') )
        if isnumeric(opVal)
            ptab = opVal;
            pti = size(ptab,1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'checkrate') || strcmpi(opName,'cr') )
        if isnumeric(opVal)
            checkrate = opVal(1);
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
end

%% some variables are checked and some are derived from input
K1 = floor(sqrt(K));    % faktorize K so that K == K1*K2
K2 = K/K1;
while rem(K2,1) > 0
    K1 = K1-1;
    K2 = K/K1;
end

%% set filename to store Ds into
if (saveDs == 1) || (saveDs == 2)
    if strcmpi(samet(1),'j')   % java implemented method
        if strcmpi(samet,'javaMP') || strcmpi(samet,'javaMatchingPursuit') || ...
                strcmpi(samet,'javaBMP') || strcmpi(samet,'javaBasicMatchingPursuit')
            xx = 'jM';
        elseif strcmpi(samet,'javaOMP') || strcmpi(samet,'javaOrthogonalMatchingPursuit')
            xx = 'jo';
        elseif strcmpi(samet,'javaORMP') || strcmpi(samet,'javaOrderRecursiveMatchingPursuit')
            xx = 'jO';
        elseif strcmpi(samet,'javaPS') || strcmpi(samet,'javaPartialSearch')
            xx = 'jP';
        else
            xx = 'jx';
        end
    elseif strcmpi(samet(1),'m')   % mex (SPAMS) implemented method
        if strcmpi(samet,'mexLasso') || strcmpi(samet,'LARS') || strcmpi(samet,'LASSO') 
            xx = 'mL';
        elseif strcmpi(samet,'mexOMP')
            xx = 'mO';
        else
            xx = 'mx';
        end
    else
        xx = samet(1:min(2,length(samet)));
    end
    xx = [xx,int2str(saopt.tnz)];
    if isnumeric(im) && (numel(im) == 1)
        imno = floor(im);
    end
    if ischar(vecMet)
        dictfile = sprintf('D%02i%s%s',imno,vecMet,xx);
    elseif isnumeric(vecMet)
        dictfile = sprintf('D%02iv%02i%s',imno,vecMet(1),xx);
    else
        dictfile = sprintf('D%02i%s%s',imno,'x',xx);
    end
    if (saveDs == 2)
        c = clock;
        dictfile = sprintf('%s-%02i%02i%02i',dictfile,c(3),c(4),c(5));
    end
    dictfile = [cat_D,cat_sep,dictfile,'.mat'];  
    if exist(dictfile,'file')
        disp([mfile,': dictionary ',dictfile,' exist. Load this and return.']);
        Ds = load(dictfile);
        return;
    end
elseif (saveDs == 3)
    if ((numel(dictfile) < 5) || ~strcmpi(dictfile((end-3):end),'.mat'))
        dictfile = [dictfile,'.mat']; 
    end
    if (numel(strfind(dictfile,'/')) == 0) && (numel(strfind(dictfile,'\')) == 0)
        dictfile = [cat_D,cat_sep,dictfile];   
    end
    if exist(dictfile,'file')
        disp([mfile,': dictionary ',dictfile,' exist. It will be overwritten.']);
    end
else   % saveDs == 0
    dictfile = 'not saved';
end
if verbose >= 1
    disp([mfile,': saveDs = ',int2str(saveDs),', samet = ',samet, ...
        ', dictfile = ''',dictfile,'''.']);
end

%% load training image
[trainingImage, imfile, imtext] = teximage(im, imrange);

%% initialize variables, including dictionary
StartOfLineDisplayed = false;    % used when (verbose == 1)
imrange = [min(trainingImage(:)), max(trainingImage(:))];
rows = 5:(size(trainingImage,1)-4);   % we will not use the edges of training image.
cols = 5:(size(trainingImage,2)-4);
L1 = min(50,numel(rows));
L2 = min(60,numel(cols));    
L = L1*L2; %  is number of tv from the image each time.
if ischar(vecMet)
    Nabcdefghij = [9,25,29,13,21,25,17,9,13,73];          % N for the letters
    N = Nabcdefghij(double(vecMet)-96); 
else
    N = size(makevector(ones(49),vecMet,25,25,0),1);
end
%
if (size(D,1) == N) && (size(D,2) == K)
    % D is probably ok
else
    temp = randperm(numel(rows));
    row = rows(temp(1:K1));
    temp = randperm(numel(cols));
    col = cols(temp(1:K2));
    D = makevector(trainingImage, vecMet, row, col, 0);
end
% normalization of dictionary
g = 1./sqrt(sum( D.*D ));
D = D .* repmat(g, N, 1); % normalize dictionary

if (numel(A) == 1)
    B = D/A;   % note A is scalar here
    A = eye(K)*A;
end
if ~((size(A,1) == K) && (size(A,2) == K)) 
    A = eye(K);
    B = D;
end
tvtot = tvp + sum( minibatch(:,1).*minibatch(:,2) );
if ~exist('r2avg','var')
    temp = randperm(numel(rows));
    row = rows(temp(1:K1));
    temp = randperm(numel(cols));
    col = cols(temp(1:K2));
    X = makevector(trainingImage, vecMet, row, col, 0);
    W = sparseapprox(X, D, samet, saopt);
    R = X - D*W;
    r2avg = mean( sqrt(sum(R.*R)) );
end
X = [];
if ~exist('pti','var');
    pti = 0;      % last written line in ptab
end
if exist('ptab','var') && (size(ptab,2) == numel(ptc)) ;
    ptab = [ptab; zeros(ceil(tvtot/checkrate)+1, numel(ptc))];
else
    pti = 0;      % last written line in ptab
    ptab = zeros(ceil(tvtot/checkrate)+1, numel(ptc));
end
ptno = tvp;   % tvp when last line in ptab was written
if (numel(ptc) > 0)   % report some results as learning goes by
    pti = pti + 1;
    ptab(pti,:) = dictionaryProperties(D, ptc, tvp, r2avg, trace(A), 0, rcond(A), 0);
end

%% Display info
if (verbose > 0)  % verbose
    disp(' ');
    disp([mfile,' started ',datestr(now)]);
end
if (verbose > 1)  % very verbose
    if exist('cat_D','var'); disp(['  catalog for dictionaries: ',cat_D]); end;
    if exist('cat_train','var'); disp(['  catalog for training images: ',cat_train]); end;
    if exist('cat_sep','var'); disp(['  catalog separator character: ''',cat_sep,'''.']); end;
    disp(['  imno = ',int2str(imno),' and imfile = ',imfile]);
    disp(['  imrange = [',num2str(imrange(1)),' ',num2str(imrange(2)),']']);
    disp(['  reading ',int2str(L),' training vectors each time.']);
    if ischar(vecMet)
        disp(['  method to use in ''makevector.m'' is given as ''',vecMet,''' --> N = ',int2str(N)]);
    elseif isnumeric(vecMet)
        disp(['  method to use in ''makevector.m'' is given as ',int2str(vecMet(1)),' --> N = ',int2str(N)]);
    end
    disp(['  number of dictionary atoms K = ',int2str(K)]);
    disp(['  thus size of dictionary D is ',int2str(size(D,1)),'x',int2str(size(D,2)),'.']);
    disp(['  method for sparse approximation samet = ',samet]);
    disp(['  number of atoms in each sparse representation = ',int2str(saopt.tnz)]);
    disp(['  number of mini-batches to do = ',int2str(sum( minibatch(:,1)))]);
    disp(['  millions training vectors (mtv) to process = ',sprintf('%6.2f',tvtot/1000000)]);
    disp(['  initially lambda = ',num2str(lam0),', it increases to 1 after ',...
           sprintf('%7.3f',lam1*tvtot/1000000),' mtv are processed.']);
    disp(['  outlayer limit, outlim = ',num2str(outlim)]);
    if (numel(dlim) == 2)
        disp(['  dictionary elements are restricted to be in range = [',...
            num2str(dlim(1)),' ',num2str(dlim(2)),'].']);
    end
    if (numel(ptc) > 0)
        disp(['  for each ',int2str(checkrate),' tv the following dictionary properties are reported:']);
        disp( ['  ',propertyhead(ptc, 10)] );
        disp( ['  ',propertyline(ptab(pti,:), ptc, 10)] );
    end
end

%%  minibatch variant of RLS-DLA
xno = 0;
for linje = 1:size(minibatch,1)
    for bno = 1:minibatch(linje,1)
        % if isnan(D(1,1))
        %     error(['D(1,1) is NaN, tvp=',int2str(tvp),', linje=',int2str(linje),', bno=',int2str(bno)]);
        % end
        batchsize = minibatch(linje,2);
        if (xno+batchsize) > size(X,2)        % get more training vectors
            temp = randperm(numel(rows));
            row = rows(temp(1:L1));
            temp = randperm(numel(cols));
            col = cols(temp(1:L2));
            X = makevector(trainingImage, vecMet, row, col, 0);
            xno = 0;
        end
        if (xno+batchsize) > size(X,2)        % still not enough for a minibatch ??
            error(['batchsize (',int2str(batchsize),') is too large.']);
        end
        %
        Xbatch = X(:,xno+(1:batchsize));
        xno = xno+batchsize;
        W = sparseapprox(Xbatch, D, samet, saopt);
        temp = nnz(isnan(W));
        if temp
            disp(' ');
            disp(['  (623) ',int2str(temp),' elements of W are NaN.  Should stop now.']);
            pause(60)
        end
        temp = nnz(isnan(D));
        if temp
            disp(' ');
            disp(['  (629) ',int2str(temp),' elements of D are NaN.  Should stop now.']);
            pause(60)
        end
        R = Xbatch - D*W;
        r2 = sqrt(sum(R.*R));   % 2-norm of approximation error
        if (outlim < inf)
            I = (r2 <= (outlim*r2avg));
            Xbatch = Xbatch(:,I);
            W = W(:,I);
            r2 = r2(I);
            batchsize = numel(r2);
            outnum = outnum + (numel(I)-nnz(I));
        end
        r2avg = 0.75*r2avg + 0.25*mean(r2);
        %
        lam = lambdafun(tvp, 'Q', lam1*tvtot, lam0, 1).^batchsize;
        A = lam*A + full(W*W');
        B = lam*B + full(Xbatch*W');
        try
            D = B/A;    % it may happen that A is ill conditioned
            temp = nnz(isnan(D));
            if temp
                disp(' ');
                disp(['  (612) ',int2str(temp),' elements of D are NaN.  Should stop now.']);
                pause(60);
            end
            % then a warning is displayed and no error thrown
            % the "normal" place to adjust dictionary and A is ca 80 lines below.
        catch ME
            if StartOfLineDisplayed
                disp('  '); 
                StartOfLineDisplayed = false;
            end
            disp(['  662: ',ME.message]);
            if strcmpi(samet((end-1):end),'MP') 
                % here we replace 'unused' atoms by atoms from training set
                d = diag(A);
                counter = 0;
                while (min(d) < 0.1*mean(d))  && (counter < min(20, size(Xbatch,2)))
                    counter = counter + 1;
                    [~, atom_no] = min(d);
                    D(:,atom_no) = Xbatch(:,counter); 
                    A(atom_no,:) = 0;
                    A(:,atom_no) = 0;
                    A(atom_no,atom_no) = 0.5*mean(d);
                    d = diag(A);
                end
                B = D*A;
                disp(['  replace ',int2str(counter),' atoms in D by training vectors.']);
            else % start with 'initial' values again
                temp = randperm(numel(rows));
                row = rows(temp(1:K1));
                temp = randperm(numel(cols));
                col = cols(temp(1:K2));
                D = makevector(trainingImage, vecMet, row, col, 0);
                A = eye(K);
                B = D;
                disp('  replace the whole dictionary D by training vectors.');
            end
        end
        % normalization of dictionary
        g = 1./sqrt(sum( D.*D ));
        D = D .* repmat(g, N, 1); % normalize dictionary, not A and B though

        tvp = tvp+batchsize;
        %
        if (numel(ptc) > 0) % report results (ptc: properties to check)
            if (tvp >= (ptno + checkrate))
                % 21.feb.2012: moved the report part after the check part
                pti = min(pti+1, size(ptab,1));
                ptno = tvp;   % tvp when last line in ptab was written
                %
                % the matrix A is also checked at checkrate
                d = diag(A);
                md = mean(d);
                if (min(d) < 0);   % should not happen, numerical error??
                    if StartOfLineDisplayed
                        disp('  '); 
                        StartOfLineDisplayed = false;
                    end
                    disp('  negative element in diagonal of A is an error, but will be corrected here.');
                    A = eye(K);   
                    B = D;
                    nof_replaced = K;
                    rcond_A = 1;
                elseif (min(d) < 0.01*md)
                    % small diagonal value mean an atom not much used
                    % adjust this in different ways depending on the used vector seletion method.
                    % display message only if (verbose > 1)
                    I = ((100*d) < md);
                    nof_replaced = nnz(I);
                    if (nof_replaced > 0)
                        if (verbose > 1) && numel(find(strcmpi(ptc, 'nofr')))
                            % fprintf('\n  %9.3f    replace %i atoms with training vectors.\n',tvp/1e6,nof_replaced);
                            fprintf('\n  %9.3f    replace %i atoms with other atoms.\n',tvp/1e6,nof_replaced);
                        end
                        [~,temp] = sort(d,'descend');
                        % D(:,I) = X(:, ceil(rand(1,nof_replaced)*size(X,2))); %#ok<AGROW>
                        D(:,I) = D(:,temp(1:nof_replaced)) + 0.1*randn(N,nof_replaced); 
                        A(:,I) = zeros(K,nof_replaced);
                        A(I,:) = zeros(nof_replaced,K);
                        A(I,I) = (0.5*md)*eye(nof_replaced);
                        g = 1./sqrt(sum( D.*D ));
                        D = D .* repmat(g, N, 1); % normalize dictionary
                    end
                    %  only when min(d) is small rcond(A) is checked
                    rcond_A = rcond(A);
                    if rcond_A < 1e-8
                        A = A +(0.0001*md)*eye(K);
                        if (verbose > 1) && (numel( find(strcmpi(ptc, 'rcondA')) ) == 0)
                            fprintf('  %9.3f    adjust A to reduce rcond(A) from %g to %g.\n',tvp/1e6, rcond_A, rcond(A));
                        end
                    end
                    %
                    B = D*A;
                    %
                else  % min(d) is large enough 
                    nof_replaced = 0;
                    rcond_A = rcond(A);
                end
                %
                if (numel(dlim) == 2)
                    temp = nnz(isnan(D));
                    if temp
                        disp(' ');
                        disp(['  (766) ',int2str(temp),' elements of D are NaN.  Should stop now.']);
                        pause(60)
                    end
                    % the 'new' method made 21.feb.2012
                    m = mean(dlim);
                    I = (D < dlim(1));
                    nof_oobe = nnz(I);   % oobe: out-of-bounds-elements
                    if nof_oobe > 0
                        D(I) = downscaling_factor*(dlim(1)-m) + m; 
                    end
                    I = (D > dlim(2));
                    if nnz(I) > 0
                        nof_oobe = nof_oobe + nnz(I);
                        D(I) = downscaling_factor*(dlim(2)-m) + m;
                    end
                    if nof_oobe > 0
                        B = D*A;
                    end
                    temp = nnz(isnan(D));
                    if temp
                        disp(' ');
                        disp(['  (784) ',int2str(temp),' elements of D are NaN.  Should stop now.']);
                        pause(60)
                    end
                else
                    nof_oobe = 0;
                end
                %
                % report to property table (logg),
                ptab(pti,:) = dictionaryProperties(D, ptc, tvp, r2avg, trace(A), ...
                    nof_replaced, rcond_A, nof_oobe);
                %
                if (verbose == 1)
                    % if you want to pay extra attention during learning it may
                    % be helful to see the value of rcond(A), this is often
                    % getting very close to zero (or infinity) when learning
                    % goes off in the wrong direction.
%                     if ~StartOfLineDisplayed
%                         fprintf('  Processing up to %7.3f mtv, mtvp = %7.3f (%7.3f)',...
%                             tvtot/1e6, tvp/1e6,log10(rcond_A));
%                         StartOfLineDisplayed = true;
%                     else
%                         fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%7.3f (%7.3f)',...
%                             tvp/1e6,log10(rcond_A));
%                     end
                    % but under normal condition it is just helpful to se number 
                    % of iterations  processed, to assure that learning is
                    % progressing as expected
                    if ~StartOfLineDisplayed
                        fprintf('  Processing up to %7.3f mtv, mtvp = %7.3f ',...
                            tvtot/1e6, tvp/1e6 );
                        StartOfLineDisplayed = true;
                    else
                        fprintf('\b\b\b\b\b\b\b\b%7.3f ', tvp/1e6);
                    end
                elseif (verbose > 1)
                    if (rem(pti,10)==0); disp( ['   ',propertyhead(ptc, 10)] ); end;
                    disp( ['  ',propertyline(ptab(pti,:), ptc, 10)] );
                end
                %
            end
        else  % nothing to report (or check! obs!) 
            % used when fast execution is wanted
        end
        %
    end
end
if (pti < size(ptab,1)) && (tvp > ptno)  
    pti = pti + 1;
    % ptno = tvp;    does not need this more
    ptab(pti,:) = dictionaryProperties(D, ptc, tvp, r2avg, trace(A), 0, rcond(A), 0);
end
if (pti < size(ptab,1))
    ptab = ptab(1:pti,:);
end

%% set output argument
tu = toc(tstart);
tuh = floor(tu/3600);
tum = floor(tu/60) - tuh*60;
tus = ceil(tu) - tum*60 - tuh*3600;
t1 = sprintf('%s finished %s, time used is %i:%02i:%02i (hh:mm:ss), %6.3f ms/tvp.', ...
         mfile, datestr(now), tuh, tum, tus, tu*1000/tvp );
text = char(text, t1);
Ds = struct( 'dictfile', dictfile, 'mfile', mfile, 'vecMet', vecMet, ...
             'imfile', imfile, 'imno', imno, 'imtext', imtext, 'imrange', imrange, ...
             'samet', samet, 'saopt', saopt, ...
             'lam0', lam0, 'lam1', lam1, ...
             'N', N, 'K', K, 'D', D, ...
             'r2avg', r2avg, ...
             'ptc', {ptc}, 'ptab', ptab, 'checkrate', checkrate, ...
             'minibatch', minibatch, 'tvp', tvp, ...
             'outnum', outnum, 'outlim', outlim, 'dlim', dlim, ...
             'timeused', tu, ...
             'text', text );
%
if (saveDs > 0)
    save(dictfile,'-struct','Ds');
end


%% finish and return

if (verbose > 0) 
    disp(' ');
    disp(t1);
    if (verbose > 1)  % very verbose
        disp(['  number of training vectors processed tvp = ',int2str(tvp)]);
        disp(['  number of training vectors rejected = ',int2str(outnum)]);
    end
    disp(' ');
end

return


function v = dictionaryProperties(D, propnames, tvp, r2avg, traceA, nofr, rcondA, nofo)
v = zeros(1,numel(propnames));
m = find(strcmpi(propnames, 'tvp'));
if (numel(m)==1); v(m) = tvp; end;
m = find(strcmpi(propnames, 'mtvp'));
if (numel(m)==1); v(m) = tvp/1e6; end;
m = find(strcmpi(propnames, 'r2avg'));
if (numel(m)==1); v(m) = r2avg; end;
m = find(strcmpi(propnames, 'traceA'));
if (numel(m)==1); v(m) = traceA; end;
m = find(strcmpi(propnames, 'nofr'));
if (numel(m)==1); v(m) = nofr; end;
m = find(strcmpi(propnames, 'rcondA'));
if (numel(m)==1); v(m) = rcondA; end;
m = find(strcmpi(propnames, 'nofo'));
if (numel(m)==1); v(m) = nofo; end;
if (sum(sum(isnan(D))) > 0)
    disp('  TexDictLearn.dictionaryProperties (894): Warning D is NaN.');
    return;
end
p = dictprop(dictnormalize(D), 0);
m = find(strcmpi(propnames, 'mu'));
if (numel(m)==1); v(m) = p.mu; end;
m = find(strcmpi(propnames, 'mumin'));
if (numel(m)==1); v(m) = p.mumin; end;
m = find(strcmpi(propnames, 'muavg'));
if (numel(m)==1); v(m) = p.muavg; end;
m = find(strcmpi(propnames, 'mumse'));
if (numel(m)==1); v(m) = p.mumse; end;
m = find(strcmpi(propnames, 'fbA'));
if (numel(m)==1); v(m) = p.A; end;
m = find(strcmpi(propnames, 'fbB'));
if (numel(m)==1); v(m) = p.B; end;
m = find(strcmpi(propnames, 'betamin'));
if (numel(m)==1); v(m) = p.betamin; end;
m = find(strcmpi(propnames, 'betaavg'));
if (numel(m)==1); v(m) = p.betaavg; end;
m = find(strcmpi(propnames, 'betamse'));
if (numel(m)==1); v(m) = p.betamse; end;
%
K = size(D,2);
if sum(strcmpi(propnames, 'd1min') | strcmpi(propnames, 'd1avg')) > 0
    % 1-norm distance between atoms 
    A = zeros(K);
    for k=1:K
        A(k,:) = sum(abs( repmat( D(:,k), 1, K) - D ));
    end
    m = find(strcmpi(propnames, 'd1avg'));
    if (numel(m)==1); v(m) = sum(A(:))/(K*(K-1)); end;
    A = A + max(A(:))*eye(K);
    m = find(strcmpi(propnames, 'd1min'));
    if (numel(m)==1); v(m) = min(A(:)); end;
end
if sum(strcmpi(propnames, 'd2min') | strcmpi(propnames, 'd2avg')) > 0
    % 2-norm distance between atoms 
    A = zeros(K);
    for k=1:K
        A(k,:) = sqrt(sum( ( repmat( D(:,k), 1, K) - D ).^2 ));
    end
    m = find(strcmpi(propnames, 'd2avg'));
    if (numel(m)==1); v(m) = sum(A(:))/(K*(K-1)); end;
    A = A + max(A(:))*eye(K);
    m = find(strcmpi(propnames, 'd2min'));
    if (numel(m)==1); v(m) = min(A(:)); end;
end

return

function s = propertyhead(propnames, len)
s = blanks(len*numel(propnames));
ii = 1;
for i=1:numel(propnames);
    i2 = floor((len-length(propnames{i}))/2);
    i1 = len-length(propnames{i})-i2;
    s(ii:(ii+len-1)) = [blanks(i1),propnames{i},blanks(i2)];
    ii = ii+len;
end
return

function s = propertyline(v, propnames, len)
mudes = '.4f ';
betades = '.2f ';
s = blanks(len*numel(propnames));
ls = int2str(len-1);
ls3 = int2str(len-3);
m = find(strcmpi(propnames, 'tvp'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'i '],v(m));
end
m = find(strcmpi(propnames, 'nofo'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls3,'i   '],v(m));
end
m = find(strcmpi(propnames, 'nofr'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls3,'i   '],v(m));
end
m = find(strcmpi(propnames, 'mtvp'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
m = find(strcmpi(propnames, 'r2avg'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
m = find(strcmpi(propnames, 'traceA'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.2e '],v(m));
end
m = find(strcmpi(propnames, 'rcondA'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.2e '],v(m));
end
m = find(strcmpi(propnames, 'mu'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,mudes],v(m));
end
m = find(strcmpi(propnames, 'mumin'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,mudes],v(m));
end
m = find(strcmpi(propnames, 'muavg'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,mudes],v(m));
end
m = find(strcmpi(propnames, 'mumse'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,mudes],v(m));
end
m = find(strcmpi(propnames, 'fbA'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.4f '],v(m));
end
m = find(strcmpi(propnames, 'fbB'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.1f '],v(m));
end
m = find(strcmpi(propnames, 'betamin'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,betades],v(m));
end
m = find(strcmpi(propnames, 'betaavg'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,betades],v(m));
end
m = find(strcmpi(propnames, 'betamse'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,betades],v(m));
end
m = find(strcmpi(propnames, 'd1min'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
m = find(strcmpi(propnames, 'd1avg'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
m = find(strcmpi(propnames, 'd2min'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
m = find(strcmpi(propnames, 'd2avg'));
if (numel(m)==1)
    ii = 1+(m-1)*len;
    s(ii:(ii+len-1)) = sprintf(['%',ls,'.3f '],v(m));
end
return

