
function [nX,D,y,ye,betatrue,supp] = gendata(n,p,T,sigma,cor,ratio,seednum,isnormstan,largecor)
%========================================================================
% Input:
% n:           number of features 
% p:           number of observation 
% T:           number of nonzero elemets in  target solution 
% ratio:       range of value in betatrue = 10^ratio;
% sigma:       varince of white gausian noise 
% cor:         coorelations  of X 
% seednum:     the seed number 
% isnormstan:  1 normlized and stantalized; 0 not
% largecore:   coorelation parameter for large scale problem    
% Output:
% nX:          covariate matrix normlized and standlized if isnormstan = 1
% D:           [], if isnormstan = 0; diagnal matrix if isnormstan = 1
% y:           outcome vector 
% betature:    the real solution 
% supp:        the support of betatrue 
%======================================================================== 


rand('seed',seednum);   % fix seed 
randn('seed', seednum); % fix seed 
ratio = ratio*20;
% creat betatrue
betatrue = zeros(p,1);  
q = randperm(p);
supp = q(1:T);
valbeta = ratio/20*rand(T,1);
valbeta = valbeta- min(valbeta);
valbeta = valbeta/max(valbeta)*ratio/20;
if ratio ~= 0
betatrue(supp) = 10.^valbeta.*sign(randn(T,1));
else
betatrue(supp) = sign(randn(T,1));
end

% creat X 

t = cputime;
if ~exist('largecor','var')||isempty(largecor)
        disp(' Data generation may be slow due to Cholesky factoration ...')
        Sigma = zeros(p,p);
        % Sigma = full(Sigma);
        for k = 1:p
            for l = 1:p 
                Sigma(k,l)=cor^(abs(k-l));
            end
        end
        X = randn(n,p);
        % Sigma = full(Sigma);
        X= X*chol(Sigma);
else
        X = randn(n,p);
        for j = 2:p-1
            X(:,j) = X(:,j) + largecor*X(:,j+1)+ largecor*X(:,j-1);
        end
end
t = cputime - t;
display(sprintf('generate done using time (second) = %g',t))
	
% creat noisy observations  
eta = sigma*randn(n,1);
ye = X*betatrue; 
y =  ye + eta;

% normalize and centerlize 
if isnormstan ==1
    y = y - mean(y);
   [nX,D] = normalize(X);
else
    nX = X;
    D = [];
end
end

function [sX,D] = normalize(X)

[n,p] = size(X);
X = X-repmat(mean(X),n,1);
D = spdiags((1./sqrt(ones(1,n)*X.^2))',0,p,p);
sX = X*D;
end