text:
import Distributions
@testset "$(format("MGF: construction"))" begin
    let
        @test typeof(MGF{Poisson}(Poisson()))                 == MGF{Poisson}
        @test typeof(MGF{Poisson}(Poisson();logscale = true)) == MGF{Poisson}
        @test typeof(MGF(Poisson()))                          == MGF{Poisson}
        @test typeof(MGF(Poisson(),logscale=false))           == MGF{Poisson}
        @test typeof(MGF(:Poisson))                           == MGF{Poisson}
        #@test_throws UnrecognizedSymbolException MGF(:POisson)
        # test argument with constructor
        let
            tc = MGF{Poisson}(Poisson(),logscale=true)
            @test typeof(tc) == MGF{Poisson}
            @test tc.OPTION_LOG_SCALE == true
        end
        let
            tc = MGF(Poisson(),logscale=true)
            @test typeof(tc) == MGF{Poisson}
            @test tc.OPTION_LOG_SCALE == true
        end
        let
            tc = MGF(:Poisson,logscale=true)
            @test typeof(tc) == MGF{Poisson}
            @test tc.OPTION_LOG_SCALE == true
        end
    end
end

@testset "$(format("SampleMGF: construction"))" begin
    let
        @test typeof(SampleMGF()) == SampleMGF
    end
    # test constuctor with argument
    let
        tc = SampleMGF(logscale = false)
        @test typeof(tc) == SampleMGF
        @test tc.OPTION_LOG_SCALE == false
    end
end

@testset "$(format("MGF: evaluation[Poisson]"))" begin
    let
        for i = 1:5
            t = collect(1:0.05:1.5)
            tc_λ = rand() * 3
            tc = evaluate(MGF(:Poisson),t,λ=tc_λ)
            compare = Distributions.mgf.(Distributions.Poisson(tc_λ),t)
            @test check(:l2diff,tc,compare) < 1e-2
            tc_log = evaluate(MGF(:Poisson,logscale=true),t,λ=tc_λ)
            compare_log = log.(Distributions.mgf.(Distributions.Poisson(tc_λ),t))
            @test check(:l2diff,tc_log,compare_log) < 1e-1
        end
    end
end

@testset "$(format("MGF: evaluation[Gaussian]"))" begin
    let
        for i = 1:5
            t = collect(1:0.05:1.5)
            tc_μ = rand() * 5
            tc_σ = rand() * 5
            tc = evaluate(MGF(:Gaussian),t;μ = tc_μ, σ = tc_σ)
            compare = Distributions.mgf.(Distributions.Gaussian(tc_μ,tc_σ),t)
            @test check(:l2diff,tc,compare) < 1e-2
            tc_log = evaluate(MGF(:Gaussian,logscale=true),t;μ = tc_μ,σ = tc_σ)
            compare_log = log.(Distributions.mgf.(Distributions.Gaussian(tc_μ,tc_σ),t))
            @test check(:l2diff,tc_log,compare_log) < 1e-2
        end
    end
end

@testset "$(format("MGF: evaluation[Gamma]"))" begin
    let
        for i = 1:5
            tc_α = rand() * 5
            tc_θ = rand() * 5
            # ensure support
            t = collect(0.01:0.05: (1/tc_θ))
            tc = evaluate(MGF(:Gamma),t;α = tc_α, θ = tc_θ)
            compare = Distributions.mgf.(Distributions.Gamma(tc_α,tc_θ),t)
            @test check(:l2diff,tc,compare) < 1e-2
            tc_log = evaluate(MGF(:Gamma,logscale=true),t;α = tc_α,θ = tc_θ)
            compare_log = log.(Distributions.mgf.(Distributions.Gamma(tc_α,tc_θ),t))
            @test check(:l2diff,tc_log,compare_log) < 1e-2
        end
    end
end

@testset "$(format("MGF: evaluation[Bernoulli]"))" begin
    let
        for i = 1:5
            tc_p = rand()
            # ensure support
            t = collect(0.01:0.05: 0.5)
            tc = evaluate(MGF(:Bernoulli),t;p = tc_p)
            compare = Distributions.mgf.(Distributions.Bernoulli(tc_p),t)
            @test check(:l2diff,tc,compare) < 1e-2
            tc_log = evaluate(MGF(:Bernoulli,logscale=true),t;p = tc_p)
            compare_log = log.(Distributions.mgf.(Distributions.Bernoulli(tc_p),t))
            @test check(:l2diff,tc_log,compare_log) < 1e-2
        end
    end
end

@testset "$(format("MGF: evaluation[NegativeBinomial]"))" begin
    let
        for i = 1:5
            tc_p = rand()
            tc_r = rand() * 5
            t = collect(0.01:0.05: -log(tc_p))
            tc = evaluate(MGF(:NegativeBinomial),t;
                          p = tc_p, r = tc_r)
            compare = Distributions.mgf.(Distributions.NegativeBinomial(tc_r,tc_p),t)
            @test check(:l2diff,tc,compare) < 1e-2
            tc_log = evaluate(MGF(:NegativeBinomial,logscale=true),t;p = tc_p,r=tc_r)
            compare_log = log.(Distributions.mgf.(Distributions.NegativeBinomial(tc_r, tc_p),t))
            @test check(:l2diff,tc_log,compare_log) < 1e-2
        end
    end
end

@testset "$(format("SampleMGF: evaluation"))" begin
    # there is not way to make sure its correctness. we try to a few distributions
    # and show that as order increases the approximation gets more accurate.
    let
        total_passed = 0
        for i = 1:1000
            t            = collect(0:0.0001:0.001)
            sample       = rand(Distributions.Gamma(rand(1:100),rand(1:100)),5000)
            sample_mgf   = evaluate(SampleMGF(),t,data = sample,order = 15)
            real_mgf     = Distributions.mgf.(Distributions.fit_mle(Distributions.Gamma,sample) ,t)
            not_real_mgf = Distributions.mgf.(Distributions.fit_mle(Distributions.Normal,sample),t)
            if check(:l2diff,sample_mgf,real_mgf) < check(:l2diff,sample_mgf,not_real_mgf)
                total_passed = total_passed + 1
            end
        end
        @test total_passed/1000 > 0.80
    end

    let
        total_passed = 0
        for i = 1:1000
            t            = collect(0:0.01:0.1)
            sample       = rand(Distributions.Normal(rand(10:100),rand(1:10)),5000)
            sample_mgf   = evaluate(SampleMGF(),t,data = sample,order = 15)
            real_mgf     = Distributions.mgf.(Distributions.fit_mle(Distributions.Normal,sample) ,t)
            not_real_mgf = nothing
            try
                not_real_mgf = Distributions.mgf.(Distributions.fit_mle(Distributions.Gamma,sample),t)
            catch
                continue
            end
            if check(:l2diff,sample_mgf,real_mgf) < check(:l2diff,sample_mgf,not_real_mgf)
                total_passed = total_passed + 1
            end
        end
        @test total_passed/1000 > 0.80
    end
end


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

name = "RADAR5Builder"
version = v"0.0.2"

# Collection of sources required to build RADAR5Builder
sources = [
    "http://www.unige.ch/~hairer/radar5-v2.1.tar" =>
    "23b3b0152724b120b8d588b7dd1cd18428fe246e0c6fac886226009851e5138f",

]

