plotGRMvsGroundTruth.m

function plotGRMvsGroundTruth(GRM_values, fbeta_scores, p21_ratios, bhattacharyya_coeffs, js_divergences, is_best_so_far)
%plotGRMvsGroundTruth - Plot scatter plots of GRM vs ground truth metrics with correlation coefficient
%
% Creates scatter plots showing correlation between GRM (Graph Representativeness Metric)
% and various ground truth validation metrics to evaluate GRM effectiveness.
%
% Inputs:
%   GRM_values (vector): Actual GRM values from all BHO iterations (positive integer, higher is better)
%   fbeta_scores (vector): F-beta scores for each iteration
%   p21_ratios (vector): P21 ratios (detected/truth) for each iteration
%   bhattacharyya_coeffs (vector): Bhattacharyya coefficients for each iteration (0-1, higher=better)
%   js_divergences (vector): Jensen-Shannon divergences for each iteration (0-1, lower=better)
%   is_best_so_far (logical vector): Boolean array indicating which iterations were best so far
%
% Outputs:
%   Creates a figure with 4 subplots (2x2) showing correlations
%   - All points shown in gray for context (discarded solutions)
%   - Best-so-far points highlighted with orange markers
%   - TWO correlation coefficients displayed:
%     * Main title: r for ALL points (including discarded)
%     * Subtitle (in orange): r for BEST-SO-FAR points only
%   - No regression line fitting

% Create figure
f = figure();
f.Name = 'GRM vs Ground Truth Metrics Correlation';
f.WindowState = 'maximized';

% Ensure all inputs are column vectors
GRM_values = GRM_values(:);
fbeta_scores = fbeta_scores(:);
p21_ratios = p21_ratios(:);
bhattacharyya_coeffs = bhattacharyya_coeffs(:);
js_divergences = js_divergences(:);

% Get lengths of all arrays (before handling is_best_so_far)
lengths = [length(GRM_values), length(fbeta_scores), length(p21_ratios), ...
           length(bhattacharyya_coeffs), length(js_divergences)];

% Find minimum length to ensure all arrays are the same size
min_length = min(lengths);
max_length = max(lengths);

% Check if arrays have different lengths
if min_length ~= max_length
    warning('Metric arrays have different lengths. Truncating to minimum length (%d).', min_length);
end

% Truncate all arrays to the same length (in case of size mismatch)
if min_length > 0
    GRM_values = GRM_values(1:min_length);
    fbeta_scores = fbeta_scores(1:min_length);
    p21_ratios = p21_ratios(1:min_length);
    bhattacharyya_coeffs = bhattacharyya_coeffs(1:min_length);
    js_divergences = js_divergences(1:min_length);
else
    error('All metric arrays are empty. Cannot create correlation plot.');
end

% Handle is_best_so_far (optional parameter for backward compatibility)
if nargin < 6 || isempty(is_best_so_far)
    % If not provided, assume all iterations are best so far (for backward compatibility)
    is_best_so_far = true(size(GRM_values));
else
    is_best_so_far = is_best_so_far(:);
    % Ensure it's the same length as other arrays
    if length(is_best_so_far) ~= min_length
        warning('is_best_so_far length (%d) does not match metric arrays (%d). Using all points as best-so-far.', ...
            length(is_best_so_far), min_length);
        is_best_so_far = true(size(GRM_values));
    end
    % Ensure it's a logical array
    if ~islogical(is_best_so_far)
        is_best_so_far = logical(is_best_so_far);
    end
end

% Extract best-so-far points (only if we have valid logical indices)
if any(is_best_so_far)
    best_GRM = GRM_values(is_best_so_far);
    best_fbeta = fbeta_scores(is_best_so_far);
    best_p21 = p21_ratios(is_best_so_far);
    best_bhatt = bhattacharyya_coeffs(is_best_so_far);
    best_js = js_divergences(is_best_so_far);
else
    % If no best-so-far points, use all points
    warning('No best-so-far points found. Using all points for correlation.');
    best_GRM = GRM_values;
    best_fbeta = fbeta_scores;
    best_p21 = p21_ratios;
    best_bhatt = bhattacharyya_coeffs;
    best_js = js_divergences;
end

% Calculate correlation coefficients for ALL points
corr_fbeta_all = calculateCorrelation(GRM_values, fbeta_scores);
corr_p21_all = calculateCorrelation(GRM_values, p21_ratios);
corr_bhatt_all = calculateCorrelation(GRM_values, bhattacharyya_coeffs);
corr_js_all = calculateCorrelation(GRM_values, js_divergences);

