% Get results from "largepool_simulation" folder, which is generated by the "largepool_simulation.m"
%

clc
addpath('./Tools/');

strategy_colors=[0,0.8,0;
    1,0.8,0;
    0.8,0,0];

% only get folders matching these requirement
N_sid=50;
N_spe=50;
R_sup=3;
typstr=[num2str(R_sup),'_',num2str(N_spe),'_',num2str(N_sid)];

%%%%%%%%%%%%%%%%%%% Data Loading
fprintf('Loading....\n')
tic
% Specify the base path to the directory containing subfolders with simulation data.
largefolderpath = './largepool_simulation/';
% Get all files and directories in the largefolderpath that match the pattern '*_typstr' (e.g., folders ending with the parameter string).
files = dir([largefolderpath, 'v2_largepool_simulation_nosaving_*', typstr]);
% Identify which entries are actual subfolders (exclude '.' and '..').
isFolder = [files.isdir] & ~ismember({files.name}, {'.', '..'});
subfolders = files(isFolder);

% Initialize empty matrices: simulation_data_matrix for numerical data, simulation_filestr_data for text data (e.g., labels or headers from CSVs).
simulation_data_matrix = [];
simulation_filestr_data = [];
% Loop over each matching subfolder to combine data from their CSV files.
for f = 1:length(subfolders)
    subfoldername = subfolders(f).name;

    % get the csv file in the folder
    fs=dir( [largefolderpath, subfolders(f).name, '/ResultSummary*.csv'])
    for id=1:length(fs)
        resultf=fs(id);
        if ~isempty(resultf) & ~ismember(resultf.name,'(')
        resultfilepath = [largefolderpath, subfolders(f).name, '/', resultf.name];

        % Import the CSV data: rawdata.data contains numerical values, rawdata.textdata contains text (e.g., row labels, headers).
        rawdata = importdata(resultfilepath);

        if ~isempty(rawdata)
            size(rawdata.data)
            % Append the numerical data vertically to the combined matrix.
            simulation_data_matrix = [simulation_data_matrix; rawdata.data];

            % Create a new column with 'resultstr' repeated for each row in rawdata.textdata (including header if present).
            num_rows = size(rawdata.textdata, 1);
            new_col = repmat({subfolders(f).name}, num_rows, 1);
            % Horizontally concatenate the new column with the original textdata.
            temp_textdata = [new_col, rawdata.textdata];
            % Append vertically to the combined text data matrix.
            simulation_filestr_data = [simulation_filestr_data; temp_textdata];
        end
    end
    end
end
toc
fprintf('Data Count per condition: %f\n',size(simulation_data_matrix,1)/(N_sid*6))


%
%
% The biomass distribution of surviving species (maintaining a time-averaged biomass greater than  10^-6 ) in log-log scale

% Extinction threshold for species, set to 100 times the migration-in rate (10^-8)
extinction_thresh = 1e-6;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Column indices in the simulation_data_matrix
CRcol = 4;              % Column for cheating receptors count
cheaternumcol = 3;    % Column for pure-cheater ratio
biomasscols = 11 + (1:N_spe);  % Columns for species biomasses (N_spe species)
issteadycol = 6;        % Column indicating if steady state is reached
initialcols = (12 + N_spe + N_sid) + (1:(N_spe + N_sid + 1));  % Columns for initial conditions (not used here)

% Extract biomass data and compute totals and surviving species
biomasses_mat = simulation_data_matrix(:, biomasscols);
biomass_tot = sum(biomasses_mat, 2);
surviving_spe = sum(biomasses_mat > extinction_thresh, 2);
notsteady = simulation_data_matrix(:, issteadycol) < 0.9;  % Flag for non-steady states

% test the right threshold for extinction
figure
hist(log10(biomass_tot),100);
xlabel('Total Biomass');

extinction_thresh_totalbiomass=10^-0.5;

figure
datatohist=biomasses_mat(:);
datatohist=datatohist(datatohist>10^-6);
[histbar,histloc]=hist(log10(biomasses_mat(:)),100);
histbar=histbar/sum(histbar);
scatter(histloc,histbar,'k','filled');
xlabel('Surviving Species Biomass')
set(gca,'yscale','log')
ylabel('Fraction')

%
%
%
% How different communities properties change with cheating breadth


% Get unique values for cheating receptor counts and pure-cheater ratios
actualCR=simulation_data_matrix(:, CRcol);
actualcheater=simulation_data_matrix(:, cheaternumcol);
CBint_inrawdata=round(actualCR/N_spe);
Cheaterint_inrawdata=round(actualcheater);

purecheaternum_list = 0:5;  % Sorted unique pure-cheater ratios
cheatingbreadth_list = 0:30;  % Cheating breadth as fraction of species

