function [invalid1, invalid2, t_beta, beta, projx, f_projx, intercept, fcumout, modelout] = FCUMReg(x, t_x, y, options)

regular = options.regular;
error = options.error;
ngrid = options.ngrid;
kernel = options.kernel;
verbose = options.verbose;
xx_cov = options.xcov;
FVE_threshold = options.FVE_threshold;
xx_mu = options.xmu;
kn = options.kn;
maxk = options.maxk;
xPCA = options.PCA;
K = options.K;
npoly = options.npoly;
bw = options.bw;
ls_maxit = options.ls_maxit;
bf_maxit = options.bf_maxit;
ls_tol = options.ls_tol;
family = options.family;
V = options.V;
t_V = options.t_V;
newdata = options.newdata;
model = options.model;
selection_k = options.selection_k;

if isempty(newdata)
    newdata = unique(cell2mat(t_x));
end

n = length(y);
% Run cumulative slicing
if isempty(V) == 1
    p = fcumOptions('regular', regular, 'error', error, ...
        'verbose', verbose, 'ngrid', ngrid, 'kernel', kernel, ...
        'xcov', xx_cov, 'FVE_threshold', FVE_threshold, ...
        'xmu', xx_mu, 'kn', kn, 'maxk', maxk, 'PCA', xPCA, 'newdata', ...
        newdata, 'selection_k', selection_k);
    [invalid1, fcumout] = fcum(x, t_x, y, p);
    if invalid1 == 0
        t_beta = fcumout.out1;
        maxK = size(fcumout.Vphi, 2);
        % Obtain K e.d.r. directions
        beta = fcumout.Vphi(:,1:min(K, maxK));
        % If sparse design, replace X with X_PACE
        if regular == 0
            x = FPCAeval(fcumout.pca, [], t_beta)';
            for i = 1:n
                t_x{i} = t_beta;
            end
        end
    end
else
    invalid1 = 0;
    fcumout = []; 
    t_beta = newdata;
    [tmplam, tmpphi, ~, maxK] = getEigens(V, t_beta, t_V, maxk);
    [~, I] = sort(tmplam, 'descend'); beta = tmpphi(:,I); beta = beta(:,1:min(K, maxK));
end

if invalid1 == 0
    
    if strcmp(verbose, 'on') && maxK > K
        fprintf(['Reverting to maximum K = ' num2str(maxK) '\n']);
    end
    
    % Project X onto e.d.r. space
    projx = zeros(min(K, maxK), n);
    for i = 1:min(K, maxK)
        for j = 1:n
            tmp = interp1(t_beta, beta(:,i), t_x{j}, 'splines');
            projx(i,j) = trapz(t_x{j}, tmp.*x{j});
        end
    end
    
    if strcmp(model, 'gam')
        % Run GAM
        p = gamOptions('ngrid1', ngrid, 'npoly', npoly, ...
            'kernel', kernel, 'bw', bw, 'ls_maxit', ls_maxit, ...
            'bf_maxit', bf_maxit, 'verbose', verbose, 'ls_tol', ls_tol, ...
            'family', family);
        [invalid2, gamout] = gam(projx, y, p);
        
        if invalid2 == 0
            modelout.gam = gamout;
            f_projx = gamout.eta1;
            intercept = gamout.intercept;
        end
    elseif strcmp(model, 'lm')
        % Run LM
        X = [ones(n, 1), projx'];
        betacoef = pinv(X'*X)*X'*y;
        intercept = betacoef(1);
        f_projx = zeros(min(K, maxK), n);
        for i = 1:size(f_projx, 1)
            f_projx(i,:) = projx(i,:)*beta(i+1);
        end
        invalid2 = 0;
        modelout.betacoef = betacoef(2:end);
    end
    modelout.model = model;
end

end