text:
temperature_switch(a,b,c,p::T) where {T} = p.numerics.temperature == false ? IfElse.ifelse.(a,b,c) : c

## Electrolyte functions
function D_s_eff_isothermal(c_s_avg_p, c_s_avg_n, T_p, T_n, p::AbstractModel)
    """
    D_s_eff_isothermal evaluates diffusion coefficients of the solid phase [m²/s].
    The user may modify the script to meet specific requirements
    """

    D_sp_eff = repeat([p.θ[:D_sp]], p.N.p)
    D_sn_eff = repeat([p.θ[:D_sn]], p.N.n)

    return D_sp_eff, D_sn_eff
end

function D_s_eff(c_s_avg_p, c_s_avg_n, T_p, T_n, p::AbstractModel)
    """
    D_s_eff evaluates diffusion coefficients of the solid phase [m²/s].
    The user may modify the script to meet specific requirements
    """

    R = const_Ideal_Gas

    T_ref = 298.15

    # oftentimes T == T_ref for isothermal experiments, so the switch prevents unnecessary calculations
    D_sp_eff = p.θ[:D_sp] .* temperature_switch( T_p .== T_ref, 1, exp.(-p.θ[:Ea_D_sp]/R*(1.0./T_p .- 1.0./T_ref)), p )
    D_sn_eff = p.θ[:D_sn] .* temperature_switch( T_n .== T_ref, 1, exp.(-p.θ[:Ea_D_sn]/R*(1.0./T_n .- 1.0./T_ref)), p )

    return D_sp_eff, D_sn_eff
end

function rxn_rate_isothermal(T_p, T_n, c_s_avg_p, c_s_avg_n, p::AbstractModel)
    """
    Reaction rates (k) of cathode and anode [m^2.5/(m^0.5 s)]
    """

    k_p = p.θ[:k_p]
    k_n = p.θ[:k_n]

    return k_p, k_n
end

function rxn_rate(T_p, T_n, c_s_avg_p, c_s_avg_n, p::AbstractModel)
    """
    Reaction rates (k) of cathode and anode [m^2.5/(m^0.5 s)]
    """

    R = const_Ideal_Gas
    T_ref = 298.15

    # oftentimes T == T_ref for isothermal experiments. the if statement prevents unnecessary calculations
    k_p = p.θ[:k_p] .* temperature_switch( T_p .== T_ref, 1, exp.(-(p.θ[:Ea_k_p]./R) .* (1.0./T_p .- 1.0./T_ref)), p )
    k_n = p.θ[:k_n] .* temperature_switch( T_n .== T_ref, 1, exp.(-(p.θ[:Ea_k_n]./R) .* (1.0./T_n .- 1.0./T_ref)), p )

    return k_p, k_n
end

