% Script for fitting power law models to the avalanche data generated by
% simulation
% 1. use MLE estimate to get slope parameters for discrete distributions
% 2. use KS statistic between empirical and model distribution to determine
% x_min. 
% 3. goodness of fit: using the best-fit parameters, generate simulated 
% data (that is truly power-law) and measure the KS statistic relative to 
% the true model. Compare this to the KS statistic calculated in (2) for 
% the empirical data relative to the true model. 

% ss_str.run_f20_e3010_timesweep = false;       % pull data from 20-field, time sweep
% ss_str.run_f5_e1204_timesweep = false;       % pull data from 5-field, eps = -12, eta = 4 time sweep
% ss_str.run_f1_e1204_timesweep = false;       % pull data from 1-field, eps = -12, eta = 4 time sweep
% 
% ss_str.run_f1finesweep = false;       % pull data from 1-field, fine sweep of eta, eps
% ss_str.run_f5finesweep = false;       % pull data from 5-field, fine sweep
% ss_str.run_f20finesweep = false;       % pull data from 20-field, fine sweep
% ss_str.run_morrell1field = false;       % analyze the envfield_fixedJ sims from Mia

% ss_str.run_f5ultrafinesweep = true;       % pull data from 5-field, ultrafine sweep


% choose one field to set true (include replicate for the ultrafine sweep)
% ss_str = selectSimulation('run_morrell1field');
% ss_str = selectSimulation('run_f1smallfinesweep', 'J');
ss_str = selectSimulation('run_f1smallultrafinesweep', 'F');

% run this script to get file names and structures
setLoadingFunctionNames
% %% LOAD THE HURWITZ ZETA LOOKUP TABLE
% % 
% load('avalanches_matlab_code/analysis_results/hurwitzLookup_20220118', 'spl_hurwitzZeta', ...
%     'spl_dhurwitzZeta')
% %%
% global glu_hurwitzZeta
% glu_hurwitzZeta = @(a, x) spl_hurwitzZeta(a, x);
% global glu_dhurwitzZeta
% glu_dhurwitzZeta = @(a, x) spl_dhurwitzZeta(a, x);
% 


%% Fit power laws, determine lower cutoff, and determine goodness of fit. 

% for 1-field plots, may need smaller cutoffs - avalanches are smaller.
% (1e3, 800?) 
% for 5-field, try 5e3 and 1e3
% size_max = 5e3; % 1e4 : number of neurons is 1e3, track length 1e5; max ava is 1e8
% dur_max = 1e3; % 1e3 : track length duration is 10,000; thus, set max duration to be 10% of that
% edit on 2/15/22: set the upper cutoff to the min cutoff + 2 decades? 

% maximum number of observations
load([save_dir ava_results_file], 'all_tau_pfit', 'all_alpha_pfit', ...
    'all_gamma_pfit', 'all_gamma_pred', 'avalanche_count', 'x_min_vals');

% % save exponents and bootstrapped errorbars
% all_tau_pfit = cell(length(y_var_list), length(x_var_list));
% all_alpha_pfit = cell(length(y_var_list), length(x_var_list));
% all_gamma_pfit = cell(length(y_var_list), length(x_var_list));
% 
% % for results scanning upper and lower cutoff values
% all_tau_pfit_rng = cell(length(y_var_list), length(x_var_list));
% all_alpha_pfit_rng = cell(length(y_var_list), length(x_var_list));
% 
% all_gamma_pred = cell(length(y_var_list), length(x_var_list));

% also save the number of avalanches in each simulation
% avalanche_count = zeros(length(y_var_list), length(x_var_list));
clear review_notes
review_notes(size(all_alpha_pfit, 1), size(all_alpha_pfit, 2)) = struct('size_scaling', [], 'size_scaling_reason', [], ...
    'duration_scaling', [], 'duration_scaling_reason', [], ...
    'sizedur_scaling', [], 'sizedur_scaling_reason', [], ... 
    'crackling', [], 'crackling_reason', []);
%%
review_interactive = false;
for ii_y = 1:length(y_var_list)
    for ii_x = 1:length(x_var_list)
    %%
        try
            x_ava = load(ava_name_fn(ii_y, ii_x));
            x_fr = load(fr_name_fn(ii_y, ii_x));
            % require at least one avalanche of size 500
            if max(x_ava.sizes) > 5e2
                proceed = true;
            else
                proceed = false;
            end
        catch
            proceed = false;
            disp('file missing')
        end
        %%
        if proceed
            % proceed only if we have at least 10000 avalanches
            if length(x_ava.sizes) > 10000
                
                
                
                [n_cts, size_bins, dur_bins] = calculateJointHistogram(x_ava.sizes, x_ava.durations);
                [ave_sizes, se_sizes, dur_bins, num_obs] = calculateAveSizeByDuration(x_ava.sizes, x_ava.durations);
           
                
                %% load everything
                % Edit plotting: if this failed to run, just plot teh
                % size and duration distributions
                avalanche_count(ii_y, ii_x) = length(x_ava.sizes);
                pfit_S = all_tau_pfit{ii_y, ii_x};
                pfit_D = all_alpha_pfit{ii_y, ii_x};
                % load the selected x_min for Duration and Size
                % distributions
                x_min_D = pfit_D.x_min;
                x_min_S = pfit_S.x_min;
                % load surrogate stats
                ks_stats_Surr_S = pfit_S.ks_surrogate;
                surr_a_S = pfit_S.a_est_surr;
                ks_stats_Surr_D = pfit_D.ks_surrogate;
                surr_a_D = pfit_D.a_est_surr;
                %%
                
                g_pred_all = all_gamma_pred{ii_y, ii_x}.mhat;
                g_predSE = all_gamma_pred{ii_y, ii_x}.mSE;
                pfit_G = all_gamma_pfit{ii_y, ii_x};
%                 all_tau_pfit{ii_reps, ii_time} = pfit_S;
%                 all_alpha_pfit{ii_reps, ii_time} = pfit_D;
%                 % all_tau_pfit_rng{ii_reps, ii_time} = pfit_rng_S;
%                 % all_alpha_pfit_rng{ii_reps, ii_time} = pfit_rng_D;
%                 all_gamma_pred{ii_reps, ii_time}.mhat = g_pred_all;
%                 all_gamma_pred{ii_reps, ii_time}.mSE = g_predSE;
%                 all_gamma_pfit{ii_reps, ii_time} = pfit_G;
%                 
                close all;
                %% make plots
                % TO-DO on these plots: plot the pdf, not the counts. Then,
                % get the axes cleaned up - show the MLE line, the cutoff,
                % with proper axis limits. Get the error bars for the PDFs.
                % 
                makeMyFigure(30, 25);
                
                
                h_sz1 = subplot(3, 3, 1);
                freq_cts_size = sum(n_cts, 2);
                log_size_bins = unique(round(logspace(0, 5, 51)));  % use unique(round()) to avoid having 1, 1.25, 1.58, etc. 
                size_pdf = histcounts(x_ava.sizes, bins2edges(log_size_bins), 'Normalization', 'pdf');
                size_cts = histcounts(x_ava.sizes, bins2edges(log_size_bins), 'Normalization', 'count');
                rel_err = 1./sqrt(size_cts-1);
                
                size_pdf_SE = rel_err.*size_pdf; % get counting statistics error bars on PDF estimates
