text:
using Images, FileIO, Colors, FixedPointNumbers

# global variables. as a matter of convention, each one ends with an
# underscore. this is to make it harder to use them within the main body of the
# code without explicitly doing so.
# (width, height) of system in voxels
sys_size_ = (140,210);
# (dx, dy) voxel lengths
voxel_lengths_ = (0.12,0.12);
#time steps
dt_ = 1.;
#max time of simulation
max_t_ = 10000;
#min and max values of k and F
# #goes along y axis
# k_ = (0.0, 0.07);
# #goes along x axis
# F_ = (0.04, 0.16);
#goes along y axis
F_ = (0.0, 0.04);
#goes along x axis
k_ = (0.03, 0.07);
# #(Du Dv): diffusion constants
# D_ = (2e-5, 1e-5);
#(Du Dv): diffusion constants
D_ = (2e-3, 1e-3);
#proportion of squares to initialise
random_init_ = 0.99;
#choose boundary conditions, either dirichlet or neumann
diff_mode_ = "neumann";
checkpoint_int_ = 5;
checkpoint_name_ = "./gray_scott"

function create_initial_conditions((num_x, num_y), random_init)
    u = rand(num_x, num_y)
    @. u[u < random_init] = 0
    @. u[u >= random_init] = 1
    v = 1 .- u
    (u,v)
end;

function create_reaction_rates((min_k,max_k), (min_F, max_F), (num_x, num_y))
    ks = reshape(repeat(range(min_k, stop=max_k, length=num_x), num_y), (num_x,num_y))
    Fs = transpose(reshape(repeat(range(min_F, stop=max_F, length=num_y), num_x), (num_y,num_x)))
    reaction_rates = zeros(2, num_x, num_y)
    reaction_rates[1,:,:] = ks
    reaction_rates[2,:,:] = Fs
    reaction_rates
end;

function diffusion_voxel(i,j, dx_diff, a, (dx,dy), (num_x, num_y), dt, D, diff_mode)
    if diff_mode=="neumann"
        #diff left
        if i - 1 >= 1
            dx_diff[i-1,j] += a[i,j] * dt * D / (dx^2)
            dx_diff[i,j] -= a[i,j] * dt * D / (dx^2)
        end
        #diff right
        if i + 1 <= num_x
            dx_diff[i+1,j] += a[i,j] * dt * D / (dx^2)
            dx_diff[i,j] -= a[i,j] * dt * D / (dx^2)
        end
        #diff up
        if j - 1 >= 1
            dx_diff[i,j-1] += a[i,j] * dt * D / (dy^2)
            dx_diff[i,j] -= a[i,j] * dt * D / (dy^2)
        end
        #diff down
        if j + 1 <= num_y
            dx_diff[i,j+1] += a[i,j] * dt * D / (dy^2)
            dx_diff[i,j] -= a[i,j] * dt * D / (dy^2)
        end
    elseif diff_mode=="dirichlet"
        #diff left
        if i - 1 >= 1
            dx_diff[i-1,j] += a[i,j] * dt * D / (dx^2)
        end
        #diff right
        if i + 1 <= num_x
            dx_diff[i+1,j] += a[i,j] * dt * D / (dx^2)
        end
        #diff up
        if j - 1 >= 1
            dx_diff[i,j-1] += a[i,j] * dt * D / (dy^2)
        end
        #diff down
        if j + 1 <= num_y
            dx_diff[i,j+1] += a[i,j] * dt * D / (dy^2)
        end
        dx_diff[i,j] -= 4 * a[i,j] * dt * D / (dy^2)
    end
    dx_diff
end

function diffusion_species(a, D, (dx, dy), (num_x, num_y), dt, diff_mode)
    dx_diff = zeros(num_x, num_y)
    for i in 1:num_x
        for j in 1:num_y
            dx_diff = diffusion_voxel(i,j, dx_diff, a, (dx,dy), (num_x, num_y), dt, D, diff_mode)
        end
    end
    dx_diff
end

function diffusion_step((u,v), (Du, Dv), (dx, dy), (num_x, num_y), dt, diff_mode)
    du_diff = diffusion_species(u, Du, (dx, dy), (num_x, num_y), dt, diff_mode)
    dv_diff = diffusion_species(v, Dv, (dx, dy), (num_x, num_y), dt, diff_mode)
    (du_diff, dv_diff)
end

function reaction_step((u,v), reaction_rates, (dx,dy), (num_x, num_y), dt)
    ks = reaction_rates[1,:,:]
    Fs = reaction_rates[2,:,:]
    non_linear = zeros(num_x, num_y)
    du_reac = zeros(num_x, num_y)
    dv_reac = zeros(num_x, num_y)
    @. non_linear = (v^2) * u * dt
    @. du_reac = -non_linear + (Fs*(1-u) * dt)
    @. dv_reac = non_linear - ((ks + Fs) * v * dt)
    (du_reac, dv_reac)
end;

function pde_step!((u, v), D, reaction_rates, voxel_lengths, sys_size, dt, diff_mode)
    du_reac, dv_reac = reaction_step((u,v), reaction_rates, voxel_lengths, sys_size, dt)
    du_diff, dv_diff = diffusion_step((u,v), D, voxel_lengths, sys_size, dt, diff_mode)
    @. u += du_diff + du_reac
    @. v += dv_diff + dv_reac
    (u, v)
end;

function pixel_color_transform(pixel)
    output = zeros(3)
    #red to blue colormap
    # if pixel <= 0
        # output = RGB(1,0,0)
    # elseif pixel < 0.25
        # output = RGB(1-4pixel, 4pixel, 0)
    # elseif pixel <=0.5
        # output = RGB(0, 1-4(pixel-0.25), 4(pixel-0.25))
    # else
        # output = RGB(0,0,1)
    # end

    #white to blue colormap
    if pixel <= 0
        output = RGB(1.,1.,1.)
    elseif pixel <=0.5
        output = RGB(1-2pixel,1-2pixel,1)
    else
        output = RGB(0.,0.,1.)
    end
    output
end

function render_img((u,v))
    x,y = size(u)
    canvas = zeros(RGB,x,y)
    for i in 1:x
        for j in 1:y
            canvas[i,j] = pixel_color_transform(v[i,j])
        end
    end
    canvas
end

function set_up_simulate(max_t, dt, sys_size, random_init, k, F)
    num_steps = round(max_t/dt)
    u, v = create_initial_conditions(sys_size, random_init)
    reaction_rates = create_reaction_rates(k, F, sys_size)
    (num_steps, (u,v), reaction_rates)
end

function simulate!((u,v), D, reaction_rates, voxel_lengths, sys_size, dt, num_steps, diff_mode, checkpoint_int, checkpoint_name)
    x,y = size(u)
    num_frames = floor(Int, num_steps/checkpoint_int)
    index_count = 1
    for i in 1:num_steps
        u,v = pde_step!((u,v), D, reaction_rates, voxel_lengths, sys_size, dt, diff_mode)
        if ((i%checkpoint_int) - 1) == 0
            println("Processing step ", i, " out of ", num_steps)
            img = render_img((u,v))
            save(string(checkpoint_name,"_",index_count,".jpg"), img);
            index_count += 1
        end
    end
    (u,v)
end

function main(max_t, dt, sys_size, random_init, k, F, D, voxel_lengths, diff_mode, checkpoint_int, checkpoint_name)
    num_steps, (u,v), reaction_rates = set_up_simulate(max_t, dt, sys_size, random_init, k, F)
    u, v = simulate!((u,v), D, reaction_rates, voxel_lengths, sys_size, dt, num_steps, diff_mode, checkpoint_int, checkpoint_name)
end

main(max_t_, dt_, sys_size_, random_init_, k_, F_, D_, voxel_lengths_, diff_mode_, checkpoint_int_, checkpoint_name_)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ ee9049ce-90f4-11eb-0130-db0c30bb4ab6
begin
	using CSV, DataFrames, DataFramesMeta
	using CategoricalArrays
	using Statistics, StatsBase
	using PyCall
	using MLJ
	using MLJLinearModels
	using MLJDecisionTreeInterface

	skl = pyimport("sklearn")
	pp = pyimport("sklearn.preprocessing")
end

# ╔═╡ 9cff05d2-9112-11eb-190d-23f1f040f5de
using Plots

# ╔═╡ 1122232c-90f5-11eb-0ddd-8972932ead6b
md"
# ADNI Neuroimaging Data Analysis Part 1: Alzheimer's Prediction

This dataset originates from the ADNI is the Alzheimers Disease Neuroimaging Initative. The goal in this notebook is to develop a model to predict Alzheimer's disease (and the probability) based on 579 lipidomics (from mass spectrometry) and 116 MRI neuroimaging features.

CN = Cognitively Normal

LMCI = Late Mild Cognitive Impairment

AD = Alzheimer's Disease

It is a cross sectional dataset consists of 760 total patients (213 CN, 375 LMCI, and 172 AD). For now, the LMCI patients will be removed, and models will be developed to predict the probability of Alzheimer's disease at baseline. Later on, in Part 2, the best model (as assessed by cross validation in this notebook) will be used to predict P(AD) for the LMCI patients.

The patients in this lipidomic/MRI baseline cross-sectional study were also followed for 8 years (in a separate longitudinal dataset), and we will determine whether a greater P(AD) at baseline in the LMCI patients is correlated with a worse cognitive outcome/trajectory over 8 years.

"

# ╔═╡ a312b3b0-90f9-11eb-2c9c-5f40f97dd471
rawdata = CSV.read("ADNI_cross_sect.csv",DataFrame); dropmissing!(rawdata)

# ╔═╡ d4f179a2-90f9-11eb-0e1a-2ff102a97066
md"

## Preprocessing

The first thing we do is split off the lipid features from the MRI features. The MRI features need to be normalized by the Estimated Total Intracranial Volume (ETIV) column. Following this we split off the LMCI patients.

For the classifier models in MLJ, the labels (DX_bl) must be converted to ordered categorical factors (with CN as the base level).
"

# ╔═╡ 2a4b1ffe-90fa-11eb-1e8f-8beeeb8a2085
begin

#Preprocess MRI features beforehand by dividing by ETIV (Est. Total Intracranial volume)

lipid_data = rawdata[:,2:580]
MRI_data_pre = rawdata[:,581:696]

names_lipid = names(lipid_data)
names_MRI = names(MRI_data_pre)

ETIV = rawdata.ETIV
DX_bl = rawdata[:,[:DX_bl]]

MRI_data = MRI_data_pre ./ ETIV

newdata = hcat(lipid_data,MRI_data,DX_bl)

#Split off non MCI (CN and AD) from MCI for now

MCIdata = @linq newdata |> where(:DX_bl .== "LMCI")
NoMCIdata = @linq newdata |> where(:DX_bl .!= "LMCI")

MCIdata.DX_bl = categorical(MCIdata.DX_bl)
NoMCIdata.DX_bl = categorical(NoMCIdata.DX_bl,levels=["CN","AD"],ordered=true)

end

# ╔═╡ 96c864aa-90fa-11eb-2ca9-979422c5a39e
md"

For the cross validation loop, we need to have a train-test split function. Additionally, we will be testing 5 models: 4 L2-regularized Logistic Regressions with no features log transformed, lipids log transformed, MRI volumes log transformed, and both log transformed as well as a Random Forest model (on untransformed features). Additionally, a 6th model which is the average of the 4 logistic model predictions will be tested as well

The data (after transformation, if any) will be robustly scaled (by median and IQR) to keep features on the same units. Robust scaling is used so that potential outliers in the 579+116 lipid+MRI features do not affect the scaling as much. We use scikit-learn from Python through Julia to perform the scaling. The sklearn RobustScaler fitted object is stored in a Julia struct that contains a bool indicating whether or not log transform was applied to lipids or MRI volumes.

"

# ╔═╡ 7bad01a2-90fb-11eb-183c-c371719efda8
begin

function train_test_split(data::DataFrame,response::Symbol;prop_train=0.8)
    ntot = nrow(data)
    n_train = Int(floor(prop_train*ntot))
    train_idx = sample(1:ntot,n_train,replace=false)

    train_X = data[train_idx,Not(response)]
    train_Y = data[train_idx,response]

    test_X = data[Not(train_idx),Not(response)]
    test_Y = data[Not(train_idx),response]

    return train_X,train_Y,test_X,test_Y
end

struct PreprocessParams
    lipidnames :: Vector{String}
    MRInames :: Vector{String}
    log_lipid::Bool
    log_MRI::Bool
    RS:: PyObject
end

"""Preprocess(data;log_lipid,log_MRI)

This function does further preprocessing including log transforming the
lipid or MRI features separately. And then using RobustScaler() to scale
the data for use in ML modeling.

Note: MRI features should be normalized by ETIV beforehand

"""
function preprocess(data,lipidnames,MRInames;log_lipid=false,log_MRI=false)
    lipids = Matrix(data[:,lipidnames])
    MRI = Matrix(data[:,MRInames])
    if log_lipid == true
        lipids = log.(lipids)
    end
    if log_MRI == true
        MRI = log.(MRI)
    end
    X = hcat(lipids,MRI)
    RS = pp.RobustScaler()
    RS.fit(X)
    return PreprocessParams(lipidnames,MRInames,log_lipid,log_MRI,RS)