function D_eff_linear(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    D_eff_linear evaluates the diffusion coefficients for the electrolyte phase [m^2/s]
    """

    D_eff_p = p.θ[:D_p].*p.θ[:ϵ_p].^p.θ[:brugg_p].*ones(p.N.p)
    D_eff_s = p.θ[:D_s].*p.θ[:ϵ_s].^p.θ[:brugg_s].*ones(p.N.s)
    D_eff_n = p.θ[:D_n].*p.θ[:ϵ_n].^p.θ[:brugg_n].*ones(p.N.n)

    return D_eff_p, D_eff_s, D_eff_n
end

function D_eff_linear_one_term(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    D_eff_linear evaluates the diffusion coefficients for the electrolyte phase [m^2/s]
    """

    D_eff_p = p.θ[:D_e].*p.θ[:ϵ_p].^p.θ[:brugg_p].*ones(p.N.p)
    D_eff_s = p.θ[:D_e].*p.θ[:ϵ_s].^p.θ[:brugg_s].*ones(p.N.s)
    D_eff_n = p.θ[:D_e].*p.θ[:ϵ_n].^p.θ[:brugg_n].*ones(p.N.n)

    return D_eff_p, D_eff_s, D_eff_n
end

function D_eff(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    D_eff evaluates the diffusion coefficients for the electrolyte phase [m^2/s]
    """

    func(c_e, T) = 1e-4*10.0^((-4.43 - 54.0 / (T - 229-5e-3 * c_e) - 0.22e-3 * c_e))

    D_eff_p = p.θ[:ϵ_p].^p.θ[:brugg_p].*func.(c_e_p, T_p)
    D_eff_s = p.θ[:ϵ_s].^p.θ[:brugg_s].*func.(c_e_s, T_s)
    D_eff_n = p.θ[:ϵ_n].^p.θ[:brugg_n].*func.(c_e_n, T_n)

    return D_eff_p, D_eff_s, D_eff_n
end

function K_eff(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    K_eff evaluates the conductivity coefficients for the electrolyte phase [S/m]
    """

    func(c_e, T) = @. 1e-4*c_e*((-10.5 + 0.668*1e-3*c_e + 0.494*1e-6*c_e^2) + (0.074 - 1.78*1e-5*c_e - 8.86*1e-10*c_e^2)*T + (-6.96*1e-5 + 2.8*1e-8*c_e)*T^2)^2

    K_eff_p = p.θ[:ϵ_p].^p.θ[:brugg_p] .* func(c_e_p, T_p)
    K_eff_s = p.θ[:ϵ_s].^p.θ[:brugg_s] .* func(c_e_s, T_s)
    K_eff_n = p.θ[:ϵ_n].^p.θ[:brugg_n] .* func(c_e_n, T_n)

    return K_eff_p, K_eff_s, K_eff_n
end

function K_eff_isothermal(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    K_eff_isothermal evaluates the conductivity coefficients for the electrolyte phase [S/m]
    """

    func(c_e) = @. 4.1253*1e-2 + 5.007*1e-4*c_e - 4.7212*1e-7*c_e^2 + 1.5094*1e-10*c_e^3 - 1.6018*1e-14*c_e^4

    K_eff_p = p.θ[:ϵ_p].^p.θ[:brugg_p] .* func(c_e_p)
    K_eff_s = p.θ[:ϵ_s].^p.θ[:brugg_s] .* func(c_e_s)
    K_eff_n = p.θ[:ϵ_n].^p.θ[:brugg_n] .* func(c_e_n)

    return K_eff_p, K_eff_s, K_eff_n
end

## OCV – Positive Electrodes
function OCV_LCO(θ_p, T=298.15, p=nothing)
    T_ref = 25 + 273.15

    # Define the OCV for the positive electrode
    U_p = (-4.656 .+ 88.669θ_p.^2 .- 401.119θ_p.^4 .+ 342.909θ_p.^6 .- 462.471θ_p.^8 .+ 433.434θ_p.^10)./(-1 .+ 18.933θ_p.^2 .- 79.532θ_p.^4 .+ 37.311θ_p.^6 .- 73.083θ_p.^8 .+ 95.96θ_p.^10)

    # Compute the variation of OCV with respect to temperature variations [V/K]
    ∂U∂T_p = -0.001(0.199521039 .- 0.928373822θ_p .+ 1.364550689000003θ_p.^2 .- 0.6115448939999998θ_p.^3)./(1 .- 5.661479886999997θ_p +11.47636191θ_p.^2 .- 9.82431213599998θ_p.^3 .+ 3.048755063θ_p.^4)

    # if T == T_ref exactly (which may be true for isothermal runs), don't calculate ∂U∂T
    U_p += temperature_switch( T .== T_ref, 0, ∂U∂T_p.*(T .- T_ref), p )

    return U_p, ∂U∂T_p
end

## OCV – Negative Electrodes
function OCV_LiC6(θ_n, T=298.15, p=nothing)
    T_ref = 25 + 273.15

    # Calculate the open circuit voltage of the battery in the negative electrode
    U_n = @. 0.7222 + 0.1387*θ_n + 0.029*sqrt_ReLU(θ_n) - 0.0172./θ_n + 0.0019./(sqrt_ReLU(θ_n;minval=1e-4).*θ_n) + 0.2808*exp(0.9-15*θ_n)-0.7984*exp(0.4465*θ_n - 0.4108)

    # Compute the variation of OCV with respect to temperature variations [V/K]
    ∂U∂T_n = @. 0.001*(0.005269056 +3.299265709*θ_n-91.79325798*θ_n.^2+1004.911008*θ_n.^3-5812.278127*θ_n.^4 + 19329.7549*θ_n.^5 - 37147.8947*θ_n.^6 + 38379.18127*θ_n.^7-16515.05308*θ_n.^8)
    ∂U∂T_n ./= @. (1-48.09287227*θ_n+1017.234804*θ_n.^2-10481.80419*θ_n.^3+59431.3*θ_n.^4-195881.6488*θ_n.^5 + 374577.3152*θ_n.^6 - 385821.1607*θ_n.^7 + 165705.8597*θ_n.^8)

    U_n += temperature_switch( T .== T_ref, 0, ∂U∂T_n.*(T .- T_ref), p )

    return U_n, ∂U∂T_n
end

function OCV_NMC(θ_p, T=298.15, p=nothing)
    # Define the OCV for the positive electrode
    U_p = @. -10.72θ_p^4+ 23.88θ_p^3 - 16.77θ_p^2 + 2.595θ_p + 4.563

    # Compute the variation of OCV with respect to temperature variations [V/K]
    ∂U∂T_p = zeros(length(U_p))

    return U_p, ∂U∂T_p
end

function OCV_LiC6_with_NMC(θ_n, T=298.15, p=nothing)
    # Define the OCV for the negative electrode
    U_n = @. 0.1493 + 0.8493exp(-61.79θ_n) + 0.3824exp(-665.8θ_n) -
        exp(39.42θ_n-41.92) - 0.03131atan(25.59θ_n - 4.099) -
        0.009434atan(32.49θ_n - 15.74)

    # Compute the variation of OCV with respect to temperature variations [V/K]
    ∂U∂T_n = zeros(length(U_n))

    return U_n, ∂U∂T_n
end

## Thermodynamic factors
function thermodynamic_factor_linear(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    Thermodynamic factor for the activity coefficient. The term `(1-t₊)` is included elsewhere,
    do not include the multiple in this function
    """
    func(c_e, T) = ones(length(c_e))

    ν_p = func(c_e_p, T_p)
    ν_s = func(c_e_s, T_s)
    ν_n = func(c_e_n, T_n)

    return ν_p, ν_s, ν_n
end

function thermodynamic_factor(c_e_p, c_e_s, c_e_n, T_p, T_s, T_n, p::AbstractModel)
    """
    Thermodynamic factor for the activity coefficient. The term `(1-t₊)` is included elsewhere,
    do not include the multiple in this function
    """
    func(c_e, T) = @. (0.601 - 0.24(c_e/1000)^0.5 +0.982*(1-0.0052(T-293))*(c_e/1000)^1.5)

    ν_p = func(c_e_p, T_p)
    ν_s = func(c_e_s, T_s)
    ν_n = func(c_e_n, T_n)

    return ν_p, ν_s, ν_n
end

## Reaction Rates
sqrt_ReLU(x;minval=0) = sqrt(max(minval,x))

function rxn_BV(c_s_star, c_e, T, η, k_i, λ_MHC, c_s_max, p;
    sqrt=sqrt_ReLU
    )
    # The modified sqrt function is to avoid errors when concentrations momentarily become non-physical

    # Butler-Volmer rate equation

    α = 0.5
    F = const_Faradays
    R = const_Ideal_Gas

    # simplified version which is faster to calculate
    if α == 0.5
        j_i = @. 2.0k_i*sqrt(c_e*c_s_star*(c_s_max - c_s_star))*(sinh(0.5*F*η/(R*T)))
    else
        j_i = @. k_i*(c_e)^(1 - α)*c_s_star^α*(c_s_max - c_s_star)^(1 - α)*(exp(α*F*η/(R*T)) - exp(-(α*F*η/(R*T))))
    end

    return j_i
end

function MHC_kfunc(η, λ)
    a = 1.0 + sqrt(λ)

    k = sqrt(π * λ) * (1.0 - erf((λ - sqrt((a) + η ^ 2)) / (2.0 * sqrt(λ)))) / (1.0 + exp(-η))

    return k
end

function rxn_MHC(c_s_star, c_e, T, η, k_i, λ_MHC, c_s_max, p;
    sqrt=sqrt_ReLU
    )
    # The modified sqrt function is to avoid errors when concentrations momentarily become non-physical

    # MHC rate equation
    # See Zeng, Smith, Bai, Bazant 2014
    # Convert to "MHC overpotential"

    # normalizing the concentrations
    F = const_Faradays
    R = const_Ideal_Gas

    η̂ = η.*(F./(R.*T))
    θ_i = c_s_star./c_s_max
    ĉ_e = c_e./p.θ[:c_e₀]
    α = 0.5

    k₀ = k_i ./ MHC_kfunc(0.0, λ_MHC)

    log_ReLU(x;minval=0) = log(max(minval,x))
    η_f = η̂ .+ log_ReLU.(ĉ_e./θ_i;minval=1e-4)

    if α == 0.5

        a = 1.0 + sqrt(λ_MHC)

        k₀ = @. k_i / ((1.0 - erf((λ_MHC - sqrt((a))) / (2.0 * sqrt(λ_MHC)))) / 2.0)

        coeff_rd_ox = @. k₀*((1.0 - erf((λ_MHC - sqrt((a) + η_f ^ 2)) / (2.0 * sqrt(λ_MHC)))))

        j_i = @. coeff_rd_ox * (1.0/(1.0 + exp(-η_f))*p.θ[:c_e₀]*c_s_star - 1.0/(1.0 + exp(+η_f))*c_e*c_s_max) * sqrt((1.0 - c_s_star/c_s_max)/p.θ[:c_e₀])
    elseif false

        k_rd = k₀.*MHC_kfunc.(-η_f, λ_MHC)
        k_ox = k₀.*MHC_kfunc.(+η_f, λ_MHC)

        j_i = @. (k_ox*p.θ[:c_e₀]*c_s_star - k_rd*c_e*c_s_max)*sqrt((c_s_max - c_s_star)/(p.θ[:c_e₀] * c_s_max))
    else

        k_rd = k₀.*MHC_kfunc.(-η_f, λ_MHC)
        k_ox = k₀.*MHC_kfunc.(+η_f, λ_MHC)

        act_R = θ_i./(1.0 .- θ_i)

        γ_ts = 1.0./(1.0 .- θ_i) # needed for ecd_extras

        ecd_extras = ĉ_e.^(1.0 .- α) .* act_R.^α ./ (γ_ts.*sqrt.(ĉ_e.*θ_i))

        j_i = ecd_extras.*(k_rd.*ĉ_e .- k_ox.*θ_i)
        j_i .*= (-p.θ[:c_e₀].^(1.0 .- α).*c_s_max) # need this to match dimensional eqn.

        j_i = @. (c_s_max - c_s_star)*((p.θ[:c_e₀]*c_s_star*k_ox - c_e*c_s_max*k_rd)*sqrt((c_s_star*(c_s_max - c_s_star))/(c_s_max*c_e*p.θ[:c_e₀])))
    end

    return j_i
end

function rxn_BV_γMod_01(c_s_star, c_e, T, η, k_i, λ_MHC, c_s_max, p)
    # Butler-Volmer rate equation

    # normalizing the concentrations
    F = const_Faradays
    R = const_Ideal_Gas
    θ_i = c_s_star./c_s_max
    ĉ_e = c_e./p.θ[:c_e₀]
    η̂ = η.*(F./(R.*T))
    act_R = θ_i./(1.0 .- θ_i)
    α = 0.5

    γ_ts = 1.0./(θ_i.*(1.0 .- θ_i)) # needed for ecd_extras

    ecd = @. k_i * ĉ_e^(1-α) * act_R^(α) / γ_ts
    j_i = @. ecd * (exp(-α*η̂) - exp((1.0-α)*η̂))

    j_i .*= (-p.θ[:c_e₀].^(1.0 .- α).*c_s_max)

    return j_i

end


text:
#=
Julia script to convert data of hills in Czech Republic into useable format ("N50° 44.161 E15° 44.380" is not a nice field)

Original data were downloaded from https://www.ultratisicovky.cz/ke-stazeni/
and saved in OpenOffice as CSV file and header was manually modified.

Author: Pavel Perina
Date:   2021-08-27
=#

import CSV
import DataFrames

dfIn  = CSV.File("cz-prom100.csv",select=[4,5,12],type=String) |> DataFrames.DataFrame
dfOut = DataFrames.DataFrame(Summit = String[], Elevation = Float64[], Latitude=Float64[], Longitude=Float64[])

nRows = size(dfIn)[1]
for row in 1:nRows
#    println(row)
    mountain  = replace(dfIn[row,1], r"^\s+|\s+$|\s+(?=\s)" => "")
    elevation = parse(Float64, match(r"(\d+)", dfIn[row,2])[1])
    m = match(r"^(\s*)N(?<ndeg>\d*)° (?<nmin>([0-9]*[.])?[0-9]*) E(?<edeg>\d*)° (?<emin>([0-9]*[.])?[0-9]*)", dfIn[row,3])
    north = parse(Float64, m["ndeg"]) + parse(Float64, m["nmin"])/60.0
    east  = parse(Float64, m["edeg"]) + parse(Float64, m["emin"])/60.0
    push!(dfOut, (mountain, elevation, north, east))
end
CSV.write("data-cz-prom100.tsv",dfOut; delim='\t',quotestrings=true)

ext: jl

max_stars_repo_name: pavel-perina/panorama-jl

max_issues_repo_name: pavel-perina/panorama-jl

max_forks_repo_name: pavel-perina/panorama-jl


text: Dict{String, Any}("x" => Dict{String, Any}("y" => Dict{String, Any}("z" => Dict{String, Any}("w" => Dict{String, Any}()))))


text:
function generate_next(obj::AlternatingSampler)
    return jcall(obj, "generateNext", void, ())
end

function get_current(obj::AlternatingSampler)
    return jcall(obj, "getCurrent", jlong, ())
end

function get_mid_point(obj::AlternatingSampler)
    return jcall(obj, "getMidPoint", jdouble, ())
end

function in_range(obj::AlternatingSampler)
    return jcall(obj, "inRange", jboolean, ())
end

function set_radius(obj::AlternatingSampler, arg0::jdouble)
    return jcall(obj, "setRadius", void, (jdouble,), arg0)
end


text:
# enable iteration for Multisets

import Base.iterate

# In iterating in a multiset, we should see each element repeatedly
# corresponding to its multiplicity in the multiset.

# We bulid a state as follows:
# ss = (val, k, next)
# val is the val we returned
# k is the number of times we returned val
# next is the state returned by iterate on M.data

function iterate(M::Multiset, ss)
    val, k, s = ss
    if k < M[val]
        return (val, (val, k + 1, s))
    end

    result = iterate(M.data, s)
    if isnothing(result)
        return nothing
    end
    pr, st = result
    val = pr[1]
    return (val, (val, 1, st))
end

function iterate(M::Multiset)
    clean!(M)  # get rid of all multiplicity zero items (if any)
    result = iterate(M.data)
    if isnothing(result)
        return nothing
    end
    pr, st = result   # this is what iterate(M.data) gave us
    val = pr[1]
    return (val, (val, 1, st))
end


text:
# --- Interrupt Handling --- #

is_interrupt_pending() = ccall(ut_is_interrupt_pending[], Cuchar, ()) != 0

function check_for_interrupt()
    if is_interrupt_pending()
        throw(InterruptException())
    end
end

# --- Exception Handling --- #

# MatlabException is basically identical to MxArray, but is an exception
# MATLAB is responsible for both creating and modifying the exception data
# so no garbage collection should be performed on the Julia side
mutable struct MatlabException <: Exception
    ptr::Ptr{Cvoid}
    MatlabException(ptr::Ptr{Cvoid}) = new(ptr)
end

# creates a MatlabException from an id and a message
function MatlabException(id::String, msg::String)

    # create error ID and message text
    msgID = MATLAB.mxarray("jl:"*replace(id, r"[^\w\s]"=>"_"))
    msgtext = MATLAB.mxarray(msg)

    # call MException to generate a MATLAB exception
    args = [MATLAB.mxarray(msgID), MATLAB.mxarray(msg)]
    mx = call_matlab(1, "MException", args)[1]

    # wrap as a MatlabException and return
    return MatlabException(mx.ptr)
end

# creates a MATLABException from an exception
function MatlabException(exn::Exception)

    # print formatted error message to buffer
    buf = IOBuffer()
    showerror(buf, exn)

    # move to beginning of buffer
    seek(buf, 0)

    # set id as the exception type
    id = string(typeof(exn))

    # convert to string and format for use with MATLAB
    msg = replace(read(buf, String), "\\"=>"\\\\")

    return MatlabException(id, msg)
end

# don't create a new MATLABException from a MATLABException
MatlabException(mexn::MatlabException) = mexn

# adds a cause to a MatlabException
function MatlabException(mexn::MatlabException, cause::MatlabException)

    # wrap MATLAB exceptions as MxArrays
    baseException = MATLAB.MxArray(mexn.ptr, false)
    causeException = MATLAB.MxArray(cause.ptr, false)

    # call addCause to append causeException to baseException
    args = [MATLAB.MxArray(mexn.ptr), MATLAB.MxArray(cause.ptr)]
    newException = call_matlab(1, "addCause", args)[1]

    # wrap result as an exception
    return MatlabException(newException.ptr)
end

# creates a MATLABException from an exception and a backtrace
function MatlabException(exn::Exception, bt::Array{Union{Ptr{Nothing}, Base.InterpreterIP},1})

    # create MATLABException from exception
    mexn = MatlabException(exn)

    # print formatted backtrace to buffer
    buf = IOBuffer()
    print(buf, "Julia backtrace:")
    Base.show_backtrace(buf, bt)

    # move to beginning of buffer
    seek(buf, 0)

    # convert to string and format for use with MATLAB
    msg = replace(read(buf, String), "\\"=>"\\\\")

    # create cause MatlabException
    cause = MatlabException("backtrace", msg)

    # add cause to the MatlabException and return the result
    return MatlabException(mexn, cause)
end

max_stars_repo_name: taylormcd/mexjulia

max_issues_repo_name: taylormcd/mexjulia

max_forks_repo_name: taylormcd/mexjulia


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test

# Tests for @__LINE__ inside and outside of macros
@test (@__LINE__) == 6

macro macro_caller_lineno()
    @test 9 == (@__LINE__) != __source__.line > 12
    return __source__.line
end

@test @macro_caller_lineno() == (@__LINE__) > 12

# @__LINE__ in a macro expands to the location of the macrocall in the source
# while __source__.line is the location of the macro caller
macro nested_LINE_expansion()
    return quote
        return (@emit_LINE, $(__source__.line))
    end
end
macro nested_LINE_expansion2()
    return :((@emit_LINE, $(__source__.line)))
end
macro emit_LINE()
    return quote
        (@__LINE__, $(__source__.line))
    end
end
@test (@emit_LINE) == ((@__LINE__) - 3, @__LINE__)
@test @nested_LINE_expansion() == ((@__LINE__() - 4, @__LINE__() - 12), @__LINE__())
@test @nested_LINE_expansion2() == ((@__LINE__() - 5, @__LINE__() - 9), @__LINE__())

loaded_files = String[]
push!(Base.include_callbacks, (mod::Module, fn::String) -> push!(loaded_files, fn))
include("test_sourcepath.jl")
@test length(loaded_files) == 1 && endswith(loaded_files[1], "test_sourcepath.jl")
pop!(Base.include_callbacks)
thefname = "the fname!//\\&\1*"
include_string_test_func = include_string(@__MODULE__, "include_string_test() = @__FILE__", thefname)
@test include_string_test_func() == thefname
@test include_string(@__MODULE__, "Base.source_path()", thefname) == Base.source_path()
@test basename(@__FILE__) == "loading.jl"
@test isabspath(@__FILE__)

@test isdir(@__DIR__)
@test @__DIR__() == dirname(@__FILE__)
let exename = `$(Base.julia_cmd()) --compiled-modules=yes --startup-file=no`,
    wd = sprint(show, abspath(pwd(), "")),
    s_dir = sprint(show, joinpath(realpath(tempdir()), ""))
    @test wd != s_dir
    @test readchomp(`$exename -E "@__DIR__" -i`) == wd
    @test readchomp(`$exename -E "cd(()->eval(:(@__DIR__)), $s_dir)" -i`) == s_dir
    @test readchomp(`$exename -E "@__DIR__"`) == wd # non-interactive
end

# Issue #5789 and PR #13542:
mktempdir() do dir
    cd(dir) do
        let true_filename = "cAsEtEsT.jl", lowered_filename="casetest.jl"
            touch(true_filename)
            @test Base.isfile_casesensitive(true_filename)
            @test !Base.isfile_casesensitive(lowered_filename)

            # check that case-sensitivity only applies to basename of a path:
            if isfile(lowered_filename) # case-insensitive filesystem
                mkdir("cAsEtEsT")
                touch(joinpath("cAsEtEsT", true_filename))
                @test Base.isfile_casesensitive(joinpath("casetest", true_filename))
                @test !Base.isfile_casesensitive(joinpath("casetest", lowered_filename))
            end
        end

        # Test Unicode normalization; pertinent for OS X
        let nfc_name = "\U00F4.jl"
            touch(nfc_name)
            @test Base.isfile_casesensitive(nfc_name)
        end
    end
end

import Base: UUID, SHA1, PkgId, load_path, identify_package, locate_package, version_slug

saved_load_path = copy(LOAD_PATH)
saved_depot_path = copy(DEPOT_PATH)
push!(empty!(LOAD_PATH), "project")
push!(empty!(DEPOT_PATH), "depot")

@test load_path() == [abspath("project","Project.toml")]

@testset "project & manifest identify_package & locate_package" begin
    local path
    for (names, uuid, path) in [
        ("Foo",     "767738be-2f1f-45a9-b806-0234f3164144", "project/deps/Foo1/src/Foo.jl"       ),
        ("Bar.Foo", "6f418443-bd2e-4783-b551-cdbac608adf2", "project/deps/Foo2.jl/src/Foo.jl"    ),
        ("Bar",     "2a550a13-6bab-4a91-a4ee-dff34d6b99d0", "project/deps/Bar/src/Bar.jl"        ),
        ("Foo.Baz", "6801f525-dc68-44e8-a4e8-cabd286279e7", "depot/packages/9HkB/TCSb/src/Baz.jl"),
        ("Foo.Qux", "b5ec9b9c-e354-47fd-b367-a348bdc8f909", "project/deps/Qux.jl"                ),
    ]
        n = map(String, split(names, '.'))
        pkg = identify_package(n...)
        @test pkg == PkgId(UUID(uuid), n[end])
        @test joinpath(@__DIR__, normpath(path)) == locate_package(pkg)
    end
    @test identify_package("Baz") == nothing
    @test identify_package("Qux") == nothing
    @testset "equivalent package names" begin
        local classes = [
            ["Foo"],
            ["Bar", "Foo.Bar"],
            ["Foo.Baz", "Bar.Baz", "Foo.Bar.Baz"],
            ["Bar.Foo", "Foo.Bar.Foo", "Foo.Baz.Foo", "Bar.Baz.Foo"],
            ["Foo.Qux", "Foo.Baz.Qux", "Bar.Baz.Qux", "Foo.Bar.Foo.Qux",
             "Bar.Foo.Qux", "Foo.Baz.Foo.Qux", "Bar.Baz.Foo.Qux", "Foo.Bar.Baz.Foo.Qux"],
            ["Baz", "Qux", "Bar.Qux", "Bar.Baz.Bar", "Bar.Foo.Bar", "Bar.Foo.Baz",
             "Bar.Foo.Qux.Foo", "Bar.Foo.Qux.Bar", "Bar.Foo.Qux.Baz"],
        ]
        for i = 1:length(classes)
            A = classes[i]
            for x in A
                X = identify_package(map(String, split(x, '.'))...)
                for y in A
                    Y = identify_package(map(String, split(y, '.'))...)
                    @test X == Y
                end
                for j = i+1:length(classes)
                    B = classes[j]
                    for z in B
                        Z = identify_package(map(String, split(z, '.'))...)
                        @test X != Z
                    end
                end
            end
        end
    end
end

@testset "project & manifest import" begin
    @test !@isdefined Foo
    @test !@isdefined Bar
    import Foo
    @test @isdefined Foo
    @test !@isdefined Bar
    import Bar
    @test @isdefined Foo
    @test @isdefined Bar

    @testset "module graph structure" begin
        local classes = Dict(
            "Foo1" => [Foo],
            "Bar"  => [Bar, Foo.Bar],
            "Baz"  => [Foo.Baz, Bar.Baz, Foo.Bar.Baz],
            "Foo2" => [Bar.Foo, Foo.Bar.Foo, Foo.Baz.Foo, Bar.Baz.Foo],
            "Qux"  => [Foo.Qux, Foo.Baz.Qux, Bar.Baz.Qux, Foo.Bar.Foo.Qux,
                       Bar.Foo.Qux, Foo.Baz.Foo.Qux, Bar.Baz.Foo.Qux,
                       Foo.Bar.Baz.Foo.Qux],
        )
        for (i, (this, mods)) in enumerate(classes)
            for x in mods
                @test x.this == this
                for y in mods
                    @test x === y
                end
                for (j, (that, mods′)) in enumerate(classes)
                    i == j && continue
                    for z in mods′
                        @test x !== z
                    end
                end
            end
        end
    end
    @test Foo.which == "path"
end

function gen_entry_point(entry::String, pkg::PkgId)
    mkpath(dirname(entry))
    open(entry, "w") do io
        print(io, """
        __precompile__(true)
        module $(pkg.name)
        uuid = $(pkg.uuid === nothing ? "nothing" : "Base.UUID(\"$(pkg.uuid)\")")
        name = "$(pkg.name)"
        end
        """)
    end
end

function gen_project_file(project_file::String, pkg::PkgId, deps::Pair{String,UUID}...)
    mkpath(dirname(project_file))
    open(project_file, "w") do io
        println(io, "name = $(repr(pkg.name))")
        pkg.uuid !== nothing && println(io, "uuid = $(repr(string(pkg.uuid)))")
        println(io, "\n[deps]")
        for (name, uuid) in deps
            println(io, "$name = ", repr(string(uuid)))
        end
    end
end

function gen_implicit(dir::String, pkg::PkgId, proj::Bool, deps::Pair{String,UUID}...)
    entry_point = joinpath(dir, pkg.name, "src", "$(pkg.name).jl")
    gen_entry_point(entry_point, pkg)
    proj && gen_project_file(joinpath(dir, pkg.name, "Project.toml"), pkg, deps...)
    return entry_point
end

function gen_depot_ver(depot::String, pkg::PkgId, deps::Pair{String,UUID}...)
    pkg.uuid === nothing && return nothing, nothing
    tree = SHA1(rand(UInt8, 20)) # fake tree hash
    dir = joinpath(depot, "packages", version_slug(pkg.uuid, tree))
    entry = joinpath(dir, "src", "$(pkg.name).jl")
    gen_entry_point(entry, pkg)
    gen_project_file(joinpath(dir, "Project.toml"), pkg, deps...)
    return tree, entry
end

let n = 0
    global function gen_explicit(dir::String, pkg::PkgId=PkgId("Env$(n += 1)"))
        gen_project_file(joinpath(dir, "Project.toml"), pkg)
        close(open(joinpath(dir, "Manifest.toml"), "w"))
    end
end

function gen_manifest(dir::String, name::String, uuid::UUID, tree::SHA1,
    deps::Pair{String,UUID}...; toplevel::Bool = true)
    toplevel && open(joinpath(dir, "Project.toml"), "a") do io
        println(io, "$name = \"$uuid\"")
    end
    open(joinpath(dir, "Manifest.toml"), "a") do io
        println(io, "[[$name]]")
        println(io, "uuid = \"$uuid\"")
        println(io, "git-tree-sha1 = \"$(bytes2hex(tree.bytes))\"")
        if !isempty(deps)
            println(io, "    [$name.deps]")
            for (n, u) in deps
                println(io, "    $n = \"$u\"")
            end
        end
        println(io)
    end
end

const name = "Flarp"
const uuidA = UUID("b2cb3794-8625-4058-bcde-7eeb13ac1c8b")
const uuidB = UUID("1513c021-3639-4616-a37b-ee45c9d2f773")
const uuids = [nothing, uuidA, uuidB]

ft(::UUID)    =  true:true
ft(::Nothing) = false:true

@testset "direct dependency loading: implict + implicit" begin
    for uuid1 in uuids, proj1 in ft(uuid1),
        uuid2 in uuids, proj2 in ft(uuid2)
        pkg1 = uuid1 === nothing ? PkgId(name) : PkgId(uuid1, name)
        pkg2 = uuid2 === nothing ? PkgId(name) : PkgId(uuid2, name)
        empty!(LOAD_PATH)
        mktempdir() do dir1
            push!(LOAD_PATH, dir1)
            path1 = gen_implicit(dir1, pkg1, proj1)
            @test identify_package(name) == pkg1
            @test locate_package(pkg1) == path1
            path = uuid1 == coalesce(uuid2, uuid1) ? path1 : nothing
            @test locate_package(pkg2) == path
            mktempdir() do dir2
                push!(LOAD_PATH, dir2)
                path2 = gen_implicit(dir2, pkg2, proj2)
                @test identify_package(name) == pkg1
                @test locate_package(pkg1) == path1
                path = uuid1 == coalesce(uuid2, uuid1) ? path1 : path2
                @test locate_package(pkg2) == path
            end
        end
    end
end

@testset "direct dependency loading: explicit + explicit" begin
    mktempdir() do depot
        push!(empty!(DEPOT_PATH), depot)
        pkgs  = [PkgId(uuid, name) for uuid in uuids]
        pairs = [gen_depot_ver(depot, pkg) for pkg in pkgs, _ in 1:2]
        trees = first.(pairs)
        paths = last.(pairs)
        for i = 1:length(uuids), k = 1:2,
            j = 1:length(uuids), l = 1:2
            uuids[i] !== nothing && uuids[j] !== nothing || continue
            empty!(LOAD_PATH)
            mktempdir() do dir1
                push!(LOAD_PATH, dir1)
                gen_explicit(dir1)
                gen_manifest(dir1, name, uuids[i], trees[i,k])
                @test identify_package(name) == pkgs[i]
                @test locate_package(pkgs[i]) == paths[i,k]
                path = uuids[i] == uuids[j] ? paths[i,k] : nothing
                @test locate_package(pkgs[j]) == path
                mktempdir() do dir2
                    push!(LOAD_PATH, dir2)
                    gen_explicit(dir2)
                    gen_manifest(dir2, name, uuids[j], trees[j,l])
                    @test identify_package(name) == pkgs[i]
                    @test locate_package(pkgs[i]) == paths[i,k]
                    path = uuids[i] == uuids[j] ? paths[i,k] : paths[j,l]
                    @test locate_package(pkgs[j]) == path
                end
            end
        end
    end
end

@testset "direct dependency loading: explicit + implicit" begin
    mktempdir() do depot
        push!(empty!(DEPOT_PATH), depot)
        pkgs  = [PkgId(uuid, name) for uuid in uuids]
        pairs = [gen_depot_ver(depot, pkg) for pkg in pkgs, _ in 1:2]
        trees = first.(pairs)
        paths = last.(pairs)
        for i = 1:length(uuids), k = 1:2,
            j = 1:length(uuids), l = 1:2, proj in ft(uuids[j])
            uuids[i] !== nothing || continue
            empty!(LOAD_PATH)
            mktempdir() do dir1
                push!(LOAD_PATH, dir1)
                gen_explicit(dir1)
                gen_manifest(dir1, name, uuids[i], trees[i,k])
                @test identify_package(name) == pkgs[i]
                @test locate_package(pkgs[i]) == paths[i,k]
                path = uuids[i] == coalesce(uuids[j], uuids[i]) ? paths[i,k] : nothing
                @test locate_package(pkgs[j]) == path
                mktempdir() do dir2
                    push!(LOAD_PATH, dir2)
                    path2 = gen_implicit(dir2, pkgs[j], proj)
                    @test identify_package(name) == pkgs[i]
                    @test locate_package(pkgs[i]) == paths[i,k]
                    path = uuids[i] == coalesce(uuids[j], uuids[i]) ? paths[i,k] : path2
                    @test locate_package(pkgs[j]) == path
                end
            end
        end
    end
end

const uuidT = UUID("a54bd003-d8dc-4161-b186-d5516cd448e9")

@testset "indirect dependency loading: explicit + explicit" begin
    mktempdir() do depot
        push!(empty!(DEPOT_PATH), depot)
        # generate top-level package
        top = PkgId(uuidT, "TopLevel")
        top_tree, _ = gen_depot_ver(depot, top)
        # generate dependency packages
        pkgs  = [PkgId(uuid, name) for uuid in uuids]
        pairs = [gen_depot_ver(depot, pkg) for pkg in pkgs, _ in 1:2]
        trees = first.(pairs)
        paths = last.(pairs)
        for i = 1:length(uuids), k = 1:2, s = false:true,
            j = 1:length(uuids), l = 1:2, t = false:true
            uuids[i] !== nothing && uuids[j] !== nothing || continue
            empty!(LOAD_PATH)
            mktempdir() do dir1
                push!(LOAD_PATH, dir1)
                gen_explicit(dir1)
                gen_manifest(dir1, name, uuids[i], trees[i,k], toplevel=false)
                s && gen_manifest(dir1, top.name, top.uuid, top_tree, name => uuids[i])
                @test identify_package(top, name) == (s ? pkgs[i] : nothing)
                @test locate_package(pkgs[i]) == paths[i,k]
                path = uuids[i] == uuids[j] ? paths[i,k] : nothing
                @test locate_package(pkgs[j]) == path
                mktempdir() do dir2
                    push!(LOAD_PATH, dir2)
                    gen_explicit(dir2)
                    t && gen_manifest(dir1, top.name, top.uuid, top_tree, name => uuids[j])
                    gen_manifest(dir2, name, uuids[j], trees[j,l], toplevel=false)
                    pkg = (s ? pkgs[i] : t ? pkgs[j] : nothing)
                    @test identify_package(top, name) == pkg
                    @test locate_package(pkgs[i]) == paths[i,k]
                    path = uuids[i] == uuids[j] ? paths[i,k] : paths[j,l]
                    @test locate_package(pkgs[j]) == path
                end
            end
        end
    end
end

append!(empty!(DEPOT_PATH), saved_depot_path)
append!(empty!(LOAD_PATH), saved_load_path)

max_stars_repo_name: nreigl/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
"""
    BlackScholes()

[Black-Scholes model](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model).
"""
struct BlackScholes <: Model end

"""
    evaluate(O, BlackScholes())

Evaluate option `O` using `BlackScholes` model.

# Arguments
- `O::Option`: Option
"""
function evaluate(O::EuropeanOption, m::BlackScholes)
    d1 = (log(O.s / O.k) + (O.r + O.σ * O.σ / 2) * O.t) / (O.σ * √O.t)
    d2 = d1 - O.σ * √O.t

    if O.call == -1
        return cdf(Normal(), -d2) * O.k * exp(-O.r * O.t) - cdf(Normal(), -d1) * O.s
    elseif O.call == 1
        return O.s * cdf(Normal(), d1) - exp(-O.r * O.t) * O.k * cdf(Normal(), d2)
    end
end

"""
    evaluate(o)

Evaluate option `o` using Back-Scholes model as default valuation model.
"""
evaluate(O::EuropeanOption) = evaluate(O::EuropeanOption, m::BlackScholes)


text:
function plot(plon,plat,data; cmap = "bone", coastlinecolor="red")
	ccrs = pyimport("cartopy.crs")
	ax = subplot(projection=ccrs.EqualEarth())
	#ax = subplot(projection=ccrs.Orthographic(-10, 45))
	ax.outline_patch.set_linewidth(0.1)
	#ax.set_extent([-180, 180, -90, 90])
	ax.set_global()
	pcolormesh(plon,plat,data,cmap=cmap,shading="auto",transform=ccrs.PlateCarree())
	#ax.coastlines(resolution="50m",linewidth=0.1)
	return nothing
end

"""
	makeplots(wavname,satellite_name)

Decodes the APT signal in `wavname` as recorded by gqrx using
wide FM-mono demodulation.
The file name `wavname` should  have the followng structure:
`string_date_time_frequency.wav` like `gqrx_20190811_075102_137620000.wav`.
Date and time of the file name are in UTC (not local time). `satellite_name` is
the name of the satellite (generally `"NOAA 15"`, `"NOAA 18"`, `"NOAA 19"`).

# Example:

```julia
import OpenWeatherDecoder

wavname = "gqrx_20190825_182745_137620000.wav"
OpenWeatherDecoder.makeplots(wavname,"NOAA 15", coastlinecolor = "red", cmap=)
```

"""

function decodeall(;
				   starttime = nothing,
				   eop = nothing,
				   qrange = (0.01,0.99),
				   coastlinecolor = "red",
				   cmap = "bone",
				   tles = get_tle(:weather),
				   dpi = 2200)

	NOAAFreqs = Dict(
		"1371" => "NOAA 19",
		"1379" => "NOAA 18",
		"1376" => "NOAA 15"
	)

	foreach(glob("*.wav","files")) do f
		freqmatch = match(r".*_.*_(.*).wav", f)
		prefix = replace(f,r".wav$" => "")
		freq = freqmatch[1]

		subfreq = SubString(freq, 1, 4)
		satellite_name = get(NOAAFreqs, subfreq, 0)

		OpenWeatherDecoder.decodewav(f, satellite_name,
			starttime=starttime,
			eop=eop,
			prefix=prefix,
			qrange=qrange,
			coastlinecolor=coastlinecolor,
			cmap=cmap,
			tles=tles,
			dpi=dpi)
	end

end

function plotall(;
				   starttime = nothing,
				   eop = nothing,
				   qrange = (0.01,0.99),
				   coastlinecolor = "red",
				   cmap = "bone",
				   tles = get_tle(:weather),
				   dpi = 2200)

	NOAAFreqs = Dict(
		"1371" => "NOAA 19",
		"1379" => "NOAA 18",
		"1376" => "NOAA 15"
	)

	foreach(glob("*.wav","files")) do f
		freqmatch = match(r".*_.*_(.*).wav", f)
		prefix = replace(f,r".wav$" => "")
		freq = freqmatch[1]

		subfreq = SubString(freq, 1, 4)
		satellite_name = get(NOAAFreqs, subfreq, 0)

		OpenWeatherDecoder.makeplots(f, satellite_name,
			starttime=starttime,
			eop=eop,
			prefix=prefix,
			qrange=qrange,
			coastlinecolor=coastlinecolor,
			cmap=cmap,
			tles=tles,
			dpi=dpi)
	end

end

function decodewav(wavname,satellite_name;
				   starttime = nothing,
				   eop = nothing,
				   prefix = replace(wavname,r".wav$" => ""),
				   qrange = (0.01,0.99),
				   coastlinecolor = "red",
				   cmap = "bone",
				   tles = get_tle(:weather),
				   dpi = 2200)
	println(prefix)
	if starttime == nothing
		starttime = OpenWeatherDecoder.starttimename(basename(wavname))
	end

	y,Fs,nbits,opt = FileIO.load(wavname)

	datatime,(channelA,channelB),data = OpenWeatherDecoder.decode(y,Fs)

	vmin,vmax = quantile(view(data,:),[qrange[1],qrange[2]])
	data[data .> vmax] .= vmax;
	data[data .< vmin] .= vmin;

	# save raw image
	imagenames = (
		rawname = prefix * ".png",
		channel_a = prefix * "_channel_a.png",
		channel_b = prefix * "_channel_b.png")

	img = colorview(Gray, data[:,1:3:end]./maximum(data))

	hist_adapt = adjust_histogram(img, AdaptiveEqualization(nbins = 256, rblocks = 4, cblocks = 4, clip = 0.2))

	FileIO.save(imagenames.rawname, hist_adapt)

	return imagenames
end

function makeplots(wavname,satellite_name;
				   starttime = nothing,
				   eop = nothing,
				   prefix = replace(wavname,r".wav$" => ""),
				   qrange = (0.01,0.99),
				   coastlinecolor = "red",
				   cmap = "bone",
				   tles = get_tle(:weather),
				   dpi = 2200)

	imagenames = (
		rawname = prefix * ".png",
		channel_a = prefix * "_channel_a.png",
		channel_b = prefix * "_channel_b.png")
	fig = figure()

	Alon,Alat,Adata = OpenWeatherDecoder.georeference(replace(prefix, "files/" => "") * ".png",satellite_name,'a')

	OpenWeatherDecoder.plot(Alon,Alat,Adata; coastlinecolor=coastlinecolor, cmap=cmap)
	savefig(imagenames.channel_a,dpi=dpi,pad_inches=0, bbox_inches="tight", transparent=true)
	fig.clf()

	# Blon,Blat,Bdata = OpenWeatherDecoder.georeference(
	#     channelB,satellite_name,datatime,starttime, eop = eop, tles = tles)
	# OpenWeatherDecoder.plot(Blon,Blat,Bdata; coastlinecolor=coastlinecolor, cmap=cmap)
	# plt.grid(linestyle = "--",color=grid_color)
	# savefig(imagenames.channel_b,dpi=dpi,pad_inches=0, bbox_inches="tight", transparent=false)
	# close(fig)
	return imagenames
end

max_stars_repo_name: rctngle/OpenWeatherDecoder

max_issues_repo_name: rctngle/OpenWeatherDecoder

max_forks_repo_name: rctngle/OpenWeatherDecoder


text:
## THIS IS A SCRATCH-WORK FILE  -- PAVEL ##

# matrix version?
function QuellerGoodnight2(data::PopData, ind1::String, ind2::String; alleles::T) where T <: NamedTuple

    #TODO NEED TO CHECK TO CONFIRM EQUATIONS

    numerator1 = 0.0
    numerator2 = 0.0
    denominator1 = 0.0
    denominator2 = 0.0

    geno1 = genotypes(data, ind1)
    geno2 = genotypes(data, ind2)

    # drop missing values
    loc, geno1, geno2 = collect.(skipmissings(Symbol.(loci(data)), geno1, geno2))

    # convert to 2-dimensional matrix for each individual
    gen1 = hcat(getindex.(geno1, 1), getindex.(geno1, 2))
    gen2 = hcat(getindex.(geno2, 1), getindex.(geno2, 2))
    # ((a == c) + (a == d) + (b == c) + (b == d))
    id_mtx = reduce(*, (gen1 .== gen2) + (gen1 .== reverse(gen2, dims = 2)), dims = 2)
    a_b_frqs = [alleles[loc[i]][gen1[i,j]] for i in 1:size(gen1)[1], j in 1:2]
    c_d_frqs = [alleles[loc[i]][gen2[i,j]] for i in 1:size(gen2)[1], j in 1:2]

    numerator1 = reduce(+, id_mtx - (2 * reduce(+, a_b_frqs, dims = 2)))
    numerator2 = reduce(+, id_mtx - (2 * reduce(+, c_d_frqs, dims = 2)))
    # reduce(+,)
    denominator1 = 2.0 * (1.0 + map(x ->x[1] == x[2],eachrow(gen1)) - foldl(-, eachcol(a_b_frqs)))
    denominator2 = 2.0 * (1.0 + map(x ->x[1] == x[2],eachrow(gen2)) - foldl(-, eachcol(a_b_frqs)))
    #return [alleles[loc[i]][gen1[i,j]] for i in 1:size(gen1)[1], j in 1:2]
    return numerator2
    a,b = gen1
    c,d = gen2

        ident = ((a == c) + (a == d) + (b == c) + (b == d))
        numerator1 += ident - 2.0 * (alleles[loc][a] + alleles[loc][b])
        numerator2 += ident - 2.0 * (alleles[loc][c] + alleles[loc][d])

        denominator1 += (2.0 * (1.0 + (a==b) - alleles[loc][a] - alleles[loc][b]))
        denominator2 += (2.0 * (1.0 + (c==d) - alleles[loc][c] - alleles[loc][d]))

    return (numerator1/denominator1 + numerator2/denominator2)/2.0
end

function _relatedness_noboot_fix(data::PopData, sample_names::Vector{String}; method::F, inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    n_samples = data.metadata.samples
    sample_pairs = pairwisepairs(sample_names)
    n_per_loci = DataFrames.combine(groupby(data.genodata, :locus), :genotype => nonmissing => :n)[:, :n]
    allelefrequencies = allelefreq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.genodata, :name)
    idx = Base.Threads.Atomic{Int}()
    @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
        @inbounds geno1 = popdata_idx[(ind1,)].genotype
        @inbounds @sync Base.Threads.@threads for ind2 in sample_names[sample_n+1:end]
            Base.Threads.atomic_add!(idx, 1)
            @inbounds geno2 = popdata_idx[(ind2,)].genotype

            # get index for genotype appearing missing in at least one individual in the pair
            keep_idx = nonmissings(geno1, geno2)

            # populate shared_loci array
            @inbounds shared_loci[idx.value] = length(keep_idx)
            @inbounds [relate_vecs[i][idx.value] = @inbounds mth(geno1[keep_idx], geno2[keep_idx], loci_names[keep_idx], allelefrequencies, loc_n = n_per_loci[keep_idx], n_samples = n_samples, inbreeding = inbreeding) for (i,mth) in enumerate(method)]

            update!(p, idx.value)
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    out_df = DataFrame(:sample_1 => getindex.(sample_pairs, 1), :sample_2 => getindex.(sample_pairs, 2), :n_loci => shared_loci)
    [out_df[:, mth] = relate_vecs[i] for (i, mth) in enumerate(method_colnames)]
    return out_df
end

function relatedness2(data::PopData, sample_names::Vector{String}; method::F, iterations::Int64 = 0, interval::Tuple{Float64, Float64} = (0.025, 0.975), resample::String = "all", inbreeding::Bool = false) where F
    all(data.metadata[data.metadata.sampleinfo.name .∈ Ref(sample_names), :ploidy] .== 2) == false && error("Relatedness analyses currently only support diploid samples")
    errs = ""
    all_samples = samplenames(data)
    if sample_names != all_samples
        [errs *= "$i," for i in sample_names if i ∉ all_samples]
        errs != "" && error("Samples not found in the PopData: " * errs)
    end
    if eltype(method) != Function
        method = [method]
    end
    relate_mthds = [:QuellerGoodnight, :Ritland, :Lynch, :LynchLi, :LynchRitland, :Wang, :Loiselle, :Blouin, :Moran, :LiHorvitz, :dyadicLikelihood]
    [errs *= "$i is not a valid method\n" for i in Symbol.(method) if i ∉ relate_mthds]
    errs != "" && throw(ArgumentError(errs * "Methods are case-sensitive. Please see the docstring (?relatedness) for additional help."))
    if iterations > 0
        if resample == "all"
            _relatedness_boot_all(data, sample_names, method = method, iterations = iterations, interval = interval, inbreeding = inbreeding)
        elseif resample == "nonmissing"
            _relatedness_boot_nonmissing(data, sample_names, method = method, iterations = iterations, interval = interval, inbreeding = inbreeding)
        else
            throw(ArgumentError("Please choose from resample methods \"all\" or \"nonmissing\""))
        end
    else
        _relatedness_noboot_fix2(data, sample_names, method = method, inbreeding = inbreeding)
    end
end

function relatedness2(data::PopData; method::F, iterations::Int64 = 0, interval::Tuple{Float64, Float64} = (0.025, 0.975), resample::String = "all", inbreeding::Bool = false) where F
    sample_names = samplenames(data) |> collect
    relatedness2(data, sample_names, method = method, iterations = iterations, interval = interval, resample = resample, inbreeding = inbreeding)
end

## original function -- threading does not work as expected ###
function _relatedness_noboot(data::PopData, sample_names::Vector{String}; method::F, inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    n_samples = data.metadata.samples
    sample_pairs = pairwisepairs(sample_names)
    n_per_loci = DataFrames.combine(groupby(data.genodata, :locus), :genotype => nonmissing => :n)[:, :n]
    allelefrequencies = allelefreq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.genodata, :name)
    idx = 0
    @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
        @inbounds geno1 = popdata_idx[(ind1,)].genotype
        @inbounds @sync Base.Threads.@spawn for ind2 in sample_names[sample_n+1:end]
            idx += 1
            @inbounds geno2 = popdata_idx[(ind2,)].genotype

            # get index for genotype appearing missing in at least one individual in the pair
            keep_idx = nonmissings(geno1, geno2)

            # populate shared_loci array
            @inbounds shared_loci[idx] = length(keep_idx)
            @inbounds [relate_vecs[i][idx] = @inbounds mth(geno1[keep_idx], geno2[keep_idx], loci_names[keep_idx], allelefrequencies, loc_n = n_per_loci[keep_idx], n_samples = n_samples, inbreeding = inbreeding) for (i,mth) in enumerate(method)]

            update!(p, idx)
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    out_df = DataFrame(:sample_1 => getindex.(sample_pairs, 1), :sample_2 => getindex.(sample_pairs, 2), :n_loci => shared_loci)
    [out_df[:, mth] = relate_vecs[i] for (i, mth) in enumerate(method_colnames)]
    return out_df
end

function _relatedness_boot_nonmissing(data::PopData, sample_names::Vector{String}; method::F, iterations::Int, interval::Tuple{Float64, Float64} = (0.025, 0.975), inbreeding::Bool) where F
    loci_names = Symbol.(loci(data))
    sample_pairs = pairwisepairs(sample_names)
    n_samples = data.metadata.samples
    n_per_loci = DataFrames.combine(groupby(data.genodata, :locus), :genotype => nonmissing => :n)[:, :n]
    allelefrequencies = allelefreq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.genodata, :name)
    idx = 0
    @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
        @inbounds geno1 = popdata_idx[(ind1,)].genotype
        @inbounds @sync Base.Threads.@spawn for ind2 in sample_names[sample_n+1:end]
            idx += 1
            @inbounds geno2 = popdata_idx[(ind2,)].genotype

            # get index for genotype appearing missing in at least one individual in the pair
            keep_idx = nonmissings(geno1, geno2)
            # generate nonmissing genotype data
            gen1, gen2, loc, n_per_loc = (i[keep_idx] for i in (geno1, geno2, loci_names, n_per_loci))

            # populate shared_loci array
            @inbounds shared_loci[idx] = length(loc)

            @inbounds for (i, mthd) in enumerate(method)
                @inbounds relate_vecs[i][idx] = mthd(gen1, gen2, loc, allelefrequencies, loc_n = n_per_loc, n_samples = n_samples, inbreeding = inbreeding)
                boot_out = _bootstrapgenos_nonmissing(gen1, gen2, loc, n_per_loc, allelefrequencies, method = mthd, iterations = iterations, inbreeding = inbreeding)
                @inbounds boot_means[i][idx], boot_medians[i][idx], boot_ses[i][idx], boot_CI[i][idx] = _bootstrapsummary(boot_out, interval)
            end
            update!(p, idx)
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_df = DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end

    return out_df
end

# method using kwargs and Symbols as methods. Incomplete.
function relatedness(data::PopData, sample_names::Vector{String}; kwargs...)
    #method::F, iterations::Int64 = 0, interval::Tuple{Float64, Float64} = (0.025, 0.975), resample::String = "all") where F
    kw_dict = Dict(kwargs...)
    all(data.metadata[data.metadata.sampleinfo.name .∈ Ref(sample_names), :ploidy] .== 2) == false && error("Relatedness analyses currently only support diploid samples")
    all_samples = samplenames(data)
    errs = ""
    if sample_names != all_samples
        [errs *= "$i," for i in sample_names if i ∉ all_samples]
        errs != "" && error("Samples not found in the PopData: " * errs)
    end
    method = kw_dict[:method]
    if !(typeof(method) <: AbstractVector)
        method = [method]
    end
    !(haskey(kw_dict, :method)) && throw(ArgumentError("Please supply the keyword argument: method. See the docstring (?relatedness) for additional help"))
    relate_mthds = [:QuellerGoodnight, :Ritland, :Lynch, :LynchLi, :LynchRitland, :Wang, :Loiselle, :Blouin, :Moran, :LiHorvitz, :dyadicLikelihood]
    [errs *= "$i is not a valid method\n" for i in method if i ∉ relate_mthds]
    errs != "" && throw(ArgumentError(errs * "Methods are case-sensitive. Please see the docstring (?relatedness) for additional help."))
    if haskey(kw_dict, :iterations)
        if kw_dict[:iterations] > 0
            if haskey(kw_dict, :resample)
                "all" != kw_dict[:resample] != "nonmissing" && throw(ArgumentError("Please choose from resample methods \"all\" or \"nonmissing\""))
                kw_dict[:resample] == "all" && _relatedness_boot_all(data, sample_names, method = method, iterations = iterations, interval = interval)
                kw_dict[:resample] == "nonmissing" && _relatedness_boot_nonmissing(data, sample_names, method = method, iterations = iterations, interval = interval)
            else
                _relatedness_boot_all(data, sample_names, method = method, iterations = iterations, interval = interval)
            end
    else
        _relatedness_noboot(data, sample_names, method = method)
    end
end

#= functions as Symbols for arguments

julia> function do_thing(x, operation::Symbol)
           getfield(@__MODULE__, operation)(x)
       end

do_thing (generic function with 1 method)

julia> do_thing(1, :sin)
0.8414709848078965

julia> do_thing(1, :sdfgsdfsdf)
ERROR: UndefVarError: sdfgsdfsdf not defined
Stacktrace:
 [1] do_thing(::Int64, ::Symbol) at ./REPL[26]:2
 [2] top-level scope at REPL[28]:1
=#

# extremely incomplete
function relatedness_bootstrap(data::PopData, sample_names::Vector{String}; method::F, iterations::Int = 100, interval::Tuple{Float64, Float64} = (0.025, 0.975)) where F
    loci_names = Symbol.(loci(data))
    sample_pairs = pairwisepairs(sample_names)
    n_loci = length(loci_names)
    n_samples = data.metadata.samples
    n_per_loci = DataFrames.combine(groupby(data.genodata, :locus), :genotype => nonmissing => :n)[:, :n]
    allelefrequencies = allelefreq(data)
    relate_vecs = map(i -> Vector{Union{Missing,Float64}}(undef, length(sample_pairs)), 1:length(method))
    boot_means, boot_medians, boot_ses = map(i -> deepcopy(relate_vecs), 1:3)
    boot_CI = map(i -> Vector{Union{Missing,Tuple{Float64,Float64}}}(undef, length(sample_pairs)), 1:length(method))
    shared_loci = Vector{Int}(undef, length(sample_pairs))
    p = Progress(length(sample_pairs), dt = 1, color = :blue)
    popdata_idx = groupby(data.genodata, :name)
    idx = 0
    for i in 1:iterations
        boot_idx = rand(1:n_loci, n_loci)
        @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
            geno1 = popdata_idx[(ind1,)].genotype
            @inbounds @sync Base.Threads.@spawn for ind2 in sample_names[sample_n+1:end]
                idx += 1
                geno2 = popdata_idx[(ind2,)].genotype

                ind1_boot, ind2_boot, loc_boot, n_per_loc = map(i -> getindex(i, boot_idx), [ind1, ind2, loci_names, n_per_loc])
                # faster/cheaper n counting
                #loc, gen1, gen2, n_per_loc_sm = collect.(skipmissings(loc_boot, ind1_boot, ind2_boot, n_per_loci))
                # filter out loci missing in at least one individual in the pair
                loc, gen1, gen2, n_per_loc_sm = collect.(skipmissings(loci_names[boot_idx], geno1[boot_idx], geno2[boot_idx], n_per_loci[boot_idx]))

                # populate shared_loci array
                @inbounds shared_loci[idx] = length(loc)
                @inbounds for (i, mthd) in enumerate(method)
                    @inbounds relate_vecs[i][idx] = mthd(gen1, gen2, loc, allelefrequencies, loc_n = n_per_loc, n_samples = n_samples)
                    boot_out = mthd(gen1, gen2, loc, allelefrequencies, loc_n = n_per_loc, n_samples = n_samples)
                    #boot_out = bootstrap_locus(data, mthd, ind1, ind2, iterations, allelefrequencies)
                    @inbounds boot_means[i][idx], boot_medians[i][idx], boot_ses[i][idx], boot_CI[i][idx] = _bootstrapsummary(boot_out, iterations, interval)
                end
            update!(p, idx)
        end
    end
    method_colnames = [Symbol("$i") for i in method]
    boot_mean_colnames = [Symbol("$i"*"_mean") for i in method]
    boot_median_colnames = [Symbol("$i"*"_median") for i in method]
    boot_se_colnames = [Symbol("$i"*"_SE") for i in method]
    CI_percent = convert(Int64, round(interval[2] - interval[1], digits = 2) * 100)
    boot_CI_colnames = [Symbol("$i"*"_CI_"*"$CI_percent") for i in method]

    out_df = DataFrame(:sample_1 => map(i -> i[1], sample_pairs), :sample_2 => map(i -> i[2], sample_pairs), :n_loci => shared_loci)
    @inbounds for (i, mth) in enumerate(method_colnames)
        out_df[:, mth] = relate_vecs[i]
        out_df[:, boot_mean_colnames[i]] = boot_means[i]
        out_df[:, boot_median_colnames[i]] = boot_medians[i]
        out_df[:, boot_se_colnames[i]] = boot_ses[i]
        out_df[:, boot_CI_colnames[i]] = boot_CI[i]
    end

    return out_df
end


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

mutable struct DiskEncryptionSetUpdate <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, DiskEncryptionSetUpdateProperties } # spec name: properties
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function DiskEncryptionSetUpdate(;properties=nothing, tags=nothing)
        o = new()
        validate_property(DiskEncryptionSetUpdate, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(DiskEncryptionSetUpdate, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type DiskEncryptionSetUpdate

const _property_map_DiskEncryptionSetUpdate = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"), Symbol("tags")=>Symbol("tags"))
const _property_types_DiskEncryptionSetUpdate = Dict{Symbol,String}(Symbol("properties")=>"DiskEncryptionSetUpdateProperties", Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ DiskEncryptionSetUpdate }) = collect(keys(_property_map_DiskEncryptionSetUpdate))
Swagger.property_type(::Type{ DiskEncryptionSetUpdate }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DiskEncryptionSetUpdate[name]))}
Swagger.field_name(::Type{ DiskEncryptionSetUpdate }, property_name::Symbol) =  _property_map_DiskEncryptionSetUpdate[property_name]