% Calculate correlation coefficients for BEST-SO-FAR points only
corr_fbeta_best = calculateCorrelation(best_GRM, best_fbeta);
corr_p21_best = calculateCorrelation(best_GRM, best_p21);
corr_bhatt_best = calculateCorrelation(best_GRM, best_bhatt);
corr_js_best = calculateCorrelation(best_GRM, best_js);

% Count valid points
n_all = length(GRM_values);
n_best = length(best_GRM);

% Subplot 1: GRM vs F-beta Score
subplot(2, 2, 1);
% Plot all points in darker gray for context
scatter(GRM_values, fbeta_scores, 50, 'filled', 'MarkerFaceColor', [0.5 0.5 0.5], 'MarkerEdgeColor', [0.4 0.4 0.4], 'MarkerFaceAlpha', 0.6);
hold on;
% Highlight best-so-far points with same size, using P21 Ratio color
scatter(best_GRM, best_fbeta, 50, 'filled', 'MarkerFaceColor', [0.8 0.4 0.2], 'MarkerEdgeColor', [0.6 0.2 0.1], 'LineWidth', 1.5);
hold off;
xlabel('GRM (higher is better)', 'FontSize', 12);
ylabel('F-beta Score', 'FontSize', 12);
title(sprintf('GRM vs F-beta Score\nr = %.4f (n=%d)', corr_fbeta_all, n_all), 'FontSize', 14);
% Add subtitle with best-so-far correlation
subtitle(sprintf('Best so far: r = %.4f (n=%d)', corr_fbeta_best, n_best), 'FontSize', 11, 'Color', [0.8 0.4 0.2]);
set(gca, 'FontSize', 14);
grid on;

% Subplot 2: GRM vs P21 Ratio
subplot(2, 2, 2);
% Plot all points in darker gray for context
scatter(GRM_values, p21_ratios, 50, 'filled', 'MarkerFaceColor', [0.5 0.5 0.5], 'MarkerEdgeColor', [0.4 0.4 0.4], 'MarkerFaceAlpha', 0.6);
hold on;
% Highlight best-so-far points with same size, using P21 Ratio color
scatter(best_GRM, best_p21, 50, 'filled', 'MarkerFaceColor', [0.8 0.4 0.2], 'MarkerEdgeColor', [0.6 0.2 0.1], 'LineWidth', 1.5);
hold off;
xlabel('GRM (higher is better)', 'FontSize', 12);
ylabel('P_{21} Ratio (Detected/Truth)', 'FontSize', 12);
title(sprintf('GRM vs P_{21} Ratio\nr = %.4f  (n=%d)', corr_p21_all, n_all), 'FontSize', 14);
% Add subtitle with best-so-far correlation
subtitle(sprintf('Best so far: r = %.4f (n=%d)', corr_p21_best, n_best), 'FontSize', 11, 'Color', [0.8 0.4 0.2]);
set(gca, 'FontSize', 14);
grid on;

% Subplot 3: GRM vs Bhattacharyya Coefficient
subplot(2, 2, 3);
% Check if we have valid data (not all NaN)
valid_bhatt = ~isnan(GRM_values) & ~isnan(bhattacharyya_coeffs);
valid_best_bhatt = ~isnan(best_GRM) & ~isnan(best_bhatt);
if any(valid_bhatt)
    % Plot all points in darker gray for context
    scatter(GRM_values(valid_bhatt), bhattacharyya_coeffs(valid_bhatt), 50, 'filled', 'MarkerFaceColor', [0.5 0.5 0.5], 'MarkerEdgeColor', [0.4 0.4 0.4], 'MarkerFaceAlpha', 0.6);
    hold on;
    % Highlight best-so-far points with same size, using P21 Ratio color
    if any(valid_best_bhatt)
        scatter(best_GRM(valid_best_bhatt), best_bhatt(valid_best_bhatt), 50, 'filled', 'MarkerFaceColor', [0.8 0.4 0.2], 'MarkerEdgeColor', [0.6 0.2 0.1], 'LineWidth', 1.5);
    end
    hold off;
else
    text(0.5, 0.5, 'No valid data', 'HorizontalAlignment', 'center', 'Units', 'normalized');
end
xlabel('GRM (higher is better)', 'FontSize', 12);
ylabel('Bhattacharyya Coefficient', 'FontSize', 12);
if isnan(corr_bhatt_all)
    title('GRM vs Bhattacharyya Coefficient\nr = NaN (no valid data)', 'FontSize', 14);
else
    title(sprintf('GRM vs Bhattacharyya Coefficient\nr = %.4f (n=%d)', corr_bhatt_all, sum(valid_bhatt)), 'FontSize', 14);