end

function MLJ.transform(preprocessor::PreprocessParams,data::DataFrame)
    lipids = Matrix(data[:,preprocessor.lipidnames])
    MRI = Matrix(data[:,preprocessor.MRInames])
    if preprocessor.log_lipid==true
        lipids = log.(lipids)
    end
    if preprocessor.log_MRI==true
        MRI = log.(MRI)
    end
    X = hcat(lipids,MRI)
    X_scaled = preprocessor.RS.transform(X)
    return(X_scaled)
end

end

# ╔═╡ 1469e9fa-90fc-11eb-1414-ff8d8148f92b
md"

## Evaluation Metrics

In addition to using the typical classification evaluation metrics (Accuracy, Sensitivity, Specificity) we would also like to assess the probability predictions of CN/AD in a continuous sense. To do so, we can use the **cross entropy loss**. In order to more easily interpret this loss, we normalize it by the cross entropy (logistic) loss corresponding to probability predictions that are simply the mean values of CN/AD in the test set. This is known as a test set **Pseudo-R2**.

$R^{2}_{ps} = 1 - \frac{CE(model)}{CE(mean)}$

"

# ╔═╡ 04560fc0-90fd-11eb-16b7-27cc0762a0be
begin

function PsR2(phat,Ytrue)
    @assert(classes.(phat) == classes.(Ytrue))
    nclasses = length(classes(Ytrue))
    pnull = UnivariateFinite(classes(Ytrue),repeat([1/nclasses],nclasses))
    pnull = repeat([pnull],length(Ytrue))
    nulldev = mean(cross_entropy(pnull,Ytrue))
    modeldev = mean(cross_entropy(phat,Ytrue))
    PsR2 = 1 - modeldev/nulldev
    return(PsR2)
end

end

# ╔═╡ 15b9ff7e-90fd-11eb-2201-e9254b329c88
md"

## Model Averaging

As mentioned earlier, we will be using 4 different logistic regression models corresponding to the possible combos of transforming none/either/both of the lipid and MRI features. Below is the function used to combine the results of these models. Note that we will be using weights of 1/4 for simplicity. Ideally, the model weights for model averaging itself will also be computed through cross validation, but in this case the data is very limited and a cross validation loop is already being used for model selection (and within that, the Logistic Regression L2 regularizer $\lambda$ will also be selected by CV).

"

# ╔═╡ 1658b74c-90fd-11eb-2eef-3747749fa98d
begin

function PredAverage(preds,wts)
    @assert(length(unique(classes.(preds))) == 1)
    cats = classes(preds[1])
    probs = [pdf(p,cats) for p in preds]
    wt_probs = wts .* probs
    avg_prob = sum(wt_probs)
    avg_prob = avg_prob ./ sum(avg_prob,dims=2) #normalize in case
    return(UnivariateFinite(classes(preds[1]),avg_prob))
end

end

# ╔═╡ b1b21132-90fd-11eb-2048-53056d1e8f5d
md"

## Cross-Validation for Model Assessment

Below is the cross-validation setup for this analysis. Note that Logistic Regression hyperameter $\lambda$ for L2 regularization and the min_purity_increase for a Random Forest (200 trees) will itself be selected in an inner cross validation setup. The $\lambda$ will range from 0.001 to 1000 with 10 points spaced on the log scale, and the min_purity_increase will range from 1e-8 to 1 with 10 points also spaced on the log scale.

Every iteration the accuracy, sensitivity, specificity, and pseudo R2 will be computed for each of the 6 models. Also note that scaling is performed on the training only and this is used to scale the test set in each iteration, so as to avoid data leakage.
"

# ╔═╡ 7dea0c46-90fe-11eb-15b3-d5860ce7da6b
begin

@load LogisticClassifier pkg=MLJLinearModels
@load RandomForestClassifier pkg = DecisionTree

end

# ╔═╡ 8a98f1a0-90fe-11eb-3f93-0f23cada5615
begin

@time begin

runs = 20
AccResults=DataFrame(BGLM_nolog=fill(-9.9,runs),
                    BGLM_liplog = fill(-9.9,runs),
                    BGLM_MRIlog = fill(-9.9,runs),
                    BGLM_bothlog = fill(-9.9,runs),
                    BGLM_avg = fill(-9.9,runs),
                    RF = fill(-9.9,runs))
SensResults = deepcopy(AccResults)
SpecResults = deepcopy(AccResults)
PsR2Results = deepcopy(AccResults)

for i=1:runs

    Xtrain,Ytrain,Xtest,Ytest = train_test_split(NoMCIdata,:DX_bl)

    prep_nolog = preprocess(Xtrain,names_lipid,names_MRI)
    Xtrain_nolog_sc = MLJ.transform(prep_nolog,Xtrain)
    Xtest_nolog_sc = coerce(MLJ.transform(prep_nolog,Xtest),Continuous)

    prep_liplog = preprocess(Xtrain,names_lipid,names_MRI,log_lipid=true)
    Xtrain_liplog_sc = MLJ.transform(prep_liplog,Xtrain)
    Xtest_liplog_sc = MLJ.transform(prep_liplog,Xtest)

    prep_MRIlog = preprocess(Xtrain,names_lipid,names_MRI,log_MRI=true)
    Xtrain_MRIlog_sc = MLJ.transform(prep_MRIlog,Xtrain)
    Xtest_MRIlog_sc = MLJ.transform(prep_liplog,Xtest)

    prep_bothlog = preprocess(Xtrain,names_lipid,names_MRI,log_lipid=true,log_MRI=true)
    Xtrain_bothlog_sc = MLJ.transform(prep_bothlog,Xtrain)
    Xtest_bothlog_sc = MLJ.transform(prep_bothlog,Xtest)

    BGLM = MLJLinearModels.LogisticClassifier(penalty=:l2)
    λ_range = range(BGLM,:lambda,lower=0.001,
                    upper=1000,scale=:log);

    TunedBGLM = TunedModel(model=BGLM,resampling=CV(nfolds=5),
                            tuning=Grid(resolution=10),range=λ_range,
                            measure=cross_entropy);

    RF = MLJDecisionTreeInterface.RandomForestClassifier(n_trees=200)
    cp_range = range(RF,:min_purity_increase,lower=1e-8,upper=1,scale=:log)

    TunedRF = TunedModel(model=RF,resampling=CV(nfolds=5),
                        tuning=Grid(resolution=10),range=cp_range,
                        measure=cross_entropy)

    TunedBGLM_nolog = machine(TunedBGLM,Xtrain_nolog_sc,Ytrain)
    fit!(TunedBGLM_nolog,verbosity=0)
    Ptest_BGLM_nolog = MLJ.predict(TunedBGLM_nolog,Xtest_nolog_sc)
    Ytest_BGLM_nolog = MLJ.predict_mode(TunedBGLM_nolog,Xtest_nolog_sc)

    AccResults.BGLM_nolog[i] = accuracy(Ytest_BGLM_nolog,Ytest)
    SpecResults.BGLM_nolog[i] = specificity(Ytest_BGLM_nolog,Ytest)
    SensResults.BGLM_nolog[i] = sensitivity(Ytest_BGLM_nolog,Ytest)
    PsR2Results.BGLM_nolog[i] = PsR2(Ptest_BGLM_nolog,Ytest)

    TunedBGLM_liplog = machine(TunedBGLM,Xtrain_liplog_sc,Ytrain)
    fit!(TunedBGLM_liplog,verbosity=0)
    Ptest_BGLM_liplog = MLJ.predict(TunedBGLM_liplog,Xtest_liplog_sc)
    Ytest_BGLM_liplog = MLJ.predict_mode(TunedBGLM_liplog,Xtest_liplog_sc)

    AccResults.BGLM_liplog[i] = accuracy(Ytest_BGLM_liplog,Ytest)
    SpecResults.BGLM_liplog[i] = specificity(Ytest_BGLM_liplog,Ytest)
    SensResults.BGLM_liplog[i] = sensitivity(Ytest_BGLM_liplog,Ytest)
    PsR2Results.BGLM_liplog[i] = PsR2(Ptest_BGLM_liplog,Ytest)

    TunedBGLM_MRIlog = machine(TunedBGLM,Xtrain_MRIlog_sc,Ytrain)
    fit!(TunedBGLM_MRIlog,verbosity=0)
    Ptest_BGLM_MRIlog = MLJ.predict(TunedBGLM_MRIlog,Xtest_MRIlog_sc)
    Ytest_BGLM_MRIlog = MLJ.predict_mode(TunedBGLM_MRIlog,Xtest_MRIlog_sc)

    AccResults.BGLM_MRIlog[i] = accuracy(Ytest_BGLM_MRIlog,Ytest)
    SpecResults.BGLM_MRIlog[i] = specificity(Ytest_BGLM_MRIlog,Ytest)
    SensResults.BGLM_MRIlog[i] = sensitivity(Ytest_BGLM_MRIlog,Ytest)
    PsR2Results.BGLM_MRIlog[i] = PsR2(Ptest_BGLM_MRIlog,Ytest)

    TunedBGLM_bothlog = machine(TunedBGLM,Xtrain_bothlog_sc,Ytrain)
    fit!(TunedBGLM_bothlog,verbosity=0)
    Ptest_BGLM_bothlog = MLJ.predict(TunedBGLM_bothlog,Xtest_bothlog_sc)
    Ytest_BGLM_bothlog = MLJ.predict_mode(TunedBGLM_bothlog,Xtest_bothlog_sc)

    AccResults.BGLM_bothlog[i] = accuracy(Ytest_BGLM_bothlog,Ytest)
    SpecResults.BGLM_bothlog[i] = specificity(Ytest_BGLM_bothlog,Ytest)
    SensResults.BGLM_bothlog[i] = sensitivity(Ytest_BGLM_bothlog,Ytest)
    PsR2Results.BGLM_bothlog[i] = PsR2(Ptest_BGLM_bothlog,Ytest)

    Ptest_BGLMavg = PredAverage([Ptest_BGLM_nolog,Ptest_BGLM_liplog,Ptest_BGLM_MRIlog,
                            Ptest_BGLM_bothlog],fill(0.25,4))
    Ytest_BGLMavg = mode.(Ptest_BGLMavg)

    AccResults.BGLM_avg[i] = accuracy(Ytest_BGLMavg,Ytest)
    SpecResults.BGLM_avg[i] = specificity(Ytest_BGLMavg,Ytest)
    SensResults.BGLM_avg[i] = sensitivity(Ytest_BGLMavg,Ytest)
    PsR2Results.BGLM_avg[i] = PsR2(Ptest_BGLMavg,Ytest)

    TunedRF_nolog = machine(TunedRF,Xtrain_nolog_sc,Ytrain)
    fit!(TunedRF_nolog,verbosity=0)
    Ptest_RF = MLJ.predict(TunedRF_nolog,Xtest_nolog_sc)
    Ytest_RF = MLJ.predict_mode(TunedRF_nolog,Xtest_nolog_sc)

    AccResults.RF[i] = accuracy(Ytest_RF,Ytest)
    SpecResults.RF[i] = specificity(Ytest_RF,Ytest)
    SensResults.RF[i] = sensitivity(Ytest_RF,Ytest)
    PsR2Results.RF[i] = PsR2(Ptest_RF,Ytest)

end

end

end

# ╔═╡ f9b7f8dc-90ff-11eb-1e85-b1f58eebf547
md"

## Aggregating Results
"

# ╔═╡ 09d361ac-9100-11eb-27f7-79c209bae840
begin

function aggregate_results(results::DataFrame)
    return(combine(results,nrow,
    :BGLM_nolog=>mean=>:meanBGLM_nolog,
    :BGLM_nolog=>std=>:sdBGLM_nolog,
    :BGLM_liplog=>mean=>:meanBGLM_liplog,
    :BGLM_liplog=>std=>:sdBGLM_liplog,
    :BGLM_MRIlog=>mean=>:meanBGLM_MRIlog,
    :BGLM_MRIlog=>std=>:sdBGLM_MRIlog,
    :BGLM_bothlog=>mean=>:meanBGLM_bothlog,
    :BGLM_bothlog=>std=>:sdBGLM_bothlog,
    :BGLM_avg=>mean => :meanBGLM_avg,
    :BGLM_avg=>std=> :sdBGLM_avg,
    :RF => mean => :meanRF,
    :RF => std => :sdRF))
end

end

# ╔═╡ 151a6af6-9100-11eb-27fa-a59d7c446f0b
md"

## Accuracy

All of the logistic regression models perform very close, though the models with both lipid/MRI features log transformed perform the best ($81.2\% \pm 5.5\%$) while the model with only lipids log transformed is 2nd best ($80.2\% \pm 5.6\%$). The averaged logistic regression performs has the 3rd best accuracy, but the lowest SD ($80.1\% \pm 4.9%$)

Random Forest performs the worst ($74.4\% \pm 5.4\%$), and is likely too complex of a model for the small sample size of this dataset with LMCI excluded.