%                 plot(log10(size_bins), log10(freq_cts_size),'.')
%                 plot(log10(log_size_bins), log10(size_pdf), '.')
                errorbar(log10(log_size_bins), log10(size_pdf), ...
                    log10(size_pdf) - log10(size_pdf-size_pdf_SE), ...  % lower errorbar
                    log10(size_pdf + size_pdf_SE) - log10(size_pdf), '.', ...  % upper errorbar
                    'capsize', 0, 'color', 0.2*[1 1 1])         
                
                hold on
                s_cutoff = [log10(pfit_S.x_min) max(log10(size_bins))];
                use_bins = log10(log_size_bins) >= (s_cutoff(1)) & isfinite(log10(size_pdf));
                b0 = mean((pfit_S.a_hat*log10(log_size_bins(use_bins)) + log10(size_pdf(use_bins))));
%                 b0 = log10(sum(n_cts(1, :)));
                plot(s_cutoff,b0 - pfit_S.a_hat*s_cutoff, '-', 'linewidth', 1);
                
                s_y_lims = [min(b0 - pfit_S.a_hat*s_cutoff) max(log10(size_pdf)) + .2];
                plot(log10(pfit_S.x_min)*[ 1 1], s_y_lims, 'k')
                % xlim([0 dur_bins(end-1)])
                n_tail_S = sum(freq_cts_size(size_bins >= pfit_S.x_min));
%                 ylim([-1 b0])
                ylabel('size pdf (log_{10})')
                xlabel('size bin')
                title({['\tau = ' num2str(pfit_S.a_hat, '%1.3f') ' ' num2str(pfit_S.se_a_hat, ' (%1.3f SE)')]; ...
                    ['s_{min} = ' num2str(pfit_S.x_min) ', n_{tail} = ' num2str(n_tail_S)]})
                
                set(gca, 'color', 'none')
                axis square
                
                h_dr1 = subplot(3, 3, 2);
                
                log_dur_bins = unique(round(logspace(0, 5, 51)));  % use unique(round()) to avoid having 1, 1.25, 1.58, etc. 
                dur_pdf = histcounts(x_ava.durations, bins2edges(log_dur_bins), 'Normalization', 'pdf');
                dur_cts = histcounts(x_ava.durations, bins2edges(log_dur_bins), 'Normalization', 'count');
                rel_err_dur = 1./sqrt(dur_cts-1);
                
                dur_pdf_SE = rel_err.*dur_pdf; % get counting statistics error bars on PDF estimates
%                 plot(log10(size_bins), log10(freq_cts_size),'.')
%                 plot(log10(log_size_bins), log10(size_pdf), '.')
                errorbar(log10(log_dur_bins), log10(dur_pdf), ...
                    log10(dur_pdf) - log10(dur_pdf-dur_pdf_SE), ...  % lower errorbar
                    log10(dur_pdf + dur_pdf_SE) - log10(dur_pdf), '.', ...  % upper errorbar
                    'capsize', 0, 'color', 0.2*[1 1 1])                
               
                
%                 plot(log10(dur_bins), log10(sum(n_cts, 1)),'.')
                hold on
                d_cutoff = [log10(pfit_D.x_min) max(log10(dur_bins))];
                use_bins = log10(log_dur_bins) >= (d_cutoff(1)) & isfinite(log10(dur_pdf));
                b0 = mean((pfit_D.a_hat*log10(log_dur_bins(use_bins)) + log10(dur_pdf(use_bins))));
%              
                
%                 b0 = log10(sum(n_cts(:, 1)));
                plot(d_cutoff, b0 - pfit_D.a_hat*d_cutoff, '-', 'linewidth', 1);
                
                d_y_lims = [min(b0 - pfit_D.a_hat*d_cutoff) max(log10(dur_pdf)) + .2];
                plot(log10(pfit_D.x_min)*[ 1 1], d_y_lims, 'k')
                
%                 plot(log10(pfit_D.x_min)*[ 1 1], [0 5], 'k')
                
                hold on
                % title(['\alpha = ' num2str(p_fit_A.mhat, '%1.3f') ' ' num2str(p_fit_A.mSE, '(%1.3f SE)')])
                title({['\alpha = ' num2str(pfit_D.a_hat, '%1.3f') ' ' num2str(pfit_D.se_a_hat, ' (%1.3f SE)')]; ...
                    ['d_{min} = ' num2str(pfit_D.x_min)]})
                
                ylabel('log_{10} number of observations')
                xlabel('duration bin')
                set(gca, 'color', 'none')
                axis square
                
                % start at the duration cutoff
                use_G_fit = sum(cellfun(@(x) x.x_cutoff(1) <= d_cutoff(1), pfit_G));
                h_sd1 = subplot(3, 3, 3);
                hold on
                errorbar(log10(dur_bins), log10(ave_sizes),log10(ave_sizes) - log10(ave_sizes - se_sizes), ...
                    log10(ave_sizes + se_sizes)- log10(ave_sizes), '.', ...
                    'capsize', 0, 'color', 0.2*[1 1 1])
                plot([0 d_cutoff(2)], pfit_G{use_G_fit}.bhat + pfit_G{use_G_fit}.mhat*[0 d_cutoff(2)], '-', 'linewidth', 1);
                plot(d_cutoff, pfit_G{use_G_fit}.bhat + pfit_G{use_G_fit}.mhat*d_cutoff, '-', 'linewidth', 1);
                
                % plot(log10(dur_bins(num_obs < 30)), log10(ave_sizes(num_obs < 30)), '.')
                % plot(log10(dur_bins(num_obs < 2)), log10(ave_sizes(num_obs < 2)), '.')
                title({['\gamma_{fit} = ' num2str(pfit_G{use_G_fit}.mhat, '%1.3f') ' ' num2str(pfit_G{use_G_fit}.mSE, '(%1.3f SE)')]; ...
                    ['\gamma_{pred} = ' num2str(g_pred_all(x_min_vals == x_min_S,x_min_vals == x_min_D), '%1.3f')...
                    ' ' num2str(g_predSE(x_min_vals == x_min_S,x_min_vals == x_min_D), '(%1.3f SE)')]});
                
                xlabel('duration')
                ylabel('average size')
                set(gca, 'color', 'none')
                legend({'<s> by dur', 'extrap to 0', 'fitted line'}, 'Location', 'best')
                axis square
                xlim([0 3])
                ylim([0 3])
                try
                    suptitle([par_title_str x_lab_sum ' = ' num2str(x_var_list(ii_x)) ', ' y_lab_sum num2str(y_var_list{ii_y})])
                catch
                    suptitle([par_title_str x_lab_sum ' = ' num2str(x_var_list(ii_x)) ', ' y_lab_sum num2str(y_var_list(ii_y))])
                    
                end
                subplot(3, 3, 4)
                errorbar(log10(pfit_S.x_min_vals), pfit_S.a_vals, pfit_S.se_a_vals)
                ylabel('\tau_{MLE}')
                xlabel('s_{min}')
                yyaxis right
                plot(log10(pfit_S.x_min_vals), pfit_S.ks_stats)
                frac_tail_S = cellfun(@(x) mean(x_ava.sizes > x), num2cell(x_min_vals));