# Bash recipe for building across all platforms
script = raw"""
cd "$WORKSPACE/srcdir/RADAR5-V2.1"

# On Windows move the library to `bin`.
if [[ ${target} == *-mingw32 ]]
then
    libdir="bin"
else
    libdir="lib"
fi

mkdir -p "${prefix}/${libdir}"

if [[ ${nbits} == 32 ]]
then
    echo "*** 32-bit BUILD ***"
    ${FC} ${LDFLAGS} -shared -fPIC -O3 -o "${prefix}/${libdir}/libradar5.${dlext}" radar5.f dc_decdel.f decsol.f contr5.f
else
    echo "*** 64-bit BUILD ***"
    ${FC} ${LDFLAGS} -shared -fPIC -fdefault-integer-8 -O3 -o "${prefix}/${libdir}/libradar5.${dlext}" radar5.f dc_decdel.f decsol.f contr5.f
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libradar5", :libradar5)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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

ext: jl


text:
# error analysis for ISDF in BSE

# %%
# loading packages
using BenchmarkTools, JLD2, FileIO, DelimitedFiles, LinearAlgebra, FFTW, Statistics, Arpack, HDF5
using BSE_k_ISDF

BLAS.set_num_threads(1)
FFTW.set_num_threads(1)

# DIAMOND EXAMPLE

# %% set up problem
example_path = "diamond/131313_20/"

N_1d = 20
N_rs = (N_1d, N_1d, N_1d)
N_core = 0
N_v = 4
N_c = 10
N_ks = (13, 13, 13)
@time prob = BSE_k_ISDF.BSEProblemExciting(N_core, N_v, N_c, N_ks, N_rs, example_path);

# %% error for different N_μ
N_k_samples = 20

# variable parameters
N_μ_vec = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 180, 220, 260, 300, 350, 400, 450, 500, 600]

u_v, u_c = prob.u_v, prob.u_c

N_sub = 30
F_vv = BSE_k_ISDF.qrcp(u_v, N_sub);
F_cc = BSE_k_ISDF.qrcp(u_c, N_sub);
F_vc = BSE_k_ISDF.qrcp(u_v, u_c, N_sub^2);

# for plotting
rc = BSE_k_ISDF.lattice_matrix(prob) * (mod.(BSE_k_ISDF.r_lattice(prob) .+ 0.375, 1.0) .- 0.375)
r_μ_vv = rc[:, F_vv.p];
r_μ_cc = rc[:, F_cc.p];
r_μ_vc = rc[:, F_vc.p];
save(example_path * "/real_space_grid_$(N_sub).jld2", "rc", rc, "r_μ_vv", r_μ_vv, "r_μ_cc", r_μ_cc, "r_μ_vc", r_μ_vc)

# compute errors in isdf

errors_Z_vv = []
errors_Z_cc = []
errors_Z_vc = []

@time for N_μ_vv in N_μ_vec
    r_μ_vv_indices = F_vv.p[1:N_μ_vv]
    ζ_vv = BSE_k_ISDF.assemble_ζ(u_v, r_μ_vv_indices)
    error_Z_vv = BSE_k_ISDF.isdf_error_estimate(u_v, ζ_vv, r_μ_vv_indices, N_k_samples)
    # error_Z_vv = BSE_k_ISDF.isdf_error(u_v, ζ_vv, r_μ_vv_indices)

    push!(errors_Z_vv, error_Z_vv)
end
@time for N_μ_cc in N_μ_vec
    r_μ_cc_indices = F_cc.p[1:N_μ_cc]
    ζ_cc = BSE_k_ISDF.assemble_ζ(u_c, r_μ_cc_indices)
    error_Z_cc = BSE_k_ISDF.isdf_error_estimate(u_c, ζ_cc, r_μ_cc_indices, N_k_samples)
    # error_Z_cc = BSE_k_ISDF.isdf_error(u_c, ζ_cc, r_μ_cc_indices)

    push!(errors_Z_cc, error_Z_cc)
end
@time for N_μ_vc in N_μ_vec
    r_μ_vc_indices = F_vc.p[1:N_μ_vc]
    ζ_vc = BSE_k_ISDF.assemble_ζ(u_v, u_c, r_μ_vc_indices)
    error_Z_vc = BSE_k_ISDF.isdf_error_estimate(u_v, u_c, ζ_vc, r_μ_vc_indices, N_k_samples)
    # error_Z_vc = BSE_k_ISDF.isdf_error(u_v, u_c, ζ_vc, r_μ_vc_indices)

    push!(errors_Z_vc, error_Z_vc)
end

save(example_path * "/errors_Z.jld2", "N_μ_vec", N_μ_vec, "errors_Z_vv", errors_Z_vv, "errors_Z_cc", errors_Z_cc, "errors_Z_vc", errors_Z_vc)

# %% compute spectra for different N_μ

# load erros in M
N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc =  load(example_path * "/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

# fixed parameters
N_iter = 200

#variable parameters
Z_tol_vec = [0.5, 0.2, 0.1, 0.05]

# compute spectra
for Z_tol in Z_tol_vec
    N_μ_vv = N_μ_vec[findfirst(errors_Z_vv .<= Z_tol)]
    N_μ_cc = N_μ_vec[findfirst(errors_Z_cc .<= Z_tol)]
    N_μ_vc = N_μ_vec[findfirst(errors_Z_vc .<= Z_tol)]
#     @show N_μ_vv, N_μ_cc, N_μ_vc
# end

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)

    H = BSE_k_ISDF.setup_H(prob, isdf)

    for direction in 1:3
        d = BSE_k_ISDF.optical_absorption_vector(prob, direction)

        α, β, U = BSE_k_ISDF.lanczos(H, normalize(d), N_iter)

        # save results
        save(example_path * "/optical_absorption_lanczos_$(Z_tol)_$(direction).jld2", "alpha", α, "beta", β, "norm_d", norm(d))
    end

    eigenvalues, eigenvectors = eigs(H, which=:SR, nev = 1, maxiter=1000)

    save(example_path * "/eigs_$(Z_tol).jld2", "eigenvalues", eigenvalues, "eigenvectors", eigenvectors)
end

# %% save spectra for different M_tol and N_iter

# load reference spectrum
Erange_reference = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT")[19:end, 1]
Ha_to_eV = 27.211396641308
Erange = Erange_reference ./ Ha_to_eV

absorption_reference = zeros(length(Erange_reference), 3)
absorption_reference[:, 1] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT")[19:end, 3]
absorption_reference[:, 2] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC22.OUT")[19:end, 3]
absorption_reference[:, 3] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC33.OUT")[19:end, 3]

save(example_path * "/optical_absorption_reference.jld2", "Erange", Erange, "absorption", absorption_reference)

# compute spectra from α and β
σ = 0.0055 # broadening used in reference computation in exciting
g = ω -> 1 / π * σ / (ω^2 + σ^2)

for N_iter in [50, 100, 150, 200]
    for Z_tol in Z_tol_vec
        absorption = zeros(length(Erange))
        for direction in 1:3
            α, β, norm_d = load(example_path * "/optical_absorption_lanczos_$(Z_tol)_$(direction).jld2", "alpha", "beta", "norm_d")

            absorption .+= BSE_k_ISDF.lanczos_optical_absorption(α, β, N_iter, g, Erange, norm_d^2 * 8 * pi^2 / (size(prob)[3] * BSE_k_ISDF.Ω0_volume(prob)))
        end
        absorption .*= 1 / 3

        save(example_path * "/optical_absorption_$(Z_tol)_$(N_iter).jld2", "Erange", Erange, "absorption", absorption)
    end
end

# %% compute errors

f = h5open(example_path * "/exciton_coefficients.h5")
eigenvalue_reference = read(f["eigvec-singlet-TDA-BAR-full"]["0001"]["evals"])[1]
close(f)
Erange, absorption_reference = load(example_path * "/optical_absorption_reference.jld2", "Erange", "absorption")

N_iter = 150
Z_tol_vec = [0.5, 0.2, 0.1, 0.05]

errors_optical_absorption = []
errors_ground_state_energy_abs = []
errors_ground_state_energy_rel = []
for Z_tol in Z_tol_vec
    absorption = load(example_path * "/optical_absorption_$(Z_tol)_$(N_iter).jld2", "absorption")

    error_optical_absorption = norm(absorption - absorption_reference[:, 1], 1) / norm(absorption_reference[:, 1], 1)
    push!(errors_optical_absorption, error_optical_absorption)

    eigenvalues = load(example_path * "/eigs_$(Z_tol).jld2", "eigenvalues")

    push!(errors_ground_state_energy_abs, abs(eigenvalues[1] - eigenvalue_reference))
    push!(errors_ground_state_energy_rel, abs(eigenvalues[1] - eigenvalue_reference) / eigenvalue_reference)
end

save(example_path * "/errors_spectrum.jld2", "Z_tol_vec", Z_tol_vec,  "errors_optical_absorption", errors_optical_absorption, "errors_ground_state_energy_abs", errors_ground_state_energy_abs, "errors_ground_state_energy_rel", errors_ground_state_energy_rel)

# %% GRAPHENE EXAMPLE

example_path = "graphene/"

N_1d = 20
N_rs = (15, 15, 50)
N_core = 0
N_v = 4
N_c = 5
N_ks = (42, 42, 1)
@time prob = BSE_k_ISDF.BSEProblemExciting(N_core, N_v, N_c, N_ks, N_rs, example_path);

# %% error for different N_μ
N_k_samples = 20

# variable parameters
N_μ_vec = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 180, 220, 260, 300, 350, 400, 450, 500, 600]

u_v, u_c = prob.u_v, prob.u_c

N_sub = 30
F_vv = BSE_k_ISDF.qrcp(u_v, N_sub);
F_cc = BSE_k_ISDF.qrcp(u_c, N_sub);
F_vc = BSE_k_ISDF.qrcp(u_v, u_c, N_sub^2);

# for plotting
rc = BSE_k_ISDF.lattice_matrix(prob) * (mod.(BSE_k_ISDF.r_lattice(prob) .+ 0.375, 1.0) .- 0.375)
r_μ_vv = rc[:, F_vv.p];
r_μ_cc = rc[:, F_cc.p];
r_μ_vc = rc[:, F_vc.p];
save(example_path * "/real_space_grid_$(N_sub).jld2", "rc", rc, "r_μ_vv", r_μ_vv, "r_μ_cc", r_μ_cc, "r_μ_vc", r_μ_vc)

# compute errors in isdf

errors_Z_vv = []
errors_Z_cc = []
errors_Z_vc = []

@time for N_μ_vv in N_μ_vec
    r_μ_vv_indices = F_vv.p[1:N_μ_vv]
    ζ_vv = BSE_k_ISDF.assemble_ζ(u_v, r_μ_vv_indices)
    error_Z_vv = BSE_k_ISDF.isdf_error_estimate(u_v, ζ_vv, r_μ_vv_indices, N_k_samples)
    # error_Z_vv = BSE_k_ISDF.isdf_error(u_v, ζ_vv, r_μ_vv_indices)

    push!(errors_Z_vv, error_Z_vv)
end
@time for N_μ_cc in N_μ_vec
    r_μ_cc_indices = F_cc.p[1:N_μ_cc]
    ζ_cc = BSE_k_ISDF.assemble_ζ(u_c, r_μ_cc_indices)
    error_Z_cc = BSE_k_ISDF.isdf_error_estimate(u_c, ζ_cc, r_μ_cc_indices, N_k_samples)
    # error_Z_cc = BSE_k_ISDF.isdf_error(u_c, ζ_cc, r_μ_cc_indices)

    push!(errors_Z_cc, error_Z_cc)
end
@time for N_μ_vc in N_μ_vec
    r_μ_vc_indices = F_vc.p[1:N_μ_vc]
    ζ_vc = BSE_k_ISDF.assemble_ζ(u_v, u_c, r_μ_vc_indices)
    error_Z_vc = BSE_k_ISDF.isdf_error_estimate(u_v, u_c, ζ_vc, r_μ_vc_indices, N_k_samples)
    # error_Z_vc = BSE_k_ISDF.isdf_error(u_v, u_c, ζ_vc, r_μ_vc_indices)

    push!(errors_Z_vc, error_Z_vc)
end

save(example_path * "/errors_Z.jld2", "N_μ_vec", N_μ_vec, "errors_Z_vv", errors_Z_vv, "errors_Z_cc", errors_Z_cc, "errors_Z_vc", errors_Z_vc)

# %% compute spectra for different N_μ

# load erros in M
N_μ_vec, errors_Z_vv, errors_Z_cc, errors_Z_vc =  load(example_path * "/errors_Z.jld2", "N_μ_vec", "errors_Z_vv", "errors_Z_cc", "errors_Z_vc")

# fixed parameters
N_iter = 200

#variable parameters
Z_tol_vec = [0.5, 0.2, 0.1, 0.05, 0.02, 0.01]

# compute spectra
for Z_tol in Z_tol_vec
    N_μ_vv = N_μ_vec[findfirst(errors_Z_vv .<= Z_tol)]
    N_μ_cc = N_μ_vec[findfirst(errors_Z_cc .<= Z_tol)]
    N_μ_vc = N_μ_vec[findfirst(errors_Z_vc .<= Z_tol)]
#     @show N_μ_vv, N_μ_cc, N_μ_vc
# end

    isdf = BSE_k_ISDF.ISDF(prob, N_μ_vv, N_μ_cc, N_μ_vc)

    H = BSE_k_ISDF.setup_H(prob, isdf)

    for direction in 1:3
        d = BSE_k_ISDF.optical_absorption_vector(prob, direction)

        α, β, U = BSE_k_ISDF.lanczos(H + 0.01 * I, normalize(d), N_iter)

        # save results
        save(example_path * "/optical_absorption_lanczos_$(Z_tol)_$(direction).jld2", "alpha", α, "beta", β, "norm_d", norm(d))
    end

    eigenvalues, eigenvectors = eigs(H, which=:SR, nev = 1, maxiter=1000)

    save(example_path * "/eigs_$(Z_tol).jld2", "eigenvalues", eigenvalues, "eigenvectors", eigenvectors)
end

# %% save spectra for different M_tol and N_iter

# load reference spectrum
Erange_reference = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT")[19:end, 1]
Ha_to_eV = 27.211396641308
Erange = Erange_reference ./ Ha_to_eV

absorption_reference = zeros(length(Erange_reference), 3)
absorption_reference[:, 1] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT")[19:end, 3]
absorption_reference[:, 2] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC22.OUT")[19:end, 3]
absorption_reference[:, 3] = readdlm(example_path * "/EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC33.OUT")[19:end, 3]

save(example_path * "/optical_absorption_reference.jld2", "Erange", Erange, "absorption", absorption_reference)

# compute spectra from α and β
σ = 0.0036 # broadening used in reference computation in exciting
g = ω -> 1 / π * σ / (ω^2 + σ^2)

for N_iter in [50, 100, 150, 200]
    for Z_tol in Z_tol_vec
        absorption = zeros(length(Erange))
        for direction in 1:2
            α, β, norm_d = load(example_path * "/optical_absorption_lanczos_$(Z_tol)_$(direction).jld2", "alpha", "beta", "norm_d")

            absorption .+= BSE_k_ISDF.lanczos_optical_absorption(α, β, N_iter, g, Erange, norm_d^2 * 8 * pi^2 / (size(prob)[3] * BSE_k_ISDF.Ω0_volume(prob)))
        end
        absorption .*= 1 / 2

        save(example_path * "/optical_absorption_$(Z_tol)_$(N_iter).jld2", "Erange", Erange, "absorption", absorption)
    end
end


text:
"""
    Dataset <: AbstractDataset

An AbstractDataset that stores a set of named columns

The columns are normally AbstractVectors stored in memory.

# Constructors
```julia
Dataset(pairs::Pair...; makeunique::Bool=false)
Dataset(pairs::AbstractVector{<:Pair}; makeunique::Bool=false)
Dataset(ds::AbstractDict)
Dataset(kwargs...)

Dataset(columns::AbstractVecOrMat, names::Union{AbstractVector, Symbol};
          makeunique::Bool=false)

Dataset(table)
Dataset(::DatasetRow)
```

# Keyword arguments

- `makeunique` : if `false` (the default), an error will be raised

(note that not all constructors support these keyword arguments)

# Details on behavior of different constructors

It is allowed to pass a vector of `Pair`s, a list of `Pair`s as positional
arguments, or a list of keyword arguments. In this case each pair is considered
to represent a column name to column value mapping and column name must be a
`Symbol` or string. Alternatively a dictionary can be passed to the constructor
in which case its entries are considered to define the column name and column
value pairs. If the dictionary is a `Dict` then column names will be sorted in
the returned `Dataset`.

In all the constructors described above column value can be a vector which is
consumed as is or an object of any other type (except `AbstractArray`). In the
latter case the passed value is automatically repeated to fill a new vector of
the appropriate length. As a particular rule values stored in a `Ref` or a
`0`-dimensional `AbstractArray` are unwrapped and treated in the same way.

It is also allowed to pass a vector of vectors or a matrix as as the first
argument. In this case the second argument must be
a vector of `Symbol`s or strings specifying column names, or the symbol `:auto`
to generate column names `x1`, `x2`, ... automatically.

