function res = rbreginvunder(Amat, rhs, lambda, invR, blocks, varargin)
%
% res=rbreginvunder(Amat, rhs, lambda, invR)
%
% Perform forward simulations at all sources and all wavelengths based on the input structure

% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     Amat: the left-hand-side matrices (a containers.Map object) at specified wavelengths
%     rhs: the right-hand-side vectors for all sources (independent of wavelengths)
%     lambda: the Tikhonov regularization parameter
%     invR: the inversion of the upper triangular matrix of QR-decomposed L
%
% output:
%     res: the least-square solution of the matrix equation
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

% solve an overdetermined Gauss-Newton normal equation
%  delta_mu=inv(L'L)*J'*inv(J*(inv(L'L))J' + lambda*I)*(y-phi)
%
% see Deng2015, to accelerate the computaiton, we do not take the inversion
% of L'L, instead, we perform QR decomposition of L, so that L=Q*R, and
% pass invR=inv(R) matrix here, this allows us to compute
%    Z=J*invR
%    J*(inv(L'L))J' = Z*Z' -> Hess
%    inv(L'L)*J'    = invR*Z'
%
% so we have
%    delta_mu=(invR*Z')*inv(Z*Z' + lambda*I)*(y-phi)

Alen = size(Amat, 2);
idx = find(sum(Amat) ~= 0);
if (length(idx) < size(Amat, 2))
    Amat = Amat(:, idx);
end

emptydata = find(sum(Amat') ~= 0);
if (length(emptydata) < size(Amat, 1))
    Amat = Amat(emptydata, :);
    rhs = rhs(emptydata);
end

if (nargin >= 4 && ~isempty(invR))
    nx = size(invR, 1);
    oldnx = Alen / length(fieldnames(blocks));
    if (nx > length(idx) / length(fieldnames(blocks)))
        Lidx = idx(idx <= nx);
        invR = invR(Lidx, Lidx);
        oldnx = nx;
        nx = size(invR, 1);
    end
    if (nx == size(Amat, 2))
        Amat = Amat * invR;  % Z=J*invR
    else
        if (isempty(blocks)) % assume the Hess matrix size is multiples of LTL
            for i = 1:nx:size(Amat, 1)
                Amat(:, i:i + nx - 1) = Amat(:, i:i + nx - 1) * invR;
            end
        else
            %             len=cumsum([1; structfun(@(x) x(2), blocks)]);
            len = cumsum([1; structfun(@(x) x(2), blocks(1))]);
            len(2:end) = len(2:end) - (oldnx - nx) .* [1:length(len(2:end))]';
            for i = 1:length(fieldnames(blocks))
                %                 Amat(:,(i-1)*nx+1:i*nx) = Amat(:,(i-1)*nx+1:i*nx)*invR;
                if (nx == len(i + 1) - len(i))  % if the block size match LTL
                    Amat(:, len(i):len(i + 1) - 1) = Amat(:, len(i):len(i + 1) - 1) * invR;
                end
            end
        end
    end
end

rhs = rhs(:);

Hess = Amat * Amat'; % Gauss-Hessian matrix, approximation to Hessian (2nd order)

% [Hess,Gdiag]=rbnormalizediag(Hess);
Hess(1:1 + size(Hess, 1):end) = Hess(1:1 + size(Hess, 1):end) + lambda;
[Hess, Gdiag] = rbnormalizediag(Hess);

res = Gdiag(:) .* rbfemsolve(Hess, Gdiag(:) .* rhs, varargin{:});

if (nargin >= 4 && ~isempty(invR))
    nx = size(invR, 1);
    if (nx == size(Amat, 2))
        res = invR * (Amat' * res);
    else
        if (isempty(blocks)) % assume the Hess matrix size is multiples of LTL
            for i = 1:nx:size(Amat, 1)
                res(i:i + nx - 1) = invR * (Amat(:, i:i + nx - 1)' * res(i:i + nx - 1));
            end
        else
            len = cumsum([1; structfun(@(x) x(2), blocks(1))]);
            len(2:end) = len(2:end) - (oldnx - nx) .* [1:length(len(2:end))]';
            res0 = [];
            for i = 1:length(fieldnames(blocks))
                %                 res0((i-1)*nx+1:i*nx)=invR*(Amat(:,(i-1)*nx+1:i*nx)'*res);%(len(i):len(i+1)-1));
                if (nx == len(i + 1) - len(i))  % if the block size match LTL
                    res0(len(i):len(i + 1) - 1) = invR * (Amat(:, len(i):len(i + 1) - 1)' * res); % (len(i):len(i+1)-1));
                end
            end
            res = res0';
            clear res0;
        end
    end
else
    res = Amat' * res;
end

if (length(idx) < Alen)
    res0 = zeros(Alen, size(res, 2));
    res0(idx, :) = res;
    res = res0;
end
