% 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
% 
clear;
% ss_str.run_f5ultrafinesweep = true;       % pull data from 5-field, ultrafine sweep
% ss_str = selectSimulation('run_morrell1field');
% ss_str.run_f1smallfinesweep = true;       % pull data from 1-field,small network, ultrafine sweep

% have to set the second input to the appropriate string (just one at a
% time, so 'A' (rather than a list of replicates). 
ss_str = selectSimulation('run_f1smallultrafinesweep', 'A');

% set only one of the above to true/false, then run this script
setLoadingFunctionNames
% %% LOAD THE HURWITZ ZETA LOOKUP TABLE
% 
% [glu_hurwitzZeta, glu_dhurwitzZeta] = loadhurwitzZetaLookup()



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

% for 1-field plots, may need smaller cutoffs - avalanches are smaller.
% (s = 1e3, d= 800?) 
% for 5-field, try s = 5e3 and d = 1e3
dur_max = 3e3; % 1e3 : track length duration is 10,000; set max duration to be 30% of that
size_max = 3e4; % the size scales with duration^1.2 or so; 1000^1.3 is 8e3. 

% maximum number of avalanches used in statistics (drawn randomly)
n_obs_max = 500000;
n_obs_min = 10000;

% fit power law, boot strap (for size-duration scaling only)
num_draws = 100;

% number of surrogates to draw for model determined by selected value of lower cut-off
n_samples_surr = 2000;

% require size of at least 3
x_min_vals = union(unique(round(logspace(0.5, 2, 21))), unique(round(logspace(2, 2.5, 5))));
    
% 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));
all_gamma_powerfit = 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));

% save gamma fit checks
gamma_check_dir = [plot_dir '/gamma_fit_check/'];
if ~exist(gamma_check_dir, 'dir')
    mkdir(gamma_check_dir);
end
%% 
for ii_reps = 1:length(y_var_list)
    for ii_time = 1:length(x_var_list)
         %%
        try
            x_ava = load(ava_name_fn(ii_reps, ii_time));
            x_fr = load(fr_name_fn(ii_reps, ii_time));
            % 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 1000 avalanches
            if length(x_ava.sizes) > n_obs_min
                
                % select up to 40000 of the avalanches
                sel_events = unique(ceil(length(x_ava.sizes)*rand(n_obs_max, 1)));
                
                x_ava.sizes = x_ava.sizes(sel_events);
                x_ava.durations = x_ava.durations(sel_events);
                
                
                
                [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);
                
                % determine system size cutoffs
                %%
                % use Clauset algorithm to determine x_min for power-law fit to size
                % distributio
                max_tol = 1.1; %(find minimum within max_tol times true minimum)
                [x_min_S, pfit_S] = findXMin(x_ava.sizes, x_min_vals, size_max, max_tol);
                %     [x_min_S2, pfit_S2] = findXMin_twoDecade(x_ava.sizes, x_min_vals);
                
                %%
                % %     %%
                % %     t1 = tic;
                % %     [x_minmax_S, pfit_rng_S] = findXMinXMax(x_ava.sizes, x_min_vals, x_max_vals);
                % %     disp(['Min-Max scan for sizes completes in ' num2str(toc(t1), '%1.0f') ' s.'])
                %
                %     d_cutoff = [0.5 2.25];
                [x_min_D, pfit_D] = findXMin(double(x_ava.durations), x_min_vals, dur_max, max_tol);
                %%
                % %     t1 = tic;
                % %     [x_minmax_D, pfit_rng_D] = findXMinXMax(double(x_ava.durations), x_min_vals, x_max_vals);
                % %     disp(['Min-Max scan for durations completes in ' num2str(toc(t1), '%1.0f') ' s.'])
                
                %
                %% calculate gamma_pred
                g_pred_all =bsxfun(@(x, y) (x - 1)./(y - 1), pfit_D.a_vals, pfit_S.a_vals');
                rel_err_D = pfit_D.se_a_vals./pfit_D.a_vals;
                rel_err_S = pfit_S.se_a_vals./pfit_S.a_vals;
                g_predSE = g_pred_all.*bsxfun(@(x,y) sqrt(x.^2 + y.^2), rel_err_D, rel_err_S');
                %     g_predSE = g_pred.*sqrt((p_fit_S.mSE/p_fit_S.mhat)^2 + (p_fit_A.mSE/p_fit_S.mhat)^2);
                
                % use the linear fit method for the size vs duration power law
                %     p_fit_G = fitPowerLawBootstrap(log10(dur_bins), log10(ave_sizes), log10(x_min_D) + [0 2], num_draws);
                [~, pfit_G] = fitPowerLawDecade(log10(dur_bins), log10(ave_sizes), [0 2.5], num_draws, false);
                
                
                [gamma_fit, gamma_range, minS_gamma_fit, minS_gamma_range] = ...
                    fitGammaRange(dur_bins, ave_sizes, true, ...
                    [param_str(ii_reps, ii_time) '_gammafit.pdf']);
                
                gamma_fit_info.gamma_fit = gamma_fit;
                gamma_fit_info.gamma_range = gamma_range;
                gamma_fit_info.minS_gamma_fit = minS_gamma_fit;
                gamma_fit_info.minS_gamma_range = minS_gamma_range;
                
                %% generate surrogate datasets for goodness of fit
                n_total = length(x_ava.sizes);
                n_tailS = sum(x_ava.sizes >= pfit_S.x_min);
                n_tailD = sum(x_ava.durations >= pfit_D.x_min);
                size_cts = sum(n_cts, 2);
                dur_cts = sum(n_cts, 1)';
                size_cdf = cumsum(size_cts)/sum(size_cts);
                dur_cdf = cumsum(dur_cts)/sum(dur_cts);
                %%
                t1 = tic;
                [ks_stats_Surr_S, surr_a_S] = computeKSStatSurrogateData(pfit_S.a_hat, ...
                    pfit_S.x_min, size_max, n_total, n_tailS, size_bins, n_samples_surr);
                disp(['Surrogate for sizes complete in ' num2str(toc(t1), '%1.0f') ' s.'])
                
                %%
                t1 = tic;
                [ks_stats_Surr_D, surr_a_D] = computeKSStatSurrogateData(pfit_D.a_hat, ...
                    pfit_D.x_min, dur_max, n_total, n_tailD, dur_bins, n_samples_surr);
                
                disp(['Surrogate for durations complete in ' num2str(toc(t1), '%1.0f') ' s.'])
                
                %% Save surrogate draws
                pfit_S.ks_surrogate = ks_stats_Surr_S;
                pfit_S.a_est_surr = surr_a_S;
                pfit_D.ks_surrogate = ks_stats_Surr_D;
                pfit_D.a_est_surr = surr_a_D;
                
                %% save everything
                avalanche_count(ii_reps, ii_time) = length(x_ava.sizes);
                
                all_tau_pfit{ii_reps, ii_time} = pfit_S;
                all_alpha_pfit{ii_reps, ii_time} = pfit_D;
                all_gamma_pred{ii_reps, ii_time}.mhat = g_pred_all;
                all_gamma_pred{ii_reps, ii_time}.mSE = g_predSE;
                all_gamma_pred{ii_reps, ii_time}.g_hat = (pfit_D.a_hat - 1)/(pfit_S.a_hat - 1);
                
                all_gamma_pfit{ii_reps, ii_time} = pfit_G;
                all_gamma_powerfit{ii_reps, ii_time} = gamma_fit_info;
                
                close all;
                save([save_dir 'ava_decade_analysis_' datestr(now, 'yyyymmdd')])
                
                
            end
        end
        close all
    end
    
end

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

%% 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)

    % 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



