function [finalAmat, finalrhs, nblock] = rbcreateinv(Amat, ymeas, ymodel, params, output)
%
% [finalAmat, finalrhs, nblock]=rbcreateinv(Amat, ymeas, ymodel, params, output)
%
% Convert an inverse problem to the log-amplitude and unwrapped phase form
% i.e. Amat*x=[ymeas - ymodel] =>
% [Alogamp Aphase]*x=[log10(ymeas)-log10(ymodel);angle(ymeas)-angle(ymodel)]
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     Amat: a complex-valued LHS matrix
%     ymeas: complex valued measurement data
%     ymodel: complex valued model prediction
%     param: unknown parameters for multiple wavelengths
%     output: 'complex': complex-valued system,
%                     finalAmat=Amat; finalrhs=ymeas-ymodel
%             'real': real-valued system,
%                     finalAmat=[Re(Amat) -Im(Amat);Im(Amat) Re(Amat)]
%                     finalrhs =[Re(ymeas-ymodel) ; Im(ymeas-ymodel)]
%             'logphase': real-valued system,
%                     finalAmat=see manual
%                     finalrhs =[log10(ymeas)-log10(ymodel) ; angle(ymeas)-angle(ymodel)]
%
% output:
%     finalAmat: converted LHS matrix, see above
%     finalrhs: converted RHS matrix, see above
%     nblock: after conversion, x can be converted to [Re(x);Im(x)]; if
%          this happens, nblock is set to 2, otherwise 1
%
% license:
%     GPL version 3, see LICENSE_GPLv3.txt files for details
%
% -- this function is part of Redbird-m toolbox
%

% convert an inverse problem to the log-amplitude and unwrapped phase form
%

if (nargin < 5)
    output = 'complex';
end

if (nargin < 4)
    param = [];
end

wavelengths = {''};

if (isa(Amat, 'containers.Map'))
    wavelengths = Amat.keys;
else
    Amat = containers.Map({''}, {Amat});
    ymeas = containers.Map({''}, {ymeas});
    ymodel = containers.Map({''}, {ymodel});
end

newA = containers.Map();
newrhs = containers.Map();

nblock = 1;

if (strcmp(output, 'complex'))
    newA = Amat;
    for waveid = wavelengths
        wv = waveid{1};
        newrhs(wv) = ymeas(wv) - ymodel(wv);
    end
else
    for waveid = wavelengths
        wv = waveid{1};
        rhs = ymeas(wv) - ymodel(wv);
        if (strcmp(output, 'real'))
            newA(wv) = real(Amat(wv));
            newrhs(wv) = real(rhs);

            if (~isreal(rhs) && ~isreal(Amat(wv)))
                newA(wv) = [real(Amat(wv)) -imag(Amat(wv))
                            imag(Amat(wv)) real(Amat(wv))];
                newrhs(wv) = [real(rhs); imag(rhs)];
                nblock = 2;
            end
        elseif (strcmp(output, 'logphase'))
            temp = conj(ymodel(wv)) ./ abs(ymodel(wv) .* ymodel(wv));
            temp = repmat(temp(:), 1, size(Amat(wv), 2)) .* Amat(wv);
            if (isreal(ymodel(wv)))
                newA(wv) = real(temp);
                temp = log(abs(ymeas(wv))) - log(abs(ymodel(wv)));
                newrhs(wv) = temp(:);
            else
                newA(wv) = [real(temp); imag(temp)];
                temp = log(abs(ymeas(wv))) - log(abs(ymodel(wv)));
                ptemp = angle(ymeas(wv))   - angle(ymodel(wv));
                newrhs(wv) = [temp(:); ptemp(:)];
                nblock = 2;
            end
        end
    end
end

finalAmat = [];
finalrhs = [];

if (isa(params, 'struct'))
    paramlist = intersect(fieldnames(params), {'hbo', 'hbr', 'water', 'lipids', 'aa3'});
    if (isempty(paramlist))
        error('specified parameters are not supported');
    end
    extins = rbextinction(wavelengths, paramlist);
    for i = 1:length(wavelengths)
        wv = waveid{i};
        wvAmat = [];
        for j = 1:length(paramlist)
            wvAmat = [wvAmat, newA(wv) * extins(i, j)];
        end
        finalAmat = [finalAmat; wvAmat];
        finalrhs = [finalrhs; newrhs(wv)];
    end
else
    for i = 1:length(wavelengths)
        wv = waveid{i};
        finalAmat = [finalAmat; newA(wv)];
        finalrhs = [finalrhs; newrhs(wv)];
    end
end