function check_required(o::DiskEncryptionSetUpdate)
    true
end

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


text:
using DataFrames

export solve_poisson!

"""
   solve_poisson!( efield, particle_group, kernel_smoother, maxwell_solver, rho )

Accumulate rho and solve Poisson
 - `particle_group` : Particles
 - `maxwell_solver` : Maxwell solver (FEM 1D)
 - `kernel_smoother_0` : Particle-Mesh method
 - `rho` : preallocated array for Charge density
 - `efield_dofs` : spline coefficients of electric field (1D)
"""
function solve_poisson!(
    efield_dofs::Vector{Float64},
    particle_group::AbstractParticleGroup,
    kernel_smoother_0::AbstractParticleMeshCoupling,
    maxwell_solver::Maxwell1DFEM,
    rho::Vector{Float64},
)
    fill!(rho, 0.0)

    for i_part in 1:(particle_group.n_particles)
        xi = particle_group.array[1, i_part]
        wi = get_charge(particle_group, i_part)
        add_charge!(rho, kernel_smoother_0, xi, wi)
    end

    return compute_e_from_rho!(efield_dofs, maxwell_solver, rho)
end

"""
    pic_diagnostics_transfer( particle_group, kernel_smoother_0,
                            kernel_smoother_1, efield_dofs, transfer)

Compute ``\\sum_{particles} w_p ( v_1,p e_1(x_p) + v_2,p e_2(x_p)) ``

- `particle_group`
- `kernel_smoother_0`  : Kernel smoother (order p+1)
- `kernel_smoother_1`  : Kernel smoother (order p)
- `efield_dofs` : coefficients of efield

"""
function pic_diagnostics_transfer(
    particle_group::ParticleGroup{1,2}, kernel_smoother_0, kernel_smoother_1, efield_dofs
)
    transfer = 0.0::Float64

    for i_part in 1:(particle_group.n_particles)
        xi = particle_group.array[1, i_part]
        wi = get_charge(particle_group, i_part)
        v1 = particle_group.array[2, i_part]
        v2 = particle_group.array[3, i_part]

        efield_1 = evaluate(kernel_smoother_1, xi, efield_dofs[1])
        efield_2 = evaluate(kernel_smoother_0, xi, efield_dofs[2])

        transfer += (v1 * efield_1 + v2 * efield_2) * wi
    end

    return transfer