If a single positional argument is passed to a `Dataset` constructor then it
is assumed to be of type that implements the
[Tables.jl](https://github.com/JuliaData/Tables.jl) interface using which the
returned `Dataset` is materialized.

Finally it is allowed to construct a `Dataset` from a `DatasetRow`.

# Notes

The `allowmissing` function is called on all columns passed to constructor before being added to the output data set.

By default an error will be raised if duplicates in column names are found. Pass
`makeunique=true` keyword argument (where supported) to accept duplicate names,
in which case they will be suffixed with `_i` (`i` starting at 1 for the first
duplicate).

If an `AbstractRange` is passed to a `Dataset` constructor as a column it is
always collected to a `Vector`. As a general rule
`AbstractRange` values are always materialized to a `Vector` by all functions in
InMemoryDatasets.jl before being stored in a `Dataset`.

`Dataset` can store only columns that use 1-based indexing. Attempting
to store a vector using non-standard indexing raises an error.

The `Dataset` type is designed to allow column types to vary and to be
dynamically changed also after it is constructed.

# Examples
```jldoctest
julia> Dataset((a=[1, 2], b=[3, 4])) # Tables.jl table constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         3
   2 │        2         4

julia> Dataset([(a=1, b=0), (a=2, b=0)]) # Tables.jl table constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset("a" => 1:2, "b" => 0) # Pair constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset([:a => 1:2, :b => 0]) # vector of Pairs constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset(Dict(:a => 1:2, :b => 0)) # dictionary constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset(a=1:2, b=0) # keyword argument constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset([[1, 2], [0, 0]], [:a, :b]) # vector of vectors constructor
2×2 Dataset
 Row │ a         b
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0

julia> Dataset([1 0; 2 0], :auto) # matrix constructor
2×2 Dataset
 Row │ x1        x2
     │ identity  identity
     │ Int64?    Int64?
─────┼────────────────────
   1 │        1         0
   2 │        2         0
```
"""
struct Dataset <: AbstractDataset
    columns::Vector{AbstractVector}
    colindex::Index
    attributes::Attributes
    # the inner constructor should not be used directly
    function Dataset(columns::Union{Vector{Any}, Vector{AbstractVector}},
                       colindex::Index; copycols::Bool=true)
        if length(columns) == length(colindex) == 0
            return new(AbstractVector[], Index(Dict{Symbol, Int}(), Symbol[], Dict{Int, Function}(), Int[], Int[], false, colindex.perm, colindex.starts, 1, false), Attributes())
        elseif length(columns) != length(colindex)
            throw(DimensionMismatch("Number of columns ($(length(columns))) and number of " *
                                    "column names ($(length(colindex))) are not equal"))
        end

        len = -1
        firstvec = -1
        for (i, col) in enumerate(columns)
            if col isa AbstractVector
                if len == -1
                    len = length(col)
                    firstvec = i
                elseif len != length(col)
                    n1 = _names(colindex)[firstvec]
                    n2 = _names(colindex)[i]
                    throw(DimensionMismatch("column :$n1 has length $len and column " *
                                            ":$n2 has length $(length(col))"))
                end
            end
        end
        len == -1 && (len = 1) # we got no vectors so make one row of scalars
        # it is not good idea to use threads when we have many rows (memory wise)
        if length(columns) > 100
            Threads.@threads for i in eachindex(columns)
              columns[i] = _preprocess_column(columns[i], len, copycols)
            end
        else
            for i in eachindex(columns)
              columns[i] = _preprocess_column(columns[i], len, copycols)
            end
        end

        for (i, col) in enumerate(columns)
            firstindex(col) != 1 && _onebased_check_error(i, col)
        end

        new(convert(Vector{AbstractVector}, columns), colindex, Attributes())
    end
end

function _preprocess_column(col::Any, len::Integer, copycols::Bool)
    if col isa AbstractRange
        return allowmissing(collect(col))
    elseif col isa AbstractVector
        if isa(col, BitVector)
            return convert(Vector{Union{Bool, Missing}}, col)
        else
            _res = allowmissing(col)
            if copycols
                _res === col ? copy(_res) : _res
            else
                _res
            end
        end
    elseif col isa Union{AbstractArray{<:Any, 0}, Ref}
        x = col[]
        return fill!(allocatecol(Union{Missing, typeof(x)}, len), x)
    elseif col isa AbstractArray
        throw(ArgumentError("adding AbstractArray other than AbstractVector " *
                            "as a column of a data set is not allowed"))
    else
        return fill!(allocatecol(Union{Missing, typeof(col)}, len), col)
    end
end

# Create Dataset
Dataset(df::Dataset) = copy(df)

# Create Dataset
function Dataset(pairs::Pair{Symbol, <:Any}...; makeunique::Bool=false,
                   )::Dataset
    colnames = [Symbol(k) for (k, v) in pairs]
    columns = Any[v for (k, v) in pairs]
    return Dataset(columns, Index(colnames, makeunique=makeunique)
                     )
end

# Create Dataset
function Dataset(pairs::Pair{<:AbstractString, <:Any}...; makeunique::Bool=false)::Dataset
    colnames = [Symbol(k) for (k, v) in pairs]
    columns = Any[v for (k, v) in pairs]
    return Dataset(columns, Index(colnames, makeunique=makeunique))
end

# Create Dataset
# this is needed as a workaround for Tables.jl dispatch
function Dataset(pairs::AbstractVector{<:Pair}; makeunique::Bool=false)
    if isempty(pairs)
        return Dataset()
    else
        if !(all(((k, v),) -> k isa Symbol, pairs) || all(((k, v),) -> k isa AbstractString, pairs))
            throw(ArgumentError("All column names must be either Symbols or strings (mixing is not allowed)"))
        end
        colnames = [Symbol(k) for (k, v) in pairs]
        columns = Any[v for (k, v) in pairs]
        return Dataset(columns, Index(colnames, makeunique=makeunique))
    end
end

# Create Dataset
function Dataset(d::AbstractDict)
    if all(k -> k isa Symbol, keys(d))
        colnames = collect(Symbol, keys(d))
    elseif all(k -> k isa AbstractString, keys(d))
        colnames = [Symbol(k) for k in keys(d)]
    else
        throw(ArgumentError("All column names must be either Symbols or strings (mixing is not allowed)"))
    end

    colindex = Index(colnames)
    columns = Any[v for v in values(d)]
    df = Dataset(columns, colindex)
    d isa Dict && select!(df, sort!(propertynames(df)))
    return df
end

# Create Dataset
function Dataset(; kwargs...)
    if isempty(kwargs)
        Dataset([], Index())
    else
        cnames = Symbol[]
        columns = Any[]
        copycols = true
        for (kw, val) in kwargs
            if kw === :copycols
                if val isa Bool
                    copycols = val
                else
                    throw(ArgumentError("the `copycols` keyword argument must be Boolean"))
                end
            elseif kw === :makeunique
                    throw(ArgumentError("the `makeunique` keyword argument is not allowed " *
                                        "in Dataset(; kwargs...) constructor"))
            else
              push!(cnames, kw)
              push!(columns, val)
            end
        end
        Dataset(columns, Index(cnames), copycols=copycols)
    end
end

# Create Dataset
function Dataset(columns::AbstractVector, cnames::AbstractVector{Symbol};
                   makeunique::Bool=false, copycols::Bool=true)::Dataset
    if !(eltype(columns) <: AbstractVector) && !all(col -> isa(col, AbstractVector), columns)
        throw(ArgumentError("columns argument must be a vector of AbstractVector objects"))
    end
    return Dataset(collect(AbstractVector, columns),
                     Index(convert(Vector{Symbol}, cnames), makeunique=makeunique),
                     copycols=copycols)
end

# Create Dataset
Dataset(columns::AbstractVector, cnames::AbstractVector{<:AbstractString};
          makeunique::Bool=false, copycols::Bool=true) =
    Dataset(columns, Symbol.(cnames), makeunique=makeunique, copycols=copycols)

# Create Dataset
Dataset(columns::AbstractVector{<:AbstractVector}, cnames::AbstractVector{Symbol};
          makeunique::Bool=false, copycols::Bool=true)::Dataset =
    Dataset(collect(AbstractVector, columns),
              Index(convert(Vector{Symbol}, cnames), makeunique=makeunique),
              copycols=copycols)

# Create Dataset
Dataset(columns::AbstractVector{<:AbstractVector}, cnames::AbstractVector{<:AbstractString};
          makeunique::Bool=false, copycols::Bool=true) =
    Dataset(columns, Symbol.(cnames); makeunique=makeunique, copycols=copycols)

# Create Dataset
function Dataset(columns::AbstractVector, cnames::Symbol; copycols::Bool=true)
    if cnames !== :auto
        throw(ArgumentError("if the first positional argument to Dataset " *
                            "constructor is a vector of vectors and the second " *
                            "positional argument is passed then the second " *
                            "argument must be a vector of column names or :auto"))
    end
    return Dataset(columns, gennames(length(columns)), copycols=copycols)
end

# Create Dataset
Dataset(columns::AbstractMatrix, cnames::AbstractVector{Symbol}; makeunique::Bool=false) =
    Dataset(AbstractVector[columns[:, i] for i in 1:size(columns, 2)], cnames,
              makeunique=makeunique, copycols=false)

# Create Dataset
Dataset(columns::AbstractMatrix, cnames::AbstractVector{<:AbstractString};
          makeunique::Bool=false) =
    Dataset(columns, Symbol.(cnames); makeunique=makeunique)

# Create Dataset
function Dataset(columns::AbstractMatrix, cnames::Symbol)
    if cnames !== :auto
        throw(ArgumentError("if the first positional argument to Dataset " *
                            "constructor is a matrix and a second " *
                            "positional argument is passed then the second " *
                            "argument must be a vector of column names or :auto"))
    end
    return Dataset(columns, gennames(size(columns, 2)), makeunique=false)
end

# Discontinued constructors

# Create Dataset
Dataset(matrix::Matrix) =
    throw(ArgumentError("`Dataset` constructor from a `Matrix` requires " *
                        "passing :auto as a second argument to automatically " *
                        "generate column names: `Dataset(matrix, :auto)`"))

# Create Dataset
Dataset(vecs::Vector{<:AbstractVector}) =
    throw(ArgumentError("`Dataset` constructor from a `Vector` of vectors requires " *
                        "passing :auto as a second argument to automatically " *
                        "generate column names: `Dataset(vecs, :auto)`"))

# Create Dataset
Dataset(column_eltypes::AbstractVector{T}, cnames::AbstractVector{Symbol},
          nrows::Integer=0; makeunique::Bool=false) where T<:Type =
    throw(ArgumentError("`Dataset` constructor with passed eltypes is " *
                        "not supported. Pass explicitly created columns to a " *
                        "`Dataset` constructor instead."))

# Create Dataset
Dataset(column_eltypes::AbstractVector{<:Type}, cnames::AbstractVector{<:AbstractString},
          nrows::Integer=0; makeunique::Bool=false) where T<:Type =
    throw(ArgumentError("`Dataset` constructor with passed eltypes is " *
                        "not supported. Pass explicitly created columns to a " *
                        "`Dataset` constructor instead."))

"""
    copy(ds::Dataset; copycols::Bool=true)

Copy data set `ds`.
If `copycols=true` (the default), return a new  `Dataset` holding
copies of column vectors in `ds`.
If `copycols=false`, return a new `Dataset` sharing column vectors with `ds`.

> This function uses `copy` rather than `deepcopy` internally, thus, it is not safe to use it when observations are mutable.
"""
function Base.copy(ds::Dataset)
    # TODO currently if the observation is mutable, copying data set doesn't protect it
    # Create Dataset
    newds = Dataset(copy(_columns(ds)), copy(index(ds)))
    setinfo!(newds, _attributes(ds).meta.info[])
    return newds
end


text:
function MultiDirectionalSimplex(arg0::Vector{Vector{jdouble}})
    return MultiDirectionalSimplex((Vector{Vector{jdouble}},), arg0)
end

function MultiDirectionalSimplex(arg0::Vector{Vector{jdouble}}, arg1::jdouble, arg2::jdouble)
    return MultiDirectionalSimplex((Vector{Vector{jdouble}}, jdouble, jdouble), arg0, arg1, arg2)
end

function MultiDirectionalSimplex(arg0::Vector{jdouble})
    return MultiDirectionalSimplex((Vector{jdouble},), arg0)
end

function MultiDirectionalSimplex(arg0::Vector{jdouble}, arg1::jdouble, arg2::jdouble)
    return MultiDirectionalSimplex((Vector{jdouble}, jdouble, jdouble), arg0, arg1, arg2)
end

function MultiDirectionalSimplex(arg0::jint)
    return MultiDirectionalSimplex((jint,), arg0)
end

function MultiDirectionalSimplex(arg0::jint, arg1::jdouble)
    return MultiDirectionalSimplex((jint, jdouble), arg0, arg1)
end

function MultiDirectionalSimplex(arg0::jint, arg1::jdouble, arg2::jdouble)
    return MultiDirectionalSimplex((jint, jdouble, jdouble), arg0, arg1, arg2)
end

function MultiDirectionalSimplex(arg0::jint, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return MultiDirectionalSimplex((jint, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function iterate(obj::MultiDirectionalSimplex, arg0::MultivariateFunction, arg1::Comparator)
    return jcall(obj, "iterate", void, (MultivariateFunction, Comparator), arg0, arg1)
end


text:
"""
    traverse_tree(shi::ShrikeIndex{T}, x::Array{T, 2}) where T -> leaf_idxs

Route data point `x` down to a leaf node each tree and return and array of
indexes of the data stored in each corresponding leaf node
"""
function traverse_tree(shi::ShrikeIndex{T}, x::AbstractArray{T, 2}) where T
    # Compute all needed projections
    proj::Array{T, 1} = reshape(transpose(x) * shi.random_vectors, :)
    # Initial node indexes (Root index for each tree)
    node_idxs = ones(Int, shi.ntrees)
    # Location of projection value for the current depth (Starts at
    # [1, shi.depth, 2*shi.depth, ... shi.depth*shi.ntrees].)
    projection_idxs::Array{Int, 1} = collect(1:shi.depth:shi.depth*shi.ntrees);
    # For each level of the trees
    @inbounds for d in 1:shi.depth
        # Determine if the projection is on the left or right of the split
        mask::Array{Bool, 1} = proj[projection_idxs] .>=  _get_splits(shi, node_idxs)
        # Move node indexes to the right child
        node_idxs .*= 2
        # Change to left child if the projection was less than the split value
        node_idxs[mask] .+= 1
        # Increase depth of projection locations by one
        projection_idxs .+= 1
    end
    # Subtract the number of non leaf nodes (Because we don't store indexes
    # at non leaf nodes)
    node_idxs .-= 2^shi.depth - 1
    leaf_idxs = map(i -> shi.indexes[node_idxs[i], i], 1:shi.ntrees)
    return leaf_idxs
end

traverse_tree(shi::ShrikeIndex{T}, x::AbstractArray{T, 1}) where T = traverse_tree(shi, reshape(x, :, 1))

"""
    _get_splits(shi::ShrikeIndex{T}, node_idx::Array{Int, 1}) where T -> splits

Given an array of node indexes of length `shi.ntrees` corresponding to the
current index in each tree, return the split values of each node.
"""
@inline _get_splits(shi::ShrikeIndex{T}, node_idx::Array{Int, 1}) where T = map(i -> shi.splits[i, node_idx[i]], 1:shi.ntrees)

"""
    leafsize(shi::ShrikeIndex) -> ls

Return the minimum number of points in the leaf nodes
"""
leafsize(shi::ShrikeIndex) = floor(Int, shi.npoints / 2^shi.depth)

"""
    getdepth(leafsize::Int, npoints::Int) -> d

Returns the depth needed for the tree to have approximately `leafsize` points in the
leaf nodes.
"""
getdepth(leafsize::Int, npoints::Int) = floor(Int, log(2, npoints) - log(2, leafsize))

mutable struct NeighborExplorer{T}
    k::Int
    heap::BinaryMaxHeap{Tuple{T, Int}}
    npoints::Int
    idx_set::Set{Int}
end

"""
    NeighborExplorer(idx::Int, k::Int) where T -> ne

Returns a NeighborExlorer type, a modified max heap for exploring neighbors of
neighbors during the creation of a k-nearest-neighbors graph

**Parameters:**

* `idx` (Int): The current datapoint or node. (Should not be included in the k nearest neighbors)
* `k` (Int): Maximum number of neighbors to store

"""
function NeighborExplorer{T}(idx::Int, k::Int) where T
    heap = BinaryMaxHeap{Tuple{T, Int}}()
    sizehint!(heap, k+2)
    idx_set = Set{Int}([idx])
    return NeighborExplorer{T}(k, heap, 0,  idx_set)
end

"""
    Base.push!(ne::NeighborSearcher{T}, idx::Int, dist::T) where T

Store the datapoint or node with index `idx` with distances `dist` from
the current datapoint in `ne`.
"""

@inline function Base.push!(ne::NeighborExplorer{T}, idx::Int, dist::T) where T
    if idx in ne.idx_set
        return
    end
    ne.npoints += 1
    push!(ne.idx_set, idx)
    push!(ne.heap, (dist, idx))
    if ne.npoints == ne.k + 1
        pop!(ne)
    end
end

"""
    Base.pop!(ne::NeighborExplorer{T}) where T -> idx, dist

Remove the node or datapoint with the largest distance from the neighbor
explorer.
"""
@inline function Base.pop!(ne::NeighborExplorer{T}) where T
    # Keep seen indexes in the index set so we don't consider them again
    maxdist = pop!(ne.heap)
    ne.npoints -= 1
    return maxdist
end

"""
    get_dists(ne::NeighborExplorer) -> neighbor_idxs

Return the current distances of `k` nearest neighbors.
"""
@inline get_dists(ne::NeighborExplorer) = [dist for (dist, idx) in ne.heap.valtree]

"""
    get_idxs(ne::NeighborExplorer) -> neighbor_idxs

Return the current indicies to `k` nearest neighbors.
"""
@inline get_idxs(ne::NeighborExplorer) = [idx for (dist, idx) in ne.heap.valtree]

function Base.show(io::IO, ne::NeighborExplorer)
    descr = "NeighborExplorer: \n    Currently holds $(ne.npoints)/$(ne.k) points \n    $(length(ne.idx_set)) points already explored "
    print(io, descr)
end

""" Helper function for RP-Tree voting """
@inline vote!(v::Array{UInt16, 1}, idx::Array{Int, 1}) = v[idx] .+= 0x0001

"""
    candidate_idxs(shi::ShrikeIndex{T}, q::Array{T, N}, k::Int; vote_cutoff=1) where {T, N} -> cand_idx

Each RP-Tree in the `ShrikeIndex` "votes" for nearest neighbor candidates. All points that receive
more than `vote_cutoff` votes are returned. By default `vote_cutoff=1` and this returns
the set union of all RP-Tree leaf nodes in the ensemble. This is the standard RP-Tree algorithm.
"""
@inline function candidate_idxs(shi::ShrikeIndex{T}, q::AbstractArray{T, N}, k::Int; vote_cutoff=1) where {T, N}
    leaf_idxs = traverse_tree(shi, q)
    # Empty array of votes
    votes = zeros(UInt16, shi.npoints)
    # Each tree votes on candiate nearest neighbors
    ThreadsX.map(idxs -> vote!(votes, idxs), leaf_idxs)
    # Find candiate points with enough votes
    vote_mask = votes .>= UInt16(vote_cutoff)
	cand_idx::Array{Int, 1} = findall(vote_mask)
    return cand_idx
end

"""
    approx_knn(shi::ShrikeIndex{T}, q::Array{T, 2}, k::Int; vote_cutoff=1) where T -> knn_idx

For a query point `q`, find the approximate `k` nearest neighbors from the data stored in the the
ShrikeIndex. The `vote_cutoff` parameter signifies how many "votes" a point needs in order to be included
in a linear search. Increasing `vote_cutoff` speeds up the algorithm but may reduce accuracy.

"""
function ann(shi::ShrikeIndex{T}, q::AbstractArray{T, N}, k::Int; vote_cutoff::Int=1) where {T, N}
	safe_vote_cutoff = min(shi.ntrees, vote_cutoff)
    cand_idx = candidate_idxs(shi, q, k, vote_cutoff=safe_vote_cutoff)
    ncand::Int = length(cand_idx)
    cand_dist = ThreadsX.map(i -> disttopoint(i, q, shi.data), cand_idx)
    sp = sortperm(cand_dist, alg=PartialQuickSort(k))
	@inbounds knn_idx = [cand_idx[sp[i]] for i=1:min(ncand,k)]
    return knn_idx
end

""" Distance between points in a set (Helper function for threading)
"""
@inline function disttopoint(i::Int, q::AbstractArray{T, N}, X::AbstractArray{T, 2}) where {T, N}
    @inbounds x = @view X[:, i]
    return sqeuclidean(x, q)
end

"""
    vote!(vs::Array{Dict{Int,Int}, 1}, idxs::Array{Int, 1})

Helper function for creating a knn-graph from a ShrikeIndex
"""
@inline function vote!(vs::Array{Dict{Int,UInt16}, 1}, idxs::Array{Int, 1})
    @inbounds for i in idxs
         for j in idxs
            v = get(vs[i], j, 0x123)
            vs[i][j] = v + 0x0001
            u = get(vs[j], i, 0x123)
            vs[j][i] = u + 0x0001
        end
    end
end

"""
    collect_votes(shi::ShrikeIndex{T}, k::Int) -> votes

Iterate through each leaf and collect predicted neighbors. Count how often
two nodes appear together in leaf nodes. (These are the votes.) We can
filter candidate neighbors by number of votes later.
"""
function collect_votes(shi::ShrikeIndex{T}, k::Int) where T
	votes = [Dict{Int,UInt16}() for i in 1:shi.npoints]
    nleafs = 2^shi.depth
    for i in 1:shi.ntrees
        for j in 1:nleafs
            vote!(votes, shi.indexes[j, i])
        end
    end
	return votes
end

"""
    build_neighbor_explorer(i::Int, votes::Dict{Int, UInt16}, shi::ShrikeIndex{T}, k::Int, vote_cutoff::Int) -> approx_nn

Linear search through candidate nearest neighbors. Store in a heap.
"""
function build_neighbor_explorer(
	i::Int,
	votes::Dict{Int, UInt16},
	shi::ShrikeIndex{T},
	k::Int,
	vote_cutoff::Int
) where T
	approx_nn = NeighborExplorer{T}(i, k)
	nodes = collect(keys(votes))
	@inbounds x = @view shi.data[:, i]
	@inbounds for node in nodes
		if votes[node] >= vote_cutoff
			y = @view shi.data[:, node]
			push!(approx_nn, node, sqeuclidean(x,y))
		end
	end
	return approx_nn
end

"""
    _allknn(shi::ShrikeIndex{T}, k::Int; vote_cutoff=1) where T -> knn

Returns a neighbor explorer for each point to allow for neighbors of neighbors exploration
and the creation of a knn-graph.

The `vote_cutoff` parameter signifies how many "votes" a point needs in order to be
included in a linear search. Increasing `vote_cutoff` speeds up the algorithm but may reduce accuracy.
"""
function _allknn(shi::ShrikeIndex{T}, k::Int, vote_cutoff::Int) where T
	votes = collect_votes(shi, k)
    knn = ThreadsX.map(
		(i,v) -> build_neighbor_explorer(i, v, shi, k, vote_cutoff),
		1:shi.npoints,
		votes
	)
    return knn
end

"""
    _allcandidates(shi::ShrikeIndex{T}, k::Int) -> allneighbors

Optimizations for the case when `vote_cutoff=1`. Use a set union of the leaf
indexes instead of a dictionary of votes.
"""
function _allcandidates(shi::ShrikeIndex{T}, k::Int) where T
	nleafs = 2^shi.depth
    allneighbors = [Set{Int64}() for i=1:shi.npoints]
	for i in 1:shi.ntrees
		for j in 1:nleafs
			ind = shi.indexes[j, i]
			for k in ind
				indset = allneighbors[k]
				for l in ind
					push!(indset, l)
				end
			end
		end
	end
	return allneighbors
end

"""
    build_neighbor_explorer(i::Int, candidates::Set{Int}, shi::ShrikeIndex{T}, k::Int) -> approx_nn

Optimizations for the case when `vote_cutoff=1`. Don't check number of votes.
Iterate though set of candidates.
"""
function build_neighbor_explorer(
	i::Int,
	candidates::Set{Int},
	shi::ShrikeIndex{T},
	k::Int
) where T
	approx_nn = NeighborExplorer{T}(i, k)
	@inbounds x = @view shi.data[:, i]
	 for node in candidates
		@inbounds y = @view shi.data[:, node]
		push!(approx_nn, node, sqeuclidean(x,y))
	end
	return approx_nn
end

"""
	_allknn(shi::ShrikeIndex{T}, k::Int) -> approx_nns

Optimizations for finding the nearest neighbors of every point in
the case `vote_cutoff=1`
"""
function _allknn(shi::ShrikeIndex{T}, k::Int) where T
	allneighbors = _allcandidates(shi, k)
	approx_nns = ThreadsX.map(
		(i,v) -> build_neighbor_explorer(i, v, shi, k),
		1:shi.npoints,
		allneighbors
	)
	return approx_nns
end

"""
    explore!(data::Array{T, 2}, approx_nn::Array{NeighborExplorer{T}, 1}) where T

Explore neighbors of neighbors to improve accuracy of the approximate
k-nearest-neigbors stored in `knn`.

**Parameters**
1. `i` is the current datapoint
2. `data` is the array of data. Each column is a data point
3. `approx_nn` is an array of `NeighborExplorer`s where `approx_nn[i]` contains
the current best approximation to the k-nearest-neigbors of point `i`.

"""
function explore(i::Int, data::AbstractArray{T}, approx_nn::Array{NeighborExplorer{T}, 1}) where T
    x = @view data[:, i]
    new_approx_nn = deepcopy(approx_nn[i])
    for j in get_idxs(approx_nn[i])
            for l in get_idxs(approx_nn[j])
                y = @view data[:, l:l]
                dist = sqeuclidean(x, y)
                push!(new_approx_nn, l, dist)
        end
    end
    return new_approx_nn
end

"""
    allknn(shi::ShrikeIndex{T}, k::Int; vote_cutoff::Int=1, ne_iters::Int=0) where T -> approxnn_array

Returns a `shi.npoints` by `k` array of approximate nearest neighbor indexes.
That is, `approxnn_array[i,:]` contains the indexes of the k nearest neighbors of
`shi.data[:, i]`.

**Parameters**

1. The `ne_iters` assigns the number of iterations of neighbor exploration to use.
Neighbor exploration is an inexpensive way to increase accuracy.

2, The `vote_cutoff` parameter signifies how many "votes" a point needs in order to be included
in a linear search. Increasing `vote_cutoff` speeds up the algorithm but may reduce accuracy.
Passing too large of a `vote_cutoff` results in the algorithm resetting `vote_cutoff` to equal the number of trees.
"""
function allknn(shi::ShrikeIndex{T}, k::Int; vote_cutoff::Int=1, ne_iters::Int=0) where T
	safe_vote_cutoff = min(shi.ntrees, vote_cutoff)
    # Search tree via neighbor explorers
	if safe_vote_cutoff == 1
		approx_nns = _allknn(shi, k)
	else
    	approx_nns = _allknn(shi, k, safe_vote_cutoff)
	end
    # Neighbor exploration
    for i in 1:ne_iters
        approx_nns = ThreadsX.map(i -> explore(i, shi.data, approx_nns), 1:shi.npoints)
    end
    # Load neighbor explorers into a matrix
    approxnn_array = zeros(Int, shi.npoints, k)
    for (i, ne) in enumerate(approx_nns)
        @inbounds approxnn_array[i, :] = get_idxs(ne)
    end
    return approxnn_array
end

"""
    add_edges!(j::Int, approx_nnj::NeighborExplorer{T}, A::AbstractArray{U, 2}) where {T, U}

Helper function for multithreading knn-graph adjacency matrix creation
"""
function add_edges!(j::Int, approx_nnj::NeighborExplorer{T}, A::AbstractArray{U, 2}) where {T, U}
     for i in get_idxs(approx_nnj)
        @inbounds A[i, j] = 1
    end
end

"""
    knngraph(shi::ShrikeIndex{T}, k::Int, vote_cutoff; vote_cutoff::Int=1, ne_iters::Int=0, gtype::G) where {T, G} -> g

Returns a graph with `shi.npoints` node and `k * shi.npoints` edges datapoints conneceted to nearest neighbors

**Parameters**

1. `shi`: random forest of the desired data
2. `k`: the desired number of nearest neighbors
3. `vote_cutoff`: signifies how many "votes" a point needs in order to be included
in a linear search through leaf nodes. Increasing `vote_cutoff` speeds up the algorithm but may reduce accuracy.
Passing too large of a `vote_cutoff` results in the algorithm resetting `vote_cutoff` to equal the number of trees.
4. `ne_iters`: assigns the number of iterations of neighbor exploration to use. Defaults to zero.
Neighbor exploration is a way to increse knn-graph accuracy.
5. `gtype` is the type of graph to construct. Defaults to `SimpleDiGraph`. `gtype=identity` returns a sparse adjacency matrix.

"""
function knngraph(
	shi::ShrikeIndex{T},
	k::Int;
	vote_cutoff::Int=1,
	ne_iters::Int=0,
	gtype::G=SimpleDiGraph
) where {T, G}
    # Search tree via neighbor explorers
	if vote_cutoff == 1
		approx_nns = _allknn(shi, k)
	else
    	approx_nns = _allknn(shi, k, vote_cutoff)
	end
    # Neighbor exploration (In parallel)
    for i in 1:ne_iters
        approx_nns = ThreadsX.map(i -> explore(i, shi.data, approx_nns), 1:shi.npoints)
    end
    # Load neighbor explorers into a sparse adj matrix (TODO: In parallel)
    A = spzeros(Int, shi.npoints, shi.npoints)
    for j in 1:shi.npoints
        add_edges!(j, approx_nns[j], A)
    end
    # Construct graph
    g = gtype(A)
end


text:
using CUDA
using CUDA.CUSPARSE
using KernelAbstractions
using LightGraphs
using LinearAlgebra
using Metis
using SparseArrays
using TimerOutputs

"""
    AbstractPreconditioner

Preconditioners for the iterative solvers mostly focused on GPUs

"""
abstract type AbstractPreconditioner end

"""
    BlockJacobiPreconditioner

Creates an object for the block-Jacobi preconditioner

* `npart::Int64`: Number of partitions or blocks
* `nJs::Int64`: Size of the blocks. For the GPUs these all have to be of equal size.
* `partitions::Vector{Vector{Int64}}``: `npart` partitions stored as lists
* `cupartitions`: `partitions` transfered to the GPU
* `Js`: Dense blocks of the block-Jacobi
* `cuJs`: `Js` transfered to the GPU
* `map`: The partitions as a mapping to construct views
* `cumap`: `cumap` transferred to the GPU`
* `part`: Partitioning as output by Metis
* `cupart`: `part` transferred to the GPU
* `P`: The sparse precondition matrix whose values are updated at each iteration
"""
mutable struct BlockJacobiPreconditioner <: AbstractPreconditioner
    npart::Int64
    nJs::Int64
    partitions::Vector{Vector{Int64}}
    cupartitions
    Js
    cuJs
    map
    cumap
    part
    cupart
    P
    id
    function BlockJacobiPreconditioner(J, npart, device=CPU())
        if isa(J, CuSparseMatrixCSR)
            J = SparseMatrixCSC(J)
        end
        m, n = size(J)
        if npart < 2
            error("Number of partitions `npart` should be at" *
                  "least 2 for partitioning in Metis")
        end
        adj = build_adjmatrix(J)
        g = Graph(adj)
        part = Metis.partition(g, npart)
        partitions = Vector{Vector{Int64}}()
        for i in 1:npart
            push!(partitions, [])
        end
        for (i,v) in enumerate(part)
            push!(partitions[v], i)
        end
        Js = Vector{Matrix{Float64}}(undef, npart)
        nJs = maximum(length.(partitions))
        id = Matrix{Float64}(I, nJs, nJs)
        for i in 1:npart
            Js[i] = Matrix{Float64}(I, nJs, nJs)
        end
        nmap = 0
        for b in partitions
            nmap += length(b)
        end
        map = Vector{Int64}(undef, nmap)
        part = Vector{Int64}(undef, nmap)
        for b in 1:npart
            for (i,el) in enumerate(partitions[b])
                map[el] = i
                part[el] = b
            end
        end
        row = Vector{Float64}()
        col = Vector{Float64}()
        nzval = Vector{Float64}()

        for b in 1:npart
            for x in partitions[b]
                for y in partitions[b]
                    push!(row, x)
                    push!(col, y)
                    push!(nzval, 1.0)
                end
            end
        end
        P = sparse(row, col, nzval)
        if isa(device, CUDADevice)
            id = CuMatrix{Float64}(I, nJs, nJs)
            cupartitions = Vector{CuVector{Int64}}(undef, npart)
            for i in 1:npart
                cupartitions[i] = CuVector{Int64}(partitions[i])
            end
            cuJs = Vector{CuMatrix{Float64}}(undef, length(partitions))
            for i in 1:length(partitions)
                cuJs[i] = CuMatrix{Float64}(I, nJs, nJs)
            end
            cumap = cu(map)
            cupart = cu(part)
            P = CuSparseMatrixCSR(P)
        else
            cuJs = nothing
            cupartitions = nothing
            cumap = nothing
            cupart = nothing
            id = nothing
        end
        return new(npart, nJs, partitions, cupartitions, Js, cuJs, map, cumap, part, cupart, P, id)
    end
end

"""
    build_adjmatrix

Build the adjacency matrix of a matrix A corresponding to the undirected graph

"""
function build_adjmatrix(A)
    rows = Int64[]
    cols = Int64[]
    vals = Float64[]
    rowsA = rowvals(A)
    m, n = size(A)
    for i = 1:n
        for j in nzrange(A, i)
            push!(rows, rowsA[j])
            push!(cols, i)
            push!(vals, 1.0)

            push!(rows, i)
            push!(cols, rowsA[j])
            push!(vals, 1.0)
        end
    end
    return sparse(rows,cols,vals,size(A,1),size(A,2))
end

"""
    fillblock_gpu

Fill the dense blocks of the preconditioner from the sparse CSC matrix arrays

"""
function fillblock_gpu!(cuJs, partition, map, rowPtr, colVal, nzVal, part, b)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(partition)
        for j in rowPtr[partition[i]]:rowPtr[partition[i]+1]-1
            if b == part[colVal[j]]
                @inbounds cuJs[map[partition[i]], map[colVal[j]]] = nzVal[j]
            end
        end
    end
    return nothing
end

"""
    fillblock_gpu

Update the values of the preconditioner matrix from the dense Jacobi blocks

"""
function fillP_gpu!(cuJs, partition, map, rowPtr, colVal, nzVal, part, b)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i in index:stride:length(partition)
        for j in rowPtr[partition[i]]:rowPtr[partition[i]+1]-1
            if b == part[colVal[j]]
                @inbounds nzVal[j] += cuJs[map[partition[i]], map[colVal[j]]]
            end
        end
    end
    return nothing
end

function _update_gpu(j_rowptr, j_colval, j_nzval, p)
    nblocks = length(p.partitions)
    for el in p.cuJs
        el .= p.id
    end
    # Fill Block Jacobi" begin
    CUDA.@sync begin
        for b in 1:nblocks
            @cuda threads=16 blocks=16 fillblock_gpu!(p.cuJs[b], p.cupartitions[b], p.cumap, j_rowptr, j_colval, j_nzval, p.cupart, b)
        end
    end
    # Invert blocks" begin
    CUDA.@sync pivot, info = CUDA.CUBLAS.getrf_batched!(p.cuJs, true)
    CUDA.@sync pivot, info, p.cuJs = CUDA.CUBLAS.getri_batched(p.cuJs, pivot)
    p.P.nzVal .= 0.0
    # Move blocks to P" begin
    CUDA.@sync begin
        for b in 1:nblocks
            @cuda threads=16 blocks=16 fillP_gpu!(p.cuJs[b], p.cupartitions[b], p.cumap, p.P.rowPtr, p.P.colVal, p.P.nzVal, p.cupart, b)
        end
    end
    return p.P
end

"""
    function update(J::CuSparseMatrixCSR, p)

Update the preconditioner `p` from the sparse Jacobian `J` in CSR format for the GPU

1) The dense blocks `cuJs` are filled from the sparse Jacobian `J`
2) To a batch inversion of the dense blocks using CUBLAS
3) Extract the preconditioner matrix `p.P` from the dense blocks `cuJs`