Given that lipidomic and MRI biomarkers are positively skewed, it is not surprising that log transforming appears to help. A distribution of the features closer to multivariate normal may be contributing to improving linear separability of the classes

"

# ╔═╡ 4511b8d6-9100-11eb-24bc-ad6d92627180
MeanAccRes = aggregate_results(AccResults)

# ╔═╡ 5b1720f8-9100-11eb-1a15-bbf968b34dce
md"

## Sensitivity

This follows a similar trend to accuracy. The best sensitivity comes from the model with both lipids/MRI features log transformed at $76.8\% \pm 9.1\%$ followed by the model with only lipids log transformed at $76.4\% \pm 8.9\%$. Averaging logistic models again appears to worsen the mean results but improve the SD with sensitivity at $73.8 \pm 7.7\%$.

RF has the worst sensitivity at $62.0\% \pm 9.4\%$

"

# ╔═╡ ddb90aee-9100-11eb-38c1-5d37680c533e
MeanSensRes = aggregate_results(SensResults)

# ╔═╡ 0ad91074-9102-11eb-2c4a-af34d3cc5de2
md"

## Specificity

In the case of specificity, the results are similar across the logistic regression models but the MRI log transformed model has the highest specificity at $86.2\% \pm 5.9\%$ while the averaged model has the 2nd highest specificity at $85.6\% \pm 5.1\%$. The model with both lipid/MRI log transformed had a specificity of $85.6\% \pm 4.9\%$, lipid only log transformed had a specificity of $84.9\% \pm 5.3\%$ and the model with none transformed had a specificity of $83.3\% \pm 6.0\%$

RF actually has a similar specificity as well at $85.4\% \pm 5.9\%$

"

# ╔═╡ 123a91d8-9102-11eb-0127-87e3a6d887a8
MeanSpecRes = aggregate_results(SpecResults)

# ╔═╡ 78d67062-9102-11eb-37bd-090d789fe0fa
md"

## Pseudo $R^{2}$

While the classification results are decent, the Pseudo $R^{2}$ based on binary cross entropy loss is not so stellar. This indicates that obtaining exact probabilities conditional on the MRI and lipidomic features is difficult, and the data is perhaps too noisy.

We find that the logistic model with both lipid/MRI features transformed provides the highest value at $0.371 \pm 0.125$ with the model that has lipid only coming a very close 2nd at $0.370 \pm 0.124$

Thus, instead of using the exact probability predictions on the LMCI patients in the next section (after the best model is fitted to the full data), we will bin the probability predictions into 5 categories (0-20, 20-40, etc). This accounts for a greater degree of noise in the probability predictions and is justified by the fact that the classification into 2 categories performed decently well. Thus we anticipate (or rather, make an educated assumption) that having 5 catefories/bins of LMCI based on the probabilities will also be representative. Then we can examine the correlation between the binned predicted probability and longitudinal outcome of the LMCI patients.

"

# ╔═╡ 3e4b7b1c-9109-11eb-2ca6-1fef208ae521
MeanPsR2 = aggregate_results(PsR2Results)

# ╔═╡ 5167360a-910b-11eb-0c23-cf23feecd6c9
md"

## Overall Assessment and Model Refitting

From the above, it appears that while the logistic models performed very close, log transforming both the lipidomic and MRI volumetric features is the way to go, considering it scores highest in each of the metrics (except specificity, though it comes close). We will now refit this model to the full dataset and examine the coefficient values
"

# ╔═╡ 988fbdf4-910b-11eb-0581-93555296cc33
begin

YNoMCI = NoMCIdata.DX_bl

prepNoMCI_bothlog = preprocess(NoMCIdata,names_lipid,names_MRI,log_lipid=true,log_MRI=true)

XNoMCI_bothlog_sc = MLJ.transform(prepNoMCI_bothlog,NoMCIdata)

BGLM = MLJLinearModels.LogisticClassifier(penalty=:l2)
λ_range = range(BGLM,:lambda,lower=0.001,
                    upper=1000,scale=:log)

TunedBGLM = TunedModel(model=BGLM,resampling=CV(nfolds=5),
                            tuning=Grid(resolution=10),range=λ_range,
                            measure=cross_entropy);

TunedBGLM_BothLogFinal = machine(TunedBGLM,XNoMCI_bothlog_sc,YNoMCI)
fit!(TunedBGLM_BothLogFinal,verbosity=0)

end

# ╔═╡ 4beb6a82-910d-11eb-3359-61c2add2e743
bestmod=fitted_params(TunedBGLM_BothLogFinal).best_model

# ╔═╡ 964f8e5c-9110-11eb-103d-bf52c18816fd
md"

## Interpretation

"

# ╔═╡ 9658af62-910d-11eb-3cd0-4df202793260
coefs = fitted_params(TunedBGLM_BothLogFinal).best_fitted_params.coefs

# ╔═╡ e79bd23c-910d-11eb-09f5-cb095235d389
begin

p_lipid = length(names_lipid)
p_MRI = length(names_MRI)

coef_data = DataFrame(feature_type=vcat(fill("Lipid",p_lipid),fill("MRI",p_MRI)),
	feature_name = vcat(names_lipid,names_MRI),
	coefvalue = coefs,
	abscoefvalue = abs.(coefs))

sort!(coef_data,:abscoefvalue,rev=true)
first(coef_data,10)

end

# ╔═╡ 49e63d30-910f-11eb-22ac-db7a4939e314
md"

According to the above, the best $\lambda = 10$ found by cross validation of fitting the logistic regression on log transformed (and scaled) lipidomic/MRI features. Additonally, in the coefficients we see that the ones corresponding to the left/right hippocampus (HIPPL/HIPPR) volume from MRI are most predictive of Alzheimer's disease. This is not surprising considering that the hippocampus is the area of the brain most involved in memory. The next most important features are those corresponding to volumes of the medial temporal lobe (TEMPMIDL/TEMPMIDR), which is also involved in cognitive/emotional processing. An important lipidomic feature which shows up in the top 10 is related to GM3, a ganglioside. Gangliosides have been shown to enhance the formation of $A\beta$ plaque: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3025365/.

Despite there being fewer MRI features in the dataset, the top 10 most important features seem to correspond to MRI rather than lipidomic. This may be due to the fact that MRI measurements are measurements directly of the brain, whereas lipidomic measurements may only indirectly correlate once they have already impacted these brain regions. More studies would be needed to determine the mechanisms involved.

"

# ╔═╡ af5dc1e6-9110-11eb-31fb-6d8e8ccb937a
md"

## Obtain Predictions

We now obtain the predictions on the LMCI patients using the fitted logistic model on the log scale

"

# ╔═╡ 7410fe86-9111-11eb-0310-b7bf4fb1a3af
XMCI_bothlog_sc = MLJ.transform(prepNoMCI_bothlog,MCIdata)

# ╔═╡ 9b9e60bc-9111-11eb-1ddf-1b5c97254034
ProbAD_MCI = pdf(MLJ.predict(TunedBGLM_BothLogFinal,XMCI_bothlog_sc),classes(YNoMCI))[:,2]

# ╔═╡ 3b6d1282-9112-11eb-1ab8-d370b62cabbb
begin
BinAD_MCI = cut(ProbAD_MCI,[0,0.2,0.4,0.6,0.8,1])
unique(BinAD_MCI)
end

# ╔═╡ 6df06150-9112-11eb-25be-2138ea528f50
begin
MCIdata.ProbAD = ProbAD_MCI
MCIdata.BinAD = BinAD_MCI;
end

# ╔═╡ d43c8c36-9112-11eb-149f-e1186b42d138
md"

Below is a histogram of the probability predictions. Most MCI patients (125) are in the 0-0.2 bin.
"

# ╔═╡ aba6879a-9112-11eb-2a4d-ef0d645c07a8
histogram(MCIdata.ProbAD)

# ╔═╡ d160d58a-9112-11eb-2207-afb811f16b0b
combine(groupby(MCIdata,:BinAD),nrow)

# ╔═╡ 51108c58-9113-11eb-2955-99f6010fb546
md"

## Save Predictions on LMCI and Model
"

# ╔═╡ 3d12f272-9113-11eb-384b-15297a11b028
CSV.write("MCIdata_withADprobs.csv",MCIdata)

# ╔═╡ 5fbf4a82-9113-11eb-19b0-f1422fe74c04
MLJ.save("AD_TunedLogistic_BothLipidMRILog_Scaled.jlso",TunedBGLM_BothLogFinal)

# ╔═╡ Cell order:
# ╠═1122232c-90f5-11eb-0ddd-8972932ead6b
# ╠═ee9049ce-90f4-11eb-0130-db0c30bb4ab6
# ╠═a312b3b0-90f9-11eb-2c9c-5f40f97dd471
# ╟─d4f179a2-90f9-11eb-0e1a-2ff102a97066
# ╠═2a4b1ffe-90fa-11eb-1e8f-8beeeb8a2085
# ╟─96c864aa-90fa-11eb-2ca9-979422c5a39e
# ╠═7bad01a2-90fb-11eb-183c-c371719efda8
# ╟─1469e9fa-90fc-11eb-1414-ff8d8148f92b
# ╠═04560fc0-90fd-11eb-16b7-27cc0762a0be
# ╟─15b9ff7e-90fd-11eb-2201-e9254b329c88
# ╠═1658b74c-90fd-11eb-2eef-3747749fa98d
# ╟─b1b21132-90fd-11eb-2048-53056d1e8f5d
# ╠═7dea0c46-90fe-11eb-15b3-d5860ce7da6b
# ╠═8a98f1a0-90fe-11eb-3f93-0f23cada5615
# ╟─f9b7f8dc-90ff-11eb-1e85-b1f58eebf547
# ╠═09d361ac-9100-11eb-27f7-79c209bae840
# ╟─151a6af6-9100-11eb-27fa-a59d7c446f0b
# ╠═4511b8d6-9100-11eb-24bc-ad6d92627180
# ╟─5b1720f8-9100-11eb-1a15-bbf968b34dce
# ╠═ddb90aee-9100-11eb-38c1-5d37680c533e
# ╟─0ad91074-9102-11eb-2c4a-af34d3cc5de2
# ╠═123a91d8-9102-11eb-0127-87e3a6d887a8
# ╟─78d67062-9102-11eb-37bd-090d789fe0fa
# ╠═3e4b7b1c-9109-11eb-2ca6-1fef208ae521
# ╟─5167360a-910b-11eb-0c23-cf23feecd6c9
# ╠═988fbdf4-910b-11eb-0581-93555296cc33
# ╠═4beb6a82-910d-11eb-3359-61c2add2e743
# ╟─964f8e5c-9110-11eb-103d-bf52c18816fd
# ╠═9658af62-910d-11eb-3cd0-4df202793260
# ╠═e79bd23c-910d-11eb-09f5-cb095235d389
# ╟─49e63d30-910f-11eb-22ac-db7a4939e314
# ╟─af5dc1e6-9110-11eb-31fb-6d8e8ccb937a
# ╠═7410fe86-9111-11eb-0310-b7bf4fb1a3af
# ╠═9b9e60bc-9111-11eb-1ddf-1b5c97254034
# ╠═3b6d1282-9112-11eb-1ab8-d370b62cabbb
# ╠═6df06150-9112-11eb-25be-2138ea528f50
# ╠═9cff05d2-9112-11eb-190d-23f1f040f5de
# ╟─d43c8c36-9112-11eb-149f-e1186b42d138
# ╠═aba6879a-9112-11eb-2a4d-ef0d645c07a8
# ╠═d160d58a-9112-11eb-2207-afb811f16b0b
# ╟─51108c58-9113-11eb-2955-99f6010fb546
# ╠═3d12f272-9113-11eb-384b-15297a11b028
# ╠═5fbf4a82-9113-11eb-19b0-f1422fe74c04

ext: jl


text:
# TODO: seed, __getitem__, copy

"""
A tuple (i.e., product) of simpler spaces

Example usage:
tuple_obj.observation_space = spaces.Tuple((spaces.Discrete(2), spaces.Discrete(3)))
"""
mutable struct TupleSpace <: AbstractSpace
    spaces::NTuple{N, AbstractSpace} where N
    shape::Int
    TupleSpace(space_array::NTuple{N, AbstractSpace}) where N = new(space_array, length(space_array))
    TupleSpace(space_array::Array{<:AbstractSpace, 1}) = new(Tuple(space_array), length(space_array))
end

sample(tuple_obj::TupleSpace) = Tuple(sample(space) for space in tuple_obj.spaces)

function contains(x, tuple_obj::TupleSpace)
    if isa(x, Array)
        x = Tuple(x)
    end
    return isa(x, Tuple) && Base.length(x) == Base.length(tuple_obj.spaces) &&
        all(part ∈ space for (space, part) in zip(tuple_obj.spaces, x))
end

Base.length(tuple_obj::TupleSpace) = length(tuple_obj.spaces)

