function subselection(varargin)
    %% subselection(varargin)
    % Function used to select a uniformly spaced subset of a given set of
    % surfaces.
    %
    % This function will read the adjacency matrix located on the relative path
    % './result/intermediate/indicator_response.csv' and use it to compute which
    % surfaces should be included in the optimal set.
    %
    %% Dependencies
    % Here we list the dependencies needed for the subselection
    % - CVX                     Version: 2.2, Build 1148
    % - Mosek                   Version: 9.2
    %
    %% Adjacency input matrix
    % The input matrix contain rows corresponding to each surface and columns
    % corresponding to each probe point. Probes are in the same order as in the
    % './result/intermediate/probes.csv' file.
    %
    %% Optional parameter values.
    % - TargetDensity 1.0           Target density to reach by summed kernels.
    % - Debug         false         Control of the outputs.
    % - probeExclude  true          Remove probes in low/high density
    % - OptionsFile   'options.cfg' Name of the options file
    %
    %% Output Files
    % The function returns after exporting the results.
    % - weights_computed.csv    List of the computed weights.

    %% ====================================================================== %%
    %% Ensure cvx is loaded
    if (~exist('cvx_begin', 'file'))

        if (exist('load_lib', 'file'))
            load_lib('cvx');
        elseif (exist('cvx_startup', 'file'))
            cvx_startup;
        end

        if (~exist('cvx_begin', 'file'))
            error('    CVX have not been loaded.');
        end

    end

    if (~exist('Options', 'class'))
        [path, ~, ~] = fileparts(mfilename('fullpath'));
        addpath([path, '/matlab']);
    end

    %% ====================================================================== %%
    %% Input handling

    % Startup the input parser structure
    p = inputParser;

    % List of optional parameters
    addParameter(p, 'Debug', false, @validateBool);
    addParameter(p, 'OptionFile', 'options.cfg');
    addParameter(p, 'probeExclude', true, @validateBool);
    addParameter(p, 'TargetDensity', 1.0, @validateDensity);

    parse(p, varargin{:});

    Debug = p.Results.Debug;
    OptionFile = p.Results.OptionFile;
    probeExclude = p.Results.probeExclude;
    TargetDensity = p.Results.TargetDensity;

    % Read Options file and extract what is needed
    if exist(OptionFile, 'file')
        opt = Options(OptionFile);
        % Debug = opt.get('debug', false);
        probeMinDen = opt.get('threshold_void', 0.05);
        probeMaxDen = opt.get('threshold_solid', 0.95);
    end

    clear p;

    %% ====================================================================== %%
    %% Load csv-files and prepare matrices
    t_algorithm = tic;

    Probes = csvread("result/intermediate/probes.csv");
    in = csvread("result/intermediate/indicator_response.csv");
    kernel_matrix_cpp = sparse( ...
        in(:, 1) + 1, in(:, 2) + 1, in(:, 3), ... Insert the sparse elements
        max(in(:, 1)) + 1, size(Probes, 1) ... Size of the full array
    );

    fixed_idx = [];

    if (exist("result/intermediate/fixed.csv", 'file'))
        fixed_idx = csvread("result/intermediate/fixed.csv");

        % Ensure fixed_idx have correct orientation and correct for 0 indexing.
        fixed_idx = fixed_idx +1;

        if (size(fixed_idx, 1) > 1)
            fixed_idx = fixed_idx';
        end

    end

    % Create matrix for optimization
    if (size(kernel_matrix_cpp, 2) == 1)
        kernel_matrix_cpp = reshape(kernel_matrix_cpp, ...
            size(Probes, 1), (size(kernel_matrix_cpp, 1) / size(Probes, 1)));
    elseif (size(kernel_matrix_cpp, 2) == size(Probes, 1))
        kernel_matrix_cpp = kernel_matrix_cpp';
    end

    dimS = max(kernel_matrix_cpp(:));

    kernelmtx_optim = [];

    for i = 1:dimS
        kernelmtx_optim = [kernelmtx_optim; double(kernel_matrix_cpp == i)];
    end

    %% ====================================================================== %%
    %% Exclude probe points
    used_columns = 1:size(kernelmtx_optim, 2);
    columns = 1:size(kernelmtx_optim, 2);

    if (probeExclude && exist("data/width.csv", 'file'))

        information = csvread("data/information.csv");
        widths = csvread("data/width.csv");

        n_dim = information(1:3); nelx = n_dim(1); nely = n_dim(2); nelz = n_dim(3);

        % Sets the probe points to zero, i.e. the energy disappears
        kernelmtx_optim = exclude_probe_points3D(kernelmtx_optim, ...
        widths, Probes, nelx, nely, nelz, probeMaxDen, probeMinDen);

        % Excludes probepoints, where the energy is actually zero)
        kernelmtx_optim(sum(kernelmtx_optim, 2) == 0, :) = [];
        used_columns(sum(kernelmtx_optim, 1) == 0) = [];

        % Exclude streamsurfaces, which do not contribute
        kernelmtx_optim = kernelmtx_optim(:, sum(kernelmtx_optim, 1) ~= 0);
    end

    %% ====================================================================== %%
    %% Least absolute deviations
    weight_vector = ones(size(kernelmtx_optim, 1), 1);
    w = TargetDensity * weight_vector;
    [~, m] = size(kernelmtx_optim);

    %% ====================================================================== %%
    %% Setup fixed surfaces

    % Compute energy contribution of fixed surfaces
    b_boundary = kernelmtx_optim(:, fixed_idx) * ones(size(fixed_idx, 2), 1);
    % Get number of remaining surfaces
    m_relax = m - size(fixed_idx, 2);
    % Find surfaces that are not being fixed
    notFixedIndices = (1:m);
    notFixedIndices = notFixedIndices(setdiff(1:numel(notFixedIndices), fixed_idx));
    % Retain only non-fixed surfaces
    kernelmtx_optim_relax = kernelmtx_optim(:, notFixedIndices);

    % Solve relaxed optimization problem
    if (Debug); fprintf("-- Solving relaxed problem \n"); end
    t = tic;

    % Solves min norm(Ax-b), where b = w, A = kernel_mtx..., x = weights_...
    cvx_begin
    cvx_quiet (~Debug)
    cvx_solver mosek
    variable weights_relax(m_relax, 1)
    minimize norm((kernelmtx_optim_relax * weights_relax) + b_boundary - w, 1)
    subject to
    0 <= weights_relax <= 1
    cvx_end

    % Round upto one decimal
    weightsr = round(weights_relax, 1);

    % Compute the energy of the streamlines that are chosen in the relaxed
    % problem (i.e. where weightsr == 1)
    b_fixed = kernelmtx_optim_relax(:, weightsr == 1) * weightsr(weightsr == 1);

    % Add fixed index surfaces too integer problem
    b_fixed = b_fixed + b_boundary;

    % zweights holds the information about streamlines that are chosen or
    % discarded in the relaxed problem
    zweights = double((weightsr == 1)) + double((weightsr == 0));

    % Number of intermediate density streamlines
    m_inter = m_relax - sum(zweights);

    if (Debug); fprintf("-- Number of remaining variables: %i \n", m_inter); end

    % Create matrix for integer optimization problem
    kernelmtx_optim_inter = kernelmtx_optim_relax(:, zweights ~= 1);

    % Solve integer problem
    if (Debug); fprintf("-- Solving integer problem \n"); end
    cvx_begin
    cvx_quiet (Debug == false)
    cvx_solver mosek
    variable weights_inter(m_inter, 1) binary
    minimize norm((kernelmtx_optim_inter * weights_inter) + b_fixed - w, 1)
    subject to
    0 <= weights_inter <= 1
    cvx_end

    % Create final weights
    weights = weightsr;
    weights(zweights ~= 1) = weights_inter;

    % Ensure empty surfaces are not selected
    weights(sum(kernel_matrix_cpp, 1) < 1) = 0;

    % Add in fixed surfaces
    weightsFixed = zeros(size(weights, 1) + size(fixed_idx, 2), 1);
    weightsFixed(fixed_idx) = 1;
    weightsFixed(notFixedIndices) = weights;

    % Overwrite weights
    weights = weightsFixed;

    % Create boolean
    weights = round(weights);

    if (Debug); fprintf("-- Optimization concluded in %.1f seconds \n", toc(t)); end

    %% ====================================================================== %%
    %% NEW: correct for excluded columns
    used_columns_activated = used_columns(weights >= 0.5);
    weights_corrected = zeros(size(columns, 2), 1);
    weights_corrected(used_columns_activated) = 1;

    %% ====================================================================== %%
    %% Save data and exit
    csvwrite("result/intermediate/weights_computed.csv", find(weights == 1) - 1);
    fprintf("Number of surfaces selected: %i \n", sum(weights));
    fprintf("Algorithm concluded in %.1f seconds \n", toc(t_algorithm));

end

%% ========================================================================== %%
%% Validation functions
function output = validateFolder(input)

    if (~exist(input, 'dir'))
        error('Folder does not exist');
    end

    if (input(end) ~= '/')
        error('Please end paths with "/"');
    end

    output = true;
end

function output = validateDensity(input)

    if (~isnumeric(input))
        error('Must be a number.');
    end

    if (input < 0.0 || input > 1.0)
        error('Must be in interval [0,1].')
    end

    output = true;
end

function output = validateBool(input)

    if (input ~= true && input ~= false)
        error('Must be a boolean value.');
    end

    output = true;
end