end

"""
    pic_diagnostics_vvb( particle_group, kernel_smoother_1, bfield_dofs )

Compute ``\\sum_{particles} ( w_p v_1, p b(x_p) v_2, p )``

- `particle_group`    : particle group object
- `kernel_smoother_1` : Kernel smoother (order p)
- `bfield_dofs` : coefficients of bfield

"""
function pic_diagnostics_vvb(particle_group, kernel_smoother_1, bfield_dofs)
    vvb = 0.0::Float64

    for i_part in 1:(particle_group.n_particles)
        xi = particle_group.array[1, i_part]
        v1 = particle_group.array[2, i_part]
        v2 = particle_group.array[3, i_part]
        wi = particle_group.array[4, i_part]
        wi *= particle_group.charge
        wi *= particle_group.common_weight

        bfield = evaluate(kernel_smoother_1, xi[1], bfield_dofs)

        vvb += wi * v1 * v2 * bfield
    end

    return vvb
end

"""
    pic_diagnostics_poynting( maxwell_solver, degree, efield_dofs, bfield_dofs,
                              rho )

Compute ``e^T M_0^{-1}  R^T b``

- `maxwell_solver` : maxwell solver object
- `degree` : degree of finite element
- `efield_dofs` : coefficients of `efield`
- `bfield_dofs` : coefficients of `bfield`

"""
function pic_diagnostics_poynting(maxwell_solver, degree, efield_dofs, bfield_dofs)
    scratch = similar(bfield_dofs)
    # Multiply B by M_0^{-1}  R^T
    compute_e_from_b!(scratch, maxwell_solver, 1.0, bfield_dofs)

    return inner_product(maxwell_solver, efield_dofs, scratch, degree)
