% BUILD_TUNABLE_PARAM_MATRIX Builds a *transposable* tunable matrix.
%
% Arguments
% sz: Size of the matrix. If provided as a scalar, produces a symmetric matrix with given size.
% prefix: The variable name prefix
% type: Either of "ltrig", "utrig", "sym", "asym", "full".
%
% SPDX-FileCopyrightText: 2023 Deutsches Elektronen-Synchrotron DESY, MSK, LbSync <lbsync-expert@desy.de>
% SPDX-License-Identifier: GPL-3.0-or-later
%
% Author: Maximilian Schütte <maximilian.schuette@desy.de>

function X = build_tunable_param_matrix(sz, prefix, type)
    arguments
        sz (1,:) double {mustBeInteger,mustBePositive}
        prefix (1,1) string
        type (1,1) string {matlab.system.mustBeMember(type, ["ltrig", "utrig", "sym", "skew", "full"])} = "ltrig"
    end
    type = lower(type);

    assert(size(sz, 2) <= 2, "Tensors not supported!");
    if isscalar(sz)
        if type == "full"
            warning("Provided scalar size but requested full matrix. Extending size to column vector.");
            sz = [sz, 1];
        else
            sz = [sz, sz];
        end
    end

    if ismember(type, ["sym", "skew"])
        assert(sz(1) == sz(2), "Incompatible size for symmetric/skew-symmetric matrix!");
    end

    X = genmat(zeros(sz));
    switch type
        case "ltrig"
            for i = 1:sz(1)
                for j = 1:i
                    if j > sz(2), continue; end  % handle tall matrices
                    X(i,j) = realp(prefix + i + j, 0);
                end
            end
        case "utrig"
            for j = 1:sz(1)
                for i = 1:j
                    if i > sz(1), continue; end  % handle tall matrices
                    X(i,j) = realp(prefix + i + j, 0);
                end
            end
            % handle fat matrix
            if sz(2) > sz(1)
                X = [zeros(sz(1), sz(2)-sz(1)), X(:,1:sz(1))];
            end
        case "sym"
            for i = 1:sz(1)
                for j = 1:i
                    X(i,j) = realp(prefix + i + j, 0);
                    X(j,i) = X(i,j);
                end
            end
        case "skew"
            for i = 1:sz(1)
                for j = 1:i
                    if i == j, continue; end
                    X(i,j) = realp(prefix + i + j, 0);
                    X(j,i) = -X(i,j);
                end
            end
        case "full"
            for i = 1:sz(1)
                for j = 1:sz(2)
                    X(i,j) = realp(prefix + i +j, 0);
                end
            end
    end
end
