function y = ctrlScalarExpansion(u, n, isStrictPosDef, outputSqrt)
% ctrlScalarExpansion Construct covariance matrices from scalar, vector or matrix inputs
%
%    y  = ctrlScalarExpansion(u, n, isStrictPosDef, outputSqrt)
%
%   This is used for scalar expanding user inputs. There is minimal error
%   checking and corrections if the input does not correspond to a positive
%   semi-definite matrix.
%
%   An n-by-n matrix y is created from u:
%   * If u is a scalar, y has u on its diagonals. If u is negative, it is
%   replaced with 0 or eps(class(u)) based on isStrictPosDef
%   * If u is a vector, y has the elements of u on its diagonals. Negative
%   elements are treated akin to the scalar case explained above
%   * If u is a matrix, symmetricity is enforced via y = (u+u.')/2. There
%   are no positive semi-definiteness checks.
%
%   If outputSqrt is true(), then the output y is the Cholesky factor the
%   covariance matrix (L in P=L*L.'). If u is a matrix and it is not
%   positive semi-definite, then u and y*y' can differ significantly.
%
%   Inputs:
%     u              - (double or single) Scalar, vector or matrix input
%     n              - size(y,1)
%     isStrictPosDef - (bool) If true, and u is a scalar or vector,
%                      non-negative values are set to eps. Otherwise,
%                      they are set to 0.
%     outputSqrtY    - (bool) If true, the block outputs the lower
%                      triangular Cholesky factor of y
%
%   Output:
%     y              - Covariance matrix (or its Cholesky factor)

%   Copyright 2013-2018 The MathWorks, Inc.

%#codegen
coder.internal.prefer_const(n, isStrictPosDef, outputSqrt);

% Helper variables
ZERO = cast(0, 'like', u);

% If enforcing strict positive definiteness, map non-negative values to
% eps. Otherwise to 0.
if isStrictPosDef
    uLowerBound = eps(class(u));
else
    uLowerBound = ZERO;
end

y = zeros(n, n, 'like', u);
if isscalar(u)
    if u(1) <= ZERO
        u(1) = uLowerBound;
    end
    if outputSqrt
        u(1) = sqrt(u(1));
    end
    for kk = 1:n
        y(kk,kk) = u(1);
    end
elseif isvector(u)
    for kk = 1:n
        if u(kk) <= ZERO
            u(kk) = uLowerBound;
        end
        if outputSqrt
            u(kk) = sqrt(u(kk));
        end
        y(kk,kk) = u(kk);
    end
elseif ismatrix(u)
    if outputSqrt
        y = controllib.internal.util.cholSD(u);
        if isStrictPosDef
            uLowerBoundSqrt = sqrt(uLowerBound);
            for kk = 1:n
                if y(kk,kk) < ZERO
                    y(kk,kk) = uLowerBoundSqrt;
                end
            end
        end
    else
        y = (u + u.')/2;
    end
end
end