Base.:(==)(tuple_obj::TupleSpace, other::TupleSpace) = tuple_obj.spaces == other.spaces
# Base.getindex(::Box, index...)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "bsdiff_classic"
version = v"4.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.daemonology.net/bsdiff/bsdiff-4.3.tar.gz", "18821588b2dc5bf159aa37d3bcb7b885d85ffd1e19f23a0c57a58723fea85f48"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/bsdiff-4.3
perl -i -ple '$_ = "#include <sys/types.h>\n" . $_ if $. == 31' bspatch.c
cc -O3 -lbz2 -I"${prefix}/include" bsdiff.c -o bsdiff
cc -O3 -lbz2 -I"${prefix}/include" bspatch.c -o bspatch
install bsdiff bspatch "${bindir}"
"""

# Disable Windows for now, as there are many BSD-isms in the source code
# that we don't want to bother to patch out.  Things like err.h and whatnot.
platforms = filter(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("bspatch", :bspatch),
    ExecutableProduct("bsdiff", :bsdiff)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Bzip2_jll", uuid="6e34b625-4abd-537c-b88f-471c36dfa7a0")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
export convertBool, compareCones, extremalCones, distinguishedAndIntPoint, minimalByLex, minimalByDist
export BerghA

"""
    convertBool(::AbstractVector)

Takes a column vector of boolean values and converts it to a vector of indices marked 'true'.

#Examples
```jldoctest makeSmoothWithDependencies
julia> B=[true,true,false,true]

julia> convertBool(B)
[0, 1, 3]
```
"""
function convertBool(B::AbstractVector)
    out=Int64[]
    for i in 1:size(B,1)
        if B[i]==true
           append!(out,i-1)
        end
    end
    return out
end

"""
    compareCones(::Array{Int64,1},::Array{Int64,1},::Array{Int64,2},::Array{Int64,1})

    Takes in two cones (in index vector notation), a ray matrix, and a incidence vector of distinguished rays. If the cones do not have an equal number of distinguished rays, returns the difference between the two values. Otherwise, returns the difference in the cone multiplicities.

# Examples
```jldoctest AlgA
julia> compareCones([1,2],[2,3],[1 0 0; 0 1 0; 0 0 1],[1,1,0])
1

julia> compareCones([1,2],[1,3],[1 0;1 2;1 -1],[1,1,1])
1
```
"""
function compareCones(cone1::Array{Int64,1}, cone2::Array{Int64,1}, rayMatrix::Array{Int64,2}, distinguished::Array{Int64,1})
    l=size(rayMatrix,1)
    c1=convertToIncidence(cone1,l)
    c2=convertToIncidence(cone2,l)
    # Calculate the number of non-distinguished rays
    nondist1 = size(cone1,1) - dot(c1, distinguished)
    nondist2 = size(cone2,1) - dot(c2, distinguished)
    if (nondist1 - nondist2 != 0)
        return nondist1 - nondist2
    else
        # Need to use the method for calculating multiplicity of cone
        mult1 = coneMultiplicity(coneConvert(cone1,rayMatrix))
        mult2 = coneMultiplicity(coneConvert(cone2,rayMatrix))
        return mult1 - mult2
    end
end

"""
    extremalCones(::Array{Array{Int64,1},1},::Array{Int64,2},::Array{Int64,1})

    Takes a list of vectors representing cones in a fan, a ray matrix, and a vector representing the distinguished rays as 0 or 1 values, and calculates the cones that are maximal with respect to (first) the number of non-distinguished rays and (second) the multiplicity of the cone. In Bergh's algorithm A (where this ordering is used), the input S will consist only of those cones containing at least one distinguished ray and at least one interior point.

#Examples
```jldoctest AlgA
julia> extremalCones([[1,2],[2,3],[3,4]],[1 0;1 2; 1 5; 1 8],[0,1,1,0])
[[ 3 ,  4 ]]
```
"""
function extremalCones(S::Array{Array{Int64,1},1}, rayMatrix::Array{Int64,2}, distinguished::Array{Int64,1})
    # Finds the extremal cones according to # distinguished rays and multiplicity
    # distinguished is a boolean vector whose size is equal to the number of rays
    # The i-th index is 1 if the i-th ray (in rayMatrix) is distinguished
    maxCones = [S[1]]
    for i in 2:size(S,1)
        cone = S[i]
        # Compare the cone with the first element of the maximal cone list
        comp = compareCones(cone, maxCones[1], rayMatrix, distinguished)
        if comp > 0
            maxCones = [cone]
        elseif comp == 0
            push!(maxCones, cone)
        end
    end
    return maxCones
end

"""
    distinguishedAndIntPoint(::Array{Int64,1},::Array{Int64,2},::Array{Int64,1})

    Calculates if the cone formed by a subset of rays in rayMatrix indexed by the entries of cone, and with a distinguished structure given by the incidence vector dist, both contains at least one distinguished ray and contains a proper interior point.

# Examples
```jldoctest AlgA
julia> distinguishedAndMultiplicity([1,2,4],[1 0 0; 1 2 0;2 1 3; 1 0 3],[1,0,0,0])
true
```
"""
function distinguishedAndIntPoint(cone::Array{Int64,1},rayMatrix::Array{Int64,2},dist::Array{Int64,1})
    l=size(rayMatrix,1)
    if dot(convertToIncidence(cone,l),dist) > 0 #check distinguished
        C=coneConvert(cone,rayMatrix)
        if interiorPoints(C)!=nothing #check interior point
            return true
        else
            return false
        end
    else
        return false
    end
end

"""
    minimalByLex(::Array{Array{Int64,1},1})

    Given a list of vectors of equal length, returns the minimal vector with respect to lexicographic ordering.

# Examples
```jldoctest AlgA
julia> A=[[1,1,1],[2,1,3],[0,5,4]];

julia> minimalByLex(A)
[ 0 ,  5 ,  4 ]
```
"""
function minimalByLex(A::Array{Array{Int64,1},1})
    l=size(A,1)
    minimal=A[1]
    d=size(minimal,1)
    for i in 2:l
        test=A[i]
        for j in 1:d
            if minimal[j]<test[j]
                break
            elseif minimal[j]>test[j]
                minimal=test
                break
            end
        end
    end
    return minimal
end

"""
    minimalByDist(::Array{Array{Int64,1},1},::Array{Int64,1})

    Given a list of vectors (representing rays as weighted sums of other rays) and a vector of 0's and 1's representing non-distinguished and distinguished indices, returns a vector from the list such that the sum of the entries corresponding to distinguished indices is minimized.

#Examples
```jldoctest AlgA
julia> minimalByDist([[0,1,5,7],[3,3,2,2],[8,5,3,6],[2,1,1,10]],[0,1,1,0])
[ 3 , 3 , 2 , 2 ]
```
"""
function minimalByDist(A::Array{Array{Int64,1},1},D::Array{Int64,1})
    l=size(A,1)
    minimal=A[1]
    d=size(minimal,1)
    for i in 2:l
        test=A[i]
        if dot(test,D)<dot(minimal,D)
            minimal=test
        end
    end
    return minimal
end

"""
    BerghA(F::StackyFan,D::Array{Int64,1})

Given a stacky fan F and a vector of booleans D representing the distinguished structure, returns a smooth stacky fan where the distinguished rays are independent.

The algorithm is adapted from Daniel Bergh's [paper on destackification](https://arxiv.org/abs/1409.5713). In brief, it identifies non-smooth cones containing at least one distinguished ray, finds interior points in those cones, and subdivides at those points through a series of stacky barycentric subdivisions.