"""
function update(J::CuSparseMatrixCSR, p)
    _update_gpu(J.rowPtr, J.colVal, J.nzVal, p)
end
function update(J::Transpose{T, CuSparseMatrixCSR{T}}, p) where T
    Jt = CuSparseMatrixCSC(J.parent)
    _update_gpu(Jt.colPtr, Jt.rowVal, Jt.nzVal, p)
end

function _update_cpu(colptr, rowval, nzval, p)
    nblocks = length(p.partitions)
    # Fill Block Jacobi
    @inbounds for b in 1:nblocks
        for i in p.partitions[b]
            for j in colptr[i]:colptr[i+1]-1
                if b == p.part[rowval[j]]
                    p.Js[b][p.map[rowval[j]], p.map[i]] = nzval[j]
                end
            end
        end
    end
    # Invert blocks
    for b in 1:nblocks
        p.Js[b] = inv(p.Js[b])
    end
    p.P.nzval .= 0.0
    # Move blocks to P
    @inbounds for b in 1:nblocks
        for i in p.partitions[b]
            for j in p.P.colptr[i]:p.P.colptr[i+1]-1
                if b == p.part[p.P.rowval[j]]
                    p.P.nzval[j] += p.Js[b][p.map[p.P.rowval[j]], p.map[i]]
                end
            end
        end
    end
    return p.P
end

"""
    function update(J::SparseMatrixCSC, p)

