% MGGD_GENERATION
%
% This Matlab script contains three functions that can be used to generate 
% multivariate generalized Gaussian (MGGD) distributed sources. MGGD
% depends on two parameters, the scatter matrix (which is symmetric
% positive definite) and the shape parameter. If shape parameter is less
% than 1 the distribution of the marginals is more peaky, with heavier
% tails, and if it is greater than 1, it is less peaky with lighter tails.
% If shape parameter is equal to one then we generate Gaussian sources.
%
% EXAMPLE:
% 
% X = MGGD_generation(1000, 3, 0.5, 1,0);
%
% Coded by Lionel Bombrun (lionel.bombrun at u-bordeaux.fr);
%          Zois Boukouvalas (zb1 at umbc.edu);
%
% Machine Learning for Signal Processing Lab (MLSP-Lab)
% http://mlsp.umbc.edu
%
% Reference: E. Gomez, M. Gomez-Villegas, and J. Marin, "A multivariate generalization of the power exponential family of distributions," Communications in Statistics-Theory and Methods, vol. 27, no. 3, pp. 589{600}, 1998.
 

function [X,COV] = MGGD_generation(N, p, rho, beta,R_struct)
% GENERATE_MGGD
%
% Main function that generates MGGD distributed sources.
%
% INPUT: 
%   - N : Number of realizations(# of samples)
%   - p : Dimension (# of variables)
%   - rho : Correlation parameter 
%   - beta : Shape parameter (1 refers to Gaussian)
%   - R_struct : Which structure of correlation matrix R to use, 
%                0 for uniform correlation structure, 
%                1 for AR correlation structure, 
%                2 for self-defined correlation structure. When R_struct=2, rho
%                  should be a structure with two variables: rho.val and rho.idx. rho.val has
%                  two entries of correlation. rho.val(1) provides the basic correlation 
%                  and rho.idx gives specific locations of rho.val(2)
%                3 for QQ' case
%                4 for another type of self-designed correlation structure
%                  (user can change the code to design your correlation structure)
% 
% OUTPUT:
%   - X : pxN matrix containing MGGD sources 

%---------Modified by Q. Long 08/2018--------%%%%%%%%
%------------ qunfang1@umbc.edu -------------%%%%%%%%
betah = 2*beta;
etanln = 1/beta*log(2) + gammaln((p+2)/betah) -gammaln(p/betah) - log(p);
etan = exp(etanln);

X = transpose(RandSphere(N,p));

COV = abs(cov_sigma(p,rho,R_struct)); %This is the covariance matrix we want
Sigma_ = COV./etan; % Use this as the Sigma in Gomez's paper
flag = 0;
while flag == 0 % the COV might not be positive-definite so repeat generating until get a PD COV
    try 
        % X = (Sigma_)^(0.5)*X; % use the scatter matrix Sigma in Gomez's paper, rather the covariance matrix COV 
        X = chol(Sigma_,'lower')*X; % Use Cholesky decomposition instead of the square root of Sigma
        flag = 1;
    catch
        COV = abs(cov_sigma(p,rho,R_struct)); %This is the desired covariance matrix
        Sigma_ = COV./etan; % get the scatter matrix Sigma 
    end
end
tau = (gamrnd(p/(2*beta),2,1,N)).^(1/(2*beta));
tau = repmat(tau, p, 1);

X = tau.*X;

end


function Sigma = cov_sigma(p, rho, R_struct)
% COV_SIGMA
%
% Function that generates scatter matrix
%
% INPUT: 
%   - p : Dimension 
%   - rho : Correlation parameter 
% 
% OUTPUT:
%   - Sigma : pxp covariance matrix hence should be scaled by eta to get
%   the scatter matrix that will be used to generate MGGD samples in the
%   main function

Sigma = eye(p);
if R_struct == 3 % case for QQ' 
    Q = normr(rand(p)-0.2);
    Sigma = Q*Q';
elseif R_struct == 4 %self-designed correlation structure
    Sigma = eye(p);
    for i = 1:p
        for j = i+1:p
            if ismember(i,rho.idx) && ismember(j,rho.idx)
                Sigma(i,j) = rho.val(2);
            else
                Sigma(i,j) = 1e-3+0.2*rand(1);
            end
            Sigma(j,i) = Sigma(i,j);
        end
    end
else
    for i=1:p
        for j=1:p
            if R_struct == 0 % uniform correlation structure
                if i == j
                    Sigma(i,j) = 1;
                else
                    Sigma(i,j) = rho;
                end
            elseif R_struct == 1 % AR correlation structure
                Sigma(i,j) = rho^(abs(i-j));
            elseif R_struct == 2 % self-designed corr structure
                if i == j
                    Sigma(i,j) = 1;
                else
                    Sigma(i,j) = rho.val(1);
                end
            end
        end
    end
    if R_struct == 2
        for i = rho.idx
            for j = rho.idx
                if i~=j
                    Sigma(i,j) = rho.val(2);
                end
            end
        end
    end
    
end

end 


function X=RandSphere(N,dim)
% RANDSPHERE
%
% RandSphere generates uniform random points on the surface of a unit radius
% N-dim sphere centered in the origin . This script uses differents algorithms
% according to the dimensions of points:
%
%    -2D:  random generation of theta [0 2*pi]
%    -3D:  the "trig method".
%    -nD:  Gaussian distribution
%
%
% INPUT:
%
%    N: integer number representing the number of points to be generated
%    dim: dimension of points, if omitted 3D is assumed as default
%
% OUTPUT:
%
%   X: Nxdim double matrix representing the coordinates of random points
%   generated
%
% EXAMPLE:
% 
%   N=1000;
%   X=RandSphere(N);
%   hold on
%   title('RandSphere')
%   plot3(X(:,1),X(:,2),X(:,3),'.k');
%   axis equal
%
%Authors: Luigi Giaccari,Ed Hoyle

switch dim
    case 3 %3D
        
        %trig method
        X=zeros(N,dim);%preallocate
        X(:,3)=rand(N,1)*2-1;%z
        t=rand(N,1)*2*pi;
        r=sqrt(1-X(:,3).^2);
        X(:,1)=r.*cos(t);%x
        X(:,2)=r.*sin(t);%y
    case 2 %2D
        
        %just a random generation of theta
        X(:,2)=rand(N,1)*2*pi;%theta: use y as temp value
        X(:,1)=cos(X(:,2));%x
        X(:,2)=sin(X(:,2));%y
    
    otherwise %nD
        
        %use gaussian distribution
        X=randn(N,dim);
        X=bsxfun(@rdivide,X,sqrt(sum(X.^2,2)));
        
end

end