% Initialize result matrices
num_ratios = length(purecheaternum_list);
num_breadths = length(cheatingbreadth_list);
ExtinctedFrac_mean_mat = zeros(num_ratios, num_breadths);  % Mean extinction fraction
NonExtinction_CommunitySize_meanvar_mat = zeros(num_ratios, num_breadths, 2);  % Mean and SEM for community size (non-extinct)
NonExtinction_CommunityBiomass_meanvar_mat = zeros(num_ratios, num_breadths, 2);  % Mean and SEM for biomass (non-extinct)
NonsteadyFrac_mean_mat = zeros(num_ratios, num_breadths);  % Fraction of non-steady states (non-extinct)
AlphaDiversityIndex_meanvar_mat = zeros(num_ratios, num_breadths, 2);  % Mean and SEM for alpha diversity

% Loop over pure-cheater ratios and cheating breadths
for k = 1:num_ratios
    for i = 1:num_breadths
        % Find rows matching current ratio and receptor count
        row_indices = find(CBint_inrawdata == cheatingbreadth_list(i) & ...
            Cheaterint_inrawdata == purecheaternum_list(k));

        % Extract relevant data for these rows
        biomasses = biomasses_mat(row_indices, :);
        totbiomass = biomass_tot(row_indices);
        survivingspeciescount = surviving_spe(row_indices);

        % Identify non-extinct and non-steady communities
        nonextinctionlist = find(survivingspeciescount > 0);
        nonsteadylist = find(notsteady(row_indices) == 1);

        alphadiversityvalues = Shannon_Alpha_Index(biomasses, extinction_thresh);



        % Compute extinction fraction
        ExtinctedFrac_mean_mat(k, i) = 1 - length(nonextinctionlist) / length(row_indices);

        % Community size for non-extinct communities
        datalist = survivingspeciescount(nonextinctionlist);
        NonExtinction_CommunitySize_meanvar_mat(k, i, 1) = mean(datalist);
        NonExtinction_CommunitySize_meanvar_mat(k, i, 2) = std(datalist) / mean(datalist);

        % Biomass for non-extinct communities
        datalist = totbiomass(nonextinctionlist);
        NonExtinction_CommunityBiomass_meanvar_mat(k, i, 1) = mean(datalist);
        NonExtinction_CommunityBiomass_meanvar_mat(k, i, 2) = std(datalist) / mean(datalist);

        % Alpha diversity (across all, including extinct?)
        datalist = alphadiversityvalues;
        AlphaDiversityIndex_meanvar_mat(k, i, 1) = mean(datalist);
        AlphaDiversityIndex_meanvar_mat(k, i, 2) = std(datalist) / mean(datalist);

        % Non-steady fraction among non-extinct communities
        NonsteadyFrac_mean_mat(k, i) = length(intersect(nonsteadylist, nonextinctionlist)) / ...
            length(nonextinctionlist);
    end
end



xtick_values = 0:5:50;
[~, xtick_indices] = ismember(xtick_values, cheatingbreadth_list);
xtick_indices = xtick_indices(xtick_indices > 0);
xtick_labels = num2cell(cheatingbreadth_list(xtick_indices));

ytick_indices = 1:length(purecheaternum_list);
ytick_labels = num2cell(100*purecheaternum_list/N_spe);

plot_data = {
    100*ExtinctedFrac_mean_mat, 'Extinction Fraq %';
    NonExtinction_CommunitySize_meanvar_mat, 'Community Size (non extinction)';
    NonExtinction_CommunityBiomass_meanvar_mat, 'Community Biomass(non extinction)';
    100*NonsteadyFrac_mean_mat, 'NonSteady Frac %(non extinction)';
    AlphaDiversityIndex_meanvar_mat, 'Alpha Diversity'
    };

plotbyimagesc = 1;
if plotbyimagesc
    for i = 1:size(plot_data, 1)
        figure;
        valuemat=plot_data{i,1};

        imagesc(valuemat(:,:,1));
        title(plot_data{i,2});
        colorbar;
        xlabel('Cheating Breadth');
        ylabel('Pure-Cheater %');

        set(gca, 'XTick', xtick_indices);
        set(gca, 'XTickLabel', xtick_labels);
        set(gca, 'YTick', ytick_indices);
        set(gca, 'YTickLabel', ytick_labels);
    end
end


