classdef quadratic_bilevel
    %{
    Upper-level
    ===========
    minimise    1/2 (x' P_upper x) + 1/2 (y' Q_upper y) + (x' R_upper y) + c_upper' x + d_upper' y
    subject to  A_upper * x + B_upper * y  >= rhs_upper
                x >= x_lower_bound
                x <= x_upper_bound
                y solves lower-level
    
    Lower-level
    ===========
    minimise    1/2 (x' P_lower x) + 1/2 (y' Q_lower y) + 1/2 (x' R_lower y) + c_lower' x + d_lower' y
    subject to  A_lower * x + B_lower * y >= rhs_lower
                y >= y_lower_bound
                y <= y_upper_bound
    %}
    properties(Constant)
        name = 'quadratic_bilevel';
        category = 'archetype';
        subcategory = '';
        datasets = {
            'aiyoshi_shimizu_1984_ex2.json'; 'an_et_al_2009.json'; 'bard_1988_ex1.json';
            'bard_1988_ex2.json'; 'bard_1991_ex1.json'; 'bard_book_1998.json';
            'calamai_vicente_1994a.json'; 'calamai_vicente_1994b.json'; 'calamai_vicente_1994c.json';
            'clark_westerberg_1990a.json'; 'dempe_etal_2012.json'; 'dempe_franke_2011_ex41.json';
            'dempe_franke_2011_ex42.json'; 'dempe_franke_2014_ex38.json'; 'dempe_lohse_2011_ex31a.json';
            'dempe_lohse_2011_ex31b.json'; 'desilva_1978.json'; 'falk_liu_1995.json';
            'floudas_etal_2013.json'; 'gumus_floudas_2001_ex.json'; 'hatz_etal_2013.json';
            'henderson_quandt_1958.json'; 'henrion_surowiec_2011.json'; 'lampariello_sagratella_2017_ex31.json';
            'lampariello_sagratella_2017_ex32.json'; 'lampariello_sagratella_2017_ex33.json';
            'lampariello_sagratella_2017_ex35.json'; 'lucchetti_etal_1987.json'; 'macal_hurter_1997.json';
            'morgan_patrone_2006a.json'; 'muu_quy_2003_ex1.json'; 'muu_quy_2003_ex2.json';
            'outrata_1990_ex1_a.json'; 'outrata_1990_ex1_b.json'; 'outrata_1990_ex1_c.json';
            'outrata_1990_ex1_d.json'; 'outrata_1990_ex1_e.json'; 'outrata_cervinka_2009.json';
            'sahin_ciric_1998_ex2.json'; 'shimizu_aiyoshi_1981_ex1.json'; 'shimizu_aiyoshi_1981_ex2.json';
            'shimizu_etal_1997a.json'; 'sinha_malo_deb_2014_tp6.json'; 'tuy_etal_2007.json';
            'wan_wang_lv_2011.json'; 'yezza_1996_ex3.json'; 'yezza_1996_ex4.json';
            % Toll settings
            'toll_setting_p1.json'; 'toll_setting_p2.json'; 'toll_setting_p3.json';
            'toll_setting_p4.json'; 'toll_setting_p5.json';
        };
        paths = fullfile('bolib3', 'data', 'quadratic_bilevel', quadratic_bilevel.datasets);
        x0 = [0.2; 2.0];
        y0 = [4.0; 4.6];
        required_data_fields = {
            'n_x', 'n_y', ...
            'A_upper', 'B_upper', 'c_upper', 'd_upper', 'rhs_upper', 'P_upper', 'Q_upper', 'R_upper', ...
            'A_lower', 'B_lower', 'c_lower', 'd_lower', 'rhs_lower', 'P_lower', 'Q_lower', 'R_lower', ...
            'x_upper_bound', 'x_lower_bound', ...
            'y_upper_bound', 'y_lower_bound'
        };
    end


    methods(Static)

        % Upper-level objective function
        % 1/2 (x' P_upper x) + 1/2 (y' Q_upper y) + (x' R_upper y) + c_upper' x + d_upper' y
        function obj = F(x, y, data)
            obj = 0.5 * x' * data.P_upper * x + ...
                  0.5 * y' * data.Q_upper * y + ...
                  1.0 * x' * data.R_upper * y + ...
                  dot(data.c_upper, x) + ...
                  dot(data.d_upper, y);
        end
        
        % Upper-level inequality constraints
        % Ax + By  >= rhs
        % x        >= x_lower_bound
        % x        <= x_upper_bound
        function ineq = G(x, y, data)
            ineq = [];
            if ~isempty(data.rhs_upper)
                ineq = [ineq; data.A_upper * x + data.B_upper * y - data.rhs_upper];
            end
            if ~isempty(data.x_lower_bound)
                ineq = [ineq; x - data.x_lower_bound];
            end
            if ~isempty(data.x_upper_bound)
                ineq = [ineq; data.x_upper_bound - x];
            end
        end
        
        % Upper-level equality constraints
        function val = H(~, ~, ~)
            val = [];
        end 
        
        % Lower-level objective function
        % minimise f(x,y) := c'x + d'y
        function obj = f(x, y, data)
            obj = 0.5 * x' * data.P_lower * x +...
                  0.5 * y' * data.Q_lower * y +...
                  1.0 * x' * data.R_lower * y +...
                  dot(data.c_lower, x) +... 
                  dot(data.d_lower, y);
        end
        
        % Lower-level inequality constraints
        function ineq = g(x, y, data)
            ineq = [];
            if ~isempty(data.rhs_lower)
                ineq = [ineq; data.A_lower * x + data.B_lower * y - data.rhs_lower];
            end
            if ~isempty(data.y_lower_bound)
                ineq = [ineq; y - data.y_lower_bound];
            end
            if ~isempty(data.y_upper_bound)
                ineq = [ineq; data.y_upper_bound - y];
            end
        end
        
        % Lower-level equality constraints
        function val = h(~, ~, ~)
            val = [];
        end

        % If the bilevel program is parameterized by data, this function should
        %  provide code to read data file and return an appropriate structure.
        function data = read_data(filepath)
            % Read JSON file
            file = fopen(filepath, 'r');
            raw = fread(file, inf, 'uint8=>char')';
            fclose(file);
            data = jsondecode(raw);

            % Add a field for the filepath
            data.filepath = filepath;

            % Bounds are optional.  If they are not provided we replace them with an empty array
            bounds = {'x_upper_bound', 'x_lower_bound', 'y_upper_bound', 'y_lower_bound'};
            for i = 1:length(bounds)
                if ~isfield(data, bounds{i})
                    data.(bounds{i}) = [];
                end
            end

            % Validate required fields
            for i = 1:length(quadratic_bilevel.required_data_fields)
                key = quadratic_bilevel.required_data_fields{i};
                assert(isfield(data, key), ...
                    sprintf("Missing required field '%s' in '%s'.", key, filepath));
            end
        end


        % Key are the function/variable names
        % Values are their dimention
        function n = dimention(key, data)
            n = dictionary( ...
                'x', data.n_x, ...
                'y', data.n_y, ...
                'F', 1, ...
                'G', length(data.rhs_upper) + length(data.x_upper_bound) + length(data.x_lower_bound), ...
                'H', 0, ...
                'f', 1, ...
                'g', length(data.rhs_lower) + length(data.y_upper_bound) + length(data.y_lower_bound), ...
                'h', 0 ...
            );
            if isKey(n,key)
                n = n(key);
            end
        end
    end
end