%{
Area Under Curve (AUC) Analysis

This script analyzes and visualizes area under the curve metrics for nerve response data.
It plots relationships between frequency, ramp rate, and various AUC metrics.

Dependencies:
- './external_dependencies/gramm' - Gramm plotting library
- 'results/AUC_data.mat' - Contains preprocessed data and AUC calculations
%}

% Add dependency path
addpath('./external_dependencies/gramm');

% Load data
load('results/AUC_data.mat','preprocessed_data_table_sorted','all_AUC','upper_bounds_for_AUC','lower_bounds_for_AUC','all_AUC_noise');

% Select metric to analyze (0-4)
flag_metric = 3;
switch flag_metric
    case 0
        % Raw AUC 
        metric = real(all_AUC);
        y_scale = 'linear';
        y_label = 'AUC (N*s)';
    case 1
        % Normalize AUC by onset response duration
        % Results indicate that onset response strength is generally
        % smaller for slower ramp rates
        metric = all_AUC./(upper_bounds_for_AUC-lower_bounds_for_AUC);
        y_scale = 'log';
        y_label = 'AUC per onset duration (N*s/s)';
    case 2
        % Ratio of onset AUC to noise AUC
        metric = real(all_AUC./all_AUC_noise);
        y_scale = 'linear';
        y_label = 'AUC onset / AUC noise';
    case 3
        % Difference between onset AUC and noise AUC
        metric = real(all_AUC-all_AUC_noise);
        y_scale = 'linear';
        y_label = 'AUC onset - noise (N*s)';
    case 4
        % Absolute difference between onset AUC and noise AUC
        metric = abs(all_AUC-all_AUC_noise);
        y_scale = 'log';
        y_label = 'AUC onset - noise (N*s)';
    otherwise
        error('Invalid metric value');
end

%% Plot AUC vs. frequency at a fixed ramp rate (0.05 mA/sec)
figure('Position', [680, 345, 718, 533], 'Color', [1 1 1]);
g1 = gramm('x', preprocessed_data_table_sorted.frequency_kHz, ...
    'y', metric, ...
    'color', preprocessed_data_table_sorted.nerve_ID, ...
    'subset', preprocessed_data_table_sorted.ramp_rate_mA_per_sec == 0.05);
g1.set_names('x', 'frequency (kHz)', 'y', y_label, 'color', 'nerve');
g1.axe_property('FontSize', 14, 'YScale', y_scale);
g1.set_text_options("base_size", 14, "legend_title_scaling", 1, "legend_scaling", 1, "big_title_scaling", 1);
g1.set_order_options('x', 1);
g1.geom_line();
g1.geom_point();
g1.draw();

%% Plot AUC vs. ramp rate for different frequencies
% Define x-tick values (rounded to first decimal place of μA for display)
x_tick = round(0.2 ./ [64, 32, 16, 8, 4, 2, 1, 0.5, 0.25]*1e3, 1); 

figure('Position', [79.5, 559, 1302.5, 319]);
clearvars g2;
g2 = gramm('x', preprocessed_data_table_sorted.ramp_rate_mA_per_sec*1e3, ...
    'y', metric, ...
    'color', preprocessed_data_table_sorted.nerve_ID,...
    'subset', ~isinf(preprocessed_data_table_sorted.ramp_rate_mA_per_sec));
    
% Create facet grid by frequency (rounded to nearest 10 kHz)
g2.facet_grid([], round(preprocessed_data_table_sorted.frequency_kHz / 10) * 10,...
    'scale', 'independent');
g2.set_names('x', 'ramp rate (\muA/s)', 'y', y_label, 'column', 'kHz', 'color', 'nerve');
g2.axe_property('FontSize', 14, 'XScale', 'log', 'YScale', y_scale,...
    'XTick', x_tick, 'XLim', [min(x_tick), max(x_tick)]);
g2.set_text_options("base_size", 14, "legend_title_scaling", 1, "legend_scaling", 1, "big_title_scaling", 1, "Interpreter", "tex");
g2.geom_line();
g2.geom_point();
g2.draw();

% Ensure the zero point appears on the y axis and set consistent y-limits
for i = 1:length(g2.facet_axes_handles)
    axes(g2.facet_axes_handles(i));
    ylim_i = get(gca, 'YLim');
    % Set limits based on manually inspecting approximate negative and positive limits
    ylim(max(ylim_i)*[-0.05, 1.05]); 
    yline(0, 'k--');
end

%% Plot with transformed y-axis to better visualize data across orders of magnitude
% Transform the metric using logarithmic scaling above a threshold for better visualization
threshold_for_lin_to_log = 25; 
transformed_metric = threshold_for_lin_to_log*log(1 + metric/threshold_for_lin_to_log);

% Plot original vs transformed metrics for reference
figure; 
scatter(metric, transformed_metric); 
axis equal; 
hold on; 
% Plot diagonal reference line
x = linspace(min(metric), max(metric), 100);
y = x;
plot(x, y, 'k--');
title('Original vs Transformed Metric');
xlabel('Original Metric');
ylabel('Transformed Metric');

% Plot AUC vs. ramp rate with frequencies as colors, faceted by nerve
figure('Position', [79.5, 559, 1302.5, 319]);
clearvars g2;
g2 = gramm('x', preprocessed_data_table_sorted.ramp_rate_mA_per_sec*1e3, ...
    'y', transformed_metric, ...
    'color', round(preprocessed_data_table_sorted.frequency_kHz / 10) * 10,...
    'subset', ~isinf(preprocessed_data_table_sorted.ramp_rate_mA_per_sec));
g2.facet_grid([], preprocessed_data_table_sorted.nerve_ID,...
    'scale', 'independent');
g2.set_names('x', 'ramp rate (\muA/s)', 'y', y_label, 'column', 'nerve', 'color', 'kHz');
g2.axe_property('FontSize', 14, 'XScale', 'log', 'YScale', y_scale,...
    'XTick', x_tick, 'XLim', [min(x_tick), max(x_tick)]);
g2.set_text_options("base_size", 14, "legend_title_scaling", 1, "legend_scaling", 1, "big_title_scaling", 1, "Interpreter", "tex");
g2.geom_line();
g2.geom_point();
g2.set_color_options('map', 'matlab');
g2.draw();

% Configure y-axis display to show original values despite using transformed scale
for i = 1:length(g2.facet_axes_handles)
    axes(g2.facet_axes_handles(i));
    yline(0, 'k--');
    ylim_i = get(gca, 'YLim');
    ylim(max(ylim_i)*[-0.03, 1]);
    
    % Set custom tick values that correspond to meaningful original values
    desired_ytick_i = [-25, 0, 5.^[2:5]];
    ytick_i = sign(desired_ytick_i).*(threshold_for_lin_to_log*log(1 + abs(desired_ytick_i)/threshold_for_lin_to_log));
    set(gca, 'YTick', ytick_i);
    
    % Set y-axis limits based on desired original value range
    desired_bounds = [-40, 3125];
    ylim(sign(desired_bounds).*(threshold_for_lin_to_log*log(1 + abs(desired_bounds)/threshold_for_lin_to_log)));
    
    % Set tick labels to show original values instead of transformed values
    yticklabel_i = arrayfun(@(g) num2str(g, '%d'), desired_ytick_i, 'UniformOutput', false);
    set(gca, 'YTickLabel', yticklabel_i);
end