function [invalid1, t_beta, beta, projx, fcumout, newprojx, newy] = FCUMClass(x, t_x, y, newx, newt_x, 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;
V = options.V;
t_V = options.t_V;
newdata = options.newdata;
selection_k = options.selection_k;
method = options.method;

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

n = length(y); 
class = unique(y);
numClass = length(class);
catProp = arrayfun(@(x) sum(x == y)/n, class);
% 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);
    if strcmpi(method, 'CUVE')
        [invalid1, fcumout] = fcuve(x, t_x, y, p);
    elseif strcmpi(method, 'SAVE')
        [invalid1, fcumout] = fsave(x, t_x, y, p);
    elseif strcmpi(method, 'CUME')
        [invalid1, fcumout] = fcum(x, t_x, y, p);
    elseif strcmpi(method, 'SIR')
        invalid1 = 0;
        fcumout = fsir(x, t_x, y, numClass, kn);
    end
    if invalid1 == 0
        t_beta = fcumout.out1;
        maxK = size(fcumout.Vphi, 2);
        % Obtain K e.d.r. directions
        if strcmp(method, 'CUME') || strcmp(method, 'SIR')
            beta = fcumout.Vphi(:,1:min(numClass-1,K));
        else            
            beta = fcumout.Vphi(:,1:min(K, maxK));
        end
        % 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(n, min(K, maxK));
    for j = 1:min(K, maxK)
        for i = 1:n
            tmp = interp1(t_beta, beta(:,j), t_x{i}, 'splines');
            projx(i,j) = trapz(t_x{i}, tmp.*x{i});
        end
    end
end

if ~isempty(newx)
    newn = length(newx);
    newprojx = zeros(newn, min(K, maxK));
    for j = 1:min(K, maxK)
        for i = 1:newn
            tmp = interp1(t_beta, beta(:,j), newt_x{i}, 'splines');
            newprojx(i,j) = trapz(newt_x{i}, tmp.*newx{i});
        end
    end
%     % Discrimination Methods
    S = cov(projx); grpS = cell(1, numClass);
    mu = zeros(numClass, min(K, maxK));
%     newy.lda = zeros(newn,1); newy.qda = zeros(newn,1);
%     lda = zeros(newn, numClass); qda = zeros(newn, numClass);
    for k = 1:numClass
        mu(k,:) = mean(projx(y==class(k),:),1);
        grpS{k} = cov(projx(y==class(k),:));
%         lda(:,k) = newprojx*S^-1*mu(k,:)' - 0.5*mu(k,:)*S^-1*mu(k,:)' + log(catProp(k));
%         qda(:,k) = -0.5*log(det(grpS{k})) - 0.5*arrayfun(@(x) (newprojx(x,:)-mu(k,:))*grpS{k}^-1*(newprojx(x,:)-mu(k,:))', 1:newn)' + log(catProp(k));
    end
%     for i = 1:newn
%         [~, I] = max(lda, [], 2);
%         newy.lda(i) = class(I(i));
%         [~, I] = max(qda, [], 2);
%         newy.qda(i) = class(I(i));
%     end
    % Mahalanobis distance
    madist = zeros(newn, numClass);
    for k = 1:numClass
        madist(:,k) = arrayfun(@(x) (newprojx(x,:)-mu(k,:))*grpS{k}^(-1)*(newprojx(x,:)-mu(k,:))', 1:newn)';
    end
    newy.ma = zeros(newn, 1);
    for i = 1:newn
        [~, I] = min(madist, [], 2);
        newy.ma(i) = class(I(i));
    end
else
    newy = [];
end

end