plotbyline=1;
if plotbyline

    linecolors = [
        0.0  0.45 0.74; % Deep Blue (low)
        0.0  0.75 0.75; % Cyan
        0.0  0.9  0.0;  % Green
        1.0  0.75 0.0;  % Yellow
        1.0  0.5  0.0;  % Orange
        1.0  0.0  0.0   % Red (high)
        ]; % for different cheater ratio
    fmt_str = sprintf('p_{cheater}=%0.0f%%,', 100 * purecheaternum_list/N_spe);fmt_str = fmt_str(1:end-1);  % 移除最后一个逗号
    legendstr = strsplit(fmt_str, ',');


    for i = 1:size(plot_data, 1)
        valuemat=plot_data{i,1};
        figure;
        hold on;
        axis square

        bh=[];
        for k=1:length(purecheaternum_list)
            mean_y=valuemat(k,:,1);
            lineColor=linecolors(k,:);

            ms=30;lw=1;
            lastp=30;
            bh(k)=scatter(cheatingbreadth_list(1:lastp), mean_y(1:lastp), ms, lineColor, 'LineWidth', lw);
            plot(cheatingbreadth_list(1:lastp), mean_y(1:lastp), 'Color', lineColor)
            hold on;
            % Define the upper and lower bounds for the shaded standard deviation
            if size(valuemat,3)==2
                std_y=valuemat(k,:,2);
                upperBound = mean_y + std_y;
                lowerBound = mean_y - std_y;
                % Create the shaded standard deviation region
                % Use fill to create a shaded area with transparency
                %%% fill([cheatingbreadth_list', fliplr(cheatingbreadth_list')], [upperBound, fliplr(lowerBound)], lineColor, 'FaceAlpha', 0.3, 'EdgeColor', 'none'); % Shading with 30% opacity
            end
        end

        ylabel(plot_data{i,2});
        xlabel('Cheating Breadth');
        if i==1
            legend(bh,legendstr)
            legend('Location', 'northwest');
            axis square
        end
    end
end

%
%
%
% Violin plots showing non-extinct community size distributions for steady (yellow) versus dynamic (green) outcomes , under different cheating breadth. Dynamic communities consistently support higher biodiversity
%

maxspecies=20;

figure;
toplotlist=1:6;
bh=[];
for k=1:length(toplotlist)
    figure; hold on; axis square;
    for i = 1:length(cheatingbreadth_list)

        touse_x=cheatingbreadth_list(i);

        % Find the indices where the 5th column matches the current unique value.

        row_indices = find(CBint_inrawdata == cheatingbreadth_list(i)  & Cheaterint_inrawdata == purecheaternum_list(toplotlist(k)));
        survivingspeciescount=surviving_spe(row_indices,:);
        notexintct_thispart=(survivingspeciescount>0);
        notsteady_thispart=(notsteady(row_indices)==1);

        exclusion_stableco_dynamicalco_colors = [0, 0, 0;          % 纯黑
    1, 0.843, 0;     % 金黄色
    0, 0.6275    0.9137]; % 蓝
        
        % oscillation and non-oscillation
        for o=0:1
            if o==0 %steady state
                surviving_count_thiscase=survivingspeciescount(notsteady_thispart==0 & notexintct_thispart);
                colorcode=exclusion_stableco_dynamicalco_colors(2,:);
            else % oscillatory
                surviving_count_thiscase=survivingspeciescount(notsteady_thispart==1 & notexintct_thispart);
                colorcode=exclusion_stableco_dynamicalco_colors(3,:);
            end

            if ~isempty(surviving_count_thiscase)
                data = min(maxspecies,surviving_count_thiscase(:));
                mean_val = mean(data);
                halfwidth = 0.46;  % Half-width for the violin and mean line; adjust based on x-range if needed
                if length(data) > 1
                    tabfreq=tabulate(data);
                    f=tabfreq(:,3);
                    xi=tabfreq(:,1);
                    if ~isempty(f)
                        maxf = max(f);
                        scale = halfwidth / maxf;
                        x_left = touse_x - scale * f;
                        x_right = touse_x + scale * f;
                        x_poly = [x_left; flipud(x_right)];
                        y_poly = [xi;flipud(xi)];
                        fill(x_poly, y_poly, colorcode, 'FaceAlpha', 0.7, 'EdgeColor', 'none');
                    end
                end
                if ~isnan(mean_val) && isfinite(mean_val)
                    bh(1+o)=plot([touse_x - halfwidth, touse_x + halfwidth], [mean_val, mean_val], ...
                        'Color', colorcode, 'LineWidth', 1);
                end
            end %if ~isempty(surviving_count_thiscase)

        end
    end
    ylim([0.9 maxspecies]);
    title(legendstr{toplotlist(k)})
    xlim([0 30])

    xlabel('Cheating Breadth');
    ylabel('Community Size (non extinction)')

    legend(bh,{'Steady','Oscillatory'})
    legend('Location', 'northwest');
end