end
% Add subtitle with best-so-far correlation
if ~isnan(corr_bhatt_best)
    subtitle(sprintf('Best so far: r = %.4f (n=%d)', corr_bhatt_best, sum(valid_best_bhatt)), 'FontSize', 11, 'Color', [0.8 0.4 0.2]);
end
set(gca, 'FontSize', 14);
grid on;

% Subplot 4: GRM vs JS Divergence
% Note: JS Divergence is "lower is better", so we might want to show negative correlation
% or convert to similarity (1 - JS_divergence). For now, showing raw JS divergence.
subplot(2, 2, 4);
% Check if we have valid data (not all NaN)
valid_js = ~isnan(GRM_values) & ~isnan(js_divergences);
valid_best_js = ~isnan(best_GRM) & ~isnan(best_js);
if any(valid_js)
    % Plot all points in darker gray for context
    scatter(GRM_values(valid_js), js_divergences(valid_js), 50, 'filled', 'MarkerFaceColor', [0.5 0.5 0.5], 'MarkerEdgeColor', [0.4 0.4 0.4], 'MarkerFaceAlpha', 0.6);
    hold on;
    % Highlight best-so-far points with same size, using P21 Ratio color
    if any(valid_best_js)
        scatter(best_GRM(valid_best_js), best_js(valid_best_js), 50, 'filled', 'MarkerFaceColor', [0.8 0.4 0.2], 'MarkerEdgeColor', [0.6 0.2 0.1], 'LineWidth', 1.5);
    end
    hold off;
else
    text(0.5, 0.5, 'No valid data', 'HorizontalAlignment', 'center', 'Units', 'normalized');
end
xlabel('GRM (higher is better)', 'FontSize', 12);
ylabel('JS Divergence (lower is better)', 'FontSize', 12);
if isnan(corr_js_all)
    title('GRM vs JS Divergence\nr = NaN (no valid data)', 'FontSize', 14);
else
    title(sprintf('GRM vs JS Divergence\nr = %.4f (n=%d)', corr_js_all, sum(valid_js)), 'FontSize', 14);
end
% Add subtitle with best-so-far correlation
if ~isnan(corr_js_best)
    subtitle(sprintf('Best so far: r = %.4f (n=%d)', corr_js_best, sum(valid_best_js)), 'FontSize', 11, 'Color', [0.8 0.4 0.2]);
end
set(gca, 'FontSize', 14);
grid on;
% Note: Negative correlation is expected (higher GRM → lower JS divergence = better)

% Add overall title
sgtitle('GRM Correlation Analysis: Evaluating Objective Function Effectiveness', ...
    'FontSize', 14, 'FontWeight', 'bold');

% Add legend to the right of the figure (moved further away from subplots)
% Create dummy scatter plots for legend (invisible, outside plot area)
ax_legend = axes('Position', [0.95 0.3 0.04 0.4], 'Visible', 'off');
hold(ax_legend, 'on');
% Best so far Solutions
h1 = scatter(ax_legend, NaN, NaN, 50, 'filled', 'MarkerFaceColor', [0.8 0.4 0.2], ...
    'MarkerEdgeColor', [0.6 0.2 0.1], 'LineWidth', 1.5);
% Discarded Solutions
h2 = scatter(ax_legend, NaN, NaN, 50, 'filled', 'MarkerFaceColor', [0.5 0.5 0.5], ...
    'MarkerEdgeColor', [0.4 0.4 0.4], 'MarkerFaceAlpha', 0.6);
hold(ax_legend, 'off');
% Create legend
legend(ax_legend, [h1, h2], {'Best so far Solutions', 'Discarded Solutions'}, ...
    'Location', 'north', 'FontSize', 12, 'Box', 'off');

end

function r = calculateCorrelation(x, y)
%calculateCorrelation - Calculate Pearson correlation coefficient
%
% Returns the Pearson correlation coefficient between x and y
% r ranges from -1 to 1, where:
%   r = 1: perfect positive correlation
%   r = -1: perfect negative correlation
%   r = 0: no correlation

% Ensure x and y are column vectors and have the same length
x = x(:);
y = y(:);

% Check if arrays have compatible sizes
if length(x) ~= length(y)
    r = NaN;
    return;
end

% Remove NaN values
valid = ~isnan(x) & ~isnan(y);
x = x(valid);
y = y(valid);

if length(x) < 2
    r = NaN;
    return;
end

% Calculate Pearson correlation coefficient using corrcoef
corr_matrix = corrcoef(x, y);
r = corr_matrix(1, 2);

end