# Examples
```jldoctest AlgA

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,1]);

julia> stackyWeights(BerghA(F,[1,1]))
[ 5 ,  2 ,  5 ,  10 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,1]);

julia> stackyWeights(BerghA(F,[1,0]));
[ 5 ,  5 ,  1 ,  2 ,  5 ,  10 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 2 5],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,1]);

julia> stackyWeights(BerghA(F,[0,1]))
[ 1 ,  5 ,  5 ,  2 ,  5 ,  1 ,  5 ,  10 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[4 1; 7 9],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,1]);

julia> stackyWeights(BerghA(F,[1,0]))
[ 609 ,  29 ,  1 ,  174 ,  29 ,  1740 ,  1218 ,  58 ,  145 ,  1044 ,  1044 ,  290 ,  290 ,  145 ,  406 ,  348 ,  261 ,  14616 ,  14616 ,  609 ,  3480 ,  870 ,  609 ,  174 ,  609 ,  9744 ,  14616 ,  1218 ,  58 ,  145 ,  435 ,  725 ,  1305 ,  1740 ,  6960 ,  3480 ,  870 ,  3480 ,  58464 ,  6090 ,  3480 ,  1392 ,  696 ,  1044 ,  2088 ,  261 ,  174 ,  261 ,  609 ,  406 ,  609 ,  609 ,  406 ,  609 ,  1218 ,  812 ,  1218 ,  2088 ,  261 ,  1044 ,  1160 ,  1740 ,  1740 ,  870 ,  1305 ,  1305 ,  14616 ,  10440 ,  145 ,  435 ,  609 ,  116 ,  580 ,  290 ,  580 ,  1740 ,  3480 ,  3480 ,  261 ,  522 ,  261 ,  522 ,  522 ,  1218 ,  1218 ,  1218 ,  1218 ,  2436 ,  2436 ,  4176 ,  2088 ,  3480 ,  2610 ,  29232 ,  6090 ,  1218 ,  290 ,  145 ,  290 ,  12180 ,  261 ,  522 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[323 179; 44 135],INPUT_CONES=[[0,1]]);

julia> F=addStackStructure(X,[1,1]);

julia> stackyWeights(BerghA(F,[1,1]))
[ 491602456800 ,  49160245680 ,  294961474080 ,  468192816 ,  12173013216 ,  73038079296 ,  12173013216 ,  97384105728 ,  73038079296 ,  196640982720 ,  131093988480 ,  156064272 ,  936385632 ,  468192816 ,  196640982720 ,  49160245680 ,  9832049136 ,  292152317184 ,  5899229481600 ,  3932819654400 ,  589922948160 ,  393281965440 ,  12173013216 ,  8115342144 ,  12173013216 ,  5899229481600 ,  589922948160 ,  737403685200 ,  51126655507200 ,  292152317184 ,  51126655507200 ,  5899229481600 ,  5899229481600 ,  589922948160 ,  589922948160 ,  196640982720 ,  2949614740800 ,  5899229481600 ,  294961474080 ,  589922948160 ,  196640982720 ,  589922948160 ,  393281965440 ,  936385632 ,  24346026432 ,  24346026432 ,  11798458963200 ,  1179845896320 ,  737403685200 ,  1474807370400 ,  2949614740800 ,  5899229481600 ,  589922948160 ,  11798458963200 ,  1179845896320 ,  2949614740800 ,  1474807370400 ,  5899229481600 ,  102253311014400 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 3; 4 5 6; 2 3 1],INPUT_CONES=[[0,1,2]]);

julia> F=addStackStructure(X,[1,1,1]);

julia> stackyWeights(BerghA(F,[1,1,1]))
[ 28 ,  21 ,  84 ,  28 ,  84 ,  84 ,  42 ,  84 ,  168 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0 2; 2 1 1; 5 3 9],INPUT_CONES=[[0,1,2]]);

julia> F=addStackStructure(X,[1,1,1]);

julia> stackyWeights(BerghA(F,[1,1,1]))
[ 4 ,  4 ,  8 ,  4 ,  4 ,  8 ]

julia> X=Polymake.fulton.NormalToricVariety(INPUT_RAYS=[1 0; 4 3; 1 5],INPUT_CONES=[[0,1],[1,2]]);

julia> F=addStackStructure(X,[1,1,1]);

julia> stackyWeights(BerghA(F,[1,1,1]))
[ 4 ,  34 ,  6 ,  68 ,  34 ,  68 ,  34 ,  6 ,  12 ]
```
"""
function BerghA(F::StackyFan,D::Array{Int64,1};verbose::Bool=false)
    if verbose==true
        println("==algorithm is running in verbose mode==")
        println(" ")
        println("=======")
    end
    X=deepcopy(F)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
    coneList=getCones(X.fan)
    dim=size(rayMatrix,2)
    numRays=size(rayMatrix,1)

    #check if the vector D has length equal to the number of rays in F
    if numRays != size(D,1)
        error("length of vector representing distinguished structure does not agree with number of rays in stacky fan.")
    end

    #A0: initialization
    i=0
    while(true)
        rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
        numRays=size(rayMatrix,1)
        coneList=getCones(X.fan)

       #Debugging
        #if verbose==true
            #println("Ray matrix: $rayMatrix")
            #println("Cone list: $coneList")
            #sW=stackyWeights(X)
            #println("Stacky weights: $sW")
            #println("Distinguished rays: $D")
        #end
        #coneMultiplicities=Int64[]
        #for cone in coneList
        #    C=coneConvert(cone,rayMatrix)
        #    push!(coneMultiplicities,coneMultiplicity(C))
        #end
        #if verbose==true
            #println("Cone multiplicities: $coneMultiplicities")
        #end
       #End debugging

        #A1: Find S the set of cones that contain a distinguised ray and an interior lattice point
        #Note: cones in S are 1-indexed.
        S=filter(cone->distinguishedAndIntPoint(cone,rayMatrix,D),coneList)
        # If S is empty, the program terminates.
        if S==[]
            break
        end

        #A2 - find extremal cones
        Smax=extremalCones(S,rayMatrix,D)

        #Print information on the number of extremal cones, their number of non-distinguished rays, and their multiplicity
        #The algorithm is structured to first reduce the number of non-distinguished rays in extremal cones, and then reduce the multiplicity of said cones,
            #so this information can be used to track the algorithm's progress
        if verbose==true
            Smaxcount=size(Smax,1)
            println("Number of extremal cones: $Smaxcount")
            testCone=Smax[1]
            c1=convertToIncidence(testCone,numRays)
            nonDist=size(testCone,1)-dot(c1,D)
            mult=coneMultiplicity(coneConvert(testCone,rayMatrix))
            println("Maximal non-distinguished rays and multiplicity: $nonDist, $mult")
        end

        #A2 - find interior points in Smax
        intPoints=[]
        for cone in Smax
            #C=rowMinors(rayMatrix,cone)
            C=coneConvert(cone,rayMatrix)
            coneIntPoints=interiorPoints(C)
            for point in coneIntPoints
               push!(intPoints,(point,cone)) #the point is stored as a tuple along with the cone containing it
            end
        end

        #A2 - find stacky points (in terms of coefficients) derived from interior points
        P=Array{Int64,1}[]
        for (point,cone) in intPoints
            stackyPoint=coneRayDecomposition(cone,rayMatrix,point,stackyWeights(X)) #each interior point is rewritten as a string of coefficients
                #corresponding to its representation as a sum of stacky rays
            push!(P,stackyPoint)
        end

        #A2 - find element of P such that the sum of the coefficients corresponding to distinguished rays is minimal.
            #This invariant does not produce a unique ray, so there is a degree of arbitrary selection.
        psi=minimalByDist(P,D)
        #if verbose==true
            #println("Psi: $psi")
        #end

        #A3 - perform root construction
        X=rootConstructionDistinguishedIndices(X,D,psi)

        #A3 - modify psi with respect to root construction
        for i in 1:length(psi)
            if D[i]==1 && psi[i]>0
                psi[i]=1
            end
        end

        #A5 - perform repeated stacky barycentric star subdivision with respect to psi.
        while(count(x->x>0,psi)>1)
            #A4 - perform stacky star subdivision
            # Get the indices of the non-zero coefficients in psi - this is used to determine the cone
                #containing the support of psi, which will be subdivided

            supportCone=findall(x->x!=0,psi)
            #find the stacky barycenter of that cone, which becomes the exceptional (blowup) ray
            exceptional=findStackyBarycenter(supportCone,X)

            #Performing a blowup may cause the rays defined in the StackyFan struct to be reordered - this property is inherited from
                #the Oscar/Polymake fan object. Since D (distinguished rays) and psi are lists, they must be reordered to match the
                #new order of the rays of X. This reordering is accomplished through defining dictionaries before the blowup is performed.

            code_rays = mapslices(encode, Polymake.common.primitive(X.fan.RAYS), dims=2)
            # Track the indices of distinguished rays
            D_pairs = map((x,y) -> (x,y), code_rays, D)
            D_Dict = Dict(D_pairs)
            # Track psi as a linear combination of the generators
            psiPairs = map((x,y) -> (x,y), code_rays,psi)
            psiDict = Dict(psiPairs)

            #perform the blowup
            X=stackyBlowup(X,[x-1 for x in supportCone],exceptional)

            G=gcd(exceptional) #since the blowup ray may not be primitive, it is made primitive and then assigned a stacky value so its stacky form is unchanged.
            primExcep=Polymake.common.primitive(exceptional)

            # Update the dictionaries storing fan information
            D_Dict[encode(primExcep)]=1
            psiDict[encode(primExcep)]=1

            #create new lists
            newRays=slicematrix(convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS))))
            newD=Int64[]
            newpsi=Int64[]
            for ray in newRays
                E=encode(ray)
                excepCode=encode(primExcep)
                push!(newD,D_Dict[E])
                #A4 - modify psi
                if E==excepCode
                    push!(newpsi,1)
                elseif psiDict[E]>1
                    push!(newpsi,psiDict[E]-1)
                else
                    push!(newpsi,0)
                end
            end
            psi=newpsi
            D=newD
        end
        if verbose==true
            println("=======")
        end
        i+=1
    end
    if verbose==true
        println("Number of steps: $i")
    end
    return X
end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
module MLJText

import MLJModelInterface
using ScientificTypesBase
import ScientificTypes: DefaultConvention
import CorpusLoaders
using SparseArrays
using TextAnalysis
using Statistics

const MMI = MLJModelInterface
const STB = ScientificTypesBase
const CL = CorpusLoaders

const PKG = "MLJText"          # substitute model-providing package name

const ScientificNGram{N} = NTuple{<:Any,STB.Textual}
const NGram{N} = NTuple{<:Any,<:AbstractString}

include("scitypes.jl")
include("utils.jl")
include("abstract_text_transformer.jl")
include("tfidf_transformer.jl")
include("count_transformer.jl")
include("bm25_transformer.jl")

export TfidfTransformer, BM25Transformer, CountTransformer

end # module


text:
using OrdinaryDiffEq
import OrdinaryDiffEq: solve
using GeneralizedGenerated

f1(du, u, p, t) = begin
    du[1] = p[1]*1/u[1]
end
u0 = [1.0]
p = [1.0]
t = (0.0,1.0)
du = similar(u0)
f1(du, u0, p, 0)
@show du
prob = ODEProblem(f1, u0, t, p)
solve(prob, Tsit5())
println("The old fashioned way works")

fex = :((du, u, p, t) -> begin
        du[1] = p[1]*1/u[1]
        end)

f2 = eval(fex)
prob = ODEProblem(f2, u0, t, p)
soln = solve(prob, Tsit5())
println("Works fine with eval")

f3   = mk_function(fex)

f3(du, u0, p, 0)
@show du
println("We can call f3")
f4 = (du, u, p, t) -> f3(du,  u, p, t)
prob = ODEProblem(f4, u0, t, p)
soln = solve(prob, Tsit5())


text:
function parse_contents(contents, filename)
    df = DataFrame()

    isnothing(contents) && return df

    content_type, content_string = split(contents, ',')
    decoded = base64decode(content_string)

    #read CSV file
    try
        if occursin("csv", filename)
            str = String(decoded)
            df =  CSV.read(IOBuffer(str), DataFrame)
        end
    catch
    end

    return df
end

function show_table(df, filename, last_modified)
    if !isnothing(last_modified)
        last_modified = Libc.strftime(last_modified)
    end

    if isempty(df)
        filename = "An empty or invalid CSV file was loaded."
        last_modified = "Try again."
        color = "red"
    else
        color = "green"
    end
    return html_div([
        html_div("File Loaded: $filename", style = (color = color,)),
        html_div("Last Modified: $(last_modified)", style = (color = color,)),
        dash_datatable(
            data=[Dict(pairs(NamedTuple(eachrow(df)[j]))) for j in 1:nrow(df)],
            columns=[Dict("name" =>i, "id" => i) for i in names(df)],
            page_size=10
        ),
    ])
end

function build_network(lead_time_df)
    #parse
    filter!(i -> !isnothing(i.lead_time), lead_time_df) #remove any empty rows
    transform!(lead_time_df,
        :lead_time =>
            ByRow(i -> i isa String ? eval(Meta.parse(i)) : i)
            => :lead_time
    )
    transform!(lead_time_df,
        :lead_time =>
            ByRow(j -> j isa Number ? [j] : j) #put numbers in Array for lead time sampling
            => :lead_time_arr
    )

    #initialize network
    node_names = union(lead_time_df.source,lead_time_df.destination)
    num_nodes = length(node_names)
    node_dict = Dict(n => i for (i,n) in enumerate(node_names))
    net = MetaDiGraph(num_nodes)
    set_prop!(net, :node_dictionary, node_dict)

    #add materials if provided
    materials = unique(lead_time_df.material)
    set_prop!(net, :materials, materials)

    for df in groupby(lead_time_df, [:source, :destination])
        src = node_dict[df.source[1]]
        dst = node_dict[df.destination[1]]
        if src != dst #store transportation lead times
            add_edge!(net, src, dst)
            set_prop!(net, src, dst, :lead_time, Dict(df.material .=> df.lead_time_arr))
        else #store production times
            set_prop!(net, src, :production_time, Dict(df.material .=> df.lead_time))
        end
    end

    return net
end

function build_bom!(net, bom_df)
    isempty(bom_df) && return

    #extract list of materials
    mats0 = get_prop(net, :materials)
    mats = union(bom_df.input,bom_df.output,mats0)

    #build bom
    bom = zeros(length(mats), length(mats))
    for row in eachrow(bom_df)
        src = row.input
        dst = row.output
        src_loc = findfirst(i -> i == src, mats)
        dst_loc = findfirst(i -> i == dst, mats)
        bom[src_loc,dst_loc] = -row.value
    end

    set_prop!(net, :materials, mats)
    set_prop!(net, :bill_of_materials, bom)
end

function build_demand!(net, demand_df)
    #parse values in third colun (distributions)
    transform!(demand_df,
        "demand_distribution" => #if its a string, parse it; if the parsed value is a number then put make it a singleton array
            ByRow(
                i -> i isa String ? eval(Meta.parse(i)) : i |>
                j -> j isa Number ? [j] : j
            )
            => :demand_distribution,
        "demand_frequency" =>
            ByRow(
                i -> i isa String ? eval(Meta.parse(i)) : i
            )
            => :demand_frequency
    )

    #add demand distribution and frequency for each node
    node_dict = get_prop(net, :node_dictionary)
    demand_grp = groupby(demand_df, "node") #group by node
    for df in demand_grp
        node_name = df.node[1]
        node_id = node_dict[node_name]
        set_props!(net, node_id, Dict(
            :demand_distribution => Dict(
                df.material .=> df.demand_distribution
            ),
            :demand_frequency => Dict(
                df.material .=> df.demand_frequency
            )
        ))
    end
end

function run_policy!(net, policy_df, policy_variable, policy_type, num_periods, backlog)
    #parse initial inventory
    transform!(policy_df,
        "initial_inventory" => #initial inventory
            ByRow(
                i -> i isa String && lowercase(i) == "unlimited" ? Inf :
                     i isa String ? eval(Meta.parse(i)) : i
            )
            => :initial_inventory,
        "param1" => #parameter 1
            ByRow(
                i -> i isa String ? eval(Meta.parse(i)) : i
            )
            => :param1,
        "param2" => #parameter 2
            ByRow(
                i -> i isa String ? eval(Meta.parse(i)) : i
            )
            => :param2,
        "review_period" => #review period
            ByRow(
                i -> i isa String ? eval(Meta.parse(i)) : i
            )
            => :review_period
    )

    #update materials list
    mats = get_prop(net, :materials)
    union!(mats, policy_df.material)

    #set initial inventory levels
    node_dict = get_prop(net, :node_dictionary)
    policy_grp = groupby(policy_df, "node") #group by node
    for df in policy_grp
        node_name = df.node[1] #node id is the first column
        node_id = node_dict[node_name]
        set_prop!(net, node_id, :initial_inventory, Dict(
            df.material .=> df.initial_inventory
        ))
    end

    #build dictionaries for policy
    param1 = Dict((node_dict[i[:node]],i[:material]) => i[:param1] for i in eachrow(policy_df))
    param2 = Dict((node_dict[i[:node]],i[:material]) => i[:param2] for i in eachrow(policy_df))
    review_period = Dict((node_dict[i[:node]],i[:material]) => i[:review_period] for i in eachrow(policy_df))

    #create and run simulation
    env = SupplyChainEnv(
        net,
        num_periods;
        backlog,
        evaluate_profit = false,
        capacitated_inventory = false,
        seed = rand(0:1000)
    )
    simulate_policy!(env, param1, param2; policy_variable, policy_type, review_period)

    return env
end


text:
using BinaryBuilder

name = "OpenMPI"
version = v"4.0.2"
sources = [
    ArchiveSource("https://download.open-mpi.org/release/open-mpi/v4.0/openmpi-$(version).tar.gz",
                  "662805870e86a1471e59739b0c34c6f9004e0c7a22db068562d5388ec4421904"),
]

script = raw"""
# Enter the funzone
cd ${WORKSPACE}/srcdir/openmpi-*

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --enable-shared=yes --enable-static=no --disable-mpi-fortran --without-cs-fs

# Build the library
make "${flags[@]}" -j${nproc}

# Install the library
make "${flags[@]}" install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
#platforms = supported_platforms()
platforms = filter(!Sys.iswindows, supported_platforms())

products = [
    LibraryProduct("libmpi", :libmpi)
    ExecutableProduct("mpiexec", :mpiexec)
]

dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: bzinberg/Yggdrasil

max_issues_repo_name: bzinberg/Yggdrasil

max_forks_repo_name: bzinberg/Yggdrasil


text:
include("../src/autodiff_module.jl")