%                 ks_statS_raw = pfit_S.ks_stats./frac_tail_S;
                % hold on
                % plot(log10(pfit_S.x_min_vals),ks_statS_raw)
                ylabel('KS statistic')
                set(gca, 'color', 'none')
                
                subplot(3, 3, 5)
                errorbar(log10(pfit_D.x_min_vals), pfit_D.a_vals, pfit_D.se_a_vals)
                ylabel('\alpha_{MLE}')
                xlabel('d_{min}')
                yyaxis right
                plot(log10(pfit_D.x_min_vals), pfit_D.ks_stats)
                frac_tail_D = cellfun(@(x) mean(x_ava.durations > x), num2cell(x_min_vals));
%                 ks_statD_raw = pfit_D.ks_stats./frac_tail_D;
                % hold on
                % plot(log10(pfit_D.x_min_vals),ks_statD_raw)
                ylabel('KS statistic')
                set(gca, 'color', 'none')
                
                subplot(3, 3, 6)
                d_low_cut = cellfun(@(x) x.x_cutoff(1), pfit_G);
                m_hat_vs_x0 = cellfun(@(x) x.mhat, pfit_G);
                se_m_hat_vs_x0 = cellfun(@(x) x.mSE, pfit_G);
                x0 = cellfun(@(x) x.x_cutoff(1), pfit_G);
                hold on
                errorbar(x0, m_hat_vs_x0, se_m_hat_vs_x0, 'color', 0.2*[1 1 1])
                
                % also plot line for gamma predicted
%                 ['\gamma_{pred} = ' num2str(g_pred_all(x_min_vals == x_min_D,x_min_vals == x_min_S), '%1.3f')...
%                     ' ' num2str(g_predSE(x_min_vals == x_min_D,x_min_vals == x_min_S), '(%1.3f SE)')]
                % prior to 5/5/22, these entries were transposed (S adn D
                % switched) which selected the wrong value. 
                g_pred_ave = g_pred_all(x_min_vals == x_min_S,x_min_vals == x_min_D);
                g_pred_SE = g_predSE(x_min_vals == x_min_S,x_min_vals == x_min_D);
                plot(x0([1 end]), g_pred_ave*[ 1 1],'k', 'linewidth', 1)
                plot(x0([1 end]), g_pred_ave*[ 1 1] + g_pred_SE, 'k', 'linewidth', 0.5)
                plot(x0([1 end]), g_pred_ave*[ 1 1] - g_pred_SE, 'k', 'linewidth', 0.5)
                xlabel('minimum duration')
                ylabel('\gamma_{fit} for 1 decade')
                set(gca, 'color', 'none')
                
                subplot(3,3,7)
                cla;
                ks_bins_S = logspace(log10(min(ks_stats_Surr_S)), log10(max(ks_stats_Surr_S)), 21);
                hold on
                histogram(log10(ks_stats_Surr_S), log10(ks_bins_S), 'EdgeColor', 'none')
                y_lims = ylim;
                ks_stat_hat = pfit_S.ks_stats(pfit_S.x_min_vals == pfit_S.x_min);
                plot(log10(ks_stat_hat)*[1 1], y_lims)
                xlabel('KS statistic')
                ylabel(['pdf (' num2str(length(ks_stats_Surr_S)) ' surrogates)'])
                title(['power law for s > ' num2str(pfit_S.x_min) ...
                    ', \tau = ' num2str(pfit_S.a_hat)])
                set(gca, 'color', 'none')
                
                subplot(3,3,8)
                cla;
                ks_bins_D = logspace(log10(min(ks_stats_Surr_D)), log10(max(ks_stats_Surr_D)), 21);
                
                hold on
                histogram(log10(ks_stats_Surr_D), log10(ks_bins_D), 'EdgeColor', 'none')
                
                % histogram(ks_stats_Surr_D,  'EdgeColor', 'none')
                y_lims = ylim;
                ks_stat_hat = pfit_D.ks_stats(pfit_D.x_min_vals == pfit_D.x_min);
                plot(log10(ks_stat_hat)*[1 1], y_lims)
                xlabel('KS statistic')
                ylabel(['pdf (' num2str(length(ks_stats_Surr_D)) ' surrogates)'])
                title(['power law for s > ' num2str(pfit_D.x_min) ...
                    ', \alpha = ' num2str(pfit_D.a_hat)])
                set(gca, 'color', 'none')
                
                subplot(3, 3, 9)
                g_xy = 1:length(pfit_S.x_min_vals);
                imagesc(g_xy, g_xy, g_pred_all)
                hold on
                sm_g_pred_all = imfilter(g_pred_all, fspecial('gaussian', [5 5], 2), 'replicate');
                
                ave_fit_g = mean(m_hat_vs_x0(d_low_cut > log10(x_min_D) & d_low_cut < 2));
                se_fit_g = std(m_hat_vs_x0(d_low_cut > log10(x_min_D) & d_low_cut < 2));
                if isfinite(se_fit_g)
                    contour(g_xy, g_xy, sm_g_pred_all, ...
                        ave_fit_g + se_fit_g*[-2 2], 'color', [1 1 1])
                end
                plot(find(x_min_S == x_min_vals), find(x_min_D == x_min_vals), 'r*')
                colorbar
                set(gca, 'xtick', g_xy(1:12:end), 'XTickLabel', x_min_vals(1:12:end))
                set(gca, 'ytick', g_xy(1:8:end), 'YTickLabel', x_min_vals(1:8:end))
                set(gca, 'ydir', 'normal')
                axis square
                xlabel('d_{min}')
                ylabel('s_{min}')
                title('\gamma_{pred}')
                %% interactive review