Update the preconditioner `p` from the sparse Jacobian `J` in CSC format for the CPU

Note that this implements the same algorithm as for the GPU and becomes very slow on CPU with growing number of blocks.

"""
function update(J::SparseMatrixCSC, p)
    _update_cpu(J.colptr, J.rowval, J.nzval, p)
end
function update(J::Transpose{T, SparseMatrixCSC{T, I}}, p) where {T, I}
    ix, jx, zx = findnz(J.parent)
    update(sparse(jx, ix, zx), p)
end

is_valid(precond::BlockJacobiPreconditioner) = _check_nan(precond.P)
_check_nan(P::SparseMatrixCSC) = !any(isnan.(P.nzval))
_check_nan(P::CuSparseMatrixCSR) = !any(isnan.(P.nzVal))

function Base.show(precond::BlockJacobiPreconditioner)
    npartitions = precond.npart
    nblock = div(size(precond.P, 1), npartitions)
    println("#partitions: $npartitions, Blocksize: n = ", nblock,
            " Mbytes = ", (nblock*nblock*npartitions*8.0)/(1024.0*1024.0))
    println("Block Jacobi block size: $(precond.nJs)")
end


text:
#=
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.
=#

function calc()
  fibs = [1, 1]
  while fibs[end] <= 4_000_000
    push!(fibs, fibs[end] + fibs[end-1])
  end
  filter(x -> iseven(x) && x < 4_000_000, fibs) |> sum |> println
end
@time calc()

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using KernelAbstractions
using Statistics
using Oceananigans.Architectures: architecture, device_event
using Oceananigans.Fields: location, AbstractDataField

instantiate(X) = X()

mask_immersed_field!(field::AbstractField, value=zero(eltype(field.grid))) = mask_immersed_field!(field, field.grid, location(field), value)
mask_immersed_field!(field, grid, loc, value) = NoneEvent()

function mask_immersed_field!(field::AbstractDataField, grid::ImmersedBoundaryGrid, loc, value)
    arch = architecture(field)
    loc = instantiate.(loc)
    return launch!(arch, grid, :xyz, _mask_immersed_field!, field, loc, grid, value; dependencies = device_event(arch))
end

@kernel function _mask_immersed_field!(field, loc, grid, value)
    i, j, k = @index(Global, NTuple)
    @inbounds field[i, j, k] = scalar_mask(i, j, k, grid, grid.immersed_boundary, loc..., value, field)
end

#####
##### mask_immersed_velocities for NonhydrostaticModel
#####

mask_immersed_velocities!(U, arch, grid) = tuple(NoneEvent())


text:
function [ tb,qb,qbs,tbv,pisfc, pib, rhou, rhow,rhb, ub,um,u ] = ...
    base( profile_method,nx,nz,dz,psurf,qsurf,q4km,ztr,temptr,ttr,tsurf,p0,cp,g,...
    rd,xk,c_v,zu,rhow,rhou,tb,tbv,qb,qbs,rhb,pib,ub,um,u)
#UNTITLED5 Summary of this function goes here
#   Detailed explanation goes here

if profile_method == 1
    #########################################
    # Define a dry and neutral environment

    # Set base state theta to 300 K and mixing ratio to 0 kg/kg

    for k = 2 : nz - 1
        tb[k] = 300.0;
        qb[k] = 0.0;
        tbv[k] = tb(k)*(1 + 0.61*qb(k));

    end

    # Assign pi
    pisfc = (psurf / p0)^xk;
    pib[2] = pisfc - g*0.5*dz/(cp*tbv(2));

    for k = 3 : nz - 1

        tbvavg[k] = 0.5*(tbv(k) + tbv(k-1));
        pib[k] = pib(k-1) - g*dz/(cp*tbvavg(k));

    end

    for k = 2 : nz - 1
        #  Assign mean state density at u point

        rhou[k] = (p0*pib(k)^(c_v/rd))/(rd*tbv(k));

        # Assign mean state density at w point

        rhow[k] = 0.5*(rhou(k) + rhou(k-1));

    end
    #Ensure there is no gradient at the boundaries
    tb[1] = tb(2);
    tb[nz] = tb(nz-1);
    pib[1] = pib(2);
    pib[nz] = pib(nz-1);
    rhou[1] = rhou(2);
    rhou[nz] = rhou(nz-1);
    rhow[2] = 1.15;
    rhow[1] = rhow(2);
    rhow[nz] = rhow(nz-1);

elseif profile_method == 2
    ## fovellwk82soundingdef
    # define a sounding on the u/scalar points

    # Assign Values to vertical grid

    for k = 2 : nz - 1

        # Assign water vapor mixing ratio
        if(zu[k] <= 4000.0)
            qb[k] = qsurf - (qsurf - q4km)*zu(k)/4000.0;
        elseif (zu[k] <= 8000.0)
            qb[k] = q4km - q4km*(zu(k)-4000.0)/4000.0;
        else
            qb[k] = 0.0;
        end

        # Assign potential temperature
        if(zu[k]<= ztr)
            tb[k] = tsurf + (ttr - tsurf)*(zu(k)/ztr)^1.25;
        else
            tb[k] = ttr*exp(g*(zu(k) - ztr)/(cp*temptr));
        end

        # Assign virtual potential temperature
        tbv[k] = tb(k)*(1 + 0.61*qb(k));

    end

    # Assign pi
    pisfc = (psurf / p0)^xk;
    pib[2] = pisfc - g*0.5*dz/(cp*tbv(2));

    for k = 3 : nz - 1

        tbvavg[k] = 0.5*(tbv(k) + tbv(k-1));
        pib[k] = pib(k-1) - g*dz/(cp*tbvavg(k));

    end

    #  Assign mean state density at u point
    for k = 2 : nz-1

        rhou[k] = (p0*pib(k)^(c_v/rd))/(rd*tbv(k));

        # Assign mean water vapor saturation mixing ratio

        t[k] = tb(k)*pib(k);
        p[k] = p0*pib(k)^(cp/rd);
        qbs[k] = (380.0/p(k))*exp(17.27*(t(k)-273.0)/(t(k)-36.0));

        # Calculate and Assign RH

        rhb[k] = qb(k)/qbs(k);

    end

    #  define density at the true surface from known information
    pisfc = (psurf/p0)^xk;
    rhow[2] = p0*(pisfc^(c_v/rd))/(rd*tsurf);

    #  define density at the other w points by interpolating from u/scalar points
    for k=3:nz-1
        rhow[k] = 0.5*(rhou(k) + rhou(k-1));
    end

elseif profile_method == 3

        #  define the Weisman Klemp sounding on the u/scalar points

        for k=2 : nz-1
            #  assign potential temperatures and relative humidities
            if zu[k]<= ztr
                tb[k] = tsurf + (ttr-tsurf)*(zu[k]/ztr)^1.25;
                rhb[k] = 1 - 0.02*(zu[k]/ztr)^1.25;
            elseif zu[k] > ztr
                tb[k] = ttr*exp(g*(zu[k]-ztr)/(cp*temptr));
                rhb[k] = 0.98;
            end
            #  first guess sans moisture
            tbv[k] = tb[k];
        end

        #  compute the non-dimensional pressure profile as a first guess
        #    (needed for qvs below)
        pisfc = (psurf / p0)^xk;
        pib[2] = pisfc - g*0.5*dz/(cp*tbv[2]);

        for k = 3 : nz - 1
            tbvavg = 0.5*(tbv[k] + tbv[k-1]);
            pib[k] = pib[k-1] - g*dz/(cp*tbvavg);
        end

        #  convert to mixing ratios
        for k=2 : nz-1
            p = p0*pib[k]^(cp/rd);
            temp = tb[k]*pib[k];
            qvs = (380.0/p) * exp((17.27*(temp-273.0))/(temp-36.0));
#             if zu(k) < 20000
            qb[k] = min(qsurf,qvs*rhb[k]);
#             else
#             qb(k) = 0;
#             end
            #  define virtual potential temperatures
            tbv[k] = tb[k]*(1 + 0.61*qb[k]);

        #   if k==98 || k == 99 || k == 100
        #     display(k)
        #     display(p)
        #     display(temp)
        #     display(qvs)
        #     display('qb')
        #     display(qb(k))
        #     display('tbv')
        #     display(tbv(k))
        #     display('rhb')
        #     display(rhb(k))
        #   end
        # end

        #  go back and recompute the non-dimensional pressure profile now that
        #    we have humidity included
        pisfc = (psurf / p0)^xk;
        pib[2] = pisfc - g*0.5*dz/(cp*tbv[2]);

        for k = 3 : nz - 1

            tbvavg = 0.5*(tbv[k] + tbv[k-1]);
            pib[k] = pib[k-1] - g*dz/(cp*tbvavg);

        end

        #  compute the density profile
        pisfc = (psurf/p0)^xk;
        rhow[2] = p0*(pisfc^(c_v/rd))/(rd*tsurf);

        #  define density at the other w points by interpolating from u/scalar points
        #  define density at u/scalar points
      for k=2:nz-1
        rhou[k] = p0*(pib[k]^(c_v/rd))/(rd*tbv[k]);
      end

#  define density at the true surface from known information
      pisfc = (psurf/p0)^xk;
      rhow[2] = p0*(pisfc^(c_v/rd))/(rd*tsurf);

#  define density at the other w points by interpolating from u/scalar points
      for k=3:nz-1
        rhow[k] = 0.5*(rhou[k] + rhou[k-1]);
      end

        #  compute the profiles of saturation mixing ratio and relative humidity
       #  define saturation mixing ratio on u/scalar points
      for k=2: nz-1
        p = p0*pib[k]^(cp/rd);
        temp = tb[k]*pib[k];
#  use Teten's formula (see Pielke's text, 2nd ed., p. 257-258 for more info)
        qbs[k] = (380.0/p) * exp((17.27*(temp-273.0))/(temp-36.0));
        rhb[k] = qb[k] / qbs[k];
      end

 end # End Method If Statement

    #Ensure there is no gradient at the boundaries
    tb[1] = tb[2];
    tb[nz] = tb[nz-1];
    tbv[1] = tbv[2];
    tbv[nz] = tbv[nz-1];
    pib[1] = pib[2];
    pib[nz] = pib[nz-1];
    rhou[1] = rhou[2];
    rhou[nz] = rhou[nz-1];
    rhow[1] = rhow[2];
    rhow[nz] = rhow[nz-1];
    qb[1] = qb[2];
    qb[nz] = qb[nz-1];

## Define U base

# Uniform Flow

    ub[:] .= -10;
    um[:,:] .= -10;
    u[:,:] .= -10;

############################################################
# High Surface Wind decreasing as altitude increases

#     usfc=-10.0;  # u-wind at the surface (m/s)
#     # * for initwinds=1 (only)
#     deltau=10.0;  # u-wind change over the layer (m/s)
#     # * for initwinds=1 (only)
#     shrdepth=3000.0;  # depth of deltau layer (m)
#     #  * for initwinds=1 (only)
#
#
#     for k=2:nz-1
#         #  assign base state values
#         if zu(k)<shrdepth
#             ub(k) = usfc + zu(k)*deltau/shrdepth;
#         else
#             ub(k) = usfc + deltau;
#         end
#         for i= 1 :nx
#             #  copy base state values into prognosed u array
#             u(i,k) = ub(k);
#             um(i,k) = ub(k);
#         end
#     end

###################################################################

# Low-Level Jet at a given z height (given by shrdepth)
#
#     usfc=0;  # u-wind at the surface (m/s)
#     # * for initwinds=1 (only)
#     deltau=-10.0;  # u-wind change over the layer (m/s)
#     # * for initwinds=1 (only)
#     shrdepth=2000.0;  # depth of deltau layer (m)
#     #  * for initwinds=1 (only)
#
#     for k=2:nz-1
#         #  assign base state values
#         if zu(k)<=shrdepth
#             ub(k) = usfc + zu(k)*deltau/shrdepth;
#         elseif zu(k) <= 2*shrdepth
#             ub(k) = 2*deltau  - deltau*(zu(k)/(shrdepth));
#         else
#             ub(k) = 0;
#         end
#     for i= 1 :nx
#             #  copy base state values into prognosed u array
#             u(i,k) = ub(k);
#             um(i,k) = ub(k);
#      end
#      end
######################################################################
    #  handle boundaries
    for i=1:nx
        u[i,1] = ub[2];
        up[i,1] = ub[2];
        um[i,1] = ub[2];
        u[i,nz] = ub[nz-1];
        up[i,nz] = ub[nz-1];
        um[i,nz] = ub[nz-1];
    end

end

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
# num_vector_load_expr(N::Symbol, W) = :(divrem($N, $W))
function num_vector_load_expr(mod, N::Expr, W::Integer)
    if N.args[1] == :length
        :(($mod).VectorizationBase.length_loads($(N.args[2]), Val{$W}()))
    elseif N.args[1] == :size && length(N.args) == 3
        :(($mod).VectorizationBase.size_loads($(N.args[2]), $(N.args[3]), Val{$W}()))
    else
        Wshift = intlog2(W)
        :($N >>> $Wshift, $N & $(W-1))
    end
end
function num_vector_load_expr(mod, N, W::Integer)
    Wshift = intlog2(W)
    :($N >>> $Wshift, $N & $(W-1))
end
### Generic fallbacks
### PaddedMatrices provide methods that determinalistically provide r = 0
@generated function length_loads(A, ::Val{W}) where W
    Wshift = intlog2(W)
    quote
        $(Expr(:meta, :inline))
        N = length(A)
        N >>> $Wshift, N & $(W - 1)
    end
end
@generated function size_loads(A, dim, ::Val{W}) where W
    Wshift = intlog2(W)
    quote
        $(Expr(:meta, :inline))
        N = size(A, dim)
        N >>> $Wshift, N & $(W - 1)
    end
end


text:
# linear algebra
cd("G:\\My Drive\\18. Github Repo\\DataScience");
# load the packages
using LinearAlgebra, Images, SparseArrays

# Basic operations
A = rand(10, 10);
A = A * A';
b = rand(10);
x = A \ b; # solving linear system

@show norm(b - A * x)

# LU factor
luA = lu(A)

# QR factor
qrA = qr(A);

# cholesky
isposdef(A)
cholA = cholesky(A)

# diag
a = Diagonal([1,2,3])

# images
X1 = load(".\\data\\khiam-small.jpg");
X1
Xgray = Gray.(X1)
Xgrayvalues = Float64.(Xgray)

# svd
sSVD_V = svd(Xgrayvalues);

# using the top 10 singular values
global img = zeros(size(Xgrayvalues));
for i = 1:50
    global img += sSVD_V.S[i] * sSVD_V.U[:, i] * sSVD_V.V[:, i]';
end

Gray.(img)

@show norm(Xgrayvalues - img)

#

max_stars_repo_path: 02. Linear Algebra.jl

max_issues_repo_path: 02. Linear Algebra.jl

max_forks_repo_path: 02. Linear Algebra.jl


text:
const ONE_MINUS_COSZERO = 1.0E-12
const THRESHOLD = 0.01
const CHANCE = 0.1

function nearest_index(x, lower, step, N)
    i = Int(round((x - lower) / step))

    # Overflow conditions
    i < 1 && (return 1)
    i > N && (return N)
    return i
end

## TODO: fix this function to remove allocations
function absorb!(p::Photon, v::CatesianVolume, a)
    p.weight -= a
    indices = nearest_index.(p.pos, v.corner1, v.d_xyz, v.Nxyz)
    @inbounds v.data[indices...] += a
    nothing
end

function absorb!(p::Photon, v::RadialVolume, a)
    p.weight -= a
    pr = norm(p.pos)
    i = nearest_index.(pr, v.r1, v.dr, v.Nr)
    @inbounds v.data[i] += a
    nothing
end

function hop!(p::Photon, s)
    hop_vec = s .* p.trj
    p.pos = p.pos + hop_vec
    nothing
end

function scatterHG(rnd, g)
    # TODO: doi: 10.1088/0031-9155/51/17/N04
    g == 0.0 && return 2.0*rnd - 1.0
    temp = (1.0 - g^2)/(1.0 - g + 2*g*rnd)
    cosθ = (1.0 + g^2 - temp^2)/(2.0*g)
end

function spin!(p::Photon, g::Real)
    # Sample a random cos(θ) from Henyey-Greenstein scattering function, inverted function
    cosθ = scatterHG(rand(), g)
    sinθ = sqrt(1.0 - cosθ^2)   # sqrt() is faster than sin().

    sinψ, cosψ = sincos(2π*rand()) # Sample psi

    ux,uy,uz = p.trj
    # New Trajectory
    if (1 - abs(uz) <= ONE_MINUS_COSZERO)      # close to perpendicular.
        ux_new = sinθ * cosψ
        uy_new = sinθ * sinψ
        uz_new = cosθ * sign(uz)   # SIGN() is faster than division.
        new_trj = [ux_new, uy_new, uz_new]
        new_trj = new_trj ./ norm(new_trj)
    else # usually use this option
        temp = sqrt(1.0 - uz^2)
        ux_new =  sinθ * (ux * uz * cosψ - uy * sinψ) / temp + ux * cosθ
        uy_new =  sinθ * (uy * uz * cosψ + ux * sinψ) / temp + uy * cosθ
        uz_new = -sinθ * cosψ * temp + uz * cosθ
        new_trj = [ux_new, uy_new, uz_new]
    end

    p.trj = new_trj
    nothing
end

function roulette!(p::Photon)
    if p.weight < THRESHOLD
        if rand() <= CHANCE
            p.weight /= CHANCE
        else
            p.live = false
        end
    end
    nothing
end

# TODO: make this function parallelizable
function simulate!(volume::Absorber, e::Emitter, num_photons)
    hopdist() = -log(rand())/(volume.prop.μ_a + volume.prop.μ_s)

    for i = 1:num_photons
        phot = makephoton(e)
        # tmp1 = SVector{3, Int32}([1,1,1])

        while phot.live
            a = phot.weight*(1 - volume.prop.albedo)

            hop!(phot, hopdist())
            absorb!(phot, volume, a)
            spin!(phot, volume.prop.g)
            roulette!(phot)
        end
    end
    return volume
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    boundbox(object)

Return the minimum axis-aligned bounding rectangle of the spatial `object`.
"""
boundbox(obj) = boundbox(domain(obj))