module dtest
using PyPlot
PyPlot.pygui(true)
using ..autodiff
using BenchmarkTools

fcn2_a(x) = x[1] + x[2]
fcn2_b(x) = x[1] * x[2]
fcn2_c(x) = x[1] / x[2]
fcn3(x) = [fcn2_a(x), fcn2_b(x), fcn2_c(x)]

mode = autodiff.ReverseMode()

x0 = [1.0, 3.0]
g = autodiff.gradient(mode, fcn2_a, x0)

# x0 = [2.0, 3.0]
g = autodiff.gradient(mode, fcn2_b, x0)

# x0 = [0.0, 3.0]
g = autodiff.gradient(mode, fcn2_c, x0)
# gd = autodiff.gradient(autodiff.ForwardMode(), fcn2_c, x0)

# x0 = [1.0, 3.0]
# @btime jac = autodiff.jacobian(mode, fcn3, x0)
# @btime jacd = autodiff.jacobian(autodiff.ForwardMode(), fcn3, x0)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using SimpleSolvers
using Test

f = x -> 1 + x^2
d = x -> 2x

x = Float64(π)
y = f(x)
z = d(x)

obj1 = UnivariateObjective(f, x)
obj2 = UnivariateObjective(f, d, x)

@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 0

@test value(obj1, x) == value(obj2, x) == y
@test f_calls(obj1) == f_calls(obj2) == 1
@test d_calls(obj1) == d_calls(obj2) == 0

@test value!(obj1, x) == value!(obj2, x) == y
@test f_calls(obj1) == f_calls(obj2) == 2
@test d_calls(obj1) == d_calls(obj2) == 0

@test value!(obj1, x) == value!(obj2, x) == y
@test f_calls(obj1) == f_calls(obj2) == 2
@test d_calls(obj1) == d_calls(obj2) == 0

@test value!!(obj1, x) == value!!(obj2, x) == y
@test f_calls(obj1) == f_calls(obj2) == 3
@test d_calls(obj1) == d_calls(obj2) == 0

@test value!(obj1, 2x) == value!(obj2, 2x)
@test f_calls(obj1) == f_calls(obj2) == 4
@test d_calls(obj1) == d_calls(obj2) == 0

SimpleSolvers.clear!(obj1)
SimpleSolvers.clear!(obj2)

@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 0

@test derivative(obj1, x) == derivative(obj2, x) == z
@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 1

@test derivative!(obj1, x) == derivative!(obj2, x) == z
@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 2

@test derivative!(obj1, x) == derivative!(obj2, x) == z
@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 2

@test derivative!!(obj1, x) == derivative!!(obj2, x) == z
@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 3

@test derivative!(obj1, 2x) == derivative!(obj2, 2x)
@test f_calls(obj1) == f_calls(obj2) == 0
@test d_calls(obj1) == d_calls(obj2) == 4

function F(x)
    1 + sum(x.^2)
end

function G!(g, x)
    g .= 0
    for i in eachindex(x,g)
        g[i] = 2x[i]
    end
end

n = 2
x = rand(n)
f = F(x)
g = SimpleSolvers.alloc_g(x)
G!(g, x)

obj1 = MultivariateObjective(F, x)
obj2 = MultivariateObjective(F, G!, x)

@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 0

@test value(obj1, x) == value(obj2, x) == f
@test f_calls(obj1) == f_calls(obj2) == 1
@test g_calls(obj1) == g_calls(obj2) == 0

@test value!(obj1, x) == value!(obj2, x) == f
@test f_calls(obj1) == f_calls(obj2) == 2
@test g_calls(obj1) == g_calls(obj2) == 0

@test value!(obj1, x) == value!(obj2, x) == f
@test f_calls(obj1) == f_calls(obj2) == 2
@test g_calls(obj1) == g_calls(obj2) == 0

@test value!!(obj1, x) == value!!(obj2, x) == f
@test f_calls(obj1) == f_calls(obj2) == 3
@test g_calls(obj1) == g_calls(obj2) == 0

@test value!(obj1, 2x) == value!(obj2, 2x)
@test f_calls(obj1) == f_calls(obj2) == 4
@test g_calls(obj1) == g_calls(obj2) == 0

SimpleSolvers.clear!(obj1)
SimpleSolvers.clear!(obj2)

@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 0

@test gradient(obj1, x) == gradient(obj2, x) == g
@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 1

@test gradient!(obj1, x) == gradient!(obj2, x) == g
@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 2

@test gradient!(obj1, x) == gradient!(obj2, x) == g
@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 2

@test gradient!!(obj1, x) == gradient!!(obj2, x) == g
@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 3

@test gradient!(obj1, 2x) == gradient!(obj2, 2x)
@test f_calls(obj1) == f_calls(obj2) == 0
@test g_calls(obj1) == g_calls(obj2) == 4


text:
@testset "Constructors and basic properties" begin
    # numbers
    for T in (Float64, ComplexF64)
        ep = ErrorPropagator(T, N_args=1)
        ep2 = ErrorPropagator(T, N_args=1)

        @test length(ep) == 0
        @test ndims(ep) == 0
        @test isempty(ep)
        @test eltype(ep) == T
        @test capacity(ep) == 2^32 - 1
        @test BinningAnalysis.nlevels(ep) == 32
        @test ep == ep2
        @test ep ≈ ep2

        data = rand(rng, 1000)
        append!(ep, data)
        @test length(ep) == 1000
        @test !isempty(ep)
        @test ep != ep2
        @test !(ep ≈ ep2)
        append!(ep2, data)
        @test ep == ep2
        @test ep ≈ ep2

        empty!(ep)
        @test length(ep) == 0
        @test isempty(ep)
    end

    # arrays
    for T in (Float64, ComplexF64)
        ep = ErrorPropagator(zeros(T, 2, 3))
        ep2 = ErrorPropagator(zeros(T, 2, 3))

        @test length(ep) == 0
        @test ndims(ep) == 2
        @test isempty(ep)
        @test eltype(ep) == Array{T, 2}
        @test capacity(ep) == 2^32 - 1
        @test BinningAnalysis.nlevels(ep) == 32
        @test ep == ep2
        @test ep ≈ ep2

        data = [rand(rng, T, 2,3) for _ in 1:1000]
        append!(ep, data)
        @test length(ep) == 1000
        @test !isempty(ep)
        @test ep != ep2
        @test !(ep ≈ ep2)
        append!(ep2, data)
        @test ep == ep2
        @test ep ≈ ep2

        empty!(ep)
        @test length(ep) == 0
        @test isempty(ep)
    end

    # Constructor arguments
    ep = ErrorPropagator(capacity=12345)
    @test capacity(ep) == 16383
    @test_throws ArgumentError ErrorPropagator(capacity=0)
    @test_throws ArgumentError ErrorPropagator(capacity=-1)

    # Test error on overflow
    ep = ErrorPropagator(N_args=1, capacity=1)
    push!(ep, 1.0)
    @test_throws OverflowError push!(ep, 2.0)

    ep = ErrorPropagator(zeros(2,2), capacity=1)
    push!(ep, rand(rng, 2,2))
    @test_throws OverflowError push!(ep, rand(rng, 2,2))

    # time series constructor (#26)
    x = rand(rng, 16)
    ep = ErrorPropagator(x)
    @test length(ep) == 16
    @test mean(ep, 1) == mean(x)
    @test all(m ≈ [mean(x)] for m in all_means(ep))

    x = [rand(rng, 2,3) for _ in 1:8]
    ep = ErrorPropagator(x)
    @test length(ep) == 8
    @test mean(ep, 1) == mean(x)
    @test all(m ≈ [mean(x)] for m in all_means(ep))
end

@testset "Check variance for Float64 values" begin
    # NOTE
    # Due to the different (mathematically equivalent) versions of the variance
    # calculated here, the values are only approximately the same. (Float error)
    StableRNGs.seed!(rng, 123)
    xs = rand(rng, Float64, 1_000_000)
    ys = rand(rng, Float64, 1_000_000)
    ep = ErrorPropagator(Float64, N_args=2)

    # Test small set (off by one errors are large here)
    for (x, y) in zip(xs[1:10], ys[1:10]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs[1:10])
    @test var(ep, 2, 1) ≈ var(ys[1:10])
    @test varN(ep, 1, 1) ≈ var(xs[1:10])/10
    @test varN(ep, 2, 1) ≈ var(ys[1:10])/10
    @test vars(ep, 1) ≈ [var(xs[1:10]), var(ys[1:10])]
    @test covmat(ep, 1) ≈ [
        cov(xs[1:10], xs[1:10]) cov(xs[1:10], ys[1:10]);
        cov(ys[1:10], xs[1:10]) cov(ys[1:10], ys[1:10])
    ]

    # Test full set
    for (x, y) in zip(xs[11:end], ys[11:end]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs)
    @test var(ep, 2, 1) ≈ var(ys)
    @test varN(ep, 1, 1) ≈ var(xs)/1_000_000
    @test varN(ep, 2, 1) ≈ var(ys)/1_000_000
    @test vars(ep, 1) ≈ [var(xs), var(ys)]
    @test covmat(ep, 1) ≈ [
        cov(xs, xs) cov(xs, ys);
        cov(ys, xs) cov(ys, ys)
    ]

    # all_* methods
    @test isapprox(first.(all_vars(ep)), [0.08325406274030639, 0.041664887985400856, 0.02086378726697563, 0.010404153111327696, 0.005186509467804135, 0.002587078870339432, 0.001287770691588086, 0.0006406912824759581, 0.00031618954047712977, 0.00015796947195956257, 7.688354009283138e-5, 3.805137200274533e-5, 1.7766726465417992e-5, 8.513139774846135e-6, 3.718644130956683e-6, 2.0176056744203308e-6, 1.198040510186349e-6, 4.1842538345271407e-7, 4.374011807750655e-7], atol=1e-6)
    @test isapprox(first.(all_varNs(ep)), zero(first.(all_varNs(ep))), atol=1e-6)
    @test isapprox(first.(all_taus(ep)), [0.0, 0.0004547119263806909, 0.0012076667550950937, -0.00012514614301650795, -0.0016210095132663804, -0.0028079043475808807, -0.005025211089571158, -0.007449064127587446, -0.01383932219829298, -0.014224629670230415, -0.02690525630387247, -0.03170980256141187, -0.06269730070559826, -0.08092234827197647, -0.1338834575186591, -0.0960946513177322, -0.020327508851475773, -0.1410081813940125, 0.37563530711097126], atol=1e-6)
    @test isapprox(first.(all_std_errors(ep)), [0.00028853780123288246, 0.00028866897299640936, 0.0002888860485864669, 0.00028850168958018525, 0.0002880696990050605, 0.0002877264740180539, 0.00028708417626479785, 0.0002863803991823409, 0.0002845166131232874, 0.000284403843712177, 0.00028066727445067965, 0.0002792384696929228, 0.00026984153262076416, 0.0002641587292141796, 0.0002469035827929826, 0.00025933284111480696, 0.00028261169003733125, 0.00024448937559525804, 0.00038183817723178037], atol=1e-6)

    @test isapprox(tau(ep, 1), -0.1338834575186591)
    @test isapprox(std_error(ep, 1), 0.0002469035827929826)
end

