function [T,p] = cholcov_(Sigma,flag)
%CHOLCOV  Cholesky-like decomposition for covariance matrix.
%   T = CHOLCOV(SIGMA) computes T such that SIGMA = T'*T.  SIGMA must be
%   square, symmetric, and positive semi-definite.  If SIGMA is positive
%   definite, then T is the square, upper triangular Cholesky factor.
%
%   If SIGMA is not positive definite, T is computed from an eigenvalue
%   decomposition of SIGMA.  T is not necessarily triangular or square in
%   this case.  Any eigenvectors whose corresponding eigenvalue is close to
%   zero (within a small tolerance) are omitted.  If any remaining
%   eigenvalues are negative, T is empty.
%
%   [T,P] = CHOLCOV(SIGMA) returns the number of negative eigenvalues of
%   SIGMA, and T is empty if P>0.  If P==0, SIGMA is positive semi-definite.
%
%   If SIGMA is not square and symmetric, P is NaN and T is empty.
%
%   [T,P] = CHOLCOV(SIGMA,0) returns P==0 if SIGMA is positive definite, and
%   T is the Cholesky factor.  If SIGMA is not positive definite, P is a
%   positive integer and T is empty.  [...] = CHOLCOV(SIGMA,1) is equivalent
%   to [...] = CHOLCOV(SIGMA).
%
%   Example:
%   Factor a rank-deficient covariance matrix C.
%       C = [2 1 1 2;1 2 1 2;1 1 2 2;2 2 2 3]
%       T = cholcov(C)
%       C2 = T'*T
%   Generate data with this covariance (aside from random variation).
%       C3 = cov(randn(10000,3)*T)
%
%   See also CHOL.

%   Copyright 1993-2009 The MathWorks, Inc. 

if nargin < 2, flag = 1; end

% Test for square, symmetric
[n,m] = size(Sigma);
wassparse = issparse(Sigma);
tol = 10*eps(max(abs(diag(Sigma))));
condition2 = all(abs(all(all(Sigma - Sigma')))<tol);
if (condition2==1) 
    % abs(all(all(Sigma - Sigma')))<tol
    [T,p] = chol(Sigma);

    if p > 0
        % Test for positive definiteness
        if flag
            % Can get factors of the form Sigma==T'*T using the eigenvalue
            % decomposition of a symmetric matrix, so long as the matrix
            % is positive semi-definite.
            [U,D] = eig(full((Sigma+Sigma')/2));
            % Use real() to elemenate imaginary part - C always returns
            % complex numbers when using eval(); zsebenik
            U = real(U);
            D = diag(real(D));
            if ~issorted(D(:,1))
                [V,D] = eig(full((Sigma+Sigma')/2));
                V = real(V);
                D = real(diag(D));
                [D,I] = sort(D);
                U = V(:, I);
            end

            % Pick eigenvector direction so max abs coordinate is positive
            [ignore,maxind] = max(abs(U),[],1);
            negloc = (U(maxind + (0:n:(m-1)*n)) < 0);
            U(:,negloc) = -U(:,negloc);

            %D = flipud(real(diag(D)));
            %D = real(diag(D));
            %D = sort(real(diag(D)));
            
            tol = repmat(eps(max(D)) * length(D), max(size(D)), 1);
            t_ = (abs(D) > tol);
            D = D(t_);
            p = sum(D<0); % number of negative eigenvalues

            if (p==0)
                T = diag(sqrt(D)) * U(:,t_)';
            else
                T = [];
            end
        else
            T = [];
        end
    end

else
    T = [];
    p = NaN;
end

if wassparse
    T = sparse(T);
end