function boundbox(obj::AbstractDomain)
  N = ncoords(obj)
  T = coordtype(obj)

  lowerleft  = MVector(ntuple(i->typemax(T), N))
  upperright = MVector(ntuple(i->typemin(T), N))

  x = MVector{N,T}(undef)
  for l in 1:nelms(obj)
    coordinates!(x, obj, l)
    for d in 1:N
      x[d] < lowerleft[d]  && (lowerleft[d]  = x[d])
      x[d] > upperright[d] && (upperright[d] = x[d])
    end
  end

  Rectangle(lowerleft, upperright)
end

function boundbox(obj::RegularGrid)
  lowerleft  = origin(obj)
  upperright = origin(obj) + (size(obj) .- 1) .* spacing(obj)
  Rectangle(lowerleft, upperright)
end


text:
using LineEdit

include("keymaps.jl")


text:
using PyPlot

## Display graph with weigthed edges
graphPlot=plot(decisionTree,"Layout','layered','NodeLabel',nodeLabel,'NodeColor','k','EdgeColor','k")
# toHL=find(decisionTree.Nodes.Sum~=0)
# for HL=1:length(toHL)
# highlight[graphPlot,toHL[HL],"Marker','s','NodeColor','r','MarkerSize",decisionTree.Nodes.Sum[toHL[HL]]+5]
# end
# toHL=find(decisionTree.Nodes.Sum==maximum(decisionTree.Nodes.Sum))
highlight[graphPlot,33,"NodeColor','r','Marker','p','MarkerSize",10]

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


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

