function [T,P,U,W,Q,B] = plsda(X, Y, ncomp)
% PLS-DA - Partial Least Squares Discriminant Analysis
% X: matriz de dados (amostras x variáveis)
% Y: vetor de classes (categórico ou numérico)
% ncomp: número de componentes

    % Converter Y para dummy variables se necessário
    if iscategorical(Y) || iscell(Y) || ischar(Y)
        [Y_dummy, classes] = dummy_encode(Y);
    else
        Y_dummy = Y;
        classes = unique(Y);
    end
    
    % Centralizar os dados
    X_centrado = X - mean(X, 1);
    Y_centrado = Y_dummy - mean(Y_dummy, 1);
    
    % Inicializar variáveis
    [n, p] = size(X);
    W = zeros(p, ncomp);  % Pesos de X
    T = zeros(n, ncomp);  % Scores de X
    P = zeros(p, ncomp);  % Loadings de X
    Q = zeros(size(Y_dummy, 2), ncomp);  % Loadings de Y
    U = zeros(n, ncomp);  % Scores de Y
    
    % Algoritmo NIPALS para PLS
    for comp = 1:ncomp
        % Inicializar u com a primeira coluna de Y
        u = Y_centrado(:, 1);
        
        % Iterar até convergência
        for iter = 1:100
            % Passo 1: Calcular pesos w
            w = (X_centrado' * u) / (u' * u);
            w = w / norm(w);
            
            % Passo 2: Calcular scores t
            t = X_centrado * w;
            
            % Passo 3: Calcar pesos c para Y
            c = (Y_centrado' * t) / (t' * t);
            c = c / norm(c);
            
            % Passo 4: Calcular scores u
            u_new = Y_centrado * c;
            
            % Verificar convergência
            if norm(u_new - u) < 1e-10
                break;
            end
            u = u_new;
        end
        
        % Passo 5: Calcular loadings p
        p = (X_centrado' * t) / (t' * t);
        
        % Passo 6: Calcular loadings q
        q = (Y_centrado' * u) / (u' * u);
        
        % Armazenar resultados
        W(:, comp) = w;
        T(:, comp) = t;
        P(:, comp) = p;
        Q(:, comp) = q;
        U(:, comp) = u;
        
        % Deflacionar matrizes
        X_centrado = X_centrado - t * p';
        Y_centrado = Y_centrado - t * q';
    end
    
    % Calcular coeficientes
    B = W / (P' * W) * Q';
    
    % Armazenar resultados
    % resultado.W = W;      % Pesos de X
    % resultado.T = T;      % Scores de X
    % resultado.P = P;      % Loadings de X
    % resultado.Q = Q;      % Loadings de Y
    % resultado.U = U;      % Scores de Y
    % resultado.B = B;      % Coeficientes
    % resultado.classes = classes;
    % resultado.ncomp = ncomp;
    
    % Calcular variância explicada
    % resultado.X_var = calc_var_explicada(T);
    % resultado.Y_var = calc_var_explicada(U);
end

function [Y_dummy, classes] = dummy_encode(Y)
% Converter variável categórica para dummy variables
    if iscategorical(Y)
        classes = categories(Y);
    else
        classes = unique(Y);
    end
    
    n_classes = length(classes);
    n_samples = length(Y);
    Y_dummy = zeros(n_samples, n_classes);
    
    for i = 1:n_classes
        if iscategorical(Y) || iscell(Y) || ischar(Y)
            Y_dummy(:, i) = strcmp(Y, classes{i});
        else
            Y_dummy(:, i) = (Y == classes(i));
        end
    end
end

function var_explicada = calc_var_explicada(scores)
% Calcular variância explicada
    var_total = sum(var(scores));
    var_explicada = zeros(1, size(scores, 2));
    
    for i = 1:size(scores, 2)
        var_explicada(i) = sum(var(scores(:, 1:i))) / var_total * 100;
    end
end