%                 review_notes(size(all_alpha_pfit)) = struct('size_scaling', [], 'size_scaling_reason', [], ...
%     'duration_scaling', [], 'duration_scaling_reason', [], ...
%     'sizedur_scaling', [], 'sizedur_scaling_reason', [], ... 
%     'crackling', []);
                if review_interactive
                size_sc = input('Size distribution review: is there scaling? [1 is yes, 0.5 for maybe, 0 for no]');
                if size_sc ~= 1
                    size_sc_q = input(['What is wrong with the scaling? Enter all that apply: \n' ...
                        '   (A) KS statistic is outside surrogate distribution \n' ...
                        '   (B) KS vs cutoff has no plateau \n' ...
                        '   (C) Lower cutoff is high (>20) \n' ...
                        '   (D) other \n']);
                    
                    h_sz1.Title.String = vertcat(h_sz1.Title.String, {['NO SCALING: error ' size_sc_q]});
                    
                else
                    size_sc_q = '';
                    h_sz1.Title.String = vertcat(h_sz1.Title.String, {'SCALING'});
                end
                
                dur_sc = input('Duration distribution review: is there scaling? [1 is yes, 0 is no] \n');
                if dur_sc ~= 1
                    dur_sc_q = input(['What is wrong with the scaling? Enter all that apply: \n' ...
                        '   (A) KS statistic is outside surrogate distribution \n' ...
                        '   (B) KS vs cutoff has no plateau \n' ...
                        '   (C) Lower cutoff is high (>20) \n' ...
                        '   (D) other \n']);
                    
                    h_dr1.Title.String = vertcat(h_dr1.Title.String, {['NO SCALING: error ' dur_sc_q]});
                else
                    h_dr1.Title.String = vertcat(h_dr1.Title.String, {'SCALING'});
                    dur_sc_q = '';
                end
                
                review_notes(ii_y, ii_x).size_scaling = size_sc;
                review_notes(ii_y, ii_x).size_scaling_reason = size_sc_q;
                review_notes(ii_y, ii_x).duration_scaling = dur_sc;
                review_notes(ii_y, ii_x).duration_scaling_reason = dur_sc_q;
                %% now look at gamma fits
                
                sizedur_sc = input('Size with duration review: is there scaling? [1 is yes, 0 is no]');
                if ~sizedur_sc
                    size_dur_sc_q = input(['What is wrong with the scaling? Enter all that apply: \n' ...
                        '   (A) Gamma is not constant over at least 2 decades \n' ...
                        '   (B) Other \n']);
                    
                    h_sd1.Title.String = vertcat(h_sd1.Title.String, {['NO SCALING: error ' size_dur_sc_q]});
                    min_duration_sd = '';
                else
                    size_dur_sc_q = '';
                    min_duration_sd = input('What is the minimum cutoff (0, 0.25, 0.5, 0.75, 1)? \n');
                    h_sd1.Title.String = vertcat(h_sd1.Title.String, {['SCALING, start at ' num2str(min_duration_sd)]});
                end
                
                review_notes(ii_y, ii_x).sizedur_scaling = sizedur_sc;
                review_notes(ii_y, ii_x).sizedur_scaling_reason = size_dur_sc_q;
                review_notes(ii_y, ii_x).min_duration_sd = min_duration_sd;
                
                has_crackles = input('Do exponents follow crackling relationship? If any scaling test failed, enter 0. \n');
                review_notes(ii_y, ii_x).crackling = has_crackles;
                % 
                end
                %% print to file
                print(gcf, '-dpdf', ['avalanches_matlab_code/plots/' dir_name_str '/review_Clausetmincutoff_' param_str(ii_y, ii_x) '.pdf'])
                %         end
                


            end
        end
        close all
    end

end

% save everything once run is complete
save([save_dir 'ava_decade_review_' datestr(now, 'yyyymmdd')], 'review_notes')


%% make summary plots
% 
% % handle the non-grid-sampled data
% 
n_samples_surr = length(all_alpha_pfit{end, 1}.ks_surrogate);

alpha_vals_all = nan(size(all_alpha_pfit));
alpha_sevals_all = nan(size(all_alpha_pfit));
alpha_ks_pval = nan(size(all_alpha_pfit));
alpha_ks_zval = nan(size(all_alpha_pfit));
alpha_lowercut = nan(size(all_alpha_pfit));

tau_vals_all = nan(size(all_alpha_pfit));
tau_sevals_all = nan(size(all_alpha_pfit));
tau_ks_pval = nan(size(all_alpha_pfit));
tau_ks_zval = nan(size(all_alpha_pfit));
tau_lowercut = nan(size(all_alpha_pfit));


gamma_vals_all = nan(size(all_alpha_pfit));
gamma_sevals_all = nan(size(all_alpha_pfit));
% rsq_vals_G = nan(size(all_alpha_pfit));
% rsq_vals_A = nan(size(all_alpha_pfit));
% rsq_vals_T = nan(size(all_alpha_pfit));

has_results = cellfun(@(x) isfield(x, 'ks_surrogate'), all_alpha_pfit); %cellfun(@(x) ~isempty(x), all_alpha_pfit);

alpha_vals_all(has_results) = cellfun(@(x) x.a_hat, all_alpha_pfit(has_results));
alpha_sevals_all(has_results) = cellfun(@(x) x.se_a_hat, all_alpha_pfit(has_results));
tau_vals_all(has_results) = cellfun(@(x) x.a_hat, all_tau_pfit(has_results));
tau_sevals_all(has_results) = cellfun(@(x) x.se_a_hat, all_tau_pfit(has_results));

alpha_lowercut(has_results) = cellfun(@(x) x.x_min, all_alpha_pfit(has_results));
alpha_ks_pval(has_results) = cellfun(@(x) mean(x.ks_surrogate >= min(x.ks_stats)), ... 
    all_alpha_pfit(has_results));
alpha_ks_zval(has_results) = cellfun(@(x) abs(min(x.ks_stats)-median(x.ks_surrogate))/std((x.ks_surrogate)), ... 
    all_alpha_pfit(has_results));

tau_lowercut(has_results) = cellfun(@(x) x.x_min, all_tau_pfit(has_results));
tau_ks_pval(has_results) = cellfun(@(x) mean(x.ks_surrogate >= min(x.ks_stats)), ... 
    all_tau_pfit(has_results));
tau_ks_zval(has_results) = cellfun(@(x) abs(min(x.ks_stats)-median(x.ks_surrogate))/std((x.ks_surrogate)), ... 
    all_tau_pfit(has_results));
    
gamma_preds_all = (alpha_vals_all - 1)./(tau_vals_all - 1);

% compute gamma fits, average and SE starting with lower cutoff from
% duration distribution fit. 
ind_has_results = find(has_results);
gamma_x_cuts = cellfun(@(x) x.x_cutoff(1), all_gamma_pfit{ind_has_results(1)});
gamma_start_inds = cellfun(@(x) ceil(sum(gamma_x_cuts <= log10(x))), ... 
    num2cell(alpha_lowercut), 'UniformOutput', false);
gamma_mhat_all_vals = cell(size(has_results));
gamma_mhat_all_vals(has_results) = cellfun(@(x) cellfun(@(y) y.mhat, x), ...
    all_gamma_pfit(has_results), 'UniformOutput', false);

gamma_ave_mhat = zeros(size(has_results));
gamma_ave_mhat(has_results) = cellfun(@(x, start_ind) mean(x(start_ind:end)), ... 
    gamma_mhat_all_vals(has_results), gamma_start_inds(has_results));

gamma_se_mhat = zeros(size(has_results));
gamma_se_mhat(has_results) = cellfun(@(x, start_ind) std(x(start_ind:end)), ... 
    gamma_mhat_all_vals(has_results), gamma_start_inds(has_results));

[gamma_r_mhat_v, gamma_pr_mhat_v] = cellfun(@(x, start_ind) corrcoef(x(start_ind:end), start_ind:length(x)), ... 
    gamma_mhat_all_vals(has_results), gamma_start_inds(has_results), 'UniformOutput', false);

gamma_r_mhat = zeros(size(has_results));
gamma_pr_mhat = zeros(size(has_results));
gamma_r_mhat(has_results) = cellfun(@(x) x(1, end), gamma_r_mhat_v);
gamma_pr_mhat(has_results) = cellfun(@(x) x(1, end), gamma_pr_mhat_v);

