classdef svr_bennett2006
    %{
        Kristin P Bennett, Jing Hu, Xiaoyun Ji, Gautam Kunapuli, and Jong-Shi Pang. 
        Model selection via bilevel optimization. 
        In the 2006 IEEE International Joint Conference on Neural Network Proceedings, pages 1922–1929.
    %}
    properties(Constant)
        name = 'svr_bennett2006';
        category = 'machine_learning';
        subcategory = 'regression';
        T_folds = 3;
        datasets = {
            'regression_auto_mpg.csv';    
            'regression_combined_cycle_power_plant.csv';
            'regression_energy_efficiency.csv';
            'regression_toy.csv';
        };
        paths = fullfile('bolib3', 'data', 'regression', svr_bennett2006.datasets);
    end


    methods(Static)
        
        % //==================================================\\
        % ||                        F                          ||
        % \\==================================================//
        % Upper-level objective function (convex)
        % Computes sum over folds of sum_i | Xw - y |
        function obj = F(~, y, data)
            obj = 0;
            d = data.d;
            for t = 1:svr_bennett2006.T_folds
                w_t = y((t-1)*d + 1 : t*d);
                residual = data.val_features{t} * w_t - data.val_labels{t};
                obj = obj + (sum(abs(residual))/length(data.val_labels{t}));
            end
        end

        
        % //==================================================\\
        % ||                        G                          ||
        % \\==================================================//
        % Upper-level inequality constraints (linear)
        % C             >= 0
        % epsilon       >= 0
        % lambda        >= 0
        % w_ub - w_lb   >= 0 
        function val = G(x, ~, data)
            d = data.d;
        
            % C, epsilon, lambda >= 0
            hp_nonnegativity = x(1:3);
        
            % w_lb <= w_ub  => w_ub - w_lb >= 0
            w_lb = x(4 : 3 + d);
            w_ub = x(4 + d : 3 + 2*d);
            w_bounds = w_ub - w_lb;
        
            % Stack them into one vector
            val = [hp_nonnegativity; w_bounds];
        end
        

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

        % //==================================================\\
        % ||                        f                          ||
        % \\==================================================//
        % Lower-level objective function (quadratic)
        % obj = sum_t [ C * sum(e_t) + 0.5 * w_t'w_t + lambda * w_t'w_t ]
        function obj = f(x, y, data)        
            % Split variables
            [hp_C, ~, hp_lambda, ~, ~, w, e] = svr_bennett2006.variable_split(x, y, data);
            obj = 0;
            for t = 1:svr_bennett2006.T_folds
                obj = obj + hp_C * sum(e{t}) + 0.5 * (w{t}' * w{t})  + hp_lambda * (w{t}' * w{t});
            end
        
        end
        

        % //==================================================\\
        % ||                        g                          ||
        % \\==================================================//
        % Lower-level inequality constraints (linear)
        % For each fold t=1,...,T:
        %     e_j >= + w*x_j - y_j - epsilon
        %     e_j >= - w*x_j + y_j - epsilon
        %     e_j >= 0
        %     w >= w_lb
        %     w <= w_ub
        function constraints = g(x, y, data)
            % Split variables
            [~, hp_epsilon, ~, w_lb, w_ub, w, e] = svr_bennett2006.variable_split(x, y, data);
        
            % Store all constraints in a cell array
            cells = cell(1, 5*svr_bennett2006.T_folds);
        
            for t = 1:svr_bennett2006.T_folds
                cells{5*(t-1) + 1} = e{t} + hp_epsilon - (data.train_features{t} * w{t}) + data.train_labels{t};
                cells{5*(t-1) + 2} = e{t} + hp_epsilon + (data.train_features{t} * w{t}) - data.train_labels{t};
                cells{5*(t-1) + 3} = e{t};
                cells{5*(t-1) + 4} = w{t} - w_lb;
                cells{5*(t-1) + 5} = w_ub - w{t};
            end
        
            % Concatenate all constraints into a single vector
            constraints = vertcat(cells{:});
        end

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

        % //==================================================\\
        % ||                  Variable Split                   ||
        % \\==================================================//
        % Upper-level variables:
        %     x = [C, epsilon, lambda, w_lb, w_ub]
        % Lower-level variables:
        %     y = [w[0], ..., w[T], e[0], ..., e[T]]
        %
        function [hp_C, hp_epsilon, hp_lambda, w_lb, w_ub, w, e] = variable_split(x, y, data)        
            % Dimensions
            d = data.d;
            n_train = data.n_train;
        
            % Upper-level decision variables
            hp_C = x(1);
            hp_epsilon = x(2);
            hp_lambda = x(3);
        
            w_lb = x(4 : 3 + d);
            w_ub = x(4 + d : 3 + 2*d);
        
            % Lower-level decision variables
            w = cell(1, svr_bennett2006.T_folds);
            e = cell(1, svr_bennett2006.T_folds);
            for t = 1:svr_bennett2006.T_folds
                idx_w_start = (t-1)*d + 1;
                idx_w_end   = t*d;
                w{t} = y(idx_w_start : idx_w_end);
        
                idx_e_start = svr_bennett2006.T_folds*d + (t-1)*n_train + 1;
                idx_e_end   = svr_bennett2006.T_folds*d + t*n_train;
                e{t} = y(idx_e_start : idx_e_end);
            end
        
        end
        
        % //==================================================\\
        % ||                     Evaluate                      ||
        % \\==================================================//
        %  Calculates the R2 and RMSE scores of the regression models represented by each fold
        function evaluate(~, y, data)
            d = data.d;
        
            % Print table headers
            fprintf('| %10s | %4s | %10s | %10s | %10s |\n', '', 'fold', 'data', 'R2 score', 'RMSE');
            fprintf('| %s | %s | %s | %s | %s |\n', repmat('-',1,10), repmat('-',1,4), repmat('-',1,10), repmat('-',1,10), repmat('-',1,10));
        
            % Evaluate for training and validation folds
            tasks = {'Training', 'Validation'};
            features_list = {data.train_features, data.val_features};
            labels_list = {data.train_labels, data.val_labels};
        
            for idx = 1:2        
                for t = 1:svr_bennett2006.T_folds
                    w_t = y((t-1)*d + 1 : t*d);
                    features = features_list{idx}{t};
                    labels = labels_list{idx}{t};
        
                    predictions = features * w_t;
        
                    % RMSE
                    rmse = sqrt(mean((labels - predictions).^2));
        
                    % R2 score
                    ss_res = sum((labels - predictions).^2);
                    ss_tot = sum((labels - mean(labels)).^2);
                    r2 = 1 - ss_res / ss_tot;
        
                    % Display
                    fprintf('| %-10s | %4d | %10s | %10.4f | %10.4f |\n', ...
                        tasks{idx}, t, sprintf('%dx%d', size(features,1), size(features,2)), r2, rmse);
                end
            end
        
        end


        % //==================================================\\
        % ||                    Read Data                      ||
        % \\==================================================//
        % The CSV file should have:
        %   - First row: headers (e.g. 'label', 'feature1', 'feature2', ...)
        %   - First column: labels (independent variable)
        %   - Remaining columns: features (dependent variables)
        % 
        % Splits the data into T training and validation folds.
        % If len(data) is not divisible by T then the last (len(data)%T) examples are dropped.
        % For example consider the dataset [1,...,17] and parameter T 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 |
        % |-----|------------------|------------------|------------------|
        % | t=1 |    validation    |      train       |      train       |
        % | t=2 |      train       |    validation    |      train       |
        % | t=3 |      train       |      train       |    validation    |
        function data = read_data(filepath)
            % Read the csv
            options = detectImportOptions(filepath);
            data_table = readtable(filepath, options);
        
            % Add a column of ones for the constant term
            data_table.constant = ones(height(data_table),1);
        
            % Total number of folds
            T = svr_bennett2006.T_folds;
        
            % Extract labels and features
            labels = table2array(data_table(:,1));
            features = table2array(data_table(:,2:end));
        
            % dimention (total number must be divsible by T)
            n_val = floor(length(labels) / T);
            n_total = n_val * T;
            n_train = n_total - n_val;

            % Preallocate cells
            val_features = cell(1, T);
            val_labels = cell(1, T);
            train_features = cell(1, T);
            train_labels = cell(1, T);

            % Form validation and training sets for each fold t=1,...,T
            for t = 1:T
                val_indices_t = ((t-1)*n_val + 1):t*n_val;
                train_indices_t = setdiff(1:n_total, val_indices_t);
                val_features{t} = features(val_indices_t, :);
                val_labels{t} = labels(val_indices_t);
                train_features{t} = features(train_indices_t, :);
                train_labels{t} = labels(train_indices_t);
            end
        
            % Form output struct
            data = struct();
            data.file = filepath;
            data.d = size(features, 2);
            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

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

    properties(Constant)
        x0 = [
            942.69;  % C
            1.10;  % Epsilon
            394.61;  % Lambda
            0.00; -0.81; -0.04; -0.01; -0.38; 0.57; 1.29; -3.18;  % w lower bound
            0.02; -0.33; -0.01; -0.00; 0.09; 0.64; 1.37; -3.18  % w upper bound
        ];
        y0 = [
            0.000349; -0.594078; 0.015399; 0.004389; -0.18728; 0.591777; 1.29; -3.18;  % w fold t=0
            0.02; -0.752486; 0.01; 0.006345; 0.059886; 0.577763; 1.29; -3.18;  % w fold t=1
            0.002926; -0.33; 0.04; 0.004406; -0.362387; 0.63628; 1.37; -3.18;  % w fold t=2
            0.0; 4.296325; 1.252274; 3.325708; 0.0; 0.0; 0.646209; 2.519647; 0.0; 0.024924; 0.252751; 0.0; 1.281405; 0.0;
            0.581238; 0.0; 0.038059; 0.855709; 2.960408; 0.377168; 4.414534; 4.779002; 0.0; 0.0; 2.371249; 0.0; 0.0; 1.435884;
            0.984823; 0.0; 1.784312; 0.376163; 3.132948; 1.021712; 1.745304; 1.584691; 0.003383; 1.570463; 0.0; 0.953865; 0.0;
            1.152161; 4.249906; 2.193958; 0.53174; 0.0; 0.126176; 0.0; 2.777644; 1.579941; 3.512486; 0.0; 1.456302; 5.110235;
            1.813491; 0.0; 2.146999; 2.270538; 0.0; 0.0; 3.080307; 7.047226; 4.250518; 1.386629; 0.0; 1.9534; 0.0; 0.0;
            1.815792; 0.330653; 0.0; 1.761197; 1.340089; 3.964859; 0.022587; 4.675208; 9.190773; 0.1355; 0.0; 0.0; 0.45286;
            2.310917; 1.220063; 0.0; 0.95974; 6.973211; 2.932587; 2.493241; 4.012014; 1.297346; 0.694085; 0.0; 0.0; 0.0;
            0.131152; 3.808762; 0.228293; 0.0; 4.630662; 0.0; 9.796482; 1.862521; 0.281683; 2.812953; 0.368985; 0.0; 3.262028;
            3.473138; 13.427186; 2.609227; 0.253462; 1.157014; 0.0; 4.722483; 0.271294; 0.32009; 0.0; 0.0; 4.427369; 0.484942;
            0.0; 0.0; 0.583326; 2.277108; 1.244839; 5.557292; 0.0; 0.062551; 4.183039; 0.0; 1.174414; 1.050416; 0.0; 0.126927;
            1.497728; 1.782488; 0.0; 0.0; 0.420509; 0.58596; 9.8e-05; 0.0; 2.014015; 0.410163; 1.873457; 0.0; 1.568642;
            0.441107; 2.063639; 0.0; 3.664282; 1.529542; 0.0; 3.346912; 13.327642; 3.667618; 0.0; 0.98113; 3.969028; 0.0; 0.0;
            0.294799; 3.886231; 1.637845; 0.045364; 4.645961; 1.265545; 2.849633; 0.0; 0.0; 1.436774; 4.233925; 1.059105;
            2.92423; 1.779308; 0.468137; 0.379592; 0.003632; 2.133575; 2.626549; 0.0; 0.214933; 2.053737; 0.0; 0.007523;
            0.315526; 0.0; 2.118025; 0.53898; 0.0; 0.0; 2.00817; 0.0; 0.0; 2.062082; 1.527982; 2.951778; 0.0; 0.0; 0.271359;
            0.964612; 2.607597; 2.080161; 0.0; 5.699111; 0.0; 1.825959; 5.0573; 3.484852; 0.410585; 0.0; 0.215467; 0.0; 0.0;
            0.122682; 0.345693; 0.576469; 0.406792; 1.008657; 0.914937; 1.049381; 3.838211; 0.974269; 1.220562; 1.821084;
            0.676337; 0.0; 1.200291; 1.029445; 0.0; 0.091608; 0.178226; 4.68604; 2.298553; 0.562112; 0.0; 0.0; 0.178522; 0.0;
            1.927432; 0.11813; 0.248471; 1.175662; 0.518371; 1.050619; 2.040495; 0.0; 7.192875; 0.0; 0.420858; 12.263202;
            0.380994; 2.078746; 0.39583; 5.609143; 0.0; 0.456673; 12.755954; 0.019168; 1.635657; 0.609954; 1.962734; 0.0;
            5.555875; 0.060674; 2.56366; 3.279727; 3.354162; 0.0; 6.494664; 0.152595; 1.265725; 0.0; 0.0; 0.0; 0.0; 3.231914;
            0.0; 0.0; 0.419252; 6.684132; 1.342993; 3.38836; 0.0; 9.311952; 0.425559; 1.437698; 0.0; 1.936583; 0.240557;
            2.150241; 0.0; 1.695011; 1.427797; 0.0; 0.0; 0.54037; 0.0; 0.0; 0.572981; 0.0; 9.360465; 0.0; 0.85105; 0.497131;
            0.0; 0.0; 0.501531; 0.0; 3.514277; 0.202377; 0.385678; 0.0; 0.960654; 2.365251; 4.652084; 3.822686; 0.0; 0.0;
            3.043464; 0.898241; 2.647794; 0.0; 1.710612; 0.0; 0.159314; 0.56644; 0.100085; 1.641978; 6.986741; 1.473045;
            0.139476; 0.951288; 0.0; 0.585494; 0.0; 0.0; 0.565002; 0.0; 0.566856; 0.0; 4.124855; 0.261632; 1.963646; 0.993579;
            1.136538; 2.192689; 1.646623; 0.0; 4.370751; 0.0; 0.778449; 3.867727; 0.625887; 9.650021; 1.577135; 0.0; 0.0; 0.0;
            6.574404; 0.764196; 1.595786; 0.0; 0.0; 3.392697; 1.939825; 0.700495; 0.052874; 9.754728; 0.0; 0.0; 3.48768;
            0.802873; 0.0; 2.043519; 0.0; 0.0; 3.921929; 1.959594; 0.0; 0.0; 0.3127; 0.0; 0.0; 2.983815; 0.0; 1.83189; 4.96362;
            5.83562; 0.0; 0.941767; 2.441303; 0.0; 0.334927; 1.322773; 2.651246; 0.0; 0.0; 0.0; 0.099974; 0.512102; 0.214042;
            0.711576; 1.192063; 2.142727; 0.0; 1.159921; 2.395009; 2.678757; 0.375006; 2.600225; 0.117835; 0.0; 2.957041;
            11.764845; 4.778022; 0.0; 0.064488; 2.119516; 0.0; 0.0; 0.202364; 4.43261; 1.040206; 0.081288; 4.686763; 2.584367;
            2.381216; 0.0; 0.031196; 0.537794; 4.35308; 2.031881; 1.740447; 1.296947; 0.0; 0.510216; 0.354055; 2.283879;
            3.607174; 0.0; 0.322203; 0.986396; 0.025445; 0.277725; 0.395995; 0.0; 2.941097; 0.98901; 0.181441; 0.0; 1.41994;
            0.0; 0.0; 1.872985; 1.565112; 2.497949; 0.0; 0.0; 0.312136; 1.431995; 2.096221; 1.882741; 0.0; 4.752209; 0.0;
            1.173437; 3.60508; 2.644296; 0.0; 0.0; 0.132417; 0.0; 0.0; 0.788642; 0.0; 0.0; 0.613296; 2.222431; 0.0; 1.664929;
            3.181359; 1.867732; 1.063342; 1.870465; 1.469997; 0.0; 0.245035; 2.612102; 0.0; 0.642298; 0.0; 4.092816; 1.507542;
            0.483124; 0.0; 0.0; 0.0; 0.0; 1.618722; 0.0; 1.212049; 1.579186; 0.0; 0.0; 2.884045; 0.0; 6.386427; 0.0; 0.227493;
            10.531423; 0.873409; 2.440361; 0.0; 3.157825; 0.0; 1.182512; 11.198405; 0.0; 2.27279; 0.298399; 1.624341; 0.0;
            4.756753; 0.556705; 2.043668; 4.110234; 4.228567; 0.0; 4.796603; 0.0; 1.350014; 0.0; 0.0; 0.0; 0.0; 2.920823; 0.0;
            0.0; 0.883152; 6.584653; 0.025556; 0.716609; 0.0; 8.783497; 0.641825; 1.524132; 0.0; 0.873738; 0.0; 2.6162;
            0.200578; 0.680832; 1.510734; 0.0; 0.0; 2.202807; 0.0; 0.0; 2.725498; 0.0; 8.552593; 0.0; 0.298605; 0.495948;
            0.0; 0.0; 0.197676; 0.0; 2.514259; 0.595785; 0.0; 0.626791; 0.0; 2.258893; 4.121462; 3.835965; 0.0; 0.672519;
            1.712922; 1.469087; 2.584109; 0.0; 3.22375; 0.0; 0.0; 2.200228; 0.0; 2.392011; 7.79248; 1.767081; 0.306815;
            1.572941; 0.0; 1.054235; 0.0; 1.429689; 0.099582; 0.0; 1.004675; 0.0; 2.93243; 0.635162; 1.19989; 3.597754;
            2.935238; 1.326419; 1.187039; 0.0; 5.292264; 0.0; 0.71673; 4.494548; 0.511313; 10.591884; 1.498288; 0.0; 0.0; 0.0;
            7.471386; 1.550287; 1.511363; 0.0; 0.0; 2.562011; 0.845157; 1.171172; 0.0; 12.579813; 0.0; 0.0; 4.535857; 0.030698;
            0.0; 2.078368; 1.421113; 0.0; 2.404634; 2.800321; 0.515842; 0.001106; 0.0; 0.0; 0.0; 3.005475; 0.0; 0.0; 4.466236;
            6.003285; 0.0; 0.0; 4.703087; 1.093062; 3.635723; 0.0; 0.0; 0.554857; 2.906585; 0.0; 0.0; 0.26434; 0.0; 1.868009;
            0.0; 0.767943; 0.0; 0.0; 0.311107; 3.068584; 0.039162; 4.032188; 4.318037; 0.0; 0.0; 2.085508; 0.0; 0.0; 1.32544;
            1.44314; 0.098778; 1.691968; 0.910304; 3.66684; 0.769006; 2.497978; 1.103729; 0.0; 0.770869; 0.0; 1.275256; 0.0;
            1.212942; 3.999437; 2.936755; 0.0; 0.0; 0.216068; 0.0; 2.618757; 2.350003; 2.97163; 1.477464; 1.003837; 4.87839;
            1.769495; 0.0; 2.109468; 2.59619; 0.0; 0.0; 3.734434; 7.18856; 4.173406; 1.484233; 0.0; 2.099632; 0.0; 0.0;
            1.082505; 0.699181; 0.0; 2.296291; 2.060906; 3.294362; 0.0; 4.71305; 9.021503; 0.0; 0.0; 0.075124; 0.149576;
            2.21331; 0.96611; 0.0; 0.101291; 7.055904; 3.12131; 1.621537; 3.487494; 0.867704; 0.862871; 0.0; 0.0; 0.0; 0.0;
            3.833386; 0.0; 0.046158; 4.818383; 0.0; 9.643603; 1.633465; 0.0; 2.952384; 0.513767; 0.0; 2.650039; 3.054194;
            13.118452; 2.230619; 0.484316; 1.363891; 0.419163; 3.972292; 0.0; 0.287555; 0.0; 0.0; 4.473218; 0.439954; 0.0; 0.0;
            0.0; 1.991387; 1.142356; 4.898578; 0.0; 0.0; 3.62821; 0.0
        ];
    end
end