mutable struct IoK8sApiMetricsV1beta1NodeMetrics <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    timestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: timestamp
    window::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Duration } # spec name: window
    usage::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: usage

    function IoK8sApiMetricsV1beta1NodeMetrics(;metadata=nothing, timestamp=nothing, window=nothing, usage=nothing)
        o = new()
        validate_property(IoK8sApiMetricsV1beta1NodeMetrics, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiMetricsV1beta1NodeMetrics, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(IoK8sApiMetricsV1beta1NodeMetrics, Symbol("window"), window)
        setfield!(o, Symbol("window"), window)
        validate_property(IoK8sApiMetricsV1beta1NodeMetrics, Symbol("usage"), usage)
        setfield!(o, Symbol("usage"), usage)
        o
    end
end # type IoK8sApiMetricsV1beta1NodeMetrics

const _property_map_IoK8sApiMetricsV1beta1NodeMetrics = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("window")=>Symbol("window"), Symbol("usage")=>Symbol("usage"))
const _property_types_IoK8sApiMetricsV1beta1NodeMetrics = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("timestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("window")=>"IoK8sApimachineryPkgApisMetaV1Duration", Symbol("usage")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}")
Base.propertynames(::Type{ IoK8sApiMetricsV1beta1NodeMetrics }) = collect(keys(_property_map_IoK8sApiMetricsV1beta1NodeMetrics))
Swagger.property_type(::Type{ IoK8sApiMetricsV1beta1NodeMetrics }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiMetricsV1beta1NodeMetrics[name]))}
Swagger.field_name(::Type{ IoK8sApiMetricsV1beta1NodeMetrics }, property_name::Symbol) =  _property_map_IoK8sApiMetricsV1beta1NodeMetrics[property_name]

function check_required(o::IoK8sApiMetricsV1beta1NodeMetrics)
    true
end

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


text:
@testset "Information Measures" begin

    s = res"ARNDCQEGHILKMFPSTWYV-"
    r = reverse(s)
    g = res"GGGGGGGGGGGGGGGGGGGGG"
    Pg = probabilities(g)
    Pgg = probabilities(g,g)
    Pggg = probabilities(g,g,g)
    Ps = probabilities(s)
    Pss = probabilities(s,s)
    Psr = probabilities(s,r)
    Psss = probabilities(s,s,s)
    count_random = count(Residue[], Residue[], pseudocounts=AdditiveSmoothing(1.0))
    prob_random = probabilities(Residue[], Residue[], pseudocounts=AdditiveSmoothing(1.0))

    @testset "Entropy" begin

        @testset "H(X)" begin

            @test entropy(Pg) ≈ 0.0
            @test entropy(Ps) ≈ log(20.0)
            @test entropy(Ps, 2) ≈ log(2, 20.0)
        end

        @testset "Joint Entropy: H(X,Y)" begin

            @test entropy(Pgg) ≈ 0.0
            @test entropy(Pss) ≈ log(20.0)
            @test entropy(Psr) ≈ log(19.0)
            @test entropy(Pss, 2) ≈ log(2, 20.0)
            @test entropy(count_random) ≈ log(400.0)
            @test entropy(prob_random) ≈ log(400.0)
        end

        @testset "Joint Entropy: H(X,Y,Z)" begin

            @test entropy(Pggg) ≈ 0.0
            @test entropy(Psss) ≈ log(20)
        end

        @testset "Using counts" begin

            @test entropy(count(g,g)) ≈ entropy(Pgg)
            @test entropy(count(s,s)) ≈ entropy(Pss)
            @test entropy(count(s,r)) ≈ entropy(Psr)
        end
    end

    @testset "Kullback-Leibler" begin

        @test kullback_leibler(Ps, [1.0/20.0 for i in 1:20]) ≈ 0.0
        @test kullback_leibler(Ps, Ps) ≈ 0.0
        @test kullback_leibler(Ps, BLOSUM62_Pi) ≈ kullback_leibler(Ps)
        @test kullback_leibler(Ps, BLOSUM62_Pi) ≈ mapreduce(i->0.05*log(0.05/BLOSUM62_Pi[i]),+,1:20)
        @test kullback_leibler(Psr, Psr) ≈ 0.0
    end

    @testset "Mutual Information" begin

        @test mutual_information(Pgg) ≈ 0.0

        P = gettablearray(prob_random)
        Q = getmarginalsarray(prob_random)[:,1] * getmarginalsarray(prob_random)[:,2]' # 0.002500000000000001
        @test mutual_information(prob_random) ≈ sum(P .* log.(P./Q)) # So, it's != 0.0
        @test mutual_information(count_random) ≈ 0.0 # It's more accurate

        Pgs = probabilities(g,s)
        P = gettablearray(Pgs)
        Q = getmarginalsarray(Pgs)[:,1] * getmarginalsarray(Pgs)[:,2]' # 0.05000000000000002
        @test mutual_information(Pgs) ≈ mapreduce(x -> isnan(x) ? 0.0 : x, +, P.*log.(P./Q))
        @test mutual_information(count(g,s)) ≈ 0.0 # It's more accurate

        # MI(X,Y)=H(X)+H(Y)-H(X,Y)

        @test mutual_information(Psr) ≈ (marginal_entropy(Psr,1) +
                                         marginal_entropy(Psr,2) - entropy(Psr))

        @test mutual_information(Pss) ≈ (marginal_entropy(Pss,1) +
                                         marginal_entropy(Pss,2) - entropy(Pss))

        @test mutual_information(Psr,2) ≈ (marginal_entropy(Psr,1,2) +
                                           marginal_entropy(Psr,2,2) - entropy(Psr,2))

        @test mutual_information(Pss,20) ≈ (marginal_entropy(Pss,1,20) +
                                            marginal_entropy(Pss,2,20) - entropy(Pss,20))

        @test isapprox(mutual_information(prob_random), 0.0, atol=1e-15)
        @test mutual_information(count_random) ≈ 0.0
        @test isapprox(mutual_information(count_random), marginal_entropy(count_random,1) +
                       marginal_entropy(count_random,2) - entropy(count_random), atol=1e-13)

        @testset "Using counts" begin

            @test mutual_information(count_random) ≈ 0.0
            @test mutual_information(count(g,g)) ≈ mutual_information(Pgg)
            @test mutual_information(count(s,s)) ≈ mutual_information(Pss)
            @test mutual_information(count(s,r)) ≈ mutual_information(Psr)
        end

        @testset "MI(X,Y,Z)" begin

            @test mutual_information(Pggg) ≈ 0.0
            @test mutual_information(probabilities(g,s,r)) ≈ 0.0
            @test mutual_information(probabilities(s,s,r)) ≈ mutual_information(count(s,s,r))

            # MI(X,Y,Z) = H(X) + H(Y) + H(Z) - H(X,Y) - H(X,Z) - H(Y,Z) + H(X,Y,Z)
            @test mutual_information(Psss) ≈ mutual_information(count(s,s,s))
            @test mutual_information(Psss) ≈ (marginal_entropy(Psss, 1) +
                                              marginal_entropy(Psss, 2) +
                                              marginal_entropy(Psss, 3) -
                                              entropy(Pss) - entropy(Pss) -
                                              entropy(Pss) + entropy(Psss))

            # MI(X,Y,Z) <= min{ H(X,Y), H(X,Z), H(Y,Z) }
            @test mutual_information(Psss) ≈ mutual_information(Pss)
        end

        @testset "Pairwise Gap Percentage" begin

            @test gap_union_percentage(count(res"AA--", res"--AA", alphabet=GappedAlphabet())) ≈ 100.0
            @test gap_intersection_percentage(count(res"AA--", res"--AA", alphabet=GappedAlphabet())) ≈ 0.0

            @test gap_union_percentage(count(res"AA--", res"--AA", alphabet=GappedAlphabet(), weights=Weights([.25, .25, .25, .25]))) ≈ 100.0
            @test gap_intersection_percentage(count(res"AA--", res"--AA", alphabet=GappedAlphabet(), weights=Weights([.25, .25, .25, .25]))) ≈ 0.0

            @test gap_union_percentage(count(res"AAA-", res"AA--", alphabet=GappedAlphabet())) ≈ 50.0
            @test gap_intersection_percentage(count(res"AAA-", res"AA--", alphabet=GappedAlphabet())) ≈ 25.0

            @test gap_union_percentage(count(res"AAA-", res"AA--", alphabet=GappedAlphabet(), weights=Weights([.2, .2, .2, .4]))) ≈ 60.0
            @test gap_intersection_percentage(count(res"AAA-", res"AA--", alphabet=GappedAlphabet(), weights=Weights([.2, .2, .2, .4]))) ≈ 40.0
        end
    end
end


text:
using Luna

a = 13e-6
gas = :Ar
pres = 5

τfwhm = 30e-15
λ0 = 800e-9
flength = 0.5e-2
energy = 1e-6