end

export TimeHistoryDiagnostics

"""
    TimeHistoryDiagnostics( particle_group, maxwell_solver,
                            kernel_smoother_0, kernel_smoother_1 )

Context to save and plot diagnostics

- `particle_group` : Particles data
- `maxwell_solver` : Maxwell solver
- `kernel_smoother_0` : Mesh coupling operator
- `kernel_smoother_1` : Mesh coupling operator
- `data` : DataFrame containing time history values
"""
struct TimeHistoryDiagnostics
    particle_group::ParticleGroup
    maxwell_solver::Maxwell1DFEM
    kernel_smoother_0::ParticleMeshCoupling1D
    kernel_smoother_1::ParticleMeshCoupling1D
    data::DataFrame
    diagnostics::Vector{Float64}
    potential_energy::Vector{Float64}

    function TimeHistoryDiagnostics(
        particle_group::ParticleGroup,
        maxwell_solver::Maxwell1DFEM,
        kernel_smoother_0::ParticleMeshCoupling1D,
        kernel_smoother_1::ParticleMeshCoupling1D,
    )
        data = DataFrame(;
            Time=Float64[],
            KineticEnergy=Float64[],
            Momentum1=Float64[],
            Momentum2=Float64[],
            PotentialEnergyE1=Float64[],
            PotentialEnergyE2=Float64[],
            PotentialEnergyB3=Float64[],
            Transfer=Float64[],
            VVB=Float64[],
            Poynting=Float64[],
            ErrorPoisson=Float64[],
        )

        diagnostics = zeros(Float64, 3)
        potential_energy = zeros(Float64, 3)

        return new(
            particle_group,
            maxwell_solver,
            kernel_smoother_0,
            kernel_smoother_1,
            data,
            diagnostics,
            potential_energy,
        )
    end
