close all;
clear; clc;

% Get all results files
files = dir("../../../results/*.out");
tw = load("../tracy-widom-approx/tracy_widom_normal.txt");
x = tw(:,1);
tw1 = tw(:,2);
nrm = tw(:,3);

N_list = []

%% Make a histogram for every results file

for f = 1:length(files)

    % Get filename
    filename = files(f).name;

    % Extract N and L
    tokens = regexp(filename,'eigs_(\d+)_(\d+)\.out', 'tokens');
    N = str2double(tokens{1}{1});
    N_list = [N_list, N];
    L = str2double(tokens{1}{2});

    % Import actual data
    data = load(strcat("../../../results/",filename));
    % And clip to just the eigenvalues
    data = data(:,2);

    % Remove outliers
    data_clean = data(data > sqrt(N));
    
    % Scale and such
    s_data = (data_clean - 2*sqrt(N))*N^(1/6);

    % Plot them
    figure(f);
    title([num2str(N)," matrix size, ",num2str(L)," data points"])
    hold on
    [counts, bins] = hist(s_data(s_data > -6), linspace(-6,4,40));
    bin_width = bins(2)-bins(1);
    pdf = counts/(sum(counts) * bin_width);
    bar(bins, pdf);
    plot(x,tw1);
    A = [bins', pdf']
    save("-ascii",["N",num2str(N),"_L",num2str(L),".txt"],"A")
end

%% Plot largest sample size for each matrix size N

skews = [];
kurts = [];
i = 0;
% Plot totals
for n = unique(N_list)
    i = i + 1;

    % Plot all histogram values for this N
    n_files = dir(strcat("../../../results/eigs_",num2str(n),"_*.out"));
    n_vals = [];
    SSE = [];
    
    % Get cumulative values
    for f = 1:length(n_files)
        filename = n_files(f).name;
        data = load(strcat("../../../results/",filename));
        data = data(:,2);
        data_clean = data(data > sqrt(n));
        s_data = (data_clean -2*sqrt(n)) * n^(1/6);
        n_vals = [n_vals, s_data'];
    end

    % Plot a histogram showing the Gaussian, TW1, and sampled distribution
    figure(n)
    title(num2str(length(n_vals)))
    hold on

    % Generate the histogram
    [counts, bins] = hist(n_vals(n_vals > -8), linspace(-8,6,70));
    bin_width = bins(2)-bins(1);
    pdf = counts/(sum(counts) * bin_width);
    
    % And plot it
    bar(bins,pdf)
    plot(bins, pdf);
    plot(x,nrm,'--')
    plot(x,tw1);

    % Save the histogram to file
    A = [bins', pdf']
    save("-ascii",[num2str(n),"_max.txt"], "A")

    % Calculate some statistics 
    [skews(i),kurts(i)]=pdf_stats(bins,pdf);
end

figure(n+1)
plot(unique(N_list),skews)
hold on
plot(unique(N_list),kurts)
set(gca, "xtick", unique(N_list))
A = [unique(N_list)', skews', kurts']
save -ascii props.txt A