grid = Grid.EnvGrid(flength, λ0, (160e-9, 3000e-9), 1e-12, thg=true)

m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff(z) = Modes.Aeff(m, z=z)

energyfun, energyfunω = Fields.energyfuncs(grid)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

linop, βfun!, frame_vel, αfun = LinearOps.make_const_linop(grid, m, λ0, thg=true)

responses = (Nonlinear.Kerr_env_thg(PhysData.γ3_gas(gas), 2π*PhysData.c/λ0, grid.to),)

    inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.collect_stats(grid, Eω,
                               Stats.ω0(grid),
                               Stats.energy(grid, energyfunω),
                               Stats.energy_λ(grid, energyfunω, (150e-9, 300e-9), label="RDW"),
                               Stats.peakpower(grid),
                               Stats.fwhm_t(grid),
                               Stats.density(densityfun))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

import FFTW
import PyPlot: pygui, plt

ω = grid.ω
t = grid.t

zout = output.data["z"]
Eout = output.data["Eω"]

Etout = FFTW.ifft(Eout, 1)

Ilog = log10.(Maths.normbymax(abs2.(Eout)))

idcs = @. (t < 30e-15) & (t >-30e-15)
It = abs2.(Etout)
Itlog = log10.(Maths.normbymax(It))
zpeak = argmax(dropdims(maximum(It, dims=1), dims=1))

Et = Etout
energy = zeros(length(zout))
for ii = 1:size(Etout, 2)
    energy[ii] = energyfun(Etout[:, ii])
end

pygui(true)
plt.figure()
plt.pcolormesh(FFTW.fftshift(ω, 1)./2π.*1e-15, zout, transpose(FFTW.fftshift(Ilog, 1)))
plt.clim(-15, 0)
plt.xlim(0.19, 1.9)
plt.colorbar()

plt.figure()
plt.pcolormesh(t*1e15, zout, transpose(It))
plt.colorbar()
plt.xlim(-30, 30)

plt.figure()
plt.plot(zout.*1e2, energy.*1e6)
plt.xlabel("Distance [cm]")
plt.ylabel("Energy [μJ]")

plt.figure()
plt.plot(t*1e15, abs2.(Et[:, 121]))
plt.xlim(-20, 20)


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

"""
    PseudocostBranching()

Branching strategy that uses historical changes in objective value to estimate
strong branching scores at lower computational cost.
"""
struct PseudocostBranching <: VariableBranchingRule end

function find_branching_var(rule::PseudocostBranching, node::Node, pool::NodePool)::Variable
    scores = [
        _pseudocost_score(
            node,
            pool,
            node.fractional_variables[j],
            node.fractional_values[j],
        ) for j = 1:length(node.fractional_variables)
    ]
    _, max_offset = findmax(scores)
    return node.fractional_variables[max_offset]
end

function _pseudocost_score(
    node::Node,
    pool::NodePool,
    var::Variable,
    x::Float64,
)::Tuple{Float64,Int}
    f_up = x - floor(x)
    f_down = ceil(x) - x
    pseudo_up = pool.history.avg_pseudocost_up * f_up
    pseudo_down = pool.history.avg_pseudocost_down * f_down
    if var in keys(pool.var_history)
        if isfinite(pool.var_history[var].pseudocost_up)
            pseudo_up = pool.var_history[var].pseudocost_up * f_up
        end
        if isfinite(pool.var_history[var].pseudocost_down)
            pseudo_down = pool.var_history[var].pseudocost_down * f_down
        end
    end
    return (pseudo_up * f_up * pseudo_down * f_down, var.index)
end


text:
# Compare Flux model from Metalhead to PyTorch model
# for a sample image

# PyTorch need to be installed

# Tested on ResNet and VGG models

using Flux
import Metalhead
using DataStructures
using Statistics
using BSON
using PyCall
using Images
using Test

using MLUtils
using Random

torchvision = pyimport("torchvision")
torch = pyimport("torch")

modellib = [
    ("vgg11", () -> Metalhead.VGG(11), torchvision.models.vgg11),
    ("vgg13", () -> Metalhead.VGG(13), torchvision.models.vgg13),
    ("vgg16", () -> Metalhead.VGG(16), torchvision.models.vgg16),
    ("vgg19", () -> Metalhead.VGG(19), torchvision.models.vgg19),
    ("resnet18", () -> Metalhead.ResNet(18), torchvision.models.resnet18),
    ("resnet34", () -> Metalhead.ResNet(34), torchvision.models.resnet34),
    ("resnet50", () -> Metalhead.ResNet(50), torchvision.models.resnet50),
    ("resnet101",() -> Metalhead.ResNet(101),torchvision.models.resnet101),
    ("resnet152",() -> Metalhead.ResNet(152),torchvision.models.resnet152),
]

tr(tmp) = permutedims(tmp,ndims(tmp):-1:1)

function normalize(data)
    cmean = reshape(Float32[0.485, 0.456, 0.406],(1,1,3,1))
    cstd = reshape(Float32[0.229, 0.224, 0.225],(1,1,3,1))
    return (data .- cmean) ./ cstd
end

# test image
guitar_path = download("https://cdn.pixabay.com/photo/2015/05/07/11/02/guitar-756326_960_720.jpg")

# image net labels
labels = readlines(download("https://raw.githubusercontent.com/pytorch/hub/master/imagenet_classes.txt"))

weightsdir = joinpath(@__DIR__,"..","weights")

for (modelname,jlmodel,pymodel) in modellib
    println(modelname)

    model = jlmodel()

    saved_model = BSON.load(joinpath(weightsdir,"$(modelname).bson"))
    Flux.loadmodel!(model,saved_model[:model])

    pytorchmodel = pymodel(pretrained=true)

    Flux.testmode!(model)

    sz = (224, 224)
    img = Images.load(guitar_path);
    img = imresize(img, sz);
    # CHW -> WHC
    data = permutedims(convert(Array{Float32}, channelview(img)), (3,2,1))
    data = normalize(data[:,:,:,1:1])

    out = model(data) |> softmax;
    out = out[:,1]

    println("  Flux:")

    for i in sortperm(out,rev=true)[1:5]
        println("    $(labels[i]): $(out[i])")
    end

    pytorchmodel.eval()
    output =  pytorchmodel(torch.Tensor(tr(data)));
    probabilities = torch.nn.functional.softmax(output[0], dim=0).detach().numpy();

    println("  PyTorch:")

    for i in sortperm(probabilities[:,1],rev=true)[1:5]
        println("    $(labels[i]): $(probabilities[i])")
    end

    @test maximum(out) ≈ maximum(probabilities)
    @test argmax(out) ≈ argmax(probabilities)

    println()
end

max_stars_repo_name: Alexander-Barth/MetalheadWeights

max_issues_repo_name: Alexander-Barth/MetalheadWeights

max_forks_repo_name: Alexander-Barth/MetalheadWeights


text:
type ServerConfig
  port
  responseConfigurations
end

@enum HTTP_METHOD GET POST PUT PATCH DELETE

type ResponseConfig
  method::HTTP_METHOD
  path::String
  responseHandler
end

type WhipRestServer
  port
  getMappings::Dict{String,Any}
  postMappings::Dict{String,Any}
  putMappings::Dict{String,Any}
  patchMappings::Dict{String,Any}
  deleteMappings::Dict{String,Any}
end

type RestMessage
  uriParams::Dict{String, String}
  body::String
end

max_stars_repo_name: WhippyCleric/WhipRest

max_issues_repo_name: WhippyCleric/WhipRest

max_forks_repo_name: WhippyCleric/WhipRest


text:
function findatomtype(alist, atomtype, typei)

atype = atomtype[alist]; # atom types of alist
ilist = alist[atype .== typei];    # return atoms of type i

return ilist

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# Uniform heatmap on the left
subplot(1,2,1)
x = LinRange(-2, 2, 40)
y = LinRange(0, pi, 20)
z = sin.(x') .+ cos.(y)
heatmap(z)
# Non uniform heatmap on the right
subplot(1,2,2)
x = [0, 2, 3, 4.5, 5]
y = [2, 3, 4.5, 5, 6, 8]
z = rand(5, 4)
heatmap(x, y, z)


text:
# reconnect logic

# Need to define get_application_rate_limit_status in order to make it available to get_endpoint_allocation

function get_application_rate_limit_status(;kwargs...)
    endp = "application/rate_limit_status.json"
    #Take kwargs array of tuples, create dict needed by OAuth process
    options = Dict{String, Any}()
    for arg in kwargs
        options[string(arg[1])] = string(arg[2])
    end

    # defines the functions
    r = get_oauth("https://api.twitter.com/1.1/$endp", options)

    #If successful API call, return JSON as Julia data structure, otherwise return error
    if r.status == 200
        success = JSON.parse(String(r.body))
    else
        error("Twitter API returned $(r.status) status")
    end

end

# retrieve endpoint remaining calls
get_endpoint_allocation  = function(endp, endp_deduction=0)

    if !@isdefined(API_INFO) || API_INFO==nothing
        global API_INFO = get_application_rate_limit_status()
    end
    #logic to acquire the correct endpoint
    endpoint_match = match( r"^.*/", endp).match # get everything before the slash
    base_endpoint =  strip(endpoint_match, '/') #remove the slash
    final_endpoint = replace(endp, ".json" => "") # remove the .json
    base_keys = keys(API_INFO["resources"][base_endpoint])
    endp_array = [match(Regex("/$(final_endpoint)(.*)"), x) for x in String.(base_keys)]
    no_limit = sum([x!=nothing for x in endp_array])==0 #endp not in list, end eval
    no_limit && return Dict("remaining" => -1, "reset" => 0, "limit"=> -1)
    endp_name = [x.match for x in endp_array if x != nothing][1]
    # subtract the deduction if necessary
    API_INFO["resources"][base_endpoint][endp_name]["remaining"] -= endp_deduction
    return API_INFO["resources"][base_endpoint][endp_name]
end

# back off loop - this will reconnect when the API says it's OK
function reconnect(endp, reconnects=0)
    while get_endpoint_allocation(endp)["remaining"]==0
        reconnects += 1
        global API_INFO = nothing # reset the API INFO because
        alloc = get_endpoint_allocation(eval(endp))
        cur_time = round(Int64, time())
        sleeptime = abs(alloc["reset"] - cur_time)
        println("Endpoint exhausted, sleeping for $sleeptime seconds..")
        sleep(sleeptime^reconnects)
    end
    return get_endpoint_allocation(endp, 1)
end


text:
#=
    contents

hi, lo = two_sum(a, b)
hi, lo = two_sum(a, b, c)
hi, lo = two_hilo_sum(a, b)
hi, lo = two_lohi_sum(a, b)

hi, lo = two_diff(a, b)
hi, lo = two_hilo_diffa, b)
hi, lo = two_lohi_diffa, b)

hi, lo = two_prod(a, b)
hi, lo = two_prod(a, b, c)

=#

"""
    two_sum(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_sum(a::T, b::T) where {T}
    hi = a + b
    db = b - (hi - a)
    da = a - (hi - b)
    lo = da + db
    return hi, lo
end

"""
    two_sum_ordered(a, b)

Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_sum_ordered(a::T, b::T) where {T}
    if abs(a) < abs(b)
        a, b = b, a
    end
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
   two_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
function two_sum(a::T, b::T, c::T) where {T}
    t0, t1 = two_sum(b, c)
    hi, t0 = two_sum(a, t0)
    t0 += t1
    hi, lo = two_hilo_sum(hi, t0)
    return hi, lo
end

# arguments sorted by magnitude

"""
    two_hilo_sum(a, b)
*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T}
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

"""
    two_lohi_sum(a, b)

*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_lohi_sum(a::T, b::T) where {T}
    hi = b + a
    lo = a - (hi - b)
    return hi, lo
end

"""
    two_diff(a, b)
Computes `s = fl(a-b)` and `e = err(a-b)`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_diff(a::T, b::T) where {T}
    hi = a - b
    db = (a - hi) - b
    da = a - (hi + b)
    lo = da + db
    return hi, lo
end

"""
    two_hilo_diff(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_hilo_diff(a::T, b::T) where {T}
    hi = a - b
    lo = (a - hi) - b
    hi, lo
end

"""
    two_lohi_diff(a, b)

*unchecked* requirement `|b| ≥ |a|`
Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_lohi_diff(a::T, b::T) where {T}
    hi = b - a
    lo = (b - hi) - a
    hi, lo
end

"""
    two_prod(a, b)

Computes `hi = fl(a*b)` and `lo = fl(err(a*b))`.
- Unchecked Precondition: !(isinf(a) | isinf(b))
"""
@inline function two_prod(a::T, b::T) where {T}
    hi = a * b
    lo = fma(a, b, -hi)
    hi, lo
end

"""
   two_prod(a, b, c)

Computes `hi = fl(a*b*c)` and `lo = err(a*b*c)`.
"""
function two_prod(a,b,c)
    a, b, c = maxminabs(a, b, c)
    cbhi, cblo = two_prod(c, b)
    hi, cba = two_prod(cbhi, a)
    lo = fma(a, cblo, cba)
    return hi, lo
end