end

export write_step!

"""
    write_step!( thdiag, time, degree, efield_dofs, bfield_dofs,
                 efield_dofs_n, efield_poisson)

write diagnostics for PIC
- `time` : Time
- `efield_dofs` : Electric field
- `efield_dofs_n` : Electric field at half step
- `efield_poisson` : Electric field compute from Poisson equation
- `bfield_dofs` : Magnetic field
- `degree` : Spline degree
"""
function write_step!(
    thdiag::TimeHistoryDiagnostics,
    time,
    degree,
    efield_dofs,
    bfield_dofs,
    efield_dofs_n,
    efield_poisson,
)
    fill!(thdiag.diagnostics, 0.0)
    fill!(thdiag.potential_energy, 0.0)

    for i_part in 1:(thdiag.particle_group.n_particles)
        v1 = thdiag.particle_group.array[2, i_part]
        v2 = thdiag.particle_group.array[3, i_part]
        wi = thdiag.particle_group.array[4, i_part]
        wi *= thdiag.particle_group.mass
        wi *= thdiag.particle_group.common_weight

        # Kinetic energy
        thdiag.diagnostics[1] += (v1^2 + v2^2) * wi
        # Momentum 1
        thdiag.diagnostics[2] += v1 * wi
        # Momentum 2
        thdiag.diagnostics[3] += v2 * wi
    end

    transfer = pic_diagnostics_transfer(
        thdiag.particle_group,
        thdiag.kernel_smoother_0,
        thdiag.kernel_smoother_1,
        efield_dofs,
    )

    vvb = pic_diagnostics_vvb(thdiag.particle_group, thdiag.kernel_smoother_1, bfield_dofs)

    poynting = pic_diagnostics_poynting(
        thdiag.maxwell_solver, degree, efield_dofs[2], bfield_dofs
    )

    thdiag.potential_energy[1] = inner_product(
        thdiag.maxwell_solver, efield_dofs[1], efield_dofs_n[1], degree - 1
    )

    thdiag.potential_energy[2] = inner_product(
        thdiag.maxwell_solver, efield_dofs[2], efield_dofs_n[2], degree
    )

    thdiag.potential_energy[3] = l2norm_squared(
        thdiag.maxwell_solver, bfield_dofs, degree - 1
    )

    return push!(
        thdiag.data,
        (
            time,
            thdiag.diagnostics...,
            thdiag.potential_energy...,
            transfer,
            vvb,
            poynting,
            maximum(abs.(efield_dofs[1] .- efield_poisson)),
        ),
    )
end

"""
    evaluate( kernel_smoother, field_dofs,  xi, n_dofs )

Evaluate the field at points xi

- `field_dofs` : field value on dofs
- `xi` : positions where the field is evaluated
"""
function evaluate(
    kernel_smoother::ParticleMeshCoupling1D, field_dofs::AbstractArray, x::AbstractArray
)
    field_grid = similar(x)
    for j in eachindex(x)
        field_grid[j] = evaluate(kernel_smoother, x[j], field_dofs)
    end
    return field_grid
end

"""
    pic_diagnostics_hpi( particle_group,  index, kinetic )

compute v(index)-part of kinetic energy

- `particle_group`
- `index` : velocity component
- `kinetic` : value of index part of kinetic energy

"""
function pic_diagnostics_hpi(particle_group, index, kinetic)
    kinetic = 0.0
    for i_part in 1:(particle_group.n_particles)
        vi = get_v(particle_group, i_part)
        wi = get_mass(particle_group, i_part)
        kinetic += kinetic_local + (vi[index]^2) * wi[1]
    end

    return kinetic
end

"""
    eval_derivative_spline( position, xmin, delta_x, n_grid,
                            field_dofs, degree, derivative )

Compute the spline coefficient of the derivative of some given spline expansion

- `position` : particle position
- `xmin` : lower boundary of the domain
- `delta_x` : step
- `n_grid` : number of grid points
- `field_dofs` : coefficients of spline representation of the field
- `degree` : degree of spline
- `derivative` : value of the derivative
"""
function eval_derivative_spline(position, xmin, delta_x, n_grid, field_dofs, degree)
    der_degree = degree - 1

    xi = (position[1] - xmin) / delta_x
    index = ceil(xi)
    xi = xi - (index - 1)
    index = index - der_degree

    spline_val = uniform_bsplines_eval_basis(der_degree, xi)

    derivative = 0.0

    for i1 in 1:degree
        ind = mod(index + i1 - 2, n_grid) + 1
        derivative +=
            spline_val[i1] * (field_dofs[ind] - field_dofs[mod(ind - 2, n_grid) + 1])
    end

    return derivative / delta_x
end


text:
@kwdef mutable struct FxpMultHandler <: AbstractHandler
    id = 0
end

@register(FxpMultHandler, *, [Number, Number] => [Number])