@testset "Check variance for Float64 vectors" begin
    StableRNGs.seed!(rng, 123)
    xs = [rand(rng, Float64, 3) for _ in 1:1_000_000]
    ys = [rand(rng, Float64, 3) for _ in 1:1_000_000]
    ep = ErrorPropagator(zeros(Float64, 3), zeros(Float64, 3))

    # Test small set (off by one errors are large here)
    for (x, y) in zip(xs[1:10], ys[1:10]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs[1:10])
    @test var(ep, 2, 1) ≈ var(ys[1:10])
    @test varN(ep, 1, 1) ≈ var(xs[1:10])/10
    @test varN(ep, 2, 1) ≈ var(ys[1:10])/10
    @test vars(ep, 1) ≈ [var(xs[1:10]), var(ys[1:10])]
    @test covmat(ep, 1) ≈ reshape([
        [cov([xs[i][j] for i in 1:10], [xs[i][j] for i in 1:10]) for j in 1:3],
        [cov([ys[i][j] for i in 1:10], [xs[i][j] for i in 1:10]) for j in 1:3],
        [cov([xs[i][j] for i in 1:10], [ys[i][j] for i in 1:10]) for j in 1:3],
        [cov([ys[i][j] for i in 1:10], [ys[i][j] for i in 1:10]) for j in 1:3]
    ], (2, 2))

    # Test full set
    for (x, y) in zip(xs[11:end], ys[11:end]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs)
    @test var(ep, 2, 1) ≈ var(ys)
    @test varN(ep, 1, 1) ≈ var(xs)/1_000_000
    @test varN(ep, 2, 1) ≈ var(ys)/1_000_000
    @test vars(ep, 1) ≈ [var(xs), var(ys)]
    @test covmat(ep, 1) ≈ reshape([
        [cov([xs[i][j] for i in 1:1_000_000], [xs[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([ys[i][j] for i in 1:1_000_000], [xs[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([xs[i][j] for i in 1:1_000_000], [ys[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([ys[i][j] for i in 1:1_000_000], [ys[i][j] for i in 1:1_000_000]) for j in 1:3]
    ], (2, 2))

    # all_std_errors for <:AbstractArray
    @test all(isapprox.(first.(all_std_errors(ep)), Ref(zeros(3)), atol=1e-2))

    @test all(isapprox.(tau(ep, 1), [-0.06939313732460628, -0.05220987781564013, 0.13728303267264075], atol=1e-6))
    @test all(isapprox.(std_error(ep, 1), [0.0002677468110272373, 0.00027324026280709543, 0.00032601860649145024], atol=1e-6))
end

@testset "Check variance for complex values" begin
    # NOTE
    # Due to the different (mathematically equivalent) versions of the variance
    # calculated here, the values are onyl approximately the same. (Float error)
    StableRNGs.seed!(rng, 123)
    xs = rand(rng, ComplexF64, 1_000_000)
    ys = rand(rng, ComplexF64, 1_000_000)
    ep = ErrorPropagator(ComplexF64, N_args=2)

    # Test small set (off by one errors are large here)
    for (x, y) in zip(xs[1:10], ys[1:10]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs[1:10])
    @test var(ep, 2, 1) ≈ var(ys[1:10])
    @test varN(ep, 1, 1) ≈ var(xs[1:10])/10
    @test varN(ep, 2, 1) ≈ var(ys[1:10])/10
    @test vars(ep, 1) ≈ [var(xs[1:10]), var(ys[1:10])]
    @test covmat(ep, 1) ≈ [
        cov(xs[1:10], xs[1:10]) cov(xs[1:10], ys[1:10]);
        cov(ys[1:10], xs[1:10]) cov(ys[1:10], ys[1:10])
    ]

    # Test full set
    for (x, y) in zip(xs[11:end], ys[11:end]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs)
    @test var(ep, 2, 1) ≈ var(ys)
    @test varN(ep, 1, 1) ≈ var(xs)/1_000_000
    @test varN(ep, 2, 1) ≈ var(ys)/1_000_000
    @test vars(ep, 1) ≈ [var(xs), var(ys)]
    @test covmat(ep, 1) ≈ [
        cov(xs, xs) cov(xs, ys);
        cov(ys, xs) cov(ys, ys)
    ]

    # all_* methods
    @test isapprox(first.(all_vars(ep)), [0.16658475406392415, 0.08321852969182086, 0.04163861706031047, 0.020801210327666575, 0.010404232627790222, 0.005208087488829505, 0.002589371281786268, 0.0012979310682935674, 0.0006498645790018465, 0.00033049882070657066, 0.00016135937775918485, 8.244734583795488e-5, 4.190438084472614e-5, 2.1161425558080538e-5, 1.0185807653328993e-5, 5.794287254556885e-6, 2.412065927770435e-6, 2.1323503403003485e-6, 2.657209370315883e-7])
    @test isapprox(first.(all_varNs(ep)), zero(first.(all_varNs(ep))), atol=1e-6)
    @test isapprox(first.(all_taus(ep)), [0.0, -0.00044330191292812904, -9.090214423418397e-5, -0.0005254725847371189, -0.0003512686978415225, 0.0002222459643324015, -0.002596252083402084, -0.0013174794805482781, -0.0006266293031622627, 0.00792831442226083, -0.0037744232968683344, 0.007097663665263987, 0.015471138766913417, 0.02061879500256958, 0.0011876815245235317, 0.07971364057455732, -0.017349603544782877, 0.41431379449307415, -0.23414840339108273])
    @test isapprox(first.(all_std_errors(ep)), [0.0004081479560942626, 0.00040796698320285887, 0.0004081108528834315, 0.00040793342915398906, 0.0004080045613037231, 0.00040823865525271385, 0.0004070869219642424, 0.0004076098748477781, 0.00040789211844213303, 0.0004113711544095814, 0.0004066045146215672, 0.0004110346447395122, 0.0004144143648054135, 0.00041647846025109046, 0.0004086324183581892, 0.00043948027087143456, 0.0004010042333750302, 0.0005519252460123237, 0.00029761324849071505])

    @test isapprox(tau(ep, 1), 0.0011876815245235317)
    @test isapprox(std_error(ep, 1), 0.0004086324183581892)
end

@testset "Check variance for complex vectors" begin
    StableRNGs.seed!(rng, 123)
    xs = [rand(rng, ComplexF64, 3) for _ in 1:1_000_000]
    ys = [rand(rng, ComplexF64, 3) for _ in 1:1_000_000]
    ep = ErrorPropagator(zeros(ComplexF64, 3), zeros(ComplexF64, 3))

    # Test small set (off by one errors are large here)
    for (x, y) in zip(xs[1:10], ys[1:10]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs[1:10])
    @test var(ep, 2, 1) ≈ var(ys[1:10])
    @test varN(ep, 1, 1) ≈ var(xs[1:10])/10
    @test varN(ep, 2, 1) ≈ var(ys[1:10])/10
    @test vars(ep, 1) ≈ [var(xs[1:10]), var(ys[1:10])]
    @test covmat(ep, 1) ≈ reshape([
        [cov([xs[i][j] for i in 1:10], [xs[i][j] for i in 1:10]) for j in 1:3],
        [cov([ys[i][j] for i in 1:10], [xs[i][j] for i in 1:10]) for j in 1:3],
        [cov([xs[i][j] for i in 1:10], [ys[i][j] for i in 1:10]) for j in 1:3],
        [cov([ys[i][j] for i in 1:10], [ys[i][j] for i in 1:10]) for j in 1:3]
    ], (2, 2))

    # Test full set
    for (x, y) in zip(xs[11:end], ys[11:end]); push!(ep, x, y) end
    @test var(ep, 1, 1) ≈ var(xs)
    @test var(ep, 2, 1) ≈ var(ys)
    @test varN(ep, 1, 1) ≈ var(xs)/1_000_000
    @test varN(ep, 2, 1) ≈ var(ys)/1_000_000
    @test vars(ep, 1) ≈ [var(xs), var(ys)]
    @test covmat(ep, 1) ≈ reshape([
        [cov([xs[i][j] for i in 1:1_000_000], [xs[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([ys[i][j] for i in 1:1_000_000], [xs[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([xs[i][j] for i in 1:1_000_000], [ys[i][j] for i in 1:1_000_000]) for j in 1:3],
        [cov([ys[i][j] for i in 1:1_000_000], [ys[i][j] for i in 1:1_000_000]) for j in 1:3]
    ], (2, 2))

    # all_std_errors for <:AbstractArray
    @test all(isapprox.(first.(all_std_errors(ep)), Ref(zeros(3)), atol=1e-2))

    @test all(isapprox.(tau(ep, 1), [0.02509990727281808, 0.0756740457843339, 0.021586917851548426], atol=1e-6))
    @test all(isapprox.(std_error(ep, 1), [0.00041845464561767837, 0.0004382193124492061, 0.00041696433722574546], atol=1e-6))
end

@testset "Type promotion" begin
    epf = ErrorPropagator(zero(1.)) # Float64 ErrorPropagator
    epc = ErrorPropagator(zero(im)) # Float64 ErrorPropagator

    # Check that this doesn't throw (TODO: is there a better way?)
    @test (append!(epf, rand(rng, 1:10, 10000)); true)
    @test (append!(epc, rand(rng, 10000)); true)
end

@testset "Sum-type heuristic" begin
    # numbers
    @test typeof(ErrorPropagator(zero(Int64))) == ErrorPropagator{Float64, 32}
    @test typeof(ErrorPropagator(zero(ComplexF16))) == ErrorPropagator{ComplexF64, 32}

    # arrays
    @test typeof(ErrorPropagator(zeros(Int64, 2,2))) == ErrorPropagator{Matrix{Float64}, 32}
    @test typeof(ErrorPropagator(zeros(ComplexF16, 2,2))) == ErrorPropagator{Matrix{ComplexF64}, 32}
end

@testset "Indexing Bounds" begin
    ep = ErrorPropagator(zero(Float64); capacity=1)
    for func in [:var, :varN, :mean, :tau]
        # check if func(ep, 0) throws BoundsError
        # It should as level 1 is now the initial level
        @test_throws BoundsError @eval $func($ep, 1, 0)
        # Check that level 1 exists
        @test (@eval $func($ep, 1, 1); true)
        # Check that level 2 throws BoundsError
        @test_throws BoundsError @eval $func($ep, 1, 2)
    end
end

@testset "_reliable_level" begin
    ep = ErrorPropagator(N_args=1)
    # Empty Binner
    @test BinningAnalysis._reliable_level(ep) == 1
    @test isnan(std_error(ep, 1, BinningAnalysis._reliable_level(ep)))

    # One Element should still return NaN (due to 1/(n-1))
    push!(ep, rand(rng, ))
    @test BinningAnalysis._reliable_level(ep) == 1
    @test isnan(std_error(ep, 1, BinningAnalysis._reliable_level(ep)))

    # Two elements should return some value
    push!(ep, rand(rng, ))
    @test BinningAnalysis._reliable_level(ep) == 1
    @test !isnan(std_error(ep, 1, BinningAnalysis._reliable_level(ep)))

    # same behavior up to (including) 63 values (31 binned in first binned lvl)
    append!(ep, rand(rng, 61))
    @test BinningAnalysis._reliable_level(ep) == 1
    @test !isnan(std_error(ep, 1, BinningAnalysis._reliable_level(ep)))

    # at 64 or more values, the lvl should be increasing
    push!(ep, rand(rng, ))
    @test BinningAnalysis._reliable_level(ep) == 2
    @test !isnan(std_error(ep, 1, BinningAnalysis._reliable_level(ep)))
end

@testset "Cosmetics (show, print, etc.)" begin
    ep = ErrorPropagator(N_args=1);
    # empty binner
    io = IOBuffer()
    println(io, ep) # compact
    show(io, MIME"text/plain"(), ep) # full

    # compact
    l = String(take!(io))
    @test l == "ErrorPropagator{Float64,32}()\nErrorPropagator{Float64,32}\n| Count: 0"
    @test length(readlines(io)) == 0

    # filled binner
    StableRNGs.seed!(rng, 123)
    append!(ep, rand(rng, 1000))
    show(io, MIME"text/plain"(), ep)

    l = String(take!(io))
    @test l == "ErrorPropagator{Float64,32}\n| Count: 1000\n| Means: [0.49413]\n| StdErrors: [0.00952]"
    @test length(readlines(io)) == 0
    close(io);
end

@testset "Error Propagation" begin
    # These tests are taken from Jackknife.jl
    StableRNGs.seed!(rng, 123)

    g(x1, x2) = x1^2 - x2
    g(v) = v[1]^2 - v[2]
    grad_g(x1, x2) = [2x1, -1.0]
    grad_g(v) = [2v[1], -1.0]

    # derivative of identity is x -> 1
    grad_identity(v) = 1

    # Real
    # comparing to results from mean(ts), std(ts)/sqrt(10) (etc)
    ts = rand(rng, Float64, 1000)
    ts2 = rand(rng, Float64, 1000)

    ep = ErrorPropagator(ts)
    @test isapprox(mean(ep, 1), 0.4938701707368938, atol=1e-12)
    @test isapprox(varN(ep, grad_identity, 1), 8.54933681103257e-5, atol=1e-12)
    @test isapprox(std_error(ep, grad_identity, 1), 0.009246262385976601, atol=1e-12)
    # check consistency with Julia's Statistics
    @test varN(ep, grad_identity, 1) ≈ std(ts)^2.0/length(ts)
    @test std_error(ep, grad_identity, 1) ≈ std(ts)/sqrt(length(ts))

    ep = ErrorPropagator(1 ./ ts)
    @test isapprox(mean(ep, 1), 8.100918442539456, atol=1e-12)
    @test isapprox(varN(ep, grad_identity, 1), 1.2769418115479338, atol=1e-12)
    @test isapprox(std_error(ep, grad_identity, 1), 1.1300185005334797, atol=1e-12)

    ep = ErrorPropagator(ts, ts2.^2)
    @test isapprox(mean(ep, g), -0.0932215528431069, atol=1e-3)
    @test isapprox(varN(ep, grad_g, 1), 0.00017206056051197917, atol=1e-4)
    @test isapprox(std_error(ep, grad_g, 1), 0.013117185693279606, atol=1e-4)

    # Complex
    ts = rand(rng, ComplexF64, 1000)
    ts2 = rand(rng, ComplexF64, 1000)

    ep = ErrorPropagator(ts)
    @test isapprox(mean(ep, 1), 0.5179888756352957 + 0.4818028227938359im, atol=1e-12)
    @test isapprox(varN(ep, grad_identity, 1), 0.00016572004644474712, atol=1e-12)
    @test isapprox(std_error(ep, grad_identity, 1), 0.012873229837330922, atol=1e-12)
    # check consistency with Julia's Statistics
    @test varN(ep, grad_identity, 1) ≈ std(ts)^2.0/length(ts)
    @test std_error(ep, grad_identity, 1) ≈ std(ts)/sqrt(length(ts))

    ep = ErrorPropagator(1 ./ ts)
    @test isapprox(mean(ep, 1), 1.1073219254991422 - 1.0192129144004423im, atol=1e-11)
    @test isapprox(varN(ep, grad_identity, 1), 0.0020457719406840333, atol=1e-12)
    @test isapprox(std_error(ep, grad_identity, 1), 0.04523021048684201, atol=1e-12)

    ep = ErrorPropagator(ts, ts2.^2)
    @test isapprox(mean(ep, g), 0.01459184115054761 + 0.015352631230090064im, atol=1e-4)
    @test isapprox(varN(ep, grad_g, 1), 0.0007026624982701395, atol=1e-5)
    @test isapprox(std_error(ep, grad_g, 1), 0.026538016921358373, atol=1e-5)
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
function operator_pass(x, state)
    if CSTParser.isbinarycall(x) || (CSTParser.isbinarysyntax(x) && headof(x).fullspan > 0)
        opprec = CSTParser.get_prec(CSTParser.isdotted(x[2]) ? valof(x[2])[2:end] : valof(x[2]))
        if (opprec in (CSTParser.ColonOp, CSTParser.PowerOp, CSTParser.DeclarationOp, CSTParser.DotOp) && valof(x[2]) != "->") || x[2].fullspan == 0
            ensure_no_space_after(x[1], state, state.offset)
            ensure_no_space_after(x[2], state, state.offset + x[1].fullspan)
        else
            ensure_single_space_after(x[1], state, state.offset)
            ensure_single_space_after(x[2], state, state.offset + x[1].fullspan)
        end
    elseif CSTParser.iswhere(x)
        ensure_single_space_after(x[2], state, state.offset + x[1].fullspan)
        n = length(x)
        offset = state.offset + x[1].fullspan + x[2].fullspan
        for i = 3:n
            a = x[i]
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif CSTParser.iscall(x) && CSTParser.isoperator(x[2]) && valof(x[2]) == ":"
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif (headof(x) === :call && valof(x.args[1]) in ("+", "*") && length(x) > 3) || headof(x) == :comparison
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_single_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function tuple_pass(x, state)
    if headof(x) === :tuple
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            i == n && continue
            wsrange = nextind(state.text, offset + a.span):offset + a.fullspan
            ws = state.text[wsrange]
            if all(in((' ', '\t')), ws)
                if headof(a) === :COMMA && !(CSTParser.ispunctuation(x[i + 1]))
                    ensure_single_space_after(a, state, offset)
                elseif !(headof(x[i + 1]) === :parameters)
                    ensure_no_space_after(a, state, offset)
                end
            end
            offset += a.fullspan
        end
    end
end

function curly_pass(x, state)
    if headof(x) === :curly
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function call_pass(x, state)
    if headof(x) === :call && !(CSTParser.isbinarycall(x) || CSTParser.ischainedcall(x) || (CSTParser.isoperator(x.args[1]) && valof(x.args[1]) == ":"))
        CSTParser.iscolon
        if issameline(state.offset, state.offset + x.span, state.lines)
            offset = state.offset + x[1].fullspan
            n = length(x)
            for (i, a) in enumerate(x)
                i == 1 && continue
                if headof(a) === :COMMA
                    ensure_single_space_after(a, state, offset)
                elseif i != n && !(headof(x[i + 1]) === :parameters)
                    ensure_no_space_after(a, state, offset)
                end
                offset += a.fullspan
            end
        else
            # space holder for splitting calls across lines
        end
    elseif headof(x) === :kw
        if state.opts.kwarg === "none"
            ensure_no_space_after(x[1], state, state.offset)
            ensure_no_space_after(x[2], state, state.offset + x[1].fullspan)
        elseif state.opts.kwarg === "single"
            ensure_exactly_single_space_after(x[1], state, state.offset)
            ensure_exactly_single_space_after(x[2], state, state.offset + x[1].fullspan)
        end
    end
end

function forloop_pass(x, state)
    if headof(x) === :for
        offset = state.offset + x[1].fullspan
        for a in x[2]
            # convert iter = I into iter in I
            if CSTParser.isassignment(a) && (a.trivia === nothing || isempty(a.trivia))
                offset += a[1].fullspan
                push!(state.edits, Edit(offset + 1:offset + 2, "in "))
                offset += a[2].fullspan
                offset += a[3].fullspan
            else
                offset += a.fullspan
            end
        end
    end
end

# TODO: move this to CSTParser?
function str_value(x)
    if CSTParser.ispunctuation(x)
        headof(x) === :LPAREN && return "("
        headof(x) === :LBRACE && return "{"
        headof(x) === :LSQUARE && return "["
        headof(x) === :RPAREN && return ")"
        headof(x) === :RBRACE && return "}"
        headof(x) === :RSQUARE && return "]"
        headof(x) === :COMMA && return ","
        headof(x) === :SEMICOLON && return ";"
        headof(x) === :AT_SIGN && return "@"
        headof(x) === :DOT && return "."
        return ""
    elseif headof(x) === :IDENTIFIER || CSTParser.isliteral(x) || CSTParser.isoperator(x) || CSTParser.iskeyword(x)
        return CSTParser.str_value(x)
    else
        s = ""
        for a in x
            s *= str_value(a)
        end
        return s
    end
end

function doc_pass(x, state)
    return
    if headof(x) === :macrocall && headof(x[1]) === :globalrefdoc
        # Align global docstring to:
        #
        # """
        # doc
        # """
        #
        # If the doc is single quoted i.e. "doc", they will be replaced with triple quotes.
        offset = state.offset + x[1].fullspan
        doc = x[2]

        val = str_value(doc)

        # s = escape_string(strip(val, ['\n']), "\$")
        s = strip(val, ['\n'])
        ds = string("\"\"\"\n", s, "\n", "\"\"\"\n")

        # Check if docstring needs to be edited
        if length(ds) != doc.fullspan || s != val
            # Remove previous docstring
            push!(state.edits, Edit(offset + 1:offset + doc.fullspan, ""))
            # Append newly formatted docstring
            push!(state.edits, Edit(offset, ds))
        end
    end
end

function kw_pass(x, state)
    if CSTParser.iskeyword(x) &&
        headof(x) in (:ABSTRACT,
                      :BAREMODULE,
                      :CONST,
                      :ELSEIF,
                      :EXPORT,
                      :FOR,
                      :FUNCTION,
                      :GLOBAL,
                      :IF,
                      :IMPORT,
                      :LOCAL,
                      :MACRO,
                      :MODULE,
                      :MUTABLE,
                      :OUTER,
                      :PRIMITIVE,
                      :STRUCT,
                      :TYPE,
                      :USING,
                      :WHILE) || (headof(x) === :DO && x.args !== nothing && length(x.args) > 1 && x.args[2].args !== nothing && x.args[2].args[1].fullspan > 0)
        ensure_exactly_single_space_after(x, state, state.offset)
    end
end

function comments_pass(text, state)
    ts = tokenize(text)
    while !Tokenize.Lexers.eof(ts)
        t = Tokenize.Lexers.next_token(ts)
        if Tokens.kind(t) == Tokens.COMMENT
            val = Vector{UInt8}(t.val)
            if length(val) > 1 && val[2] == 0x3d # multiline
                if !(val[3] in (0x20, 0x09)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 2, " "))
                end
                if length(val) > 5 # ensure single space at
                    for i = length(val) - 2:-1:3
                        if val[i] in (0x20, 0x09, 0x0a, 0x0d)
                            continue
                        else
                            push!(state.edits, Edit(t.startbyte .+ (i + 1:length(val) - 2), " "))
                            break
                        end
                    end
                end

            elseif length(val) > 1
                t.startpos == (1, 1) && val[2] == 0x21 && continue # don't mess with the shebang
                if !(val[2] in (0x20, 0x09, 0x23)) # ensure single space at start
                    push!(state.edits, Edit(t.startbyte + 1, " "))
                end
            end
        end
    end
end

function lineends_pass(text, x, state)
    n = lastindex(text)
    io = IOBuffer(reverse(text))
    while !eof(io)
        c = read(io, Char)
        if c === '\n' && !eof(io)
            Base.peek(io) == 0x0d && read(io, Char) # crlf
            i1 = i2 = position(io)
            pc = Base.peek(io)
            while !eof(io) && pc in (0x20, 0x09)
                i2 = position(io)
                pc = read(io, UInt8)
            end
            if i1 != i2 && (y = get_expr(x, n - i1); y isa CSTParser.EXPR ?
                !(headof(y) in (:STRING, :TRIPLE_STRING, :CMD, :TRIPLE_CMD)) : true)
                push!(state.edits, Edit((n - i2) + 1:(n - i1), ""))
            end
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ErrorResponse <: SwaggerModel
    error::Any # spec type: Union{ Nothing, ErrorDetails } # spec name: error

    function ErrorResponse(;error=nothing)
        o = new()
        validate_property(ErrorResponse, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type ErrorResponse

const _property_map_ErrorResponse = Dict{Symbol,Symbol}(Symbol("error")=>Symbol("error"))
const _property_types_ErrorResponse = Dict{Symbol,String}(Symbol("error")=>"ErrorDetails")
Base.propertynames(::Type{ ErrorResponse }) = collect(keys(_property_map_ErrorResponse))
Swagger.property_type(::Type{ ErrorResponse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ErrorResponse[name]))}
Swagger.field_name(::Type{ ErrorResponse }, property_name::Symbol) =  _property_map_ErrorResponse[property_name]

function check_required(o::ErrorResponse)
    true
end

function validate_property(::Type{ ErrorResponse }, name::Symbol, val)
end


text:
export ConstantBoundary

abstract ConstantBoundary{V} <: AbstractBoundary
ConstantBoundary(V) = ConstantBoundary{V}

import Base.getindex

@generated function getindex{AT, N, CB<:ConstantBoundary, IT<:Number, IM<:Tuple}(arr::AbstractArray{AT, N}, ::Type{IM}, ::Type{CB}, x::Vararg{IT, N})
	if typeof(CB.parameters[1])==TypeVar
		v = zero(AT)
	else
		v = AT(CB.parameters[1])
	end

	quote
		$(Expr(:meta, :inline))

		@boundscheck if !checkbounds(Bool, arr, x...); return $v; end
		@inbounds return arr[$(IM.parameters...), x...]
	end
end

@inline getindex{AT, N, CB<:ConstantBoundary, IT<:Number}(arr::AbstractArray{AT, N}, ::Type{CB}, x::Vararg{IT, N}) = getindex(arr, Tuple{}, CB, x...)

@generated function inplaceadd!{AT, N, CB<:ConstantBoundary, VT, IT<:Number, IM<:Tuple}(arr::AbstractArray{AT, N}, ::Type{IM}, ::Type{CB}, v::VT, x::Vararg{IT, N})
	quote
		$(Expr(:meta, :inline))

		@boundscheck if !checkbounds(Bool, arr, x...); return end
		@inbounds inplaceadd!(arr, $(IM.parameters...), v, x...)
	end
end

@inline inplaceadd!{AT, N, CB<:ConstantBoundary, VT, IT<:Number}(arr::AbstractArray{AT, N}, ::Type{CB}, v::VT, x::Vararg{IT, N}) = inplaceadd!(arr, Tuple{}, CB, v, x...)

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
function print_header(method::Union{Brent})
    @printf "Iter     Function value      Lower bound       Upper bound       Best bound\n"
end

function Base.show(io::IO, trace::OptimizationTrace{Brent})
    @printf io "Iter     Function value      Lower bound       Upper bound       Best bound\n"
    @printf io "------   --------------      -----------       -----------       ----------\n"
    for state in trace.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState{Brent})
    @printf io "%6d   %14e    %14e    %14e      %s\n" t.iteration t.value t.metadata["x_lower"] t.metadata["x_upper"] t.metadata["best bound"]

    return
end

function print_header(method::Union{GoldenSection})
    @printf "Iter     Function value      Lower bound       Upper bound\n"
end

function Base.show(io::IO, trace::OptimizationTrace{GoldenSection})
    @printf io "Iter     Function value      Lower bound       Upper bound"
    @printf io "------   --------------      -----------       -----------"
    for state in trace.states
        show(io, state)
    end
    return
end

function Base.show(io::IO, t::OptimizationState{GoldenSection})
    @printf io "%6d   %14e    %14e    %14e\n" t.iteration t.value t.metadata["x_lower"] t.metadata["x_upper"]

    return
end

function Base.show(io::IO, r::UnivariateOptimizationResults)
    @printf io "Results of Optimization Algorithm\n"
    @printf io " * Algorithm: %s\n" summary(r)
    @printf io " * Search Interval: [%f, %f]\n" lower_bound(r) upper_bound(r)
    @printf io " * Minimizer: %e\n" minimizer(r)
    @printf io " * Minimum: %e\n" minimum(r)
    @printf io " * Iterations: %d\n" iterations(r)
    @printf io " * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(%.1e*|x|+%.1e): %s\n" rel_tol(r) abs_tol(r) converged(r)
    @printf io " * Objective Function Calls: %d" f_calls(r)
    return
end