min_p = 5/n_samples_surr;
% gamma_vals_all(has_results) = cellfun(@(x) x.mhat, all_gamma_pfit(has_results));
% gamma_preds_all(has_results) = cellfun(@(x) x.mhat, all_gamma_pred(has_results));
% rsq_vals_G(has_results) = cellfun(@(x) mean(x.r2), all_gamma_pfit(has_results));
% rsq_vals_A(has_results) = cellfun(@(x) mean(x.r2), all_alpha_pfit(has_results));
% rsq_vals_T(has_results) = cellfun(@(x) mean(x.r2), all_tau_pfit(has_results));
if strcmp(x_lab_sum, '\eta')
    tc_lims = x_var_list([1 end]) + [-1 1];
    x_full = repmat(x_var_list, length(y_var_list), 1);
else
    tc_lims = [-1.1 1.1];
    x_full = repmat(log10(x_var_list), size(alpha_vals_all, 1), 1);
end
% 
% % figure(),imagesc(gamma_preds_all)
% % figure(),imagesc(gamma_preds_all, [0 2])

%  Version 1: plot lines / dots 
makeMyFigure(28, 20);

% 
subplot(3, 3, 1)
errorbar(x_full', tau_vals_all', tau_sevals_all', 'o-', 'CapSize', 0)
xlabel(x_lab_sum)
ylabel('fitted \tau')
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)
% 

subplot(3, 3, 2)
errorbar(x_full', alpha_vals_all', alpha_sevals_all', '.-', 'CapSize', 0)
xlabel(x_lab_sum)
ylabel('fitted \alpha')
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

subplot(3, 3, 3)
errorbar(x_full', gamma_ave_mhat', gamma_se_mhat', '.-', 'capsize', 0)
xlabel(x_lab_sum)
ylabel('fitted \gamma')
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

suptitle([par_title_str ' summary'])

subplot(3, 3, 4),
plot(x_full', log10(min_p + tau_ks_pval)', 'o-')
xlabel(x_lab_sum)
xlim(tc_lims)
ylabel('\tau fit KS stat p-value (log_{10})')
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)


% 
subplot(3, 3, 5),
min_p = 4/5000;
plot(x_full', log10(min_p + alpha_ks_pval)', 'o-')
xlabel(x_lab_sum)
xlim(tc_lims)
ylabel('\alpha fit KS stat p-value (log_{10})')
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

% 
subplot(3, 3, 6)
plot(x_full', gamma_preds_all', 'o-')
xlabel(x_lab_sum)
xlim(tc_lims)
ylabel('predicted \gamma')
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

% 

subplot(3, 3, 7)
plot(x_full', log10(tau_lowercut'), 'o-')
ylabel('lower cutoff - size (log_{10})')
xlabel(x_lab_sum)
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

% 
subplot(3, 3, 8)
plot(x_full', log10(alpha_lowercut'), 'o-')
ylabel('lower cutoff - dur. (log_{10})')
xlabel(x_lab_sum)
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)

subplot(3, 3, 9)
errorbar(x_full', gamma_ave_mhat' - gamma_preds_all', gamma_se_mhat', 'o-', 'capsize', 0)
xlabel(x_lab_sum)
ylabel('difference (fit - pred)')
xlim(tc_lims)
set(gca, 'color', 'none', 'box', 'off', 'fontsize', 11)


print(gcf, '-dpdf',  ['avalanches_matlab_code/plots/' dir_name_str '/summary_exponents_lines' name_str '.pdf'])