function (handler::FxpMultHandler)(netlist::Netlist,
                                   inputs::Vector{Variable},
                                   outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("FXP multiply does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN f_mult$(handler.id)
        fxp_mult #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) f_mult$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END f_mult$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
import QDLDL
using LightGraphs: SimpleGraph, maximal_cliques
import AMD, Metis, CuthillMcKee

function preprocess!(mats::AbstractVector{SparseMatrixCSC{T,S}}) where {T <: Number, S <: Integer}
    for mat in mats
        preprocess!(mat)
    end
end

"""
	preprocess!(mat::SparseMatrixCSC{T, <:Integer}) where {T}

Checks that `mat` is symmetric and drops numerical zeros.
"""
function preprocess!(mat::SparseMatrixCSC{T, <:Integer}) where {T}
    n = size(mat, 1)
    @assert size(mat) == (n, n)
	for i in 1:n, j in i+1:n
		@assert mat[i,j] .== mat[j,i]
	end
    dropzeros!(mat)
end

# Returns 0/1 sparsity pattern of a set of matrices
"""
	sparsity_pattern(mats::AbstractVector{SparseMatrixCSC{T,S}}) where {T, S <: Integer}

Returns the 0/1 aggregate sparsity pattern of the matrices in `mats` as a sparse matrix.
"""
function sparsity_pattern(mats::AbstractVector{SparseMatrixCSC{T,S}}) where {T, S <: Integer}
	n = size(mats[1], 1)
	for mat in mats
		!(n == size(mat, 1) && issymmetric(mat)) && error(ArgumentError(
			"Matrices must be symmetric and have the same dimensions."
		))
	end

	# Build 0/1 sparsity pattern matrix
	ret = spzeros(n,n)
	for mat in mats
	    @inbounds for col in 1:n, k in nzrange(mat, col)
	        ret[rowvals(mat)[k], col] = 1.0
		end
	end

	# Always assume diagonal can be nonzero
	ret[diagind(ret)] .= 1

    return ret
end
sparsity_pattern(mat::SparseMatrixCSC) = sparsity_pattern([mat])

# Gets chordal extension: a reordering + associated graph from cholesky
"""
	get_chordal_extension(sp_pattern::SparseMatrixCSC; perm="amd", verbose=false)

Returns `p, ip, L`, where `L` is a lower-triangular matrix representing the
chordal extension of the undirected graph `sp_pattern` after fill-reducing permutation
`p` (with inverse permutation `ip`).

Options for perm include
	- `"nothing"` (no permutation)
	- `amd` (approximate minimum degree)
	- `metis` (uses [METIS](http://glaros.dtc.umn.edu/gkhome/metis/metis/overview))
	- `cuthill-mckee` (Cuthill-McKee bandwidth reducing algorithm)

"""
function get_chordal_extension(sp_pattern::SparseMatrixCSC; perm="amd", cost=nothing, verbose=false)
	!issymmetric(sp_pattern) && error(ArgumentError("Matrix must be symmetric"))
	n = size(sp_pattern)[1]

	# Permutations: nothing, amd, vertex separators, TODO: add more options?
	# See survey Minimal triangulations of graphs: A survey
	# https://www.sciencedirect.com/science/article/pii/S0012365X05006060
	if perm == "none" || perm == "nothing" || isnothing(perm)
		p = nothing
	elseif perm == "amd"
		p = AMD.amd(sp_pattern)
	elseif perm == "metis"
		# Convert to Int64 for QDLDL (ret as Int32)
		# first arg = perm (second = iperm)
		p = Int.(Metis.permutation(sp_pattern)[1])
    elseif perm == "cuthill-mckee"
		p = cuthill_mckee(sp_pattern)
	elseif perm == "greedy" && !isnothing(cost)
		# NOTE: This part isn't really working well yet
		# Choose v according to criterion X
		# Let v be the eith vertex in ordering p
		# Let H be the graph obstained by eliminating v from H
		#  - make neighborhood of v a clique and then remove v
		# Return H, p
		# Criteria:
		# Greedy degree: 				v = argminᵤ d(u)
		# Greedy fill in				argminᵤ f(u) = num pairs of non-adj neighbors
		# Greedy degree + fill in		d(u) + f(u)
		# Greedy sparsest subgraph		f(u) - d(u)
		# Greedy fill in degree			d(u) + 1/n^2 * f(u)
		# Greedy degree fill in			f(u) + 1/n * d(u)
		n = size(sp_pattern, 1)
		p = Vector{Int}(undef, n)
		active_inds = Vector(1:n)
		for i in 1:n
			H = @view(sp_pattern[active_inds, active_inds])
			min_cost_vert = get_min_cost_vert(H, cost)
			p[i] = active_inds[min_cost_vert]
			deleteat!(active_inds, [min_cost_vert])
		end

	else
		error(ArgumentError("Invalid permutation"))
	end
	# TODO: add greedy
	# Algorithm: (summarized in Bodlaender & Koster, Treewith computations I. Upper bounds)
	# H = copy(G)
	# for i in 1:n

	F = QDLDL.qdldl(sp_pattern, perm=p, logical=true)

	if verbose
		num_nonzero = 2*nnz(F.L) + n
		num_nonzero_added = (num_nonzero - nnz(sp_pattern)) ÷ 2

		@info "Chordal Extension added $num_nonzero_added nonzeros."
		@info "Density: $(round(num_nonzero/n^2; digits=5))"
	end

    return F.perm, F.iperm, F.L
end

"""
	get_cliques(L::SparseMatrixCSC)

Returns the (maximal) cliques of the undirected graph represented by lower
triangular matrix `L`.
"""
function get_cliques(L::SparseMatrixCSC)
	G = SimpleGraph(L + L')
	return maximal_cliques(G)
end

function cuthill_mckee(sp::SparseMatrixCSC)
	# TODO: Custom implementation for efficiency?
	# Cuthill-McKee (https://en.wikipedia.org/wiki/Cuthill–McKee_algorithm)
	return CuthillMcKee.symrcm(sp)
end

function get_min_cost_vert(H, cost::Function)
	# NOTE: Removed check for performance
	# !issymmetric(H) && error("adj matrix must be symmetric")
	n = size(H, 1)

	min_cost = Inf
	amin_cost = 0
	for i in 1:n
		v_cost = cost(H, i)
		if v_cost < min_cost
			min_cost = v_cost
			amin_cost = i
		end
	end

	return amin_cost
end

# Some cost functions
# H is a (sparse) adj matrix of a graph (CSC format)
# i is the vertex to eval cost on
# https://www.sciencedirect.com/science/article/pii/S0890540109000947
function cost_deg(H::SparseMatrixCSC, i)
	# NOTE: assumes 0/1 matrix
	return length(H.colptr[i]:(H.colptr[i+1]-1))
end

function cost_fill_in(H, i)
	nna_neighbors = 0
	display(fieldnames(H))
	neighbors = @view(H.rowval[H.colptr[i]:(H.colptr[i+1]-1)])
	for i in neighbors, j in neighbors
		i == j && continue
		if H[i,j] == 0
			nna_neighbors += 1
		end
	end
	return nna_neighbors
end

cost_deg_plus_fill_in(H, i) = cost_deg(H, i) + cost_fill_in(H, i)
cost_sparsest_subset(H, i) 	= cost_fill_in(H, i) - cost_deg(H, i)
cost_fill_in_degree(H, i) 	= cost_deg(H, i) + 1/size(H,1)^2 * cost_fill_in(H, i)
cost_deg_fill_in(H, i) 		= 1/size(H,1) * cost_deg(H, i) + cost_fill_in(H, i)


text:
export
    LazyInference,
    LazyState

"""
    mutable struct LazyState

Maintains the state of a lazy algorithm

# Fields

- `previous_algorithm`: The last instant algorithm used, if any
- `evidence`: The evidence supplied in `prepare`
- `interventions`: The interventions supplied in `prepare`
- `placeholder_beliefs`: The placeholder beliefs supplied in `prepare`
- `next_size`: The range size to use in the next call to `refine`
- `next_depth`: The depth to use in the next call to `refine`
- `next_iteration`: The number of the next iteration
- `is_complete`: A flag indicating whether the netwowk has been fully expanded
- `order`: The order of nodes used in computations
"""
mutable struct LazyState
    previous_algorithm :: Union{InstantAlgorithm, Nothing}
    evidence :: Dict{Symbol, Score}
    interventions :: Dict{Symbol, Dist}
    placeholder_beliefs :: Dict{Symbol, Dist}
    next_size :: Int
    next_depth :: Int
    next_iteration :: Int
    is_complete :: Bool
    order :: Vector{Node}
end

"""
    LazyState(ns, nd, ni, ic)

Intantiate `LazyState` with `next_size`, `next_depth`, `next_iterator`, and `is_complete`.
"""
LazyState(ns, nd, ni, ic) = LazyState(nothing, Dict{Symbol,Score}(), Dict{Symbol,Dist}(), Dict{Symbol,Dist}(), ns, nd, ni, ic, Node[])

"""
    struct LazyInference <: IterativeAlgorithm

An iterative algorithm that expands recursively and increases the ranges of instances on every iteration.
"""
struct LazyInference <: IterativeAlgorithm
    algorithm_maker :: Function # A function that takes the current depth and range size and returns the instant algorithm to use
    increment :: Int
    start_size :: Int
    max_iterations :: Int
    start_depth :: Int
    state :: LazyState

    """
    function LazyInference(maker::Function; increment = 10, start_size = increment, max_iterations = 100, start_depth = 1)

    # Arguments
    - `maker``:  A function that takes a range size and expansion depth and returns an `InstantAlgorithm`
    - `increment`:  The increment to range size on every iteration
    - `start_size`:  The starting range size
    - `max_iterations`:  The maximum number of refinement steps
    - `start_depth`:  The depth of recursive expansion in the first iteration
    """
    function LazyInference(maker::Function;
        increment = 10, start_size = increment, max_iterations = 100, start_depth = 1)
        new(maker, increment, start_size, max_iterations, start_depth,
            LazyState(start_size, start_depth, 1, false))
    end
end

function answer(query::Query, lazyalg::LazyInference, runtime::InstantRuntime, target::VariableInstance)
    state = lazyalg.state
    return answer(query, state.previous_algorithm, runtime, target)
end

function prepare(alg::LazyInference, runtime::InstantRuntime,
    evidence::Dict{Symbol, <:Score} = Dict{Symbol, Score}(),
    interventions::Dict{Symbol, <:Dist} = Dict{Symbol, Dist}(),
    placeholder_beliefs = get_placeholder_beliefs(runtime))
    ensure_all!(runtime)
    net = get_network(runtime)
    # The evidence, interventions, and placeholder_beliefs are punted to refine to pass to the underlying algorithm
    state = alg.state
    state.evidence = evidence
    state.interventions = interventions
    state.placeholder_beliefs = placeholder_beliefs
    state.order = topsort(get_initial_graph(net))
    clear_analysis!()
end

function refine(lazyalg::LazyInference, runtime::InstantRuntime)
    state = lazyalg.state
    if !(state.is_complete || state.next_iteration > lazyalg.max_iterations)
        set_ranges!(runtime, state.evidence, state.next_size, state.next_depth, state.order)
        inferencealg = lazyalg.algorithm_maker(state.next_size, state.next_depth)
        state.previous_algorithm = inferencealg
        infer(inferencealg, runtime,
            state.evidence, state.interventions, state.placeholder_beliefs)
        state.is_complete = _complete(runtime)
        state.next_iteration += 1
        state.next_depth += 1
        state.next_size += lazyalg.increment
    end
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

import SparseArrays
import LinearAlgebra

"""
    SensitivityReport

See [`lp_sensitivity_report`](@ref).
"""
struct SensitivityReport
    rhs::Dict{ConstraintRef,Tuple{Float64,Float64}}
    objective::Dict{VariableRef,Tuple{Float64,Float64}}
end

Base.getindex(s::SensitivityReport, c::ConstraintRef) = s.rhs[c]
Base.getindex(s::SensitivityReport, x::VariableRef) = s.objective[x]

"""
    lp_sensitivity_report(model::Model; atol::Float64 = 1e-8)::SensitivityReport

Given a linear program `model` with a current optimal basis, return a
[`SensitivityReport`](@ref) object, which maps:

 - Every variable reference to a tuple `(d_lo, d_hi)::Tuple{Float64,Float64}`,
   explaining how much the objective coefficient of the corresponding variable
   can change by, such that the original basis remains optimal.
 - Every constraint reference to a tuple `(d_lo, d_hi)::Tuple{Float64,Float64}`,
   explaining how much the right-hand side of the corresponding constraint can
   change by, such that the basis remains optimal.

Both tuples are relative, rather than absolute. So given a objective coefficient
of `1.0` and a tuple `(-0.5, 0.5)`, the objective coefficient can range
between `1.0 - 0.5` an `1.0 + 0.5`.

`atol` is the primal/dual optimality tolerance, and should match the tolerance
of the solver used to compute the basis.

Note: interval constraints are NOT supported.

# Example

    model = Model(GLPK.Optimizer)
    @variable(model, -1 <= x <= 2)
    @objective(model, Min, x)
    optimize!(model)
    report = lp_sensitivity_report(model; atol = 1e-7)
    dx_lo, dx_hi = report[x]
    println(
        "The objective coefficient of `x` can decrease by \$dx_lo or " *
        "increase by \$dx_hi."
    )
    c = LowerBoundRef(x)
    dRHS_lo, dRHS_hi = report[c]
    println(
        "The lower bound of `x` can decrease by \$dRHS_lo or increase " *
        "by \$dRHS_hi."
    )
"""
function lp_sensitivity_report(model::Model; atol::Float64 = 1e-8)
    if !_is_lp(model)
        error(
            "Unable to compute LP sensitivity because model is not a linear " *
            "program (or it contains interval constraints).",
        )
    elseif !has_values(model)
        error("Unable to compute LP sensitivity: no primal solution available.")
    elseif !has_duals(model)
        error("Unable to compute LP sensitivity: no dual solution available.")
    end

    std_form = _standard_form_matrix(model)
    basis = _standard_form_basis(model, std_form)
    B = std_form.A[:, basis.basic_cols]
    @assert size(B, 1) == size(B, 2)

    n = length(std_form.columns)
    is_min = objective_sense(model) == MOI.MIN_SENSE

    x = vcat(value.(all_variables(model)), value.(std_form.constraints))
    x_B = @view x[basis.basic_cols]
    l_B = @view std_form.lower[basis.basic_cols]
    u_B = @view std_form.upper[basis.basic_cols]

    B_fact = LinearAlgebra.factorize(B)
    d = Dict{Int,Vector{Float64}}(
        # We call `collect` here because some Julia versions are missing sparse
        # matrix \ sparse vector fallbacks.
        j => B_fact \ collect(std_form.A[:, j]) for
        j = 1:length(basis.basic_cols) if basis.basic_cols[j] == false
    )

    report = SensitivityReport(
        Dict{ConstraintRef,Tuple{Float64,Float64}}(),
        Dict{VariableRef,Tuple{Float64,Float64}}(),
    )

    ###
    ### Compute RHS sensitivity
    ###

    # There is an easy case to consider: a constraint is basic, so we can just
    # take the distance between the value of the constraint and the
    # corresponding bound. Otherwise, we need to compute a search direction as
    # in `_compute_rhs_range`. This is just the negative of the search direction
    # computed above. Moreover, we have to be careful with doubly-bounded
    # variables, because our computed range doesn't take into account the
    # inactive bound.

    for (i, con) in enumerate(std_form.constraints)
        if basis.constraints[i] == MOI.BASIC
            report.rhs[con] = _basic_range(con, constraint_object(con).set)
        else
            report.rhs[con] = _compute_rhs_range(-d[i+n], x_B, l_B, u_B, atol)
        end
    end
    for (i, con) in enumerate(std_form.bounds)
        con_obj = constraint_object(con)
        if basis.bounds[i] == MOI.BASIC
            report.rhs[con] = _basic_range(con, con_obj.set)
        else
            col = std_form.columns[con_obj.func]
            t_lo, t_hi = _compute_rhs_range(-d[col], x_B, l_B, u_B, atol)
            if basis.bounds[i] == MOI.NONBASIC_AT_UPPER
                t_lo = max(t_lo, std_form.lower[col] - x[col])
            elseif basis.bounds[i] == MOI.NONBASIC_AT_LOWER
                t_hi = min(t_hi, std_form.upper[col] - x[col])
            end
            report.rhs[con] = (t_lo, t_hi)
        end
    end

    ###
    ### Compute objective sensitivity
    ###

    π = Dict{Int,Float64}(
        i => reduced_cost(var) for
        (var, i) in std_form.columns if basis.variables[i] != MOI.BASIC
    )
    for (i, c) in enumerate(std_form.constraints)
        if basis.constraints[i] != MOI.BASIC
            π[n+i] = is_min ? dual(c) : -dual(c)
        end
    end

    for (var, i) in std_form.columns
        if basis.variables[i] == MOI.NONBASIC_AT_LOWER && is_min
            @assert π[i] > -atol
            # We are minimizing and variable `i` is nonbasic at lower bound.
            # (δ⁻, δ⁺) = (-πᵢ, ∞) because increasing the objective coefficient
            # will only keep it at the bound.
            report.objective[var] = (-π[i], Inf)
        elseif basis.variables[i] == MOI.NONBASIC_AT_UPPER && !is_min
            @assert π[i] > -atol
            # We are maximizing and variable `i` is nonbasic at upper bound.
            # (δ⁻, δ⁺) = (-πᵢ, ∞) because increasing the objective coefficient
            # will only keep it at the bound.
            report.objective[var] = (-π[i], Inf)
        elseif basis.variables[i] != MOI.BASIC &&
               std_form.lower[i] < std_form.upper[i]
            @assert π[i] < atol
            # The variable is nonbasic with nonfixed bounds. This is the
            # reverse of the above two cases because the variable is at the
            # opposite bound
            report.objective[var] = (-Inf, -π[i])
        elseif basis.variables[i] != MOI.BASIC
            # The variable is nonbasic with fixed bounds. Therefore, (δ⁻, δ⁺) =
            # (-∞, ∞) because the variable can be effectively substituted out.
            # TODO(odow): is this correct?
            report.objective[var] = (-Inf, Inf)
        else
            # The variable `i` is basic. Given an optimal basis B, the reduced
            # costs are:
            #   c_bar = π = c_N - c_Bᵀ(B⁻¹N)
            # To maintain optimality, we want to find a δ such that (if
            # minimizing):
            #     c_N - (c_B + δeᵢ)ᵀ(B⁻¹N) ≥ 0
            #     c_N - c_BᵀB⁻¹N - δ(eᵢ)ᵀ(B⁻¹N) ≥ 0
            #     π_N ≥ δ * (eᵢ)ᵀ(B⁻¹N)
            # To do so, we can loop through every nonbasic variable `j`, and
            # compute
            #     dᵢⱼ = (eᵢ)ᵀB⁻¹aⱼ
            # Then, depending on the sign of dᵢⱼ, we can compute bounds on δ.
            @assert basis.variables[i] == MOI.BASIC
            t_lo, t_hi = -Inf, Inf
            e_i = sum(basis.basic_cols[ii] for ii in 1:i)
            for j in 1:length(basis.basic_cols)
                if basis.basic_cols[j]
                    continue  # Ignore basic components.
                elseif isapprox(
                    std_form.lower[j],
                    std_form.upper[j];
                    atol = atol,
                )
                    continue  # Fixed variables can be ignored.
                elseif abs(d[j][e_i]) <= atol
                    continue  # Direction is ≈0 so any value of δ is okay.
                end
                # There are three confusing sign switch opportunities. We
                # usually want to be:
                # - minimizing (switch if maximizing)
                # - with d[j][e_i] ≥ 0 (switch if ≤ 0)
                # - and nonbasic at the lower bound (switch if upper)
                # If an odd number of these things is true, then the ratio
                # forms an upper bound for δ. Otherwise, it forms a lower bound.
                stat = j <= n ? basis.variables[j] : basis.constraints[j-n]
                if xor(is_min, d[j][e_i] > atol, stat == MOI.NONBASIC_AT_LOWER)
                    t_hi = min(t_hi, π[j] / d[j][e_i])
                else
                    t_lo = max(t_lo, π[j] / d[j][e_i])
                end
            end
            report.objective[var] = (t_lo, t_hi)
        end
    end

    return report
end

_basic_range(con, set::MOI.LessThan) = (value(con) - set.upper, Inf)
_basic_range(con, set::MOI.GreaterThan) = (-Inf, value(con) - set.lower)
_basic_range(::Any, ::Any) = (0.0, 0.0)

"""
    _compute_rhs_range(d_B, x_B, l_B, u_B, atol)

Assume we start with the optimal solution `x_old`, we want to compute a step
size `t` in a direction `d` such that `x_new = x_old + t * d` is still
represented by the same optimal basis. This can be computed a la primal simplex
where we use an artificial entering variable.

    A * x_new = A * (x_old + t * d)
                = A * x_old + t * A * d
                = 0         + t * A * d  # Since A * x_old = 0
    =>  A * d = 0
    => B * d_B + N * d_N = 0
    => d_B = B \\ -(N * d_N)

Note we only have to compute the basic component of the direction vector,
because `d_N` is just zeros with a `1` in the component associated with the
artificial entering variable. Therefore, all that remains is to compute the
associated column of `N`.

If we are increasing the bounds associated with the `i`th decision variable,
then our artificial entering variable is a duplicate of the `i`th variable, and
`N * d_N = A[:, i]`.

If we are increasing the bounds associated with the `i`th affine constraint,
then our artificial entering variable is a duplicate of the slack variable
associated with the `i`th constraint, i.e., a `-1` in the `i`th row and zeros
everywhere else.

In either case:

    d_B = -(B \\ A[:, i])

Now, having computed a direction such that `x_new = x_old + t * d`. By ensuring
that `A * d = 0`, we maintained structural feasibility. Now we need to compute
bounds on `t` such that `x_new` maintains bound feasibility. That is, compute
bounds on t such that:

    l_B[j] <= x_B[j] + t * d_B[j] <= u_B[j].
"""
function _compute_rhs_range(d_B, x_B, l_B, u_B, atol)
    t_lo, t_hi = -Inf, Inf
    for j in 1:length(l_B)
        if d_B[j] > atol
            t_lo = max(t_lo, (l_B[j] - x_B[j]) / d_B[j])
            t_hi = min(t_hi, (u_B[j] - x_B[j]) / d_B[j])
        elseif d_B[j] < -atol
            t_lo = max(t_lo, (u_B[j] - x_B[j]) / d_B[j])
            t_hi = min(t_hi, (l_B[j] - x_B[j]) / d_B[j])
        else
            continue  # d_B[j] ≈ 0.0
        end
    end
    return t_lo, t_hi
end

"""
    _is_lp(model::Model)

Return `true` if `model` is a linear program.
"""
function _is_lp(model::Model)
    for (F, S) in list_of_constraint_types(model)
        # TODO(odow): support Interval constraints.
        if !(S <: Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo})
            return false
        elseif !(F <: Union{VariableRef,GenericAffExpr})
            return false
        end
    end
    return true
end

"""
    _standard_form_matrix(model::Model)

Given a problem:

    r_l <= Ax <= r_u
    c_l <=  x <= c_u

Return the standard form:

           [A -I] [x, y] = 0
    [c_l, r_l] <= [x, y] <= [c_u, r_u]

`columns` maps the variable references to column indices.
"""
function _standard_form_matrix(model::Model)
    columns = Dict(var => i for (i, var) in enumerate(all_variables(model)))
    n = length(columns)
    c_l, c_u = fill(-Inf, n), fill(Inf, n)
    r_l, r_u = Float64[], Float64[]
    I, J, V = Int[], Int[], Float64[]
    bound_constraints = ConstraintRef[]
    affine_constraints = ConstraintRef[]
    for (F, S) in list_of_constraint_types(model)
        _fill_standard_form(
            model,
            columns,
            bound_constraints,
            affine_constraints,
            F,
            S,
            c_l,
            c_u,
            r_l,
            r_u,
            I,
            J,
            V,
        )
    end
    return (
        columns = columns,
        lower = vcat(c_l, r_l),
        upper = vcat(c_u, r_u),
        A = SparseArrays.sparse(I, J, V, length(r_l), n + length(r_l)),
        bounds = bound_constraints,
        constraints = affine_constraints,
    )
end

function _fill_standard_form(
    model::Model,
    x::Dict{VariableRef,Int},
    bound_constraints::Vector{ConstraintRef},
    ::Vector{ConstraintRef},
    F::Type{VariableRef},
    S::Type,
    c_l::Vector{Float64},
    c_u::Vector{Float64},
    ::Vector{Float64},
    ::Vector{Float64},
    ::Vector{Int},
    ::Vector{Int},
    ::Vector{Float64},
)
    for c in all_constraints(model, F, S)
        push!(bound_constraints, c)
        c_obj = constraint_object(c)
        i = x[c_obj.func]
        set = MOI.Interval(c_obj.set)
        c_l[i] = max(c_l[i], set.lower)
        c_u[i] = min(c_u[i], set.upper)
    end
    return
end

function _fill_standard_form(
    model::Model,
    x::Dict{VariableRef,Int},
    ::Vector{ConstraintRef},
    affine_constraints::Vector{ConstraintRef},
    F::Type{<:GenericAffExpr},
    S::Type,
    ::Vector{Float64},
    ::Vector{Float64},
    r_l::Vector{Float64},
    r_u::Vector{Float64},
    I::Vector{Int},
    J::Vector{Int},
    V::Vector{Float64},
)
    for c in all_constraints(model, F, S)
        push!(affine_constraints, c)
        c_obj = constraint_object(c)
        @assert iszero(c_obj.func.constant)
        row = length(r_l) + 1
        set = MOI.Interval(c_obj.set)
        push!(r_l, set.lower)
        push!(r_u, set.upper)
        for (var, coef) in c_obj.func.terms
            push!(I, row)
            push!(J, x[var])
            push!(V, coef)
        end
        push!(I, row)
        push!(J, length(x) + row)
        push!(V, -1.0)
    end
    return
end

_convert_nonbasic_status(::MOI.LessThan) = MOI.NONBASIC_AT_UPPER
_convert_nonbasic_status(::MOI.GreaterThan) = MOI.NONBASIC_AT_LOWER
_convert_nonbasic_status(::Any) = MOI.NONBASIC

function _standard_form_basis(model::Model, std_form)
    variable_status = fill(MOI.BASIC, length(std_form.columns))
    bound_status = fill(MOI.BASIC, length(std_form.bounds))
    constraint_status = fill(MOI.BASIC, length(std_form.constraints))
    for (i, c) in enumerate(std_form.bounds)
        status = MOI.get(model, MOI.ConstraintBasisStatus(), c)
        c_obj = constraint_object(c)
        if status == MOI.NONBASIC
            status = _convert_nonbasic_status(c_obj.set)
        end
        if status != MOI.BASIC
            col = std_form.columns[c_obj.func]
            variable_status[col] = status
        end
        bound_status[i] = status
    end
    for (i, c) in enumerate(std_form.constraints)
        status = MOI.get(model, MOI.ConstraintBasisStatus(), c)
        if status == MOI.NONBASIC
            status = _convert_nonbasic_status(constraint_object(c).set)
        end
        constraint_status[i] = status
    end
    return (
        variables = variable_status,
        bounds = bound_status,
        constraints = constraint_status,
        basic_cols = vcat(variable_status, constraint_status) .==
                     Ref(MOI.BASIC),
    )
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
f = open("hightemp.txt")
println(countlines(f))


text:
module RealPolynomialRoots

using LinearAlgebra

include("bigfloat-utils.jl")
include("interval-arithmetic.jl")
include("descartes-bound.jl")
include("anewdsc.jl")
include("refine-roots.jl")

export ANewDsc, refine_roots

end


text:
@testset "nonmissing" begin
    U = Union{Missing,Int}
    @test nonmissing(U) == Int
end

@testset "table" begin
    T0 = Table(Continuous)
    @test T0 == Table{K} where K<:AbstractVector{<:Continuous}
    T1 = Table(Continuous, Count)
    @test T1 == Table{K} where K<:Union{AbstractVector{<:Continuous}, AbstractVector{<:Count}}
    T2 = Table(Continuous, Union{Missing,Continuous})
    @test T2 == Table{K} where K<:Union{AbstractVector{<:Union{Missing,Continuous}}}
end


text:
#####
# Tensor sets
#####

function test_tensor_sets(T)
    a = Fourier(12)
    b = Fourier(13)
    c = Fourier(11)
    d = TensorProductDict(a,b,c)

    bf = d[3,4,5]
    x1 = T(2//10)
    x2 = T(3//10)
    x3 = T(4//10)
    @test bf(x1, x2, x3) ≈ eval_element(a, 3, x1) * eval_element(b, 4, x2) * eval_element(c, 5, x3)

    # Can you iterate over the product set?
    z = zero(T)
    i = 0
    for f in d
        z += f(x1, x2, x3)
        i = i+1
    end
    @test i == length(d)
    @test abs(-0.5 - z) < 0.01

    z = zero(T)
    l = 0
    for i in eachindex(d)
        f = d[i]
        z += f(x1, x2, x3)
        l = l+1
    end
    @test l == length(d)
    @test abs(-0.5 - z) < 0.01

    # Indexing with ranges
    # @test try
    #     d[CartesianIndices(CartesianIndex(1,1,1),CartesianIndex(3,4,5))]
    #     true
    # catch
    #     false
    # end
    # Is an error thrown if you index with a range that is out of range?
end


text:
# # The exercise notebook
#
# Like the tutorial notebook, but this has solutions. You need to follow a set
# convention to mark solutions and tests for the solutions. It can be nice to
# keep exercises separate

max_stars_repo_path: Courses/Template/0500 The Jupyter notebooks-quiz.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
import Observables
import AbstractPlotting

scene = AbstractPlotting.Scene()
number = rand(5:10)
points = [Geographic(rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
top = U1(rand() * 2pi - pi)
bottom = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
color = AbstractPlotting.RGBAf0(rand(4)...)
transparency = rand(1:2) == 1 ? true : false
whirl = Whirl(scene,
              points,
              top = top,
              bottom = bottom,
              s3rotation = s3rotation,
              config = config,
              segments = segments,
              color = color,
              transparency = transparency)

## update points

points2 = [Geographic(rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
update(whirl, points2)

@test isapprox(whirl.points, points2)

## update top and bottom

top2 = U1(rand() * 2pi - pi)
bottom2 = U1(rand() * 2pi - pi)
update(whirl, top2, bottom2)

@test isapprox(whirl.top, top2)
@test isapprox(whirl.bottom, bottom2)

## update S³ rotation

s3rotation2 = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
update(whirl, s3rotation2)

@test isapprox(whirl.s3rotation, s3rotation2)

## update configuration

config2 = config * Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
value1 = getsurface(whirl.observable, whirl.segments, length(whirl.points))
update(whirl, config2)
value2 = getsurface(whirl.observable, whirl.segments, length(whirl.points))

@test isapprox(whirl.config, config2)
@test isapprox(value1, value2) == false


text:
abstract type DynamicGeneratorComponent <: DynamicComponent end
abstract type AVR <: DynamicGeneratorComponent end
abstract type Machine <: DynamicGeneratorComponent end
abstract type PSS <: DynamicGeneratorComponent end
abstract type Shaft <: DynamicGeneratorComponent end
abstract type TurbineGov <: DynamicGeneratorComponent end