classdef svm_ward2025
    %{
        Samuel Ward, Alain Zemkoho, and Selin Ahipasaoglu.
        Mathematical programs with complementarity constraints and application
        to hyperparameter tuning for nonlinear support vector machines. 2025.
    %}
    properties(Constant)
        name = 'svm_ward2025';
        category = 'machine_learning';
        subcategory = 'classification';
        K_folds = 3;
        datasets = {
            'classification_ionosphere.csv';    
            'classification_circles90.csv';
            'classification_hastie2009.csv';
            'classification_iris.csv';
            'classification_moons54.csv';
            'classification_palmer.csv'
        };
        paths = fullfile('bolib3', 'data', 'classification', svm_ward2025.datasets);
    end


    methods(Static)
        % //==================================================\\
        % ||                        F                          ||
        % \\==================================================//
        % Upper-level objective function (linear)
        % Sum of the upper-level varaibles from index 3 to end
        function val = F(x, ~, ~)
            val = sum(x(3:end)) / (length(x) - 2);
        end


        % //==================================================\\
        % ||                        G                          ||
        % \\==================================================//
        % Upper-level inequality constraints (nonconvex)
        %   - C >= 0
        %   - gamma >= 0
        %   - zeta >= 0 for all k,i
        %   - zeta^k_i - 1 + sum_j (alpha^k_j * Q^k_ij) + y^k_i * b >= 0 for each k,i
        function ineq = G(x, y, data)
            % Number of validation/training examples
            n_val = data.n_val;
            n_train = data.n_train;
        
            % Upper-level hyperparameters
            C = x(1);
            gamma = x(2);
        
            % Start with non-negativity of x
            ineq = x;
        
            % Loop over folds
            for k = 1:svm_ward2025.K_folds
                % Index the subvectors for zeta and alpha relating to the k-th fold
                idx_zeta = 2 + ((k-1)*n_val) + (1:n_val);
                idx_alpha = ((k-1)*n_train) + (1:n_train);
        
                zeta_k = x(idx_zeta);
                alpha_k = y(idx_alpha);
        
                % Extract data relating to the k-th fold
                train_labels_k = data.train_labels{k};
                train_features_k = data.train_features{k};
                val_labels_k = data.val_labels{k};
                val_features_k = data.val_features{k};
        
                % Calculate the kernel and the bias (depend on the hyperparameters C and gamma)
                Q_k = exp(-gamma * svm_ward2025.pairwise_squared_distance_matrix(val_features_k, train_features_k));
                bias_k = svm_ward2025.reconstruct_bias(alpha_k, C, gamma, train_features_k, train_labels_k);
        
                % Add the constraint Z
                constraint_Z = zeta_k - 1 + (Q_k * (alpha_k .* train_labels_k)) - (val_labels_k * bias_k);
                ineq = [ineq; constraint_Z];
            end
        end


        % //==================================================\\
        % ||                        H                          ||
        % \\==================================================//
        % Upper-level equality constraints (none)
        function val = H(~, ~, ~)
            val = [];
        end
        

        % //==================================================\\
        % ||                        f                          ||
        % \\==================================================//
        % Lower-level objective function
        function obj = f(x, y, data)
            gamma = x(2);
            obj = 0;
            for k = 1:svm_ward2025.K_folds
                idx_alpha = ((k-1)*data.n_train) + (1:data.n_train);
                alpha_k = y(idx_alpha);
        
                a_times_y = alpha_k .* data.train_labels{k};
                Q = exp(-gamma * svm_ward2025.pairwise_squared_distance_matrix(data.train_features{k}, data.train_features{k}));
        
                obj = obj + 0.5 * (a_times_y' * Q * a_times_y) - sum(alpha_k);
            end
        end
        

        % //==================================================\\
        % ||                        g                          ||
        % \\==================================================//
        % Lower-level inequality constraints
        % 0 <= alpha <= C
        function val = g(x, y, ~)
            val = [
                    y(:); 
                    y(:) - x(1)
                ];
        end
        

        % //==================================================\\
        % ||                        h                          ||
        % \\==================================================//
        % Lower-level equality constraints
        % For each fold k=1,...,K:
        %   sum_{i=1} (aplha(k,i)) * (train_labels(k,i)) == 0
        function val = h(~, y, data)
            val = zeros(svm_ward2025.K_folds, 1);

            % Iterate over the K folds
            for k = 1:(svm_ward2025.K_folds)
                idx_alpha_k = ((k-1)*data.n_train) + (1:data.n_train);
                alpha_k = y(idx_alpha_k);
                val(k) = sum(alpha_k .* data.train_labels{k});
            end
        
        end
        

        % //==================================================\\
        % ||        Pairwise Squared Distance Matrix           ||
        % \\==================================================//
        %   a: [n_a x d]
        %   b: [n_b x d]
        %   D: [n_a x n_b]  where   D(i,j) = ||a_i - b_j||^2
        function D = pairwise_squared_distance_matrix(a, b)
            % Compute squared norms
            a_norm_sq = sum(a.^2, 2);
            b_norm_sq = sum(b.^2, 2);
        
            % Pairwise squared distances
            D = bsxfun(@plus, a_norm_sq, b_norm_sq') - 2 * (a * b');
        end
        

        % //==================================================\\
        % ||                 Reconstruct Bias                  ||
        % \\==================================================//
        % Compute the SVM bias term from dual variables.
        % See Appendix B of Ward 2025.
        function b = reconstruct_bias(alpha_k, C, gamma, train_features_k, train_labels_k)
            rbf_kernel = exp(-gamma * svm_ward2025.pairwise_squared_distance_matrix(train_features_k, train_features_k));
            omega_of_alpha_k = 1 - ((2 * alpha_k / C) - 1).^4 + 1e-5;
            sum_term = train_labels_k - (rbf_kernel * (alpha_k .* train_labels_k));
            b = dot(omega_of_alpha_k' , sum_term) / sum(omega_of_alpha_k);
        end


        % //==================================================\\
        % ||                     Evaluate                      ||
        % \\==================================================//
        % Calculate mean training and validation accuracies for the SVM.
        function [mean_train_accuracy, mean_val_accuracy] = evaluate(x, y, data)        

            train_accuracies = zeros(svm_ward2025.K_folds, 1);
            val_accuracies = zeros(svm_ward2025.K_folds, 1);
            % Iterate over K folds
            for k = 1:svm_ward2025.K_folds

                % Make predictions
                training_predict = svm_ward2025.predict(x, y, data, k, data.train_features{k});
                validation_predict = svm_ward2025.predict(x, y, data, k, data.val_features{k});

                % Callcuate accuracy
                train_accuracies(k) = sum(training_predict == data.train_labels{k}) / data.n_train;
                val_accuracies(k) = sum(validation_predict == data.val_labels{k}) / data.n_val;
            end
        
            % Mean accuracies
            mean_train_accuracy = mean(train_accuracies);
            mean_val_accuracy = mean(val_accuracies);
        
            % Output
            fprintf('C:                         %.2e\n', x(1));
            fprintf('gamma:                     %.2e\n', x(2));
            fprintf('Mean training accuracy:    %.2f%%\n', mean_train_accuracy * 100);
            fprintf('Mean validation accuracy:  %.2f%%\n', mean_val_accuracy * 100);
        end

        % //==================================================\\
        % ||                     Predict                       ||
        % \\==================================================//
        % Takes the k-th lower-level SVM represented by decision variables (x, y, data) 
        % Makes classification predictions for unseen_examples
        function predictions = predict(x, y, data, k, unseen_examples)
            % Hyperparameters C and gamma
            C = x(1);
            gamma = x(2);

            % SVM dual variables for fold k
            idx_alpha = ((k-1)*data.n_train) + (1:data.n_train);
            alpha_k = y(idx_alpha);
            bias_k = svm_ward2025.reconstruct_bias(alpha_k, C, gamma, data.train_features{k}, data.train_labels{k});

            % Make predictions
            kernel = exp(-gamma * svm_ward2025.pairwise_squared_distance_matrix(unseen_examples, data.train_features{k}));
            decision = kernel * (alpha_k .* data.train_labels{k}) + bias_k;
            predictions = sign(decision);
        end



        % //==================================================\\
        % ||                    Read Data                      ||
        % \\==================================================//
        % Reads a classification dataset from a csv file.
        %     - First row should contain header names (e.g. 'label', 'feature1', 'feature2', ...)
        %     - First column should contain +1/-1 labels (e.g. [+1, -1, -1, +1, ...])
        % 
        % Splits the data into K training and validation folds.
        % If len(data) is not divisible by K then the last (len(data)%K) examples are dropped.
        % For example consider the dataset X=[1,...,17] and parameter K is set to 3.
        % The procedure can be visualized with the following table:
        % 
        % |     |   i=1,2,3,4,5    |   i=6,7,8,9,10   | i=11,12,13,14,15 |
        % |-----|------------------|------------------|------------------|
        % | k=1 |    validation    |      train       |      train       |
        % | k=2 |      train       |    validation    |      train       |
        % | k=3 |      train       |      train       |    validation    |
        function data = read_data(filepath)
            % Read CSV as table
            data_table = readtable(filepath);
            K = svm_ward2025.K_folds;
        
            % Extract labels (first column) and features (remaining columns)
            labels = table2array(data_table(:, 1));
            features = table2array(data_table(:, 2:end));
        
            % Assertions
            if ~all(ismember(labels, [-1, 1]))
                error('Data contains labels other than -1 and +1');
            end
            if size(labels, 1) ~= size(features, 1)
                error('Labels and features must have the same number of rows');
            end
        
            % Number of examples per fold
            n_val = floor(length(labels) / K);
            n_train = n_val * (K - 1);
        
            % Preallocate cells
            val_features = cell(K, 1);
            val_labels = cell(K, 1);
            train_features = cell(K, 1);
            train_labels = cell(K, 1);
        
            % Split data into folds
            for k = 1:K
                idx_val = ( (k-1)*n_val + 1 ) : (k*n_val);
        
                val_features{k} = features(idx_val, :);
                val_labels{k} = labels(idx_val, :);
        
                idx_train = setdiff(1:(n_val*K), idx_val);
        
                train_features{k} = features(idx_train, :);
                train_labels{k} = labels(idx_train, :);
            end
        

        
            % Form struct
            data = struct();
            data.file = filepath;
            data.n_val = n_val;
            data.val_features = val_features;
            data.val_labels = val_labels;
            data.n_train = n_train;
            data.train_features = train_features;
            data.train_labels = train_labels;
        end
        
        % //==================================================\\
        % ||                    Dimensions                     ||
        % \\==================================================//
        % Key are the function/variable names
        % Values are their dimention
        function n = dimention(key, data)
            n = dictionary( ...
                'x', 2 + svm_ward2025.K_folds*data.n_val, ...
                'y', svm_ward2025.K_folds*data.n_train, ...
                'F', 1, ...
                'G', 2 + 2*svm_ward2025.K_folds*data.n_val, ...
                'H', 0, ...
                'f', 1, ...
                'g', 2*svm_ward2025.K_folds*data.n_train, ...
                'h', svm_ward2025.K_folds ...
            );
            if isKey(n,key)
                n = n(key);
            end
        end
    end
    
    % Feasible point for ionosphere dataset
    properties(Constant)
        x0 = [
            10.0; 0.0316; 1.478029; 0.515214; 0.638951; 0.16882; 0.0; 0.47903; 0.0; 0.492319; 0.11637; 1.126847; 0.0; 0.0; 0.0;
            0.425634; 0.136929; 0.0; 0.0; 0.039028; 0.714666; 0.373124; 0.859488; 0.096156; 0.0; 0.0; 0.328446; 1.082523; 0.0;
            0.351459; 0.0; 0.0; 0.04637; 0.060679; 0.0; 0.183981; 0.0; 0.0; 0.152899; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.093246;
            0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.17625; 0.0; 1.995379; 0.0; 0.020577; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0;
            1.082444; 0.84025; 0.867701; 0.0; 0.0; 1.489934; 0.796183; 0.0; 0.938974; 0.8357; 0.0; 0.0; 0.0; 0.0; 0.0; 0.088078;
            0.0; 0.0; 0.319574; 0.629909; 0.0; 0.0; 0.0; 0.0; 0.737598; 0.0; 0.0; 0.0; 2.473741; 0.371078; 0.902599; 0.863943;
            0.0; 0.5875; 0.062201; 2.054707; 0.0; 0.0; 0.005878; 1.260742; 0.0; 1.77453; 0.614526; 0.607405; 0.0; 0.583617;
            0.113608; 1.866811; 1.199894; 2.051952; 0.146744; 0.0; 0.545599; 0.30468; 0.662842; 0.0
        ];
        y0 = [
            0.0; 0.0; 0.0; 1.326666; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 10.0; 0.149757; 10.0; 0.023316; 2.004749; 8.191228;
            0.0; 0.0; 0.0; 0.0; 5.154778; 0.0; 10.0; 10.0; 2.762516; 1.618138; 0.0; 2.580872; 2.239659; 0.0; 10.0; 0.695156;
            0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.965842; 0.0; 0.0; 0.402557; 0.0; 0.0; 0.0; 0.0; 1.298699; 0.0; 0.0; 10.0; 10.0;
            10.0; 10.0; 1.53898; 0.0; 10.0; 0.0; 10.0; 6.152073; 0.0; 0.0; 10.0; 0.0; 10.0; 1.958521; 0.0; 0.0; 0.03816;
            3.111662; 4.658287; 9.517782; 10.0; 0.076704; 0.0; 5.30926; 10.0; 10.0; 0.0; 10.0; 0.549232; 2.255263; 1.537917;
            0.0; 10.0; 0.0; 1.326571; 4.016173; 10.0; 0.0; 0.0; 0.0; 10.0; 3.698682; 0.0; 0.0; 10.0; 2.208426; 9.652785;
            3.525443; 3.871607; 0.0; 0.0; 0.0; 2.626525; 0.0; 0.36712; 0.0; 2.136759; 1.033481; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0;
            0.0; 4.987289; 0.0; 0.0; 0.0; 0.0; 0.364872; 0.0; 0.0; 0.0; 0.0; 10.0; 0.0; 0.0; 0.0; 10.0; 10.0; 10.0; 4.560768;
            0.0; 10.0; 0.0; 10.0; 0.0; 0.0; 0.0; 5.073862; 0.0; 7.656005; 1.36613; 4.030989; 0.0; 1.627093; 0.0; 7.948106; 10.0;
            10.0; 0.0; 0.0; 0.0; 0.0; 10.0; 0.0; 0.389106; 1.238009; 3.97614; 0.846078; 1.39715; 0.0; 0.0; 3.218744; 0.0; 10.0;
            0.0; 1.482606; 0.0; 10.0; 0.0; 3.684959; 0.301541; 0.0; 1.211874; 10.0; 0.0; 1.920627; 0.0; 0.0; 0.0; 3.462367;
            1.97356; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.389077; 2.090656; 0.0; 0.0;
            0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 10.0; 5.604287; 2.856547; 0.0; 1.282584; 0.0; 0.0; 0.0; 1.8536; 0.0; 10.0;
            8.371425; 0.0; 1.709428; 0.0; 3.178985; 4.984181; 0.0; 5.936704; 2.279125; 0.0; 0.0; 0.0; 0.0; 0.0; 1.631075
        ];
    end
end