%%  Version 2: plot colors  
at_clim = [1.7 2.3];
g_clim = [1.1 1.6];
logp_clim = [-3 -1];
xcut_clim = [0 log10(20)];
makeMyFigure(28, 20);
subplot(3, 3, 1)
imagesc(1:length(x_var_list),  1:length(y_var_list), tau_vals_all, at_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('fitted \tau')
% xlim(tc_lims)
% 

subplot(3, 3, 2)
imagesc(1:length(x_var_list),  1:length(y_var_list), alpha_vals_all, at_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('fitted \alpha')



subplot(3, 3, 3)
errorbar(x_full, gamma_ave_mhat, gamma_se_mhat, 'k.', 'capsize', 0)

imagesc(1:length(x_var_list),  1:length(y_var_list), gamma_ave_mhat, g_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('fitted \gamma')

suptitle([par_title_str ' summary'])


subplot(3, 3, 4),
min_p = 4/5000;
imagesc(1:length(x_var_list),  1:length(y_var_list), log10(min_p + tau_ks_pval), logp_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('\tau fit KS stat p-value (log_{10})')

% 
subplot(3, 3, 5),
imagesc(1:length(x_var_list),  1:length(y_var_list), log10(min_p+ alpha_ks_pval), logp_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('\alpha fit KS stat p-value (log_{10})')


subplot(3, 3, 6)
imagesc(1:length(x_var_list),  1:length(y_var_list), gamma_preds_all, g_clim)
xlabel(x_lab_sum)
ylabel(y_lab_sum)
% 
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);

colorbar
title('predicted \gamma')



subplot(3, 3, 7)
imagesc(1:length(x_var_list),  1:length(y_var_list), log10(tau_lowercut), xcut_clim)
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);

title('lower cutoff - size (log_{10})')
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar

% 
subplot(3, 3, 8)
imagesc(1:length(x_var_list),  1:length(y_var_list), log10(alpha_lowercut), xcut_clim)
title('lower cutoff - dur (log_{10})')
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);

xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar

% 
subplot(3, 3, 9)
% errorbar(x_full, gamma_ave_mhat - gamma_preds_all, gamma_se_mhat, 'k.', 'capsize', 0)

imagesc(1:length(x_var_list), 1:length(y_var_list), gamma_ave_mhat - gamma_preds_all, 0.2*[-1 1])
set(gca, 'xtick', 1:length(x_var_list), 'xticklabel', x_var_list);
set(gca, 'ytick', 1:length(y_var_list), 'yticklabel', y_var_list);
xlabel(x_lab_sum)
ylabel(y_lab_sum)
colorbar
title('difference (fit - pred)')

% 
% 
print(gcf, '-dpdf',  ['avalanches_matlab_code/plots/' dir_name_str '/summary_exponents_colors' name_str '.pdf'])

%% save analysis results for future 
% close all
% save(['avalanches_matlab_code/analysis_results/ava_analysis' name_str '_' datestr(now, 'yyyymmdd') '.mat'])

%% Function definitions

% function [ave_sizes, se_sizes, dur_bins, num_obs] = calculateAveSizeByDuration(sizes, durations)
% % Calculates the average avalanche size by duration of avalanche. Allows
% % creating duration bins larger than 1 (default dur_bin_size is 1). The
% % min_obs_cutoff drops all duration bins with fewer than the specified
% % number of observations. 
% %
%     % calculated joint counts, returns bin centers. bin width is 1. 
%     [n_cts, size_bins, dur_bins] = calculateJointHistogram(sizes, durations);
% 
%     % normalize to get p( size | duration)
%     n_cts_norm = n_cts*diag(1./sum(n_cts, 1));
%         
%     % calculate < size p(size | duration) > 
%     ave_sizes = size_bins*n_cts_norm;
%     
%     % number of observations of each duration
%     num_obs = sum(n_cts, 1);
% 
%     % calculate expectations of (s_i - mu)^2 by evaluating sum(
%     % (size - ave)^2 * p(size|duration)) (summation over duration values. 
%     var_sizes_bias = sum(bsxfun(@(x,y) (x-y).^2, size_bins', ave_sizes).*n_cts_norm, 1);
%     
%     % denominator N-1 not N
%     var_sizes = var_sizes_bias.*num_obs./(num_obs-1);
%     se_sizes = sqrt(var_sizes);
% end

% function [n_cts, size_bins, dur_bins] = calculateJointHistogram(sizes, durations)
% % This calls histcounts2 using bins for discrete data: bins are all unique
% % sizes and durations from the observed data. 
% 
%     % set bins
%     dur_bins = unique(double(durations));
%     dur_edges = [dur_bins - 0.5 dur_bins(end)+0.5];
%     size_bins = unique(sizes);
%     size_edges = [size_bins - 0.5 size_bins(end)+0.5];
%     
%     % tally observations
%     [n_cts, ~, ~] = histcounts2(sizes, double(durations), ...
%         size_edges, dur_edges, 'normalization', 'count');
% 
% end

function calc_ave_fr = numericCalculateCellFR(fr_struct)

    h1_fun = @(j_i, h) fr_struct.eta*j_i*h + fr_struct.epsilon;
    p_fire = @(j_i, h) 1./(1 + exp(-h1_fun(j_i, h)));
    
    p_h = @(h) exp(-h.^2/2)/sqrt(2*pi);
    %%
    calc_ave_fr = 0*fr_struct.J;
    for ii = 1:length(fr_struct.J)
        intgrnd = @(x) p_fire(fr_struct.J(ii), x).*p_h(x);
        calc_ave_fr(ii) = integral(intgrnd, -Inf, Inf,'AbsTol',1e-6,'RelTol',1e-3);
    end
    
%     %% simulate smapling (consistent with the analytic - sample difference)
%     sample_ave_s_i = 0*fr_struct.J;
%     num_samples = 1000;
%     h_samples = randn(num_samples, 1);
%     for ii = 1:length(fr_struct.J)
%         sample_ave_s_i(ii) = mean(p_fire(fr_struct.J(ii), h_samples));
% %         intgrnd = @(x) p_fire(fr_struct.J(ii), x).*p_h(x);
% %         ave_s_i(ii) = integral(intgrnd, -Inf, Inf,'AbsTol',1e-6,'RelTol',1e-3);
%     end
end

function p_fit = fitPowerLawBootstrap(x, y, x_cutoff, num_draws)
% fits a power-law relationship between x and y in the regions specified by
% x_cutoff and y_cutoff. points are sampled ~logarithmically spaced so that
% the unit spacing on the high end of x, y doesn't swamp R^2. 
    % reshape x, y, in case not same size
    x = x(:);
    y = y(:);
    
    if isempty(x_cutoff)
        % adaptively find reasonable values of x_cutoff
        
        x_min = min(x);
        x_max = max(x);
        n_mesh = 15;
        x_vals = linspace(x_min, x_max, n_mesh);
        min_range = 1;
        ij_off = round(min_range/(x_vals(2)-x_vals(1)));
        lenX = zeros(n_mesh);
        allRsq = zeros(n_mesh);
        xMins = zeros(n_mesh);
        xMaxs = zeros(n_mesh);
        for ii = 1:n_mesh-ij_off
            for jj = ii+ij_off:n_mesh
                try_x_cutoff = [x_vals(ii) x_vals(jj)]; 
                lenX(ii,jj) = x_vals(jj) - x_vals(ii);
                xMins(ii, jj) = x_vals(ii);
                xMaxs(ii, jj) = x_vals(jj);
                try_pfit = fitPowerLawBootstrap(x, y, try_x_cutoff, 10);
                allRsq(ii, jj) = nanmean(try_pfit.r2);
            end
        end
        
        try 
            minRsq = 0.98;
            ok_vals = find(allRsq >= minRsq);
            [largest_regions, ok_vals_ind] = max(lenX(ok_vals));
            x_cutoff = [xMins(ok_vals(ok_vals_ind)) xMaxs(ok_vals(ok_vals_ind))];
        catch
            minRsq = max(allRsq(isfinite(allRsq)));
            ok_vals = find(allRsq >= minRsq);
            [largest_regions, ok_vals_ind] = max(lenX(ok_vals));
            x_cutoff = [xMins(ok_vals(ok_vals_ind)) xMaxs(ok_vals(ok_vals_ind))];
        end
        
        if isempty(x_cutoff)
            x_cutoff = [x_min x_max];
        end
        %% now call function, with cutoff provided
        p_fit = fitPowerLawBootstrap(x, y, x_cutoff, num_draws);
    end
    

    m_pars = zeros(num_draws, 1);
    b_pars = zeros(num_draws, 1);
    r2_vals = zeros(num_draws, 1);
    
    y = y( x>x_cutoff(1) & x < x_cutoff(2));
    x = x( x>x_cutoff(1) & x < x_cutoff(2));
    
    % probabbility of keeping x is proportional to 
    prob_keep_x = diff((x));
    % scale so that ~20% of points are kept
    prob_keep_x = prob_keep_x/mean(prob_keep_x)*0.2;
    % reset values > 1 to 1, then rescale again 
    prob_keep_x(prob_keep_x > 1) = 1;
    prob_keep_x = prob_keep_x/mean(prob_keep_x)*0.5;
    
    for i_draw = 1:num_draws
        keep_inds = rand(size(prob_keep_x)) < prob_keep_x;
        N = sum(keep_inds);
        x_samp = x(keep_inds);
        y_samp = y(keep_inds);
        
%         mod = fitlm(x_samp, y_samp)
        mb = polyfit(x_samp, y_samp, 1);
        m_pars(i_draw) = mb(1);
        b_pars(i_draw) = mb(2);
        
        r2_vals(i_draw) = 1 - ((N-1)/(N-2))*sum((polyval(mb, x_samp) - y_samp).^2)/sum((y_samp -mean(y_samp)).^2);

        
    end
    p_fit.bhat = mean(b_pars);
    p_fit.mhat = mean(m_pars);
    p_fit.mSE = std(m_pars);
    p_fit.m = m_pars;
    p_fit.b = b_pars;
    p_fit.r2 = r2_vals;
    p_fit.x_cutoff = x_cutoff;
end

function [p_fit, cdf_p_fit] = fitPowerLawDecade(x, y, xcutoff, num_draws, is_dist)
% runs bootstrap powerlaw fit decade by decade  
    % reshape x, y, in case not same size
    x = x(:);
    y = y(:);
    if is_dist
        %% this has not been debugged -- using Clauset method instead for fitting distributions
        % y has log10(raw counts) in it
        y_cts = round(10.^y);
    %     bin_size = round(diff([0; 10.^x]));
        y_cts_norm = y_cts/sum(y_cts);
        % fit the complementary cdf not pdf)
        y_fit = log10( 1 - cumsum(y_cts_norm));
    else
        y_fit = y;
    end
    %%
    xcut_max = floor(max(x));
    xcuts = 1:0.25:xcut_max;
    xcut_range = length(xcuts);
    cdf_p_fit = cell(xcut_range, 1);
    for ii = 1:xcut_range
        cdf_p_fit{ii} = fitPowerLawBootstrap(x, y_fit, xcuts(ii) + [-1 0], num_draws);
    end
    %% get p_fit over specified region
    in_region = cellfun(@(x) x.x_cutoff(1) > xcutoff(1) & x.x_cutoff(2) < xcutoff(2), cdf_p_fit);
    
    p_fields = fieldnames(cdf_p_fit{1});
    for ii_f = 1:length(p_fields)
        p_fit.(p_fields{ii_f}) = cell2mat(cellfun(@(x) x.(p_fields{ii_f}), cdf_p_fit(in_region)', 'uniformoutput', false));
    
    end
    if is_dist
        % subtract 1 from the m values
        p_fit.mhat = p_fit.mhat - 1;
        p_fit.m = p_fit.m - 1;
    end
end


function [x_min, fit_res] = findXMin_twoDecade(x_i, x_min_vals)
% uses KS statistic between empirical and model distributions to determine
% x_min for power-law fitting. Fits over two decades from each value of
% x_min_vals. 

    % if not provided, try values between 1 and 100. remember x_i are
    % discrete. 
    if nargin < 2
        x_min_vals = 1:100;
    end
    
%     % if provided, drop all x_i that are greater than x_max_val
%     if nargin >= 3
%         x_i = x_i(x_i <= x_max_val);
%     end
    
    a_vals = 0*x_min_vals;
    se_a_vals = 0*x_min_vals;
    for ii = 1:length(a_vals)
        use_x_i = x_i(x_i <= 100*x_min_vals(ii));
        
        a_vals(ii) = powerLawMLEClauset(use_x_i, x_min_vals(ii));
        se_a_vals(ii) = findSEpowerlawslope(a_vals(ii), x_min_vals(ii), sum(use_x_i >= x_min_vals(ii)));
        disp(['x_min = ' num2str(x_min_vals(ii)) ' done'])
    end
    disp('done calculating MLE for slope')
    % calculate KS statistic between true CDF and model CDF. 
    ks_stats = 0*x_min_vals;
%%     ks_p_vals = 0*x_min_vals;
    for ii = 1:length(a_vals)

        ks_stats(ii) = computeKSStatPowerLaw(x_i, a_vals(ii), x_min_vals(ii), x_min_vals(ii)*100);
    end
%%    
    % the best x_min minimizes the distnace between model and empirical
    % data
    [~, ind] = min(ks_stats);

    x_min = x_min_vals(ind(1));
    a_hat = a_vals(ind);
    n_obs = sum(x_i >= x_min);
    se_a_hat = findSEpowerlawslope(a_hat, x_min, n_obs);
    fit_res.a_hat = a_hat;  % return the MLE slope at the selected xmin value
    fit_res.se_a_hat = se_a_hat;    % SE in slope
    fit_res.x_min = x_min;          % best-fit xmin value
    fit_res.x_min_vals = x_min_vals;% all candidate x_min values
    fit_res.a_vals = a_vals;        % all MLE slopes for each xmin
    fit_res.se_a_vals = se_a_vals;  % all SE of MLE slopes for each xmin.
    fit_res.ks_stats = ks_stats;    % ks_stats for model and empirical powerlaws
    
end


function [x_min, fit_res] = findXMin(x_i, x_min_vals, x_max_val, max_tol)
% uses KS statistic between empirical and model distributions to determine
% x_min for power-law fitting

    % if not provided, try values between 1 and 100. remember x_i are
    % discrete. 
    if nargin < 2
        x_min_vals = 1:100;
    end
    
    % if provided, drop all x_i that are greater than x_max_val
    if nargin >= 3
        x_i = x_i(x_i <= x_max_val);
    end
    
    a_vals = 0*x_min_vals;
    se_a_vals = 0*x_min_vals;
    for ii = 1:length(a_vals)
        a_vals(ii) = powerLawMLEClauset(x_i, x_min_vals(ii));
        se_a_vals(ii) = findSEpowerlawslope(a_vals(ii), x_min_vals(ii), sum(x_i >= x_min_vals(ii)));
        disp(['x_min = ' num2str(x_min_vals(ii)) ' done'])
    end
    disp('done calculating MLE for slope')
    % calculate KS statistic between true CDF and model CDF. 
    ks_stats = 0*x_min_vals;
%%     ks_p_vals = 0*x_min_vals;
    for ii = 1:length(a_vals)
        frac_tail = mean(x_i > x_min_vals(ii));

        ks_stats(ii) = computeKSStatPowerLaw(x_i, a_vals(ii), x_min_vals(ii), max(x_i));
    end
%%    
    % the best x_min minimizes the distnace between model and empirical
    % data
    [min_KS_stat, ind] = min(ks_stats);
    % if "max_tol" is provided, accept the first value that is less than
    % max_tol*min(ks_stats)
    if nargin >= 4
        max_tol = max(max_tol, 1);  % must be 1 or larger, return warning if not
        min_val = max_tol*min(ks_stats);
        [~, ind] = find(ks_stats < min_val);
        ind = ind(1);
        min_KS_stat = ks_stats(ind);
    end

    x_min = x_min_vals(ind);
    a_hat = a_vals(ind);
    n_obs = sum(x_i >= x_min);
    se_a_hat = findSEpowerlawslope(a_hat, x_min, n_obs);
    fit_res.a_hat = a_hat;  % return the MLE slope at the selected xmin value
    fit_res.se_a_hat = se_a_hat;    % SE in slope
    fit_res.min_KS_stat = min_KS_stat;
    fit_res.x_min = x_min;          % best-fit xmin value
    fit_res.x_min_vals = x_min_vals;% all candidate x_min values
    fit_res.a_vals = a_vals;        % all MLE slopes for each xmin
    fit_res.se_a_vals = se_a_vals;  % all SE of MLE slopes for each xmin.
    fit_res.ks_stats = ks_stats;    % ks_stats for model and empirical powerlaws
    
end

function [x_minmax, fit_res] = findXMinXMax(x_i, x_min_vals, x_max_vals)
% uses KS statistic between empirical and model distributions to determine
% lower and upper cutoffs for power-law fitting

    % if not provided, try values between 1 and 100. remember x_i are
    % discrete. 
    if nargin < 2
        x_min_vals = 1:100;
    end
    %%
    a_vals = nan(length(x_max_vals), length(x_min_vals));
    se_a_vals =nan(length(x_max_vals), length(x_min_vals));
    % require at least 500 avalanches in the tail to fit a power law
    min_num_obs = 500;
    for ii = 1:length(x_max_vals)
        % only fit power law using observations up to the upper cutoff
        use_x_i = x_i(x_i <= x_max_vals(ii));
        for jj = 1:length(x_min_vals)
            if sum(use_x_i > x_min_vals(jj)) > min_num_obs
                a_vals(ii, jj) = powerLawMLEClauset(use_x_i, x_min_vals(jj));
                se_a_vals(ii, jj) = findSEpowerlawslope(a_vals(ii, jj), ...
                    x_min_vals(jj), sum(use_x_i >= x_min_vals(jj)));
                disp(['     x_min = ' num2str(x_min_vals(jj)) ' done'])
            end
        end
        disp(['X_MAX = ' num2str(x_max_vals(ii)) ' done'])

    end
    disp('done calculating MLE for slope')
    %%
    
    % calculate KS statistic between true CDF and model CDF. 
    ks_stats = nan(length(x_max_vals), length(x_min_vals));
    n_obs_tot = nan(length(x_max_vals), length(x_min_vals));
    n_obs_tail = nan(length(x_max_vals), length(x_min_vals));
    for jj = 1:length(x_min_vals)
        n_obs_tot(ii, :) = length(use_x_i);

        if x_min_vals(jj) < max(x_max_vals)
            [ks_stats(:,jj), n_obs_tail(:, jj)] = computeKSStatPowerLaw(use_x_i, a_vals(ii, jj), x_min_vals(jj), x_max_vals);
        end
        disp(['     xmin = ' num2str(x_min_vals(jj)) ' done'])

    end
%%    
    % the best x_min minimizes the distance between model and empirical
    % data
    min_KS_stat = min(ks_stats, [], 'all');
    [i_minKS, j_minKS] = find(ks_stats == min_KS_stat);


    x_min = x_min_vals(j_minKS(1));
    x_max = x_max_vals(i_minKS(1));
    x_minmax = [x_min x_max];
    a_hat = a_vals(i_minKS(1), j_minKS(1));
    fit_res.n_obs_ptail = sum(x_i >= x_min & x_i < x_max);  % number of avalanches in power-law tail
    se_a_hat = se_a_vals(i_minKS(1), j_minKS(1));
    fit_res.a_hat = a_hat;  % return the MLE slope at the selected xmin value
    fit_res.se_a_hat = se_a_hat;    % SE in slope
    fit_res.min_KS_stat = min_KS_stat;
    fit_res.x_min = x_min;          % best-fit xmin value
    fit_res.x_max = x_max;          % best-fit xmax value
    fit_res.x_min_vals = x_min_vals;% all candidate x_min values
    fit_res.a_vals = a_vals;        % all MLE slopes for each xmin
    fit_res.se_a_vals = se_a_vals;  % all SE of MLE slopes for each xmin.
    fit_res.ks_stats = ks_stats;    % ks_stats for model and empirical powerlaws
    
end



function se_a = findSEpowerlawslope(a_hat, x_min, n)
% evaluates equation (3.6) of Clauset, Shalizi, and Newman 2009 to estimate
% the standard error of the estimate of alpha. 

    se_a = 1/sqrt(n*(hurwitzZeta(2, a_hat, x_min)/hurwitzZeta(a_hat, x_min) - ...
        (hurwitzZeta(1, a_hat, x_min)/hurwitzZeta(a_hat, x_min))^2));

end


function [ks_stats_D, surr_a_est] = computeKSStatSurrogateData(a_hat, x_min, x_max_val, n_total, n_tail, emp_bins, num_surr)
% generate surrogate datasets and measure the KS statistic between model
% and surrogate dataset
% ks_stats_D: over power law only (eqn 3.9)
% ks_stats_Surr: over entire region

% PROBLEM WITH THE EMPIRICAL BINS: they are very gappy for the large
% numbers, and we should be careful about how we draw from them. one idea
% is to interpolate the power law CDF to all possible size values, to get
% better samples... then bin with the empirical bins? hmmmmmmm
%     cdf_model = empirical_cdf;

    tail_bins = emp_bins(emp_bins >= x_min);
%     pretail_bins = emp_bins(emp_bins < x_min);

    % THis is hack to avoid out of bound problems
    tail_bins_eval = tail_bins;
%     x_max_val = 1e4;
    tail_bins_eval(tail_bins_eval > x_max_val) = x_max_val;    
    %         cdf_model = 1 - computeCDFpowerlawdistribution(x_bins_eval, a_val, x_min);

    powerLaw_ccdf = computeCDFpowerlawdistribution(tail_bins_eval, a_hat, x_min);
    pe4_val = computeCDFpowerlawdistribution(x_max_val, a_hat, x_min);
    % super dumb bad hack to interpolate from 1e4 to max(tail_bins) - should at
    % least do this in log-space to preserve power law 
    powerLaw_ccdf(tail_bins > x_max_val) = interp1([x_max_val max(tail_bins)], [pe4_val 0], tail_bins(tail_bins > x_max_val));


    % tail_bins_approx = unique(round(logspace(log10(tail_bins(1)), log10(tail_bins(end)), 201)));
    %%
%     scale_fac = 1 - empirical_cdf(emp_bins == x_min);
    %%
%     cdf_model(emp_bins >= x_min) = (1 - scale_fac) + scale_fac*(1 - powerLaw_ccdf);
    %%

%     cdf_preTail = empirical_cdf(emp_bins < x_min);
%     cdf_preTail = cdf_preTail/cdf_preTail(end);

%%
%     ks_stats_Surr = zeros(num_surr, 1);
    ks_stats_D = zeros(num_surr, 1);
    surr_a_est = zeros(num_surr, 1);
    tic
    for ii = 1:num_surr
        rng('shuffle')
        % with probability n_tail/n_total, each surrogate observation is assigned
        % to the tail
        num_tail_surr = sum(rand(n_total, 1) < n_tail/n_total);

        % those not in the tail are in the "pre-tail"
    %     num_preTail_surr = n_total - num_tail_surr;
%%

        tail_draws = rand(num_tail_surr, 1);
        tail_x_i_draws = round(interp1(1 - powerLaw_ccdf, tail_bins, tail_draws, 'previous'));
        
     %%   
        % drop tail draws that are larger than the x_max_val
        tail_x_i_draws = tail_x_i_draws(tail_x_i_draws <= x_max_val);
        
        a_s_hat = powerLawMLEClauset(tail_x_i_draws, x_min);
        surr_a_est(ii) = a_s_hat;
        ks_stats_D(ii) = computeKSStatPowerLaw(tail_x_i_draws, a_s_hat, x_min, x_max_val); 

    %     pretail_draws = rand(num_preTail_surr, 1);
    %     preTail_x_i_draws = interp1([0; cdf_preTail], [0 pretail_bins], pretail_draws, 'next');
    % 
    %     surr_cdf = histcounts([preTail_x_i_draws; tail_x_i_draws], [emp_bins Inf], ... 
    %         'Normalization', 'cdf');
    %     ks_stats_Surr(ii) = max(abs(surr_cdf - cdf_model'));
    %     
    %     % over powerlaw only; interpolating from 100 log-spaced bins in the tail
    %     % is a numerical error ov about 10^-5 (for values checked, xmin = 3, a
    %     % = 1.9)
    %     surr_a_est(ii) = a_s_hat;
    %     
    %     est_powerLaw_ccdf = 1 - computeCDFpowerlawdistribution(tail_bins_eval, a_s_hat, x_min);
    % 
    % %     pL_ccdf_approx = 1 - computeCDFpowerlawdistribution(tail_bins_approx, a_s_hat, x_min);
    % %     est_powerLaw_ccdf = interp1(tail_bins_approx, pL_ccdf_approx, tail_bins, 'spline');
    %     
    % %     powerLaw_ccdf = 1 - computeCDFpowerlawdistribution(tail_bins, a_s_hat, x_min);
    %     surr_PL_cdf = histcounts(tail_x_i_draws, [0 tail_bins-.5], 'Normalization', 'cdf');
    %     ks_stats_D(ii) = max(abs(surr_PL_cdf - est_powerLaw_ccdf));

    if mod(ii, 100) == 0
        disp(['shuffle ' num2str(ii) ' done'])
    end
    end
end

