classdef svm_linear
    %{
        Linear Support Vector Machine (SVM)
        ===================================
        
        The upper-level program seeks to choose x = [C, zeta] to minimise validation error.
        The lower-level program seeks to choose y = [b, w, xi] to minimise training error.
        Use the table below and svm_linear.pdf to reference variable names.
        
        | Variable        | Shape         | Note                                        |
        | --------------- | ------------- | ------------------------------------------- |
        | train_features  | (n_train, d)  | E.g. [[1.6,0.4], [4.7,1.5], [1.3,0.2], ...] |
        | train_labels    | (n_train,)    | E.g. [+1, -1, +1, ...]                      |
        | val_features    | (n_val, d)    |                                             |
        | val_features    | (n_val,)      |                                             |
        |                 |               |                                             |
        | w               | (d,)          | Weights                                     |
        | b               | (1,)          | Bias                                        |
        | C               | (1,)          | Regularization hyperparameter               |
        | xi              | (n_train,)    | Slack variables for training hinge loss     |
        | zeta            | (n_val,)      | Slack variables for training hinge loss     |
    %}
    properties(Constant)
        name = 'svm_linear';
        category = 'machine_learning';
        subcategory = 'classification';
        train_val_split = 0.66
        datasets = {
            'classification_ionosphere.csv';
            'classification_iris.csv';
            'classification_moons54.csv';
            'classification_palmer.csv';
        };
        paths = fullfile('bolib3', 'data', 'classification', svm_linear.datasets);
    end


    methods(Static)
        

        % //==================================================\\
        % ||                        F                          ||
        % \\==================================================//
        % Upper-level objective function (linear)
        function val = F(x, ~, ~)
            zeta = x(2:end);
            val = sum(zeta);
        end
        

        % //==================================================\\
        % ||                        G                          ||
        % \\==================================================//
        % Upper-level inequality constraints (linear)
        %   C >= 0
        %   zeta_i >= 0
        %   zeta_i >= 1 - y_i (w^T x_i + b)
        function val = G(x, y, data)       
            % Split variables
            [C, zeta, b, w, ~] = svm_linear.variable_split(x, y, data);
        
            % Compute validation hinge loss
            % val_features has size (n_val x d)
            % w has size (d x 1)
            val_hinge_loss = 1 - data.val_labels .* (data.val_features * w + b);
        
            % Combine constraints
            val = [
                C;
                zeta;
                zeta - val_hinge_loss
            ];
        end
        

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

        % //==================================================\\
        % ||                        f                          ||
        % \\==================================================//
        % Lower-level objective function (quadratic)
        % ||w||^2 + C * (1/n) * sum(xi)
        function val = f(x, y, data)
            % Reciprocal = 1 / n_train
            reciprocal = data.reciprocal;
        
            % Split variables
            [C, ~, ~, w, xi] = svm_linear.variable_split(x, y, data);
        
            % Compute objective
            val = dot(w, w) + C * reciprocal * sum(xi);
        end


        % //==================================================\\
        % ||                        g                          ||
        % \\==================================================//
        % Lower-level inequality constraints (linear)
        %   xi >= 0
        %   xi >= 1 - y_i (w^T x_i + b)
        function val = g(x, y, data)       
            % Training data
        
            % Split variables
            [~, ~, b, w, xi] = svm_linear.variable_split(x, y, data);
        
            % Compute hinge loss for training data
            train_hinge_loss = 1 - data.train_labels .* (data.train_features * w + b);
        
            % Combine constraints
            val = [
                xi;
                xi - train_hinge_loss
            ];
        end

        
        % //==================================================\\
        % ||                        h                          ||
        % \\==================================================//
        % Lower-level equality constraints
        function val = h(~, ~, ~)
            val = [];
        end

        
        % //==================================================\\
        % ||                  Variable Split                   ||
        % \\==================================================//
        % Splits the vectors x, y into their meaningful components.
        %   - Upper-level variables: X = [C; zeta_1; ...; zeta_n_val]
        %   - Lower-level variables: Y = [b; w_1; ...; w_d; xi_1; ...; xi_n_train]
        function [C, zeta, b, w, xi] = variable_split(x, y, data)
            % Upper-level variables
            x = x(:);
            C = x(1);
            zeta = x(2:end);
        
            % Lower-level variables
            y = y(:);
            b = y(1);
            w = y(2:data.d+1);
            xi = y(data.d+2:end);
        end
        
        % //==================================================\\
        % ||                     Evaluate                      ||
        % \\==================================================//
        % Calculates the training and validation accuracy scores.
        function [train_accuracy, val_accuracy] = evaluate(x, y, data)       
            % Split decision variables
            [C, ~, b, w, ~] = variable_split(x, y, data);
        
            % Make predictions
            train_predict = sign(data.train_features * w + b);
            val_predict = sign(data.val_features * w + b);
        
            % Compute accuracy
            train_accuracy = sum(data.train_labels == train_predict) / data.n_train;
            val_accuracy = sum(data.val_labels == val_predict) / data.n_val;
        
            % Display if verbose
            if verbose
                fprintf('C (hyperparameter):   %.2e\n', C);
                fprintf('Training accuracy:    %.2f%%\n', train_accuracy * 100);
                fprintf('Validation accuracy:  %.2f%%\n', val_accuracy * 100);
            end
        end


        % //==================================================\\
        % ||                    Read Data                      ||
        % \\==================================================//
        % Reads data from a CSV file and splits it into training and validation.
        %   - The CSV file must have a header row.
        %   - The first column must contain +1/-1 labels.
        %   - The remaining columns are features.
        %   - Splits the data according to the 'train_val_split' property.
        function data = read_data(filepath)
            % Read CSV
            options = detectImportOptions(filepath);
            data_table = readtable(filepath, options);
        
            % Extract labels and features
            labels = data_table{:,1};
            features = data_table{:,2:end};
        
            % Dimensions
            n_total = size(features, 1);
            d = size(features, 2);
        
            % Assertions
            if ~all(ismember(labels, [-1, 1]))
                error('Data contains labels other than -1 and +1');
            end
            if size(labels, 1) ~= n_total
                error('Labels and features must have same number of rows');
            end
        
            % Split dataset
            n_train = floor(n_total * svm_linear.train_val_split);
            n_val = ceil(n_total * (1 - svm_linear.train_val_split));
        
            train_features = features(1:n_train, :);
            train_labels = labels(1:n_train);
        
            val_features = features(n_train+1:end, :);
            val_labels = labels(n_train+1:end);
        
            % Store results in struct
            data = struct();
            data.file = filepath;
            data.d = d;
            data.n_val = n_val;
            data.n_train = n_train;
            data.reciprocal = 1 / n_train;
            data.train_features = train_features;
            data.train_labels = train_labels;
            data.val_features = val_features;
            data.val_labels = val_labels;
        end
        

        % //==================================================\\
        % ||                    Dimensions                     ||
        % \\==================================================//
        % Key are the function/variable names
        % Values are their dimention
        function n = dimention(key, data)
            n = dictionary( ...
                'x', 1 + data.n_val, ...
                'y', 1 + data.d + data.n_train, ...
                'F', 1, ...
                'G', 1 + 2*data.n_val, ...
                'H', 0, ...
                'f', 1, ...
                'g', 2*data.n_train, ...
                'h', 0 ...
            );
            if isKey(n,key)
                n = n(key);
            end
        end
    end
    

    % Feasible point for ionosphere dataset
    properties(Constant)
        x0 = [
            40.27; 2.87776; 0.0; 0.0; 0.05403; 0.625691; 0.0; 0.286515; 0.0; 0.0; 0.617832; 0.0; 0.0; 0.0; 2.161009; 0.706492;
            1.08977; 1.713132; 0.0; 0.629926; 0.547302; 1.714907; 0.0; 0.0; 0.0; 1.008483; 0.0; 1.771606; 3.186846; 0.74349;
            0.06487; 0.499154; 0.15101; 1.83496; 1.879876; 1.644198; 0.0; 0.0; 0.961719; 0.299985; 0.628953; 2.061592
        ];
        y0 = [
            -1.16497; 0.464239; 0.0; 0.397513; -0.068354; 0.747696; 0.459217; 0.461731; 0.294408; -0.01772; 0.354079; 0.183116;
            0.089647; -0.11211; 0.36672; 0.192438; -0.381318; -0.015473; 0.302549; 0.031944; 0.164399; 0.18765; -0.548938;
            0.268507; -0.141875; 0.223904; 0.096494; -0.449633; 0.283959; -0.026194; 0.082053; 0.095415; -0.070464; -0.043594;
            -0.464038; 0.0; 0.0; 0.547609; 0.0; 0.0; 0.231288; 0.0; 0.031851; 0.0; 1.108514; 0.136396; 0.0; 0.0; 0.523643;
            0.353763; 0.0; 0.0; 0.0; 0.0; 0.306526; 0.0; 0.207194; 0.0; 0.0; 0.0; 0.0; 0.594968; 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.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 1.715708;
            0.457647; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.888483; 0.844932; 0.0; 0.0; 0.0; 0.056951; 0.346067; 0.0;
            0.074463; 2.428207; 0.0; 0.0; 0.0; 0.0; 0.0
        ];
    end
end