function seq = simgauss(mu,sigma,NT,varargin)
% seq = simgauss(mu,sigma,NT,...)
%
% Simulate a Gaussian random process using eigenvalue
% decomposition (default), Cholesky decomposition or cholcov.
%
% Input arguments:
% - mu = mean (scalar or vector)
% - sigma = covariance matrix (must be symmetric)
% - NT = number of samples
%
% Optional arguments (key/value pairs):
% - 'method', can be either
%    - 'eig' use eigenvalue decomposition
%    - 'chol' use Cholesky decomposition
%    - 'cholcov' use Cholesky covariance decomposition
%    - 'cholinc' use incomplet Cholesky decomposition (default)
% - 'droptol', tolerance for incomplete Cholesky decomposition (default: 1E-15)
% - 'testpos', can be true or false; Issue a warning if resulting covariance matrix
%    is not positive semi-definite
% - 'testdist' can be true or false; apply kstest to each sequence to check if it
%    follows a Gaussian distribution

% Note:
% - sigma must be positive semi-definite. The number of
%   negative eigenvalues can be checked by [~,num] = cholcov(sigma)
% - numerical errors may lead to non positive semi-definiteness
%   - cholinc neglects these errors if smaller than droptol and
%     returns full decomposition
%   - cholcov neglects subspace with negative eigenvalue and
%     returns decomposition in complementary subspace
%   - chol throws error

args = varargin;
nargs = length(args);

% set defaults
method = 'cholinc';
droptol = 1E-15;
checkpos = false;
checkdist = false;
showcdf = false;

% parse varargs
for ii = 1:2:nargs
    switch args{ii}
      case 'method', method = args{ii+1};
      case 'droptol', droptol = args{ii+1};
      case 'testpos', checkpos = args{ii+1};
      case 'testdist', checkdist = args{ii+1};
      case 'showcdf', showcdf = args{ii+1};
      otherwise, warning(['unrecognized argument ' args{ii}])
    end
end

% noise dimension
ss = size(sigma,1);

% if empty covariance is given return empty sequence
if isempty(sigma)
    seq = [];
    return
end

% if mu is scalar, all components
% have the same mean
if length(mu) == 1
    mum = mu;
else
    % Note: the final output is [NT,ss]
    mum = diag(mu)*ones(ss,NT);
end

% check sigma>=0
% (sampling a process which does not observe this will lead to
% negative Kalman estimation variances)
[~,pmark] = cholcov(sigma);
if pmark > 0
    error(['Given covariance matrix is not ' ...
           'positive semi-definite. Smallest eigenvalue is %g'], ...
          eigs(sigma,1,'sa'))
end

% Sample uncorrelated processes and transform
switch method
  case 'eig',
    % using eigenvalue decomposition
    % (see Stengel, p.337)
    % svd would probably be better as it is
    % numerically more stable
    [V,D] = eig(sigma);
    if not(all(diag(D>=0)))
        error(['Given covariance matrix is not positive ' ...
               'semi-definite.'])
    end
    C = V*sqrt(D);
    seq = C*randn(ss,NT)+mum;
  case 'chol',
    % using Cholesky decomposition
    % (see Stengel, p.339)
    C = chol(sigma);
    seq = C'*randn(ss,NT)+mum;
  case 'cholcov',
    % Note: this neglects zero eigenvalues
    % see cholchov help for more info
    C = cholcov(sigma);
    seq = C'*randn(size(C,1),NT)+mum;
  case 'cholinc',
    % obsolete in newer matlab versions
    % use 'ichol' instead
    C = cholinc(sparse(sigma),droptol);
    alpha = max(0,max(sum(abs(sigma),2)./diag(sigma))-2);
    seq = C'*randn(size(C,1),NT)+mum;
  case 'ichol',
    % Note: this is the prefered method as
    % its most robust to numerical errors which
    % lead to non positive definiteness
    % Does NOT work for SEMI-DEFINITE matrices!
    alpha = max(0,max(sum(abs(sigma),2)./diag(sigma))-2);
    C = ichol(sparse(sigma),struct('droptol',droptol,'type','ict'));
    seq = C*randn(size(C,1),NT)+mum;
  otherwise
    error(['unrecognized method ' args{1}])
end

% output sequence as [NT,ss] array
seq = seq.';

% tests
if checkdist
    disp('Kolmogorov-Smirnov test:')
    for mm = 1:size(seq,1)
        % make sure the tabulated values span the full range of x
        xx = (min(seq(mm,:))-1/NT:1/NT:max(seq(mm,:))+1/NT);
        cdff = @(x) [x;normcdf(x,0,sqrt(sigma(mm,mm)))];
        CDF = cdff(xx);                       % tabulate cdf values
        [h,p] = kstest(seq(mm,:),CDF')
        if showcdf
            figure
            cdfplot(seq(mm,:));
            plot(CDF(1,:),CDF(2,:),'r')
        end
    end
    disp('covariance matrix:')
    cm = cov(seq)
    disp('nominal:')
    sigma
end

if checkpos
    [~,pmark] = cholcov(cov(seq));
    if pmark>0
        warning(['Covariance matrix of sample is not ' ...
                 'positive semi-definite'])
    end
end