function D = getNxK(N, K, par)
% getNxK          Get a NxK dictionary (frame) with small coherence, use L-BFGS.
%                 Start by an initial random dictionary, or a dictionary
% given as input. Then do some iterations (gradient descent) to improve
% an objective function related to coherence. Also see sedil(7).
%
%  use: D = getNxK(N, K, par);  
%       D = getNxK(D0, par);  
%       D = getNxK(N, K);  
%       D = getNxK(D0);  
%
%  examples: 
%  D0 = get6x16() + 0.1*randn(6,16);   % initial dictionary
%  D = getNxK(D0);  
%  D = getNxK(6, 16);
%  par = struct('maxIter',1000, 'memSize',5, 'gLim',1e-4, 'dfLim',0,'step',1);
%  D = getNxK(6, 16, par);
%  p = dictprop(D); dictshow(D,'plottype',3);
%  D2 = getNxK(6,16);  % D1 = get6x16();        % often find this
%  D2 = getNxK(16,144);  % D1 = get16x144();    % but not this
%  D1 = getNxK(25,100); D2 = getNxK(D1);

if ((nargin == 3) && isstruct(par))
    %  par = par;
elseif ((nargin == 2) && isstruct(K)) 
    par = K;
else
    par = struct([]);
end
if isfield(par,'maxIter'); maxIter = par.maxIter; else maxIter = 1000; end;
if isfield(par,'memSize'); memSize = par.memSize; else memSize = 5; end;
if isfield(par,'gLim'); gLim = par.gLim; else gLim = 1e-4; end;
if isfield(par,'dfLim'); dfLim = par.dfLim; else dfLim = 0; end;
if isfield(par,'step'); step = par.step; else step = 1; end;

if (nargin == 1) || ((nargin == 2) && isstruct(K))
    D = N;
    [N,K] = size(D);
elseif ((nargin == 2) && ~isstruct(K) && (numel(K) == 1)) || (nargin == 3)
    D = randn(N,K);
else
    disp('getNxK: should have 1 or 2 or 3 input arguments, see help.');
    D = [];
    return;
end

fprintf('%s: Size of D is %3ix%3i\n', mfilename(), N, K);

if memSize > 0
    s = zeros(N*K,memSize);
    y = zeros(N*K,memSize);
    rho = zeros(1,memSize);   % rho in L-BFGS as described in Wikipedia.
end

fPrev = inf;
d = zeros(N,K);
for iter = 1:maxIter
    if (memSize > 0) && (iter > 1)
        [D, f, g] = fObj( D + step*d );    
    else
        [D, f, g, h] = fObj( D + step*d );    
        if (f > fPrev)
            step = 0.8*step;
        end
    end
    df = fPrev - f;
    g1 = norm(g(:),1); 
    verbose = ((iter == 1) || (iter == maxIter) || (rem(iter,ceil(maxIter/25)) == 0));
    if verbose || (g1 < gLim) || (abs(df) < dfLim)  % display results 
        fprintf('Iter %4i: f = %10.4f, ||g|| = %10.4f,  df = %10.4f, step = %5.3f  \n', ...
            iter, f, g1, df, step );
        if ((g1 < gLim) || (abs(df) < dfLim)); 
            break; 
        end
    end
    %
    if (memSize > 0) && (iter > 1)
        s = [D(:) - xk, s(:,1:(memSize-1))];
        y = [g(:) - gk, y(:,1:(memSize-1))];
        rho = [1/(y(:,1)'*s(:,1)), rho(1:(memSize-1))];
        xk = D(:); 
        gk = g(:);       % keep previous values
        % and find (uphill) direction to next point using L-BFGS
        d = -reshape(findHgk(iter-1,memSize,gk,s,y,rho), N, K);
        d = d - repmat(sum(D.*d),N,1).*D;    % 'orthogonalize' direction
    else
        % d = -g./h;            % update direction
        d = -g/mean(h(:));   % update direction  opposite of gradient
        if (memSize > 0)
            xk = D(:); 
            gk = g(:);       % keep previous values
        end
    end
    fPrev = f; 
end

% D = dictnormalize(D,2);

return


% Function fObj(.) finds one term of f in Hawe (16), i.e. the term r(A). 
% Note that input A may be 'any point' in Euclidian space R^{NxK}.
% Normalizing the dictionary corresponds to projecting the point onto 
% a smooth Riemannian manifold M in Hawe (17), a unit ball for each column.
% The function f is found as in Hawe (7), the mask M has K(K-1)/2 non-zeros. 
% The gradient g is found by differentiating f. Tested above, sedil(2).
% Orthogonalizing g is "the orthogonal projection of the gradient onto the
% tanget space T_{Y}M", Hawe (18). Below g is first the gradient in
% Euclidian space and then the gradient in tangent space.
% The diagonal of Hessian matrix, hess, is not used in Hawe. 
% The formula is checked/verified by sedil(6).
function [A,f,g,hess] = fObj(A)
[N,K] = size(A);
A = A ./ repmat(sqrt(sum(A.^2)), N, 1);   % normalized dictionary
if nargout > 1
    C = A.' * A;                          % inner products
    M = ~tril(true(K));                   % mask or indicator matrix
    f = -sum( log(1 - C(M).^2) ); % C(M) is colum vector of relevant entries 
end
if nargout > 2
    C(1:(K+1):K^2) = 0;    % set diagonal to zeros, KEEP lower left part
    g = 2 * A * (C./(1-C.^2));
    % absg = abs(g);
    g = g - repmat(sum(A.*g),N,1).*A;     % orthogonalize g to A
    % gratio = abs(g) ./ absg;
end
if nargout > 3   
    % hess should be the diagonal elements of the Hessian (nxn matrix)
    % represent A as a vector of length n = N*K
    % then g is also a vector of length n, g(i) = df/dA(i)
    % and hess is also a vector of length n, h(i) = dg(i)/dA(i) = ddf/ddA(i)
    % hess is reshaped to NxK and returned.
    AA = A.^2;      
    CC = C.^2;   
    hess = 4*AA*(CC./((CC-1).^2)) - 2*AA*(ones(K)./(CC-1)) - 2*AA;  % ok
end 
return


% A simple implementation of L-BFGS as described in Wikipedia: 
% https://en.wikipedia.org/wiki/Limited-memory_BFGS
function z = findHgk(k,m,gk,s,y,rho)
    % k,m    are integers, k is iteration number, m is memory size
    %        first time called with k = 1, initial step gave x0 (and g0)
    %        and thus in step k=1 we only have one element in s and y
    %        k=1: s(:,1) = s0 = x1 - x0, and y(:,1) = y0 = g1 - g0
    % gk     gradient in step k, size is nx1, g_{k}
    % s,y    memory, each is nxm, s(:,i) = s_{k-i} = x_{k+1-i} - x_{k-i}
    %        y(:,i) = y_{k-i} = g_{k+1-i} - g_{k-i},  i = 1,2,..,m
    % rho    1xm, rho(i) = 1/(y(:,i)' * s(:,i))
    b = min(m,k);    % bound
    alpha = zeros(1,b);
    z = gk;   % nx1   (z is first q in wikipedia)
    for i = 1:b
        alpha(i) = rho(i)*(s(:,i)'*z);
        z = z - alpha(i)*y(:,i);
    end
    Hk0 = (s(:,1)'*y(:,1))/(y(:,1)'*y(:,1));
    z = Hk0*z;
    for i = b:(-1):1
        beta = rho(i)*(y(:,i)'*z);
        z = z + (alpha(i)-beta)*s(:,i);
    end
return
