function [t_x, x, y, t, Ex, Ey, beta] = getFunData_eg2(n, lint, ngrid, model, ...
    regular, varx, vary, evals, C1, C2)
% one-dimensional central subspace examples
% model 1 and 2 is the same as model I and II from Yao, Lei and Wu (2015)

% Create sparse and dense functional data (with measurement error) with
% sccalar response
% Input n:          Number of subjects to create data for
% Input lint:       Upper bound on time, i.e., the time of each measurement
%                   is supported on [0, lint]
% Input ngrid:      Number of time measurements for dense case
% Input model:      Specify which model of y to use
% Input basis:      Specify which orthonormal basis to use for beta
%                   1, Fourier basis
%                   2, Legendre basis
%                   3, Linear combination of fouriers
% Input regular:    0, sparse (or irregular) functional data.      
%                   2, completely balanced (regular) data.

% Output t_x(1,n):  Vector of time points (for dense case)
% Output x(1,n):    Vector of functional measurements (with meas. error)
% Output y(n,1):    Matrix of scalar responses (with regression error)
% Output Ex(1,n):   Vector of functional measurements (without meas. error)
% Output Ey(n,1):   Matrix of scalar responses (without regression error)
% Output K(scalar): True dimension of projection

t = linspace(0, lint, ngrid);
t_x = cell(1,n); 
x = cell(1,n); 
Ey = zeros(n,1); 
Ex = cell(1,n);
bx = zeros(n,3);

% FPC scores for X
pc_x = randn(n, length(evals))*diag(sqrt(evals));

% Create regression coefficient functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simplified linear single index model
%b1 = C1*[zeros(1,4) sqrt(5.^1.5) zeros(1,95)]*xeig(t, lint, length(evals)); 
b1 = C1*[1 1 1 4*((4:length(evals))-2).^-3]*xeig(t, lint, length(evals));
b2 = C2*[zeros(1,5) sqrt(6.^1.5) zeros(1,94)]*xeig(t, lint, length(evals)); 
beta = [b1; b2];

% Create data
fx = zeros(3, n);
for i = 1:n
    % Sparse case:
    if regular == 0
        % Get number of time points for subject
        nwu = floor(random('unif',5,11,1,1));
        %nwu = floor(random('unif', 5, 11, 1, 1));
        t_x{i} = sort(random('unif',0,lint,1,nwu));
    % Dense Case:
    elseif regular == 2
        % Create dense time points
        t_x{i} = t;
    end
        
    % Create eigenfunctions
    phi = xeig(t, lint, size(pc_x, 2)); 
    
    % True x
    Ex{i} = pc_x(i,:)*phi;
    tmp = interp1(t, Ex{i}, t_x{i}, 'splines');
    
    % Observed x
    x{i} = tmp + sqrt(varx)*randn(1,length(t_x{i}));
    clear tmp
    
    % Create regression functions
    for k = 1:size(beta,1)
        bx(i,k) = trapz(t, beta(k,:).*Ex{i});
    end
    
    % Create model for Y
    if model == 1
        K = 2;
        fx(1,i) = (bx(i,1));
        fx(2,i) = exp(bx(i,2)/10);
    elseif model == 2
        K = 2;
        fx(1,i) = atan(bx(i,1)*pi);
        fx(2,i) = (bx(i,2));
    elseif model == 3
        K = 2;
        fx(1,i) = 1.5*bx(i,1)/(0.5+(1.5+bx(i,2)).^2);
        fx(2,i) = 0;
    end
end


fx = fx(1:K, : );
%var(fx,0,2)
%var(bx,0,1)
Ey = sum(fx ,1)';


y = Ey + sqrt(vary)*randn(n,1);
%var(Ey)/var(sqrt(vary)*randn(n,1))  
beta = beta(1:K,:);
end
