classdef brotcorne2001
    %{
        A bilevel model for toll optimization on a multicommodity transportation network.
        Luce Brotcorne, Martine Labbe, Patrice Marcotte, and Gilles Savard. 
        Transportation Science, 35(4):345–358, 2001.
    %}
    properties(Constant)
        name = 'brotcorne2001';
        category = 'transportation';
        subcategory = '';
        datasets = {
            fullfile('bolib3', 'data', 'network_brotcorne_fig1.csv');
            fullfile('bolib3', 'data', 'network_highway.csv');
            fullfile('bolib3', 'data', 'network_tryangle345.csv');
        };
        x0 = [0.0, 0.0, 0.0, 0.0, 5.0, 0.0, -2.0];
        y0 = [0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 1.0];
    end


    methods(Static)
        % Upper-level objective function (non-convex)
        function val = F(x, y, data)
            y2d = reshape(y, [data.num_arcs, data.num_commodities]);
            val = -sum(x*y2d);
        end
        
        % Upper-level inequality constraints
        function val = G(x, ~, data)
            val = [
            data.max_toll - x.';
            data.max_toll + x.'
        ];
        end
        
        % Upper-level equality constraints
        function val = H(~, ~, ~)
            val = [];
        end
        
        % Lower-level objective function
        function val = f(x, y, data)
            y2d = reshape(y, [data.num_arcs, data.num_commodities]);
            val = sum((data.cost.' + x) * y2d);
        end
        
        % Lower-level inequality constraints
        function val = g(~, y, ~)
            val = y.';
        end
        
        % Lower-level equality constraints
        function val = h(~, y, data)
            y2d = reshape(y, [data.num_arcs, data.num_commodities]);
            val = ((data.outflow * y2d) - (data.inflow * y2d) - data.supply_demand);
            val = reshape(val, [], 1);
        end

        function data = read_data(filepath)
            % Read the entire file
            file = fopen(filepath, 'r');
            content = fread(file, '*char')';
            fclose(file);
        
            % Split content into node and arc parts at blank line
            parts = regexp(strtrim(content), '\n\s*\n', 'split');
            if length(parts) ~= 2
                error('File format incorrect: expected two tables separated by a blank line.');
            end
            node_part = parts{1};
            arc_part = parts{2};
        
            % Write each part to temporary files to read as tables
            tmp_node = [tempname, '_nodes.csv'];
            tmp_arc = [tempname, '_arcs.csv'];
            fid = fopen(tmp_node, 'w'); fprintf(fid, '%s', node_part); fclose(fid);
            fid = fopen(tmp_arc, 'w'); fprintf(fid, '%s', arc_part); fclose(fid);
        
            % Read the tables
            df_nodes = readtable(tmp_node);
            df_arcs = readtable(tmp_arc);
        
            % Delete temporary files
            delete(tmp_node);
            delete(tmp_arc);
        
            % Check expected columns
            if ~ismember('node', df_nodes.Properties.VariableNames)
                error('Node column is missing from %s.', filepath);
            end
            required_arc_cols = {'arc', 'from', 'to', 'cost', 'max_toll'};
            if ~all(ismember(required_arc_cols, df_arcs.Properties.VariableNames))
                error('One or more of the required columns are missing from %s.', filepath);
            end
        
            % Dimensions
            num_commodities = width(df_nodes) - 1;
            num_nodes = height(df_nodes);
            num_arcs = height(df_arcs);
        
            % Mapping from node names to indices
            % It is important we consider all node names as strings even if
            % they're numerical!
            node_names = string(df_nodes.node);
            node_to_index = containers.Map(node_names, 1:num_nodes);
        
            inflow = zeros(num_nodes, num_arcs);
            outflow = zeros(num_nodes, num_arcs);
        
            for arc_idx = 1:num_arcs
                i = node_to_index(string(df_arcs.from(arc_idx)));
                j = node_to_index(string(df_arcs.to(arc_idx)));
                outflow(i, arc_idx) = 1;
                inflow(j, arc_idx) = 1;
            end
        
            data = struct();
            data.num_commodities = num_commodities;
            data.num_nodes = num_nodes;
            data.num_arcs = num_arcs;
            data.inflow = inflow;
            data.outflow = outflow;
            data.cost = df_arcs.cost;
            data.max_toll = df_arcs.max_toll;
            data.supply_demand = table2array(df_nodes(:, 2:end));
        end

        % Key are the function/variable names
        % Values are their dimention
        function n = dimentions(key, data)
            n = dictionary( ...
                'x', data.num_arcs, ...
                'y', data.num_arcs * data.num_commodities, ...
                'F', 1, ...
                'G', data.num_arcs * 2, ...
                'H', 0, ...
                'f', 1, ...
                'g', data.num_arcs * data.num_commodities, ...
                'h', data.num_nodes * data.num_commodities ...
            );
            if isKey(n,key)
                n = n(key);
            end
        end
    end
end