text:
## tests of find_zero interface
using Roots
using Test

# for a user-defined method
struct Order3_Test <: Roots.AbstractSecant end

## Test the interface
@testset "find_zero interface tests" begin

    meths = [Order0(),
             Order1(), Roots.Order1B(), Roots.King(),
             Order2(), Roots.Steffensen(), Roots.Order2B(), Roots.Esser(),
             Order5(), Roots.KumarSinghAkanksha(),
             Order8(), Roots.Thukral8(),
             Order16(), Roots.Thukral16()]

    ## different types of initial values
    for m in meths
        @test find_zero(sin, 3, m) ≈ pi
        @test find_zero(sin, 3.0, m) ≈ pi
        @test find_zero(sin, big(3), m) ≈ pi
        @test find_zero(sin, big(3.0), m) ≈ pi
    end

    ## defaults for method argument
    @test find_zero(sin, 3.0) ≈ pi    # order0()
    @test find_zero(sin, (3,4)) ≈ π   # Bisection()
    @test find_zero(sin, [3,4]) ≈ π   # Bisection()

    ## test tolerance arguments
    ## xatol, xrtol, atol, rtol, maxevals, maxfneval, strict
    fn, xstar = x -> sin(x) - x + 1, 1.9345632107520243
    x0, M = 20.0, Order2()
    @test find_zero(fn, x0, M)  ≈ xstar   # needs 16 iterations, 33 fn evaluations, difference is exact

    # test of maxevals, maxfneval
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, maxevals=2)
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, maxfnevals=2)

    # tolerance on f, atol, rtol: f(x) ~ 0
    M = Order2()
    h = 1e-2
    rt = find_zero(fn, x0, M, atol=h, rtol=0.0)
    @test abs(fn(rt)) > h^2 / 100
    rt = find_zero(fn, x0, M, atol=0.0, rtol=h)
    @test abs(fn(rt)) > h^2 / 100

    ## test of tolerances xatol, xrtol with bisection
    a,b = 1.5, 2.0
    h = 1e-6
    M = Roots.Bisection()
    tracks = Roots.Tracks(Float64[], Float64[])
    find_zero(fn, (a,b), M, tracks=tracks, xatol=h, xrtol=0.0)
    u,v = tracks.xs[end-1:end]
    @test h >= abs(u-v) >= h/2

    ## test of strict
    fn, x0 = x -> cos(x) - 1, pi/4
    @test fn(find_zero(fn, x0, Order5())) <= 1e-8
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, Order5(), strict=true)

    # xn increment needs atol setting for zeros near 0.0 if strict=true
    M = Order1()
    @test_throws Roots.ConvergenceFailed find_zero(x -> x * exp(x), 1.0, M, atol=0.0, rtol=0.0, strict=true, xatol=0.0)
    @test abs(find_zero(x -> x * exp(x), 1.0, M, atol=0.0, rtol=0.0, strict=true)) <= eps()

    ## test of extreme values for fn, bisection
    c = pi
    fn = x -> Inf * sign(x - c)
    find_zero(fn, (-Inf, Inf)) ≈ c

    fn = x -> Inf * x/abs(x) # stop at NaN values
    find_zero(fn, (-Inf, Inf)) ≈ 0
end

@testset "find_zero internals" begin
    ##  init_state! method
    g1(x) = x^5 - x - 1
    x0_, xstar_ = (1.0, 2.0), 1.1673039782614187
    M = Roots.A42()
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state)
    for M in (Roots.A42(), Roots.Bisection(), Roots.FalsePosition())
        Roots.init_state!(state, M, g1, x0_)
        @test find_zero(M, g1, options, state) ≈ xstar_
    end

    ## hybrid
    g1(x) = exp(x) - x^4
    x0_, xstar_ = (5.0, 20.0), 8.613169456441398
    M = Roots.Bisection()
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state, xatol=1/2)
    find_zero(M, g1, options, state)
    M = Roots.Order1()
    Roots.init_state!(state, M, g1, state.m[1])
    Roots.init_options!(options, M)
    @test find_zero(M, g1, options, state) ≈ xstar_

    # test conversion between states/options
    Ms = vcat([Roots.FalsePosition(i) for i in 1:12], Roots.A42(), Roots.AlefeldPotraShi(), Roots.Brent(), Roots.Bisection())
    Ns = [Order1(),
          Roots.Order1B(),
          Order2(),
          Roots.Order2B(),
          Order5(),
          Order8(),
          Order16()]

    g1(x) = exp(x) - x^4
    x0_= (5.0, 20.0)
    M = Ms[1]
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state)

    for M in Ms
        @test Roots.init_state!(state, M, g1, x0_) == nothing
        @test Roots.init_options!(options, M) == nothing
        @test Roots.update_state(M, g1, state, options) == nothing
        @test Roots.assess_convergence(M, state, options) == false # none in just 1 step
    end

    x0_ = 8.0
    xstars = find_zeros(g1, -20.0, 20.0)
    for M in Ns
        @test Roots.init_state!(state, M, g1, x0_) == nothing
        @test Roots.init_options!(options, M) == nothing
        @test Roots.update_state(M, g1, state, options)  == nothing
        @test Roots.assess_convergence(M, state, options) == false  # none in 1 step
    end

    ## test creation of new methods
    ## xn - f/f' - f'' * f^2 / 2(f')^3 = xn - r1 - r1^2/r2 is third order,
    # had to previousely define:
    # struct Order3_Test <: Roots.AbstractSecant end
    function Roots.update_state(M::Order3_Test, f, o::Roots.UnivariateZeroState{T,S}, options) where {T, S}
        # xn - f/f' - f'' * f^2 / 2(f')^3 = xn - r1 - r1^2/r2 is third order
        xn_1, xn = o.xn0, o.xn1
        fxn_1, fxn = o.fxn0, o.fxn1

        f_10 = (fxn - fxn_1) / (xn - xn_1)
        xn1::T = xn - fxn / f_10
        fxn1::S = f(xn1)

        f01 = (fxn1 - fxn) /  (xn1 - xn)

        if isnan(f_10) || iszero(f_10) || isnan(f01) || iszero(f01)
            o.message = "issue with bracket. "
            o.stopped = true
            return
        end

        r1 = fxn1 / f01
        r2 = f01 / ((f01 - f_10) / (xn1 - xn_1))

        wn = xn1 - r1 - r1^2/r2
        fwn::S = f(wn)

        o.xn0, o.xn1 = xn, wn
        o.fxn0, o.fxn1 = fxn, fwn
    end

    g1(x) = exp(x) - x^4
    @test find_zero(g1, 8.3, Order3_Test()) ≈ find_zero(g1, 8.3, Order1())

end

@testset "find_zero issue tests" begin

    ## Misc tests
    Ms = [Order0(), Order1(), Roots.Order1B(), Order2(), Roots.Order2B(),
          Order5(), Order8(), Order16()]

    ## issues with starting near a maxima. Some bounce out of it, but
    ## one would expect all to have issues
    fn, xstar = x -> x^3 + 4x^2 -10,  1.365230013414097
    for M in [Order1(), Roots.Order1B(), Order2(), Roots.Order2B()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, -1.0, M)
    end
    for M in [Order0(),  Order5(), Order8(), Order16()]
        @test find_zero(fn, -1.0, M) ≈ xstar
    end

    ## non-simple root
    ## convergence can depend on relaxed convergence checked after an issue
    fn, xstar, x0 = x -> cos(x) - 1, 0.0, 0.1
    for M in Ms
        xn = find_zero(fn, x0, M)
        @test abs(fn(xn)) <= 1e-10
    end
    for M in [Order2(), Order5(), Order8(), Order16()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M, strict=true)
    end

    ## issue with large steps
    fn, x0 = x -> x^20 - 1, 0.5
    for M in Ms[2:end] # not 0, as it uses bracket
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M)
    end

    ## issue with large f''
    fn, x0 = cbrt, 1.0
    for M in [Order1(), Order2(), Order5()]
        @test_throws Roots.ConvergenceFailed find_zero(fn, x0, M)
    end
    ### these stop but only because rtol is used for checking f(xn) ~ 0
    for M in [Order8(), Order16()]
        @test abs(find_zero(fn, x0, M) -  0.0) >= 100
    end

    ## similar (http://people.sc.fsu.edu/~jburkardt/cpp_src/test_zero/test_zero.html)
    function newton_baffler(x)
        a = 1/10
        m, b = 1/4, 1/8

        if x < -a
            m*x - b
        elseif x > a
            m*x + b
        else
            (m*a+b)/a * (x+a) + (-m*a -b)
        end
    end
    for M in (Order0(), Order1(), Roots.Order1B(), Order2(), Roots.Order2B())
        @test abs(find_zero(newton_baffler, 1.0, M)) <= 1e-15
    end
    for M in (Order5(), Order8())
        @test_throws Roots.ConvergenceFailed find_zero(newton_baffler, 1.0, M)
    end

    ## Closed issues ###
    ## issue tests: put in tests to ensure closed issues don't reappear.

    ## issue #94; tolerances not matching documentation
    function test_94(;kwargs...)
        g, T = 1.62850, 14.60000
        α, t1, tf = 0.00347, 40.91375, 131.86573
        y, ya, yf = 0.0, 9000.0, 10000.0
        vy = sqrt(2g*(ya-y))
        θ0, θ1 = atan(α*tf), atan(α*(tf-t1))
        I_sintan(x) = tan(x)/2cos(x) - atanh(tan(x/2))
        I_sintan(x, y) = I_sintan(y) - I_sintan(x)
        function lhs(θ)
        tRem = (vy - T/α*(sec(θ1) - sec(θ))) / g
            val = -yf + y + vy*tRem - 0.5g*tRem^2 - T/α^2*I_sintan(θ, θ1)
            val
        end

        M = Roots.FalsePosition()
        f, x0 = lhs, [atan(α*tf), atan(α*(tf-t1))]
        F = Roots.DerivativeFree(lhs)
        state = Roots.init_state(M, F, x0)
        options = Roots.init_options(M, state)
        find_zero(M, F, options, state)

        @test state.steps <= 45 # 15
    end
    test_94()

    ## Issue with quad_step after truncated M-step PR #140
    @test find_zero(x -> tanh(x) - tan(x), 7.36842 , Order0()) ≈ 7.068582745628732

    ## Use tolerance on f, not x with bisectoin
    atol = 0.01
    u = find_zero(sin, (3, 4), atol=atol)
    @test atol >= abs(sin(u)) >= atol^2

    ## issue #159 bracket with zeros should be found
    @test find_zero(x->x+1,(-1,1)) == -1

end

struct _SampleCallableObject end
_SampleCallableObject(x) = x^5 - x - 1

mutable struct Cnt
    cnt::Int
    f
    Cnt(f) = new(0, f)
end
(f::Cnt)(x) = (f.cnt += 1; f.f(x))

@testset "find_zero with other callable types" begin

    Ms = [Order0(),
          Order1(), Roots.Order1B(),
          Order2(), Roots.Order2B(),
          Order5(),
          Order8(),
          Order16()]

    for M in Ms
        @test find_zero(_SampleCallableObject, 1.0, M) ≈ 1.1673039782614187
    end

    for M in Ms
        g = Cnt(x -> x^5 - x - 1)
        @test find_zero(g, 1.0, M) ≈ 1.1673039782614187
        @test g.cnt <= 30
    end

end

@testset "find_bracket test" begin
    fs_xs = ((x -> x^5 - x - 1, (0,2)),
             (sin, (3,4)),
             (x -> exp(x) - x^4, (5, 20))
             )

    for (f, x0) in fs_xs
        out = Roots.find_bracket(f, x0)
        @test prod(f.(out.bracket)) <= 0
    end

    ## test size of  bracket
    for M in  (Roots.BisectionExact(),  Roots.A42(), Roots.AlefeldPotraShi())
        for (fn, b) in ((x  -> x == 0. ? 0. : x/exp(1/(x*x)), (-1.0, 4.0)),
                        (x -> exp(-15*x)*(x-1) + x^15, (0., 1.)),
                        (x -> (-40.0)*x*exp(-x),  (-9, 31))
                        )
            l =  Roots.find_bracket(fn, b, M)
            @test  l.exact ||  abs(l.bracket[2] -  l.bracket[1])  <=  eps(l.xstar)
        end
    end

    ## subnormal
    x0 = nextfloat(nextfloat(0.0))
    fn = x  ->  (x-x0)
    b =  (0.0, 1.0)
    m =  Roots.__middle(b...)  # 1e-154
    M = Roots.BisectionExact()
    l = Roots.find_bracket(fn,  b,  M)
    ## Here l.exact=true, but  this checks the bracket size
    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
    @test  abs(-(l.bracket...)) <=  m

    M =  Roots.A42()
    l = Roots.find_bracket(fn,  b,  M)
    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
    @test  abs(-(l.bracket...)) <= m

    M =  Roots.AlefeldPotraShi()
    l = Roots.find_bracket(fn,  b,  M)
    @test  !(abs(-(l.bracket...)) <=  eps(l.xstar))
    @test  abs(-(l.bracket...)) <= m

    x0 = nextfloat(0.0)
    fn  = x ->  x <=  0.0 ?  x-x0  : x  +  x0
    b = (-1.0, 1.0)
    M = Roots.BisectionExact()
    l = Roots.find_bracket(fn,  b,  M)
    @test !l.exact &&   abs(-(l.bracket...))  <= maximum(eps.(l.bracket))

end


text:
module ModelerToolbox

using JLD, HDF5, MacroTools

# NOTE: exports in each file

include("jld.jl")

end # module

ext: jl


text:
using SubstitutionTilings
using SubstitutionTilings.CoreDefs
using StructEquality
using LinearAlgebra

@def_structequal struct TM <: DGroupElem
    a::Int
end

@enum TMTile begin
    A=1
    B=2
end

function Base.:*(x :: TM, y :: TM)
    return TM(x.a + y.a)
end

function SubstitutionTilings.dilate(λ, x :: TM)
    return TM(λ*x.a)
end

function SubstitutionTilings.inv(x :: TM)
    return TM(-x.a)
end

function SubstitutionTilings.id(::Type{TM})
    return TM(0)
end

function SubstitutionTilings.is_interior(tiling :: Dict, t :: TM)
    return haskey(t.a, tiling) && haskey(t.a-1, tiling) && haskey(t.a+1, tiling)
end

function SubstitutionTilings.collar_in(tiling :: Dict, t :: TM)
    collar_shape = t*[TM(x) for x in [-2,0,2]]
    collar = []
    tiling_dict = tiling
    for s in collar_shape
        if !haskey(tiling_dict, s)
            throw(UnrecognizedCollar)
        end
        push!(collar, (s => tiling_dict[s]))
    end
    return Dict(collar)
end

thue_morse = SubSystem(Dict(A => [(TM(-1) => A), (TM(1) => B)], B => [(TM(-1) => B), (TM(1) => A)]),2)

substitute(thue_morse, [(TM(0) => A)], 4)

(collars, S) = total_collaring(thue_morse, Dict([(TM(-2) => A), (TM(0) => A), (TM(2) => B)]))
length(collars)


text:
import GreyDecision.Utility: makeminmax
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Aras: aras

@testset "Aras" begin
    tol = 0.0001

    K1 = [105000.0, 120000, 150000, 115000, 135000]
    K2 = [105.0, 110, 120, 105, 115]
    K3 = [10.0, 15, 12, 20, 15]
    K4 = [4.0, 4, 3, 4, 5]
    K5 = [300.0, 500, 550, 600, 400]
    K6 = [10.0, 8, 12, 9, 9]

    decisionMat = hcat(K1, K2, K3, K4, K5, K6)
    functionlist = [minimum, maximum, minimum, maximum, maximum, minimum]
    w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40]
    row, col = size(decisionMat)
    greydec = Array{GreyNumber{Float64}, 2}(undef, row, col)
    for i in 1:row
        for j in 1:col
            greydec[i, j] = GreyNumber(Float64(decisionMat[i, j]))
        end
    end

    result = aras(greydec, w, functionlist)

    @test result[:bestIndex] == 2
    @test result[:orderings] == [3, 1, 4, 5, 2]
    @test isapprox(result[:scores][1].a, 0.81424068, atol = tol)
    @test isapprox(result[:scores][2].a, 0.89288620, atol = tol)
    @test isapprox(result[:scores][3].a, 0.76415790, atol = tol)
    @test isapprox(result[:scores][4].a, 0.84225462, atol = tol)
    @test isapprox(result[:scores][5].a, 0.86540635, atol = tol)

end


text:
struct SSHHost
    address::String
    port::Int
end

SSHHost(address) = SSHHost(address, 22)
Base.isless(a::SSHHost, b::SSHHost) = isless(a.address, b.address) && isless(a.port, b.port)

function Base.show(io::IO, h::SSHHost)
    print(io, "SSHHost(\"$(h.address):$(h.port)\")")
end

struct SSHSession
    handle::Ptr{Cvoid}
    socket::TCPSocket
end

Base.convert(::Type{Ptr{LIBSSH2_SESSION}}, sess::SSHSession) = sess.handle
Base.unsafe_convert(::Type{Ptr{LIBSSH2_SESSION}}, sess::SSHSession) = convert(Ptr{LIBSSH2_SESSION}, sess.handle)

function SSHSession(host::SSHHost)
    sock = connect(host.address, host.port)

    sess = libssh2_session_init()

    libssh2_session_handshake(sess, sock)

    return SSHSession(sess, sock)
end

isauthenticated() = isauthenticated(current_ssh_session())
isauthenticated(session::SSHSession) = libssh2_userauth_authenticated(session) != 0

authenticate!(username, password) = authenticate!(current_ssh_session(), username, password)
authenticate!(session::SSHSession, username, password) = libssh2_userauth_password(session, username, password) == 0

struct SFTPSession
    handle::Ptr{Cvoid}
    ssh_session::SSHSession
end

Base.convert(::Type{Ptr{LIBSSH2_SFTP}}, sess::SFTPSession) = sess.handle
Base.unsafe_convert(::Type{Ptr{LIBSSH2_SFTP}}, sess::SFTPSession) = convert(Ptr{LIBSSH2_SFTP}, sess.handle)

function SFTPSession(session::SSHSession=current_ssh_session())
    sftp_sess = libssh2_sftp_init(session)
    return SFTPSession(sftp_sess, session)
end

##
const _current_ssh_session = Ref{Union{SSHSession,Nothing}}(nothing)
const _current_sftp_session = Ref{Union{SFTPSession,Nothing}}(nothing)

function current_ssh_session()
	isnothing(_current_ssh_session[]) && throw("No SSH Session")
	return _current_ssh_session[]::SSHSession
end

function current_sftp_session()
	isnothing(_current_sftp_session[]) && set_current_session(SFTPSession(current_ssh_session()))
	return _current_sftp_session[]::SFTPSession
end

set_current_session(session::SSHSession) = _current_ssh_session[] = session
set_current_session(session::SFTPSession) = _current_sftp_session[] = session

libssh2_userauth_list(uname::String) = libssh2_userauth_list(current_ssh_session(), uname)


text:
its1 = Iterators.filter(x->!(occursin("started", x) || occursin("failed", x)), eachline("vanilla-raw.txt"))
its2 = Iterators.filter(x->!(occursin("started", x) || occursin("failed", x)), eachline("rr-raw.txt"))

function timings(its)
    map(its) do e
        parts = split(e, '|')
        strip(split(parts[1],'(')[1]), parse(Float64, parts[2])
    end
end

vanilla = timings(its1)
rr = timings(its2)

using DataFrames, StatsBase, Statistics, Gadfly
df1 = rename(DataFrame(vanilla), 1=>:name, 2=>:time_vanilla)
df2 = rename(DataFrame(rr), 1=>:name, 2=>:time_rr)
df = join(df1, df2, on=:name)
df[:overhead] = map(df[:time_vanilla], df[:time_rr]) do v, r
    r == 0. && return 0.
    100 * (r/v - 1)
end
df = sort(df, cols=:overhead, rev=true)
part = Iterators.partition(1:size(df, 1), 60)

overhead_mean = mean(df[:overhead])
function plot_for_range(range)
    Gadfly.plot(df[range,:],
        yintercept=[0., overhead_mean], Geom.hline(color=["black", "black"]),
        Guide.xlabel(last(range) == size(df, 1) ? "Test name" : nothing), Guide.ylabel("Overhead (%)"),
        x=:name, y=:overhead, Geom.bar, Scale.y_continuous(minvalue=-100., maxvalue=700.),
        Theme(default_color="orange"))
end

vstack((plot_for_range(range) for range in part)...)
draw(SVG("overhead.svg", 20cm ,40cm), vstack((plot_for_range(range) for range in part)...))


text:
#=
        Constants for the problem, just change nb_teams, the others depend on it
=#
const nb_teams = 10
const nb_weeks = nb_teams - 1
const nb_periods = div(nb_teams, 2)

#= Numbers of variables and number of clauses =#
nb_var = (nb_teams * nb_weeks * nb_periods)
nb_combi = (nb_teams * nb_teams * nb_weeks * nb_periods)
nb_clauses = ((nb_teams * nb_weeks) + (nb_teams * nb_weeks * sum(1:nb_periods-1)))
nb_clauses = nb_clauses + (nb_periods*nb_weeks*nb_teams) + (nb_periods*nb_weeks* binomial(nb_teams, 3))
nb_clauses = nb_clauses + (nb_teams * nb_periods) + (nb_teams * nb_periods * binomial(nb_weeks, 3))
nb_clauses = nb_clauses + binomial(nb_teams, 2) + (binomial(nb_teams, 2) * nb_weeks * nb_periods * 3)
nb_clauses = nb_clauses + nb_teams

# For better readability, we will use a tables constaining the variables
variables = reshape(1:nb_var, nb_teams, nb_weeks, nb_periods)
matchs = reshape(nb_var+1:(nb_var+nb_combi), nb_teams, nb_teams, nb_weeks, nb_periods)

"""
    function creator_sat()

    Will create a file sts.cnf which contains all the clauses of the problem
"""
function create()
    file = open("sts.cnf", "w")         # Open file in writting mode
    write(file, "c SAT version for STS problem\n")  # Comments
    write(file, "c Authors : Killian Fretaud, Rémi Garcia, Benoit Loger, Boualem Lamraoui, Benoît Le Badezet\n")  #Comments
    write(file, "p cnf ", string(nb_var + nb_combi -1), " ", string(nb_clauses), "\n") # Spécifications du problème

    #Constraint 1 : Each teams plays only once per weeks
    for t in 1:nb_teams
        for w in 1:nb_weeks
            for p in 1:nb_periods
                write(file, string(variables[t,w,p]), " ")
            end
            write(file, " 0\n")
            for p in 1:nb_periods-1
                for p2 in p+1:nb_periods
                    write(file, "-", string(variables[t,w,p]), " -", string(variables[t,w,p2]), " 0\n")
                end
            end
        end
    end

    #Constraint 2 : Only two teams by periods
    for w in 1:nb_weeks
        for p in 1:nb_periods
            for t in 1:nb_teams
                for t2 in 1:nb_teams
                    if t != t2
                        write(file, string(variables[t2,w,p], " "))
                    end
                end
                write(file, "0\n")
            end
            for t in 1:nb_teams-2
                for t2 in t+1:nb_teams-1
                    for t3 in t2+1:nb_teams
                        write(file, "-", string(variables[t,w,p]), " -", string(variables[t2,w,p]), " -", string(variables[t3,w,p]), " 0\n")
                    end
                end
            end
        end
    end

    #constraint 3 : Each teams plays at most two times in the same period
    for t in 1:nb_teams
        for p in 1:nb_periods
            for w in 1:nb_weeks
                write(file, string(variables[t,w,p]), " ")
            end
            write(file, "0\n")
            for w in 1:nb_weeks-2
                for w2 in w+1:nb_weeks-1
                    for w3 in w2+1:nb_weeks
                        write(file, "-", string(variables[t,w,p]), " -", string(variables[t,w2,p]), " -", string(variables[t,w3,p]), " 0\n")
                    end
                end
            end
        end
    end

    #Constraint 4: Each team plays angainst all other teams
    for t in 1:nb_teams-1
        for t2 in t+1:nb_teams
            s = ""
            for w in 1:nb_weeks
                for p in 1:nb_periods
                    write(file, "-", string(variables[t,w,p]), " -", string(variables[t2,w,p]), " ", string(matchs[t,t2,w,p]), " 0\n")
                    write(file, string(variables[t,w,p]), " -", string(matchs[t,t2,w,p]), " 0\n")
                    write(file, string(variables[t2,w,p]), " -", string(matchs[t,t2,w,p]), " 0\n")
                    s = s*string(matchs[t,t2,w,p])*" "
                end
            end
            write(file, s, "0\n")
        end
    end
    #Symetry breaker, we force the first week of matchs : 1v2 , 3v4 , ... , nb_teams-1 v nb_teams
    for p in 1:nb_periods
        write(file, string(variables[2*p-1,1,p]), " 0\n")
        write(file, string(variables[2*p,1,p]), " 0\n")
    end

    close(file)
end

create()

max_stars_repo_name: remi-garcia/ProjetCSP

max_issues_repo_name: remi-garcia/ProjetCSP

max_forks_repo_name: remi-garcia/ProjetCSP


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges
const MOIBC = MathOptInterface.Bridges.Constraint

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig(duals = false)

@testset "RSOCtoNonConvexQuad" begin

    @test MOIBC.RSOCtoNonConvexQuadBridge{Float64} == MOIBC.concrete_bridge_type(
        MOIBC.RSOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone)
    @test MOI.supports_constraint(MOIBC.RSOCtoNonConvexQuadBridge{Float64},
                                  MOI.VectorOfVariables,
                                  MOI.RotatedSecondOrderCone)
    @test !MOI.supports_constraint(MOIBC.RSOCtoNonConvexQuadBridge{Float64},
                                   MOI.ScalarAffineFunction{Float64},
                                   MOI.RotatedSecondOrderCone)

    bridged_mock = MOIB.Constraint.RSOCtoNonConvexQuad{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock, config,
        include = [(F, S)
                   for F in [MOI.VectorOfVariables]
                   for S in [MOI.RotatedSecondOrderCone]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.5, 1.0, 1/√2, 1/√2])
    MOIT.rotatedsoc1vtest(bridged_mock, config)

    ci = first(MOI.get(bridged_mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                    MOI.RotatedSecondOrderCone}()))

    test_delete_bridge(bridged_mock, ci, 4,
                    (
                        (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, 0),
                        (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
                    ))
end

@testset "SOCtoNonConvexQuad" begin

    @test MOIBC.SOCtoNonConvexQuadBridge{Float64} == MOIBC.concrete_bridge_type(
        MOIBC.SOCtoNonConvexQuadBridge{Float64},
        MOI.VectorOfVariables,
        MOI.SecondOrderCone)
    @test MOI.supports_constraint(MOIBC.SOCtoNonConvexQuadBridge{Float64},
                                  MOI.VectorOfVariables,
                                  MOI.SecondOrderCone)
    @test !MOI.supports_constraint(MOIBC.SOCtoNonConvexQuadBridge{Float64},
                                   MOI.ScalarAffineFunction{Float64},
                                   MOI.SecondOrderCone)

    bridged_mock = MOIB.Constraint.SOCtoNonConvexQuad{Float64}(mock)

    MOIT.basic_constraint_tests(
        bridged_mock, config,
        include = [(F, S)
                   for F in [MOI.VectorOfVariables]
                   for S in [MOI.SecondOrderCone]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2])
    MOIT.soc1vtest(bridged_mock, config)

    ci = first(MOI.get(bridged_mock,
        MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                    MOI.SecondOrderCone}()))

    test_delete_bridge(bridged_mock, ci, 3,
                       (
                           (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}, 0),
                           (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
                       ))
end


text:
"""
Test minimize function
"""

#push!(LOAD_PATH,"../src/")
using joptimise

function rosenbrock(x)
    f = (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
    return f
end

function rosenbrock!(g, x)
    # compute objective
    f = rosenbrock(x)
    # constraint
    g[1] = x[1]^2 + x[2]^2 - 1.0
    return f
end

# initial guess
x0 = [4.0; 4.0]
# bounds on variables
lx = [-5.0; -5.0]
ux = [5.0; 5.0]
# bounds on constriants
lg = [0.0]
ug = [0.0]
# number of constraints
ng = 1

## run minimizer with IPOPT
ip_options = Dict(
    "max_iter" => 2500,   # 1500 ~ 2500
    "tol" => 1e-6
)

xopt, fopt, info = minimize(rosenbrock!, x0, ng; lx=lx, ux=ux, lg=lg, ug=ug, solver="ipopt", options=ip_options, derivatives=ForwardAD());

println("Done with IPOPT!")
println(info)
println(xopt)

## run minimizer with SNOPT
sn_options = Dict(
    "Major feasibility tolerance" => 1.e-6,
    "Major optimality tolerance"  => 1.e-6,
    "Minor feasibility tolerance" => 1.e-6,
    "Major iterations limit" => 1000,
    "Major print level" => 1,
)

xopt, fopt, info = minimize(rosenbrock!, x0, ng; lx=lx, ux=ux, lg=lg, ug=ug, solver="snopt", options=sn_options, derivatives=ForwardAD());

println("Done with SNOPT!")
println(info)
println(xopt)

max_stars_repo_name: Yuricst/joptimise

max_issues_repo_name: Yuricst/joptimise

max_forks_repo_name: Yuricst/joptimise


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export exact_likelihood

function exact_likelihood(Net::NetworkConditionalHINT, X, Y)

    Zx, Zy, logdet = Net.forward(X, Y)

    loglike = sum(logpdf(0f0, 1f0, cat(Zx, Zy; dims=3)), dims=[1, 2, 3])[1, 1, 1, :]
    loglike = loglike .+ logdet

    return loglike

end

function exact_likelihood(Net::NetworkConditionalHINT, Z)

    Zx, Zy = Z[1], Z[2]

    logdet = Net.inverse(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)); logdet=true, x_lane=true)[3]

    loglike = sum(logpdf(0f0, 1f0, Zx), dims=[1, 2, 3])[1, 1, 1, :]
    loglike = loglike .- logdet

    return loglike

end

function exact_likelihood(Net::ReverseNetwork, Z)

    Zx, Zy = Z[1], Z[2]

    logdet = Net.forward(Zx, repeat(Zy, 1, 1, 1, size(Zx, 4)); logdet=true, x_lane=true)[3]

    loglike = sum(logpdf(0f0, 1f0, Zx), dims=[1, 2, 3])[1, 1, 1, :]
    loglike = loglike .- logdet

    return loglike

end


text:
"""
BVP_solve(n = 801, T = Float64; bfact=qr!)
Solve the BVP for the Chapter 2 figures and testing.
"""
function BVP_solve(n = 801, T = Float64; bfact = qr!)
    # Set it up
    bdata = bvpinit(n, T)
    #
    U0 = zeros(2n)
    FV = zeros(2n)
    # Banded matrix with the correct number of bands
    FPV = BandedMatrix{T}(Zeros(2n, 2n), (2, 4))
    #
    # Build the initial iterate
    #
    BVP_U0!(U0, n, bdata)
    #
    if bfact == qr!
        bvpout = nsol(
            Fbvp!,
            U0,
            FV,
            FPV,
            Jbvp!;
            rtol = 1.e-10,
            sham = 1,
            pdata = bdata,
            jfact = bfact,
        )
    else
        # Test for default of qr. Used for CI only.
        bvpout = nsol(Fbvp!, U0, FV, FPV, Jbvp!; rtol = 1.e-10, sham = 1, pdata = bdata)
    end
    return (bvpout = bvpout, tv = bdata.tv)
end

function BVP_U0!(U0, n, bdata)
    tv = bdata.tv
    view(U0, 1:2:2n-1) .= exp.(-.1 .* tv .* tv)
    view(U0, 2:2:2n) .= -.2 .* view(U0, 1:2:2n-1) .* tv
end


text:
import Base: rand, randn
function rand(::Type{SimpleQuaternion{T}}) where T <: Real
    args = (rand(T) for j=1:4)
    return SimpleQuaternion(args...)
end

"""
`rand(SimpleQuaternion)` returns a quaternion whose four entries
are uniform [0,1] random values.
"""
rand(::Type{SimpleQuaternion}) = rand(SimpleQuaternion{Float64})

"""
`randn(SimpleQuaternion)` returns a quaternion whose four entries are
standard normal random values.
"""
function randn(::Type{SimpleQuaternion{T}}) where T
    args = (randn(T) for j=1:4)
    return SimpleQuaternion(args...)
end

randn(::Type{SimpleQuaternion}) = randn(SimpleQuaternion{Float64})


text:
"""
Convert ASCII (RFC20) character `c` to lower case.
"""
ascii_lc(c::UInt8) = c in UInt8('A'):UInt8('Z') ? c + 0x20 : c

"""
Case insensitive ASCII character comparison.
"""
ascii_lc_isequal(a::UInt8, b::UInt8) = ascii_lc(a) == ascii_lc(b)

"""
Case insensitive ASCII string comparison.
"""
function ascii_lc_isequal(a, b)
    a = Iterators.Stateful(codeunits(a))
    b = Iterators.Stateful(codeunits(b))
    for (i, j) in zip(a, b)
        if !ascii_lc_isequal(i, j)
            return false
        end
    end
    return isempty(a) && isempty(b)
end


text:
# Quantify the error of the 6th order approximation to transverse Mercator
# projections by the invertability of points

using Geodesy
using PyPlot

z = 31
h = true
t1 = UTMfromLLA(z,h,wgs84)
t2 = LLAfromUTM(z,h,wgs84)

lats = Float64[0:1:80;]
lons = Float64[0:1:80;]
errors = Float64[distance(LLA(lats[i],lons[j]+3.0,0.0), transform(t2, transform(t1, LLA(lats[i],lons[j]+3.0,0.0)))) for i=1:length(lats), j in 1:length(lons)]
errors = convert(Matrix{Float64}, errors) # Compehension is broken (code_warntype is correct for distance but not for the array...)

figure()
imshow(log10(errors))

error_lon = [maximum(errors[:,i]) for i in 1:length(lons)]

figure()
plot(lons, error_lon)


text:
# host-level indexing

export allowscalar, @allowscalar, @disallowscalar, assertscalar

# mechanism to disallow scalar operations

@enum ScalarIndexing ScalarAllowed ScalarWarned ScalarDisallowed

const scalar_allowed = Ref(ScalarWarned)
const scalar_warned = Ref(false)

"""
    allowscalar(allow=true, warn=true)

Configure whether scalar indexing is allowed depending on the value of `allow`.

If allowed, `warn` can be set to throw a single warning instead. Calling this function will
reset the state of the warning, and throw a new warning on subsequent scalar iteration.
"""
function allowscalar(allow::Bool=true, warn::Bool=true)
    scalar_warned[] = false
    scalar_allowed[] = if allow && !warn
        ScalarAllowed
    elseif allow
        ScalarWarned
    else
        ScalarDisallowed
    end
    return
end

"""
    assertscalar(op::String)

Assert that a certain operation `op` performs scalar indexing. If this is not allowed, an
error will be thrown ([`allowscalar`](@ref)).
"""
function assertscalar(op = "operation")
    if scalar_allowed[] == ScalarDisallowed
        error("$op is disallowed")
    elseif scalar_allowed[] == ScalarWarned && !scalar_warned[]
        @warn "Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`"
        scalar_warned[] = true
    end
    return
end

"""
    @allowscalar ex...
    @disallowscalar ex...
    allowscalar(::Function, ...)

Temporarily allow or disallow scalar iteration.

Note that this functionality is intended for functionality that is known and allowed to use
scalar iteration (or not), i.e., there is no option to throw a warning. Only use this on
fine-grained expressions.
"""
macro allowscalar(ex)
    quote
        local prev = scalar_allowed[]
        scalar_allowed[] = ScalarAllowed
        local ret = $(esc(ex))
        scalar_allowed[] = prev
        ret
    end
end

@doc (@doc @allowscalar) ->
macro disallowscalar(ex)
    quote
        local prev = scalar_allowed[]
        scalar_allowed[] = ScalarDisallowed
        local ret = $(esc(ex))
        scalar_allowed[] = prev
        ret
    end
end

@doc (@doc @allowscalar) ->
function allowscalar(f::Base.Callable, allow::Bool=true, warn::Bool=false)
    prev = scalar_allowed[]
    allowscalar(allow, warn)
    ret = f()
    scalar_allowed[] = prev
    ret
end

# basic indexing

Base.IndexStyle(::Type{<:AbstractGPUArray}) = Base.IndexLinear()

function Base.getindex(xs::AbstractGPUArray{T}, i::Integer) where T
    assertscalar("scalar getindex")
    x = Array{T}(undef, 1)
    copyto!(x, 1, xs, i, 1)
    return x[1]
end

function Base.setindex!(xs::AbstractGPUArray{T}, v::T, i::Integer) where T
    assertscalar("scalar setindex!")
    x = T[v]
    copyto!(xs, i, x, 1, 1)
    return xs
end

Base.setindex!(xs::AbstractGPUArray, v, i::Integer) = xs[i] = convert(eltype(xs), v)

# Vector indexing

to_index(a, x) = x
to_index(a::A, x::Array{ET}) where {A, ET} = copyto!(similar(a, ET, size(x)...), x)
to_index(a, x::UnitRange{<: Integer}) = convert(UnitRange{Int}, x)
to_index(a, x::Base.LogicalIndex) = error("Logical indexing not implemented")

@generated function index_kernel(ctx::AbstractKernelContext, dest::AbstractArray, src::AbstractArray, idims, Is)
    N = length(Is.parameters)
    quote
        i = @linearidx dest
        is = CartesianIndices(idims)[i]
        @nexprs $N i -> @inbounds I_i = Is[i][is[i]]
        @inbounds dest[i] = @ncall $N getindex src i -> I_i
        return
    end
end

function Base._unsafe_getindex!(dest::AbstractGPUArray, src::AbstractGPUArray, Is::Union{Real, AbstractArray}...)
    if length(Is) == 1 && isa(first(Is), Array) && isempty(first(Is)) # indexing with empty array
        return dest
    end
    idims = map(length, Is)
    gpu_call(index_kernel, dest, src, idims, map(x-> to_index(dest, x), Is))
    return dest
end

# FIXME: simple broadcast getindex like function... reuse from Base
@inline bgetindex(x::AbstractArray, i) = x[i]
@inline bgetindex(x, i) = x

@generated function setindex_kernel!(ctx::AbstractKernelContext, dest::AbstractArray, src, idims, Is, len)
    N = length(Is.parameters)
    idx = ntuple(i-> :(Is[$i][is[$i]]), N)
    quote
        i = linear_index(ctx)
        i > len && return
        is = CartesianIndices(idims)[i]
        @inbounds setindex!(dest, bgetindex(src, i), $(idx...))
        return
    end
end

function Base._unsafe_setindex!(::IndexStyle, dest::T, src, Is::Union{Real, AbstractArray}...) where T <: AbstractGPUArray
    if length(Is) == 1 && isa(first(Is), Array) && isempty(first(Is)) # indexing with empty array
        return dest
    end
    idims = length.(Is)
    len = prod(idims)
    src_gpu = adapt(T, src)
    gpu_call(setindex_kernel!, dest, src_gpu, idims, map(x-> to_index(dest, x), Is), len;
             total_threads=len)
    return dest
end


text:
using Distributed
import GZip
import JLD: load, save
import FileIO: File, @format_str
import Blosc
import JSON
using Pkg
include("db_graph.jl")

# Complement a set; useful for 'compressing' large allele sets when building the scheme DB.
function complement_alleles(vector, m)
  comp_vector = Int16[]
  expected::Int16 = 1
  i::Int16 = 1
  l::Int16 = length(vector)
  while i <= l
    diff = vector[i] - expected
    if diff < 0
      i += 1
      continue
    end
    while diff > 0
      push!(comp_vector, expected)
      expected += 1
      diff -= 1
    end
    i += 1
    expected +=1
  end
  while expected <= m
    push!(comp_vector, expected)
    expected += 1
  end
  return comp_vector
end

# Combine all loci kmer classifications in a format that is easy to serialize to disk.
function combine_loci_classification(::Type{DNAKmer{k}}, results, loci) where {k}
  # Arrays to store the kmer classification votes:
  # all kmers in one big list:
  kmer_list = String[]

  # list with elements n1, l1, n2, l2, ..., where n1 is the number of kmers on
  # the kmer_list that corresponds to the the locus l1
  loci_list = Int32[]

  # weight of each kmer:
  weight_list = Int16[]

  # n1, a_1, ..., a_n1, n2, b_1, ..., b_n2, ...
  # n1 is the number of alleles corresponding to the next kmer, a_1,...,a_n1 are
  # the alleles getting the vote for the kmer;
  alleles_list = Int16[]

  # Store the allele ids for each locus; In the DB, we index the alleles from 1 to n,
  # so we have to store the original ID to translate back in the output:
  # same format as alleles list: a list size (# of alleles), then the list of alleles;
  allele_ids_per_locus = Int[]

  # locus to int.
  l2int = Dict{String,Int16}(locus => idx for (idx, locus) in enumerate(loci))
  weight::Int16 = 0

  # allele coverages: how many kmers cover each allele; same format as alleles_list;
  allele_coverage_per_locus = Int16[]

  for (locus,(kmer_class, allele_ids, kmer_weights, allele_coverages)) in zip(loci,results)
    n_alleles = length(allele_ids)
    half = n_alleles/2
    n_kmers_in_class = 0
    for (kmer, allele_list) in kmer_class
      l_as = length(allele_list)
      if l_as == n_alleles
        continue
      elseif l_as > half
        allele_list = complement_alleles(allele_list, n_alleles)
        if isempty(allele_list)
          continue
        end
        weight = -get(kmer_weights, kmer, 1)
      else
        weight = get(kmer_weights, kmer, 1)
      end
      # kmer as string:
      dnakmer = DNAKmer{k}(kmer)
      push!(kmer_list, "$dnakmer")
      push!(weight_list, weight)
      push!(alleles_list, length(allele_list))
      append!(alleles_list, allele_list)
      n_kmers_in_class += 1
    end
    push!(loci_list, n_kmers_in_class)
    push!(loci_list, l2int[locus])
    push!(allele_ids_per_locus, n_alleles)
    append!(allele_ids_per_locus, allele_ids)
    append!(allele_coverage_per_locus, allele_coverages)
    # println("MIN COV:", m_val, "MAX VAR:", maximum(allele_coverages), "VAR>0: ", [x for x in allele_coverages if x > 0 ])
    # println("MIN COV:", m_val, " MAX VAR:", maximum(allele_coverages))
  end
  return (loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus, allele_coverage_per_locus)

end

function kmer_class_for_each_locus(::Type{DNAKmer{k}}, files::Vector{String}, coverage) where {k}
  loci = [splitext(basename(file))[1] for file in files]
  if length(workers()) == 1
    results = [build_db_graph(DNAKmer{k}, file, coverage) for file in files]
  else
    results = pmap(file->build_db_graph(DNAKmer{k}, file, coverage), files)
  end
  return results, loci
end

function save_db(::Type{DNAKmer{k}}, kmer_db, loci, filename, profile, args, version) where {k}

  loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus, al_coverages = kmer_db
  d = Dict(
    "mentalist_version" => version,
    "loci_list"=> Blosc.compress(loci_list),
    "weight_list" => Blosc.compress(weight_list),
    "allele_ids_per_locus" => Blosc.compress(allele_ids_per_locus),
    "allele_coverages" => Blosc.compress(al_coverages),
    "kmer_list" => join(kmer_list,""),
    "loci"=>loci,
    "args"=>JSON.json(args)
  )
  # alleles list: potentially larger than 2G, so check limit:
  limit = floor(Int,2000000000/sizeof(alleles_list[1])) # Blosc limit is 2147483631, slightly less than 2G (2147483648). Using 2000000000 just to make it easier on the eye. :)
  if sizeof(alleles_list) <= limit
    d["alleles_list"] = Blosc.compress(alleles_list)
  else
    # we have to break this into smaller chunks:
    n_chunks = ceil(Int, length(alleles_list)/limit)
    l = length(alleles_list)
    for i = 0:n_chunks-1
      name = "alleles_list_$i"
      s = i*limit + 1
      e = min((i+1)*limit,l)
      d[name] = Blosc.compress(alleles_list[s:e])
    end
  end
  # mkdir:
  mkpath(dirname(filename))
  # database:
  save(File(format"JLD", "$filename"), d)
  # Profile:
  if profile != nothing
    cp(profile, "$filename.profile", force=true)
  end
end

max_stars_repo_name: CNRResistanceAntibiotic/MentaLiST

max_issues_repo_name: CNRResistanceAntibiotic/MentaLiST

max_forks_repo_name: CNRResistanceAntibiotic/MentaLiST


text:
module RocketFunctionActorTest

using Test
using Rocket

@testset "FunctionActor" begin

    println("Testing: actor FunctionActor")

    @testset begin
        source = from([ 1, 2, 3 ])
        values = Int[]
        actor  = (t::Int) -> push!(values, t)

        subscribe!(source, actor)

        @test values == [ 1, 2, 3 ]
    end

    @testset begin
        source = faulted(Int, "Error")
        values = Int[]
        actor  = (t::Int) -> push!(values, t)

        @test_throws ErrorException subscribe!(source, actor)
        @test values == []
    end

end

end


text:
# ref
# https://github.com/JuliaDebug/Debugger.jl/blob/c01f25c14317e6b8a4e72cab8418d6112465a369/src/repl.jl

using Yao, YaoBlocks, MLStyle
using REPL, REPL.LineEdit

export yaorepl

include("parse.jl")

function yaorepl(repl = nothing, terminal = nothing)
    if repl === nothing
        if !isdefined(Base, :active_repl)
            error("Debugger.jl needs to be run in a Julia REPL")
        end
        repl = Base.active_repl
    end
    if !isa(repl, REPL.LineEditREPL)
        error("Debugger.jl requires a LineEditREPL type of REPL")
    end

    if terminal === nothing
        terminal = Base.active_repl.t
    end
    state = YaoREPLState()

    normal_prefix = Sys.iswindows() ? "\e[33m" : "\e[38;5;166m"
    panel = LineEdit.Prompt("yao> ";
        prompt_prefix = () -> normal_prefix,
        prompt_suffix = Base.text_colors[:normal],
        complete = REPL.REPLCompletionProvider(),
        on_enter = s->true)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:debug => panel))
    REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()
    standard_keymap = Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap(standard_keymap)

    panel.on_done = (s,buf,ok)->begin
        line = String(take!(buf))
        if !ok || strip(line) == "q"
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if isempty(strip(line))
            command = panel.hist.history[end]
        else
            command = strip(line)
        end
        do_print_status = true
        cmd1 = split(command,' ')[1]
        try
            msg = yaorepl_handler(command, state)
            if msg isa AbstractString
                println(msg)
                println()
            end
        catch err
            Base.display_error(err)
            LineEdit.reset_state(s)
        end
        LineEdit.reset_state(s)
        return true
    end

    REPL.run_interface(terminal, LineEdit.ModalInterface([panel,search_prompt]))
    return state
end


text:
struct TransformedDistribution{T, U} <: Distribution{T}
    base :: Distribution{U}
    # How many more parameters does this distribution have
    # than the base distribution?
    nArgs :: Int8
    # forward is a U, arg... -> T function,
    # and backward is a T, arg... -> U function,
    # such that for any `args`, we have
    # backward(forward(u, args...), args...) == u
    # and
    # forward(backward(t, args...), args...) == t.
    # Note that if base is a continuous distribution, then
    # forward and backward must be differentiable.
    forward :: Function
    backward :: Function
    backward_grad :: Function
end

function random(d::TransformedDistribution{T, U}, args...)::T where {T, U}
    d.forward(random(d.base, args[d.nArgs+1:end]...), args[1:d.nArgs]...)
end

function logpdf_correction(d::TransformedDistribution{T, U}, x, args) where {T, U}
    log(abs(d.backward_grad(x, args...)[1]))
end

function logpdf(d::TransformedDistribution{T, U}, x::T, args...) where {T, U}
    orig_x = d.backward(x, args[1:d.nArgs]...)
    orig_logpdf = logpdf(d.base, orig_x, args[d.nArgs+1:end]...)

    if is_discrete(d.base)
        orig_logpdf
    else
        orig_logpdf + logpdf_correction(d, x, args[1:d.nArgs])
    end
end

function logpdf_grad(d::TransformedDistribution{T, U}, x::T, args...) where {T, U}
    orig_x = d.backward(x, args[1:d.nArgs]...)
    base_grad = logpdf_grad(d.base, orig_x, args[d.nArgs+1:end]...)

    if is_discrete(d.base) || !has_output_grad(d.base)
        # TODO: should this be nothing or 0?
        [base_grad[1], fill(nothing, d.nArgs)..., base_grad[2:end]...]
    else
        transformation_grad = d.backward_grad(x, args[1:d.nArgs]...)
        correction_grad = ReverseDiff.gradient(v -> logpdf_correction(d, v[1], v[2:end]), [x, args[1:d.nArgs]...])
        # TODO: Will this sort of thing work if the arguments w.r.t. which we are taking
        # gradients are themselves vector-valued?
        full_grad = (transformation_grad .* base_grad[1]) .+ correction_grad
        [full_grad..., base_grad[2:end]...]
    end
end

is_discrete(d::TransformedDistribution{T, U}) where {T, U} = is_discrete(d.base)

(d::TransformedDistribution{T, U})(args...) where {T, U} = random(d, args...)

function has_output_grad(d::TransformedDistribution{T, U}) where {T, U}
    has_output_grad(d.base)
end

function has_argument_grads(d::TransformedDistribution{T, U}) where {T, U}
    if is_discrete(d.base) || !has_output_grad(d.base)
        [fill(false, d.nArgs)..., has_argument_grads(d.base)...]
    else
        [fill(true, d.nArgs)..., has_argument_grads(d.base)...]
    end
end


text:
# """
# Function applies the transformations for loading a parameter to
# fp = (quote end).args
# sp = (quote end).args
# M == param_length == 0 => scalar
# """
function load_transformations!(
    fp, sp, b::Bounds{T}, out, shape::Vector{Int},
    partial::Bool, logjac::Bool, sptr::Union{Symbol,Nothing},
    m::Module = DistributionParameters,
    θ = Symbol("##θparameter##"), ∂θ = Symbol("##∂θparameter##"),
    exportparam::Bool = false
) where {T}
    maybe_align = x -> exportparam ? x : VectorizationBase.align(x)
    N = length(shape)
    scalar = iszero(N)
    M = prod(shape)
    if !scalar
        X = similar(shape); X[1] = 1
        for n in 2:N
            X[n] = X[n-1] * shape[n-1]
        end
    end
    use_sptr = sptr isa Symbol
    outinit = if scalar
        quote end
    elseif use_sptr
        if N == 1
            quote
                $out = $m.PtrVector{$(first(shape)),$T}(pointer($sptr, $T))
                $sptr += $(maybe_align(sizeof(T)*M))
            end
        else
            quote
                $out = $m.PtrArray{$(Tuple{shape...}),$T,$N,$(Tuple{X...}),$M,true}(pointer($sptr, $T))
                $sptr += $(maybe_align(sizeof(T)*M))
            end
        end
    else
        quote # Do we want to pad these?
            $out = $m.FixedSizeArray{$(Tuple{shape...}),$T,$N,$(Tuple{X...}),$M}(undef)
        end
    end
    outlifestart = if !scalar && !exportparam && use_sptr
        push!(outinit.args, :($m.lifetime_start!($out)))
        true
    else
        false
    end
    adjout = ReverseDiffExpressionsBase.adj(out)
    if isunbounded(b)
        if scalar
            push!(fp, :($out = $m.VectorizationBase.load($θ)))
        elseif exportparam
            isym = gensym(:i)
            loop_quote = quote
                LoopVectorization.@vvectorize_unsafe $T $m for $isym ∈ 1:$M
                    $out[$isym] = $θ[$isym]
                end
            end
            copy_q = quote
                $outinit
                $(macroexpand(m, loop_quote))
            end
            push!(fp, copy_q)
        else
            push!(fp, :($out = $m.PtrArray{$(Tuple{shape...}),$T,$N,$(Tuple{X...}),$M,true}(pointer($θ))))
        end
    elseif islowerbounded(b)
        logout = gensym(Symbol(:log_, out))
        if scalar
            outdef = ((b.lb == zero(T)) && !exportparam ) ? :(RealFloat{$(Bounds(zero(T),typemax(T))),$T,$T}(Base.exp($logout),$logout)) : :(exp($logout) + $(T(b.lb)))
            load_expr = quote
                $logout = $m.VectorizationBase.load($θ)
                $out = $outdef
            end
            logjac && push!(load_expr.args, :(target = $m.vadd(target, $logout)))
        else
            outdef = (b.lb == zero(T)) ? :(exp($logout)) : :(exp($logout) + $(T(b.lb)))
            isym = gensym(:i)
            loopbody = quote
                $logout = $θ[$isym]
                $out[$isym] = $outdef
            end
            logjac && push!(loopbody.args, :(target = $m.vadd(target, $logout)))
            loop_quote = quote
                LoopVectorization.@vvectorize_unsafe $T $m for $isym ∈ 1:$M
                    $loopbody
                end
            end
            if b.lb == zero(T) && !exportparam && sptr !== nothing
                outinit = quote
                    $out = $m.RealArray{$(Tuple{shape...}),$(Bounds(zero(T),typemax(T))),$T,$N,$(Tuple{X...}),$M,Ptr{$T}}(pointer($sptr, $T),pointer($θ))
                    $sptr += $(maybe_align(sizeof(T)*M))
                end
            end
            load_expr = quote
                $outinit
                $(macroexpand(m, loop_quote))
            end
        end
        push!(fp, load_expr)
        # (exportparam || scalar) || push!(fp, load_expr)
        if partial
            if scalar
                # push!(sp, :($m.VectorizationBase.store!($∂θ, muladd($adjout, $out, one($T)))))
                push!(sp, :($m.VectorizationBase.store!($adjout, muladd($m.VectorizationBase.load($adjout), $out, one($T)))))
            else
                isym = gensym(:i)
                storeloop = quote
                    LoopVectorization.@vvectorize_unsafe $T $m for $isym ∈ 1:$M
                        $adjout[$isym] = $m.SIMDPirates.vmuladd($adjout[$isym], $out[$isym], one($T))
                    end
                end
                push!(sp, macroexpand(m, storeloop))
            end
        end
        outlifestart && push!(sp, :($m.lifetime_end!($out)))
    elseif isupperbounded(b)
        logout = gensym(Symbol(:log_, out))
        outdef = ((b.ub == zero(T)) && !exportparam) ? :(- exp($logout)) : :($(T(b.ub)) - exp($logout))
        if scalar
            load_expr = quote
                $logout = $m.VectorizationBase.load($θ)
                $out = $outdef
            end
            logjac && push!(load_expr.args, :(target = $m.vadd(target, $logout)))
        else
            isym = gensym(:i)
            loopbody = quote
                $logout = $θ[$isym]
                $out[$isym] = $outdef
            end
            logjac && push!(loopbody.args, :(target = $m.vadd(target, $logout)))
            loop_quote = quote
                LoopVectorization.@vvectorize_unsafe $T $m for $isym ∈ 1:$M
                    $loopbody
                end
            end
            load_expr = quote
                $outinit
                $(macroexpand(m, loop_quote))
            end
        end
        push!(fp, load_expr)
        if partial
            if scalar
                push!(sp, :($m.VectorizationBase.store!($adjout, $m.SIMDPirates.vfnmadd($m.VectorizationBase.load($adjout), $out, one($T)))))
            else
                isym = gensym(:i)
                storeloop = quote
                    LoopVectorization.@vvectorize_unsafe $T $m for $isym ∈ 1:$M
                        $adjout[$isym] = $m.SIMDPirates.vfnmadd($adjout[$isym], $out[$isym], one($T))
                    end
                end
                push!(sp, macroexpand(m, storeloop))
            end
        end
        outlifestart && push!(sp, :($m.lifetime_end!($out)))
    elseif isbounded(b)
        scale = b.ub - b.lb
        invlogit = gensym(:invlogit)
        ninvlogit = gensym(:ninvlogit)
        ∂invlogit = gensym(:∂invlogit)
        if scalar
            if b.lb == 0 && b.ub == 1
                unconstrained = gensym(:unconstrained)
                q = quote
                    $unconstrained = $m.VectorizationBase.load($θ)
                    $ninvlogit = one($T) / (one($T) + exp($unconstrained))
                    $invlogit = one($T) - $ninvlogit
                    $out = RealFloat{$(Bounds(zero(T),one(T))),$T,$T}($invlogit, $unconstrained)
                    $∂invlogit = $invlogit * $ninvlogit
                end
            else
                q = quote
                    $ninvlogit = one($T) / (one($T) + exp($m.VectorizationBase.load($θ)))
                    $invlogit = one($T) - $ninvlogit
                    $out = $(b.lb == zero(T) ?
                             (scale == one(T) ? invlogit : :($scale * $invlogit)) :
                             (scale == one(T) ? :($(b.lb) + $invlogit) : :(muladd($scale, $invlogit, $(b.lb)))))
                    $∂invlogit = $invlogit * $ninvlogit
                end
            end
            logjac && push!(q.args, :( target = $m.vadd(target, $∂invlogit)))
            push!(fp, q)
            if partial
                ∂q = quote
                    VectorizationBase.store!(
                        $adjout, @fastmath one($T) - $(T(2)) * $invlogit +
                        $m.VectorizationBase.load($adjout) * $(scale == one(T) ? ∂invlogit : :($scale * $∂invlogit))
                    )
                end
                push!(sp, macroexpand(m, ∂q)) #Expand the fastmath? Why not Base.FastMath.add_fast / mul_fast directly?
            end
        else
            isym = gensym(:i)
            if b.lb == 0 && b.ub == 1
                outinit = quote
                    $out = $m.RealArray{$(Tuple{shape...}),$(Bounds(zero(T),one(T))),$T,$N,$(Tuple{X...}),$M,Ptr{$T}}(pointer($sptr, $T), pointer($θ))
                    $sptr += $(maybe_align(sizeof(T)*M))
                end
                outlifestart && push!(outinit.args, :($m.lifetime_start!($out)))
            end
            if partial
                invlogitinits = if use_sptr
                    quote
                        $outinit
                        $invlogit = $m.PtrVector{$M,$T}(pointer($sptr,$T))
                        $sptr += $(maybe_align(M*sizeof(T)))
                        $∂invlogit = $m.PtrVector{$M,$T}(pointer($sptr,$T))
                        $sptr += $(maybe_align(M*sizeof(T)))
                    end
                else
                    quote
                        $outinit
                        $invlogit = FixedSizeVector{$M,$T}(undef)
                        $∂invlogit = FixedSizeVector{$M,$T}(undef)
                    end
                end
                if outlifestart
                    push!(invlogitinits.args, :($m.lifetime_start!($invlogit)))
                    push!(invlogitinits.args, :($m.lifetime_start!($∂invlogit)))
                end
                push!(fp, invlogitinits)
                ilt = gensym(:ilt)
                ∂ilt = gensym(:∂ilt)
                loop_body = quote
                    $ninvlogit = one($T) / (one($T) + $m.SLEEFPirates.exp($θ[$isym]))
                    $ilt = one($T) - $ninvlogit
                    $invlogit[$isym] = $ilt
                    $∂ilt = $ilt * $ninvlogit
                    $∂invlogit[$isym] = $∂ilt
                    $out[$isym] = $(b.lb == zero(T) ?
                                    (scale == one(T) ? ilt : :($scale * $ilt)) :
                                    (scale == one(T) ? :($(b.lb) + $ilt) : :($m.SIMDPirates.vmuladd($scale, $ilt, $(b.lb)))))
                end
                logjac && push!(loop_body.args, :(target = $m.vadd(target, $m.SLEEFPirates.log($∂ilt))))
                loop_quote = quote
                    LoopVectorization.@vvectorize_unsafe $T $((m)) for $isym ∈ 1:$M
                        $loop_body
                    end
                end
                push!(fp, macroexpand(m, loop_quote))
                storeloop = quote
                    LoopVectorization.@vvectorize_unsafe $T $((m)) for $isym ∈ 1:$M
                        $adjout[$isym] = one($T) - $(T(2)) * $invlogit[$isym] +
                            ($adjout)[$isym] * $(scale == one(T) ? :($∂invlogit[$isym]) : :($∂invlogit[$isym] * $scale))
                    end
                end
                push!(sp, macroexpand(m, storeloop))
                if outlifestart
                    push!(sp, :($m.lifetime_end!($invlogit)))
                    push!(sp, :($m.lifetime_end!($∂invlogit)))
                end
            else
                loop_body = quote
                    $ninvlogit = one($T) / (one($T) + $m.SLEEFPirates.exp($θ[$isym]))
                    $invlogit = one($T) - $ninvlogit
                    $out[$isym] = $(b.lb == zero(T) ?
                                    (scale == one(T) ? invlogit : :($scale * $invlogit)) :
                                    (scale == one(T) ? :($(b.lb) + $invlogit) : :($m.SIMDPirates.vmuladd($scale, $invlogit, $(b.lb)))))
                end
                logjac && push!(loop_body.args, :(target = $m.vadd(target, $m.SLEEFPirates.log($invlogit * $ninvlogit))))
                loop_quote = quote
                    $outinit
                    LoopVectorization.@vvectorize_unsafe $T $((m)) for $isym ∈ 1:$M
                        $loop_body
                    end
                end
                push!(fp, macroexpand(m, loop_quote))
            end
            outlifestart && push!(sp, :($m.lifetime_end!($out)))
        end
    end
    if partial
        if scalar
            push!(fp, :($adjout = $∂θ))
        else
            push!(fp, :($adjout = $m.PtrArray{$(Tuple{shape...}),$T,$N,$(Tuple{X...}),$M,true}(pointer($∂θ))))
        end
    end
    if exportparam && scalar
        push!(fp, :($m.VectorizationBase.store!(pointer($sptr, $T), convert($T, $out)); $sptr += $(sizeof(T))))
    end
    push!(fp, :($θ += $M))
    partial && push!(fp, :($∂θ += $M))
    nothing
end


text:
using ArnoldiMethod, Random, FFTW

@testset "basic KPM" begin
    h = LatticePresets.honeycomb() |>
        hamiltonian(hopping(1.0I, range = 1/sqrt(3)), orbitals = (Val(1), Val(2))) |>
        unitcell(region = RegionPresets.circle(10))
    brange = bandrangeKPM(h)
    m1 = momentaKPM(h, order = 30, ket = randomkets(1, maporbitals = true))
    m2 = momentaKPM(h, bandrange = brange, order = 30, ket = randomkets(1, maporbitals = true))
    @test abs(m1.mulist[1]) ≈ abs(m2.mulist[1])
    dos1 = dosKPM(h, order = 10, bandrange = (-3,3), ket = randomkets(1, maporbitals = true))
    dos2 = dosKPM(h, order = 10, bandrange = (-3,3), ket = randomkets(1, r -> randn(), maporbitals = true))
    @test all(>(0), last(dos1))
    @test all(>(0), last(dos2))
    dos = dosKPM(h, order = 10, bandrange = (-3,3), ket = randomkets(2, r -> randn() * SA[1 0; 0 -1], sublats = :B))
    @test all(>(0), last(dos))
    x, energy = densityKPM(h, h, order = 100, resolution = 1, bandrange = (-3,3), ket = ketmodel(1, indices = 1, maporbitals = true))
    @test isapprox(abs(sum(energy)), 0, atol = Quantica.default_tol(eltype(energy)))
    k = ket(ketmodel(1, indices = 1, maporbitals = true), h)
    x´, energy´ = densityKPM(h, h, order = 100, resolution = 1, bandrange = (-3,3), ket = k)
    @test x´ ≈ x
    @test energy´ ≈ energy
    x´, energy´ = densityKPM(h, h[], order = 100, resolution = 1, bandrange = (-3,3), ket = k)
    @test x´ ≈ x
    @test energy´ ≈ energy
end

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
###############################################################################
#
# Calculate end pressure of a hydraulic system
#
###############################################################################
"""
    end_pressure(
                leaf::LeafHydraulics{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    end_pressure(
                root::RootHydraulics{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    end_pressure(
                stem::StemHydraulics{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    end_pressure(
                tree::TreeSimple{FT},
                flow::FT
    ) where {FT<:AbstractFloat}
    end_pressure(
                tree::TreeSimple{FT},
                f_sl::FT,
                f_sh::FT,
                r_sl::FT
    ) where {FT<:AbstractFloat}

Return the xylen end pressure(s) from flow rate(s), given
- `leaf` [`LeafHydraulics`](@ref) type struct
- `flow` Flow rate (per leaf area for [`LeafHydraulics`](@ref))
- `root` [`RootHydraulics`](@ref) type struct
- `stem` [`StemHydraulics`](@ref) type struct
- `tree` [`TreeSimple`](@ref) type struct
- `f_sl` Flow rate to sunlit leaves
- `f_sh` Flow rate to shaded leaves
- `r_sl` Fraction of sunlit leaves

Note, gravity is accounted for in root and stem; rhizosphere conductance is
    accounted for in root; extra-xylary conductance is not accounted for in
    leaf here because it calculates xylem end pressure.
"""
function end_pressure(
            leaf::LeafHydraulics{FT},
            flow::FT
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, p_history, p_ups, vc = leaf;

    p_end::FT = p_ups;

    # compute k from temperature and history, then update pressure
    for (_k, _kh, _ph) in zip(k_element, k_history, p_history)
        p_25 = p_end / f_st;
        if p_25 < _ph
            k = xylem_k_ratio(vc, p_25, f_vis) * _k;
        else
            k = _kh * _k / f_vis;
        end
        p_end -= flow / k;
    end

    return p_end
end

function end_pressure(
            root::RootHydraulics{FT},
            flow::FT
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, k_rhiz, p_gravity, p_history,
            p_osm, p_ups, sh, T_sap, vc = root;

    # make sure that p_ups is not p_25 and then convert
    p_end::FT = p_ups;
    p_25 ::FT = p_end / f_st;

    # compute pressure drop along rhizosphere, using p_25 for Θ
    _dp = flow / k_rhiz * f_vis / 10;
    for i in 1:10
        # No idea why soil_k_ratio_p25 results in unnecessary allocations
        # _f  = soil_k_ratio_p25(sh, p_25);
        _rwc  = soil_rwc(sh, p_25);
        _f    = soil_k_ratio_rwc(sh, _rwc);
        p_25 -= _dp / _f;
    end
    p_end = p_25 * f_st + p_osm * T_sap / K_25(FT);

    # compute k from temperature and history, then update pressure
    for (_k, _kh, _pg, _ph) in zip(k_element, k_history, p_gravity, p_history)
        p_25 = p_end / f_st;
        if p_25 < _ph
            k = xylem_k_ratio(vc, p_25, f_vis) * _k;
        else
            k = _kh * _k / f_vis;
        end
        p_end -= flow / k + _pg;
    end

    return p_end
end

function end_pressure(
            stem::StemHydraulics{FT},
            flow::FT
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, p_gravity, p_history, p_ups,
            vc = stem;

    p_end::FT = p_ups;

    # compute k from temperature and history, then update pressure
    for (_k, _kh, _pg, _ph) in zip(k_element, k_history, p_gravity, p_history)
        p_25 = p_end / f_st;
        if p_25 < _ph
            k = xylem_k_ratio(vc, p_25, f_vis) * _k;
        else
            k = _kh * _k / f_vis;
        end
        p_end -= flow / k + _pg;
    end

    return p_end
end

function end_pressure(
            tree::TreeSimple{FT},
            flow::FT
) where {FT<:AbstractFloat}
    # calculate the p_dos for roots
    p_dos = end_pressure(tree.root, flow);
    (tree.stem).p_ups = p_dos;

    # calculate the p_dos for stem
    p_dos = end_pressure(tree.stem, flow);
    (tree.leaf).p_ups = p_dos;

    # calculate the p_dos for leaves
    p_dos = end_pressure(tree.leaf, flow / (tree.leaf).area);

    return p_dos
end

function end_pressure(
            tree::TreeSimple{FT},
            f_sl::FT,
            f_sh::FT,
            r_sl::FT
) where {FT<:AbstractFloat}
    flow = f_sl + f_sh;

    # calculate the p_dos for roots
    p_dos = end_pressure(tree.root, flow);
    (tree.stem).p_ups = p_dos;

    # calculate the p_dos for stem
    p_dos = end_pressure(tree.stem, flow);
    (tree.leaf).p_ups = p_dos;

    # calculate the p_dos for leaves
    p_sl = end_pressure(tree.leaf, f_sl / (tree.leaf).area / r_sl);
    p_sh = end_pressure(tree.leaf, f_sh / (tree.leaf).area / (1-r_sl));

    return p_sl,p_sh
end

###############################################################################
#
# Update pressure profile along a hydraulic system
#
###############################################################################
"""
    pressure_profile!(
                leaf::LeafHydraulics{FT},
                flow::FT;
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                root::RootHydraulics{FT},
                flow::FT;
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                root::RootHydraulics{FT},
                q_in::FT,
                flow::Array{FT,1};
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                stem::StemHydraulics{FT},
                flow::FT;
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                stem::StemHydraulics{FT},
                flow::Array{FT,1};
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                tree::TreeSimple{FT},
                p_soil::FT,
                flow::FT;
                update::Bool = true
    ) where {FT<:AbstractFloat}
    pressure_profile!(
                tree::TreeSimple{FT},
                p_soil::FT,
                f_sl::FT,
                f_sh::FT,
                r_sl::FT;
                update::Bool = true
    ) where {FT<:AbstractFloat}

Update the pressure profile, given
- `leaf` [`LeafHydraulics`](@ref) type struct
- `root` [`RootHydraulics`](@ref) type struct
- `stem` [`StemHydraulics`](@ref) type struct
- `tree` [`TreeSimple`](@ref) type struct
- `p_soil` Soil water potential
- `flow` Flow rate (per leaf area for [`LeafHydraulics`](@ref))
- `f_sl` Flow rate to sunlit leaves
- `f_sh` Flow rate to shaded leaves
- `r_sl` Fraction of sunlit leaves
- `update` Optional. If true, update drought legacy. Default is `true` for
    [`TreeSimple`](@ref) but `false` for others using capacitance
"""
function pressure_profile!(
            leaf::LeafHydraulics{FT},
            flow::FT;
            update::Bool = true
) where {FT<:AbstractFloat}
    # update leaf flow
    leaf.flow = flow;

    @unpack k_element, k_history, p_history, p_ups, f_st, f_vis, vc = leaf;

    p_end::FT = p_ups;

    # compute k from temperature and history, then update pressure
    for i in eachindex(k_element)
        # update history first
        p_25 = p_end / f_st;
        if p_25 < p_history[i]
            _kr = xylem_k_ratio(vc, p_25, f_vis);
            if update
                leaf.p_history[i] = p_25;
                leaf.k_history[i] = _kr * f_vis;
            end
            k = _kr * k_element[i];
        else
            k = k_history[i] / f_vis * k_element[i];
        end

        # then use historical minimal k
        p_end -= flow / k;

        leaf.p_element[i] = p_end;
    end

    # update the leaf xylem end pressure
    leaf.p_dos  = p_end;
    leaf.p_leaf = p_end - flow / leaf.k_ox;

    return nothing
end

function pressure_profile!(
            root::RootHydraulics{FT},
            flow::FT;
            update::Bool = true
) where {FT<:AbstractFloat}
    # update root flow rate
    root.flow = flow;

    @unpack f_st, f_vis, k_element, k_history, k_rhiz, p_gravity, p_history,
            p_osm, p_ups, sh, T_sap, vc = root;

    # make sure that p_ups is not p_25 and then convert
    p_end::FT = p_ups;
    p_25 ::FT = p_end / f_st;

    # compute pressure drop along rhizosphere, using p_25 for Θ
    _dp = flow / k_rhiz * f_vis / 10;
    for i in 1:10
        # No idea why soil_k_ratio_p25 results in unnecessary allocations
        # _f  = soil_k_ratio_p25(sh, p_25);
        _rwc  = soil_rwc(sh, p_25);
        _f    = soil_k_ratio_rwc(sh, _rwc);
        p_25 -= _dp / _f;
    end
    p_end = p_25 * f_st + p_osm * T_sap / K_25(FT);
    root.p_rhiz = p_end;

    # compute k from temperature and history, then update pressure
    for i in eachindex(k_element)
        # update history first
        p_25 = p_end / f_st;
        if p_25 < p_history[i]
            _kr = xylem_k_ratio(vc, p_25, f_vis);
            if update
                root.p_history[i] = p_25;
                root.k_history[i] = _kr * f_vis;
            end
            k = _kr * k_element[i];
        else
            k = k_history[i] / f_vis * k_element[i];
        end

        # then use historical minimal k
        p_end -= flow / k + p_gravity[i];

        root.p_element[i] = p_end;
    end

    # update the leaf xylem end pressure
    root.p_dos  = p_end;

    return nothing
end

function pressure_profile!(
            root::RootHydraulics{FT},
            q_in::FT,
            flow::Array{FT,1};
            update::Bool = true
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, k_rhiz, p_gravity, p_history,
            p_osm, p_ups, sh, T_sap, vc = root;

    # make sure that p_ups is not p_25 and then convert
    p_end::FT = p_ups;
    p_25 ::FT = p_end / f_st;

    # compute pressure drop along rhizosphere, using p_25 for Θ
    _dp = q_in / k_rhiz * f_vis / 10;
    for i in 1:10
        # No idea why soil_k_ratio_p25 results in unnecessary allocations
        # _f  = soil_k_ratio_p25(sh, p_25);
        _rwc  = soil_rwc(sh, p_25);
        _f    = soil_k_ratio_rwc(sh, _rwc);
        p_25 -= _dp / _f;
    end
    p_end = p_25 * f_st + p_osm * T_sap / K_25(FT);
    root.p_rhiz = p_end;

    # compute k from temperature and history, then update pressure
    for i in eachindex(k_element)
        # update history first
        p_25 = p_end / f_st;
        if p_25 < p_history[i]
            _kr = xylem_k_ratio(vc, p_25, f_vis);
            if update
                root.p_history[i] = p_25;
                root.k_history[i] = _kr * f_vis;
            end
            k = _kr * k_element[i];
        else
            k = k_history[i] / f_vis * k_element[i];
        end

        # then use historical minimal k
        p_end -= flow[i] / k + p_gravity[i];

        root.p_element[i] = p_end;
    end

    # update the leaf xylem end pressure
    root.p_dos  = p_end;

    return nothing
end

function pressure_profile!(
            stem::StemHydraulics{FT},
            flow::FT;
            update::Bool = true
) where {FT<:AbstractFloat}
    # update stem flow
    stem.flow = flow;

    @unpack f_st, f_vis, k_element, k_history, p_gravity, p_history, p_ups,
            vc = stem;

    p_end::FT = p_ups;

    # compute k from temperature and history, then update pressure
    for i in eachindex(k_element)
        # update history first
        p_25 = p_end / f_st;
        if p_25 < p_history[i]
            _kr = xylem_k_ratio(vc, p_25, f_vis);
            if update
                stem.p_history[i] = p_25;
                stem.k_history[i] = _kr * f_vis;
            end
            k = _kr * k_element[i];
        else
            k = k_history[i] / f_vis * k_element[i];
        end

        # then use historical minimal k
        p_end -= flow / k + p_gravity[i];

        stem.p_element[i] = p_end;
    end

    # update the leaf xylem end pressure
    stem.p_dos  = p_end;

    return nothing
end

function pressure_profile!(
            stem::StemHydraulics{FT},
            flow::Array{FT,1};
            update::Bool = true
) where {FT<:AbstractFloat}
    @unpack f_st, f_vis, k_element, k_history, p_gravity, p_history, p_ups,
            vc = stem;

    p_end::FT = p_ups;

    # compute k from temperature and history, then update pressure
    for i in eachindex(k_element)
        # update history first
        p_25 = p_end / f_st;
        if p_25 < p_history[i]
            _kr = xylem_k_ratio(vc, p_25, f_vis);
            if update
                stem.p_history[i] = p_25;
                stem.k_history[i] = _kr * f_vis;
            end
            k = _kr * k_element[i];
        else
            k = k_history[i] / f_vis * k_element[i];
        end

        # then use historical minimal k
        p_end -= flow[i] / k + p_gravity[i];

        stem.p_element[i] = p_end;
    end

    # update the leaf xylem end pressure
    stem.p_dos  = p_end;

    return nothing
end

function pressure_profile!(
            tree::TreeSimple{FT},
            p_soil::FT,
            flow::FT;
            update::Bool = true
) where {FT<:AbstractFloat}
    # update the profile in root
    (tree.root).p_ups = p_soil;
    pressure_profile!(tree.root, flow; update=update);

    # update the profile in stem
    (tree.stem).p_ups = (tree.root).p_dos;
    pressure_profile!(tree.stem, flow; update=update);

    # update the profile in leaf
    (tree.leaf).p_ups = (tree.stem).p_dos;
    pressure_profile!(tree.leaf, flow; update=update);

    return nothing
end

function pressure_profile!(
            tree::TreeSimple{FT},
            p_soil::FT,
            f_sl::FT,
            f_sh::FT,
            r_sl::FT;
            update::Bool = true
) where {FT<:AbstractFloat}
    flow = f_sl + f_sh;

    # update the profile in root
    (tree.root).p_ups = p_soil;
    pressure_profile!(tree.root, flow; update=update);

    # update the profile in stem
    (tree.stem).p_ups = (tree.root).p_dos;
    pressure_profile!(tree.stem, flow; update=update);

    # update the profile in sunlit leaf
    (tree.leaf).p_ups = (tree.stem).p_dos;
    pressure_profile!(tree.leaf, f_sl / (tree.leaf).area / r_sl;
                         update=update);
    pressure_profile!(tree.leaf, f_sh / (tree.leaf).area / (1-r_sl);
                         update=update);

    return nothing
end

function pressure_profile!(
            tree::GrassLikeOrganism{FT},
            mode::NonSteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack leaves, roots = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.q_in, root.q_element; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in leaf
    for leaf in leaves
        leaf.p_ups = p_mean;
        pressure_profile!(leaf, leaf.q_in; update=update);
    end

    return nothing
end

function pressure_profile!(
            tree::GrassLikeOrganism{FT},
            mode::SteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack leaves, roots = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.flow; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in leaf
    for leaf in leaves
        leaf.p_ups = p_mean;
        pressure_profile!(leaf, leaf.flow; update=update);
    end

    return nothing
end

function pressure_profile!(
            tree::PalmLikeOrganism{FT},
            mode::NonSteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack leaves, roots, trunk = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.q_in, root.q_element; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in trunk
    trunk.p_ups = p_mean;
    pressure_profile!(trunk, trunk.q_element; update=update);

    # update the profile in leaf
    for leaf in leaves
        leaf.p_ups = trunk.p_dos;
        pressure_profile!(leaf, leaf.q_in; update=update);
    end

    return nothing
end

function pressure_profile!(
            tree::PalmLikeOrganism{FT},
            mode::SteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack leaves, roots, trunk = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.flow; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in trunk
    trunk.p_ups = p_mean;
    pressure_profile!(trunk, trunk.flow; update=update);

    # update the profile in leaf
    for leaf in leaves
        leaf.p_ups = trunk.p_dos;
        pressure_profile!(leaf, leaf.flow; update=update);
    end

    return nothing
end

function pressure_profile!(
            tree::TreeLikeOrganism{FT},
            mode::NonSteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack branch, leaves, roots, trunk = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.q_in, root.q_element; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in trunk
    trunk.p_ups = p_mean;
    pressure_profile!(trunk, trunk.q_element; update=update);

    # update the profile in leaf
    for i in eachindex(leaves)
        stem = branch[i];
        leaf = leaves[i];
        stem.p_ups = trunk.p_dos;
        pressure_profile!(stem, stem.q_element; update=update);
        leaf.p_ups = stem.p_dos;
        pressure_profile!(leaf, leaf.q_in; update=update);
    end

    return nothing
end

function pressure_profile!(
            tree::TreeLikeOrganism{FT},
            mode::SteadyStateMode;
            update::Bool = false
) where {FT<:AbstractFloat}
    @unpack branch, leaves, roots, trunk = tree;

    # update the profile in roots
    p_mean::FT = 0;
    for root in roots
        pressure_profile!(root, root.flow; update=update);
        p_mean += root.p_dos;
    end
    p_mean /= length(roots);

    # update the profile in trunk
    trunk.p_ups = p_mean;
    pressure_profile!(trunk, trunk.flow; update=update);

    # update the profile in leaf
    for i in eachindex(leaves)
        stem = branch[i];
        leaf = leaves[i];
        stem.p_ups = trunk.p_dos;
        pressure_profile!(stem, stem.flow; update=update);
        leaf.p_ups = stem.p_dos;
        pressure_profile!(leaf, leaf.flow; update=update);
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# Check that we can read the official HDF5 example files
urlbase = "http://www.hdfgroup.org/ftp/HDF5/examples/files/exbyapi/"

using HDF5

fcmp = [0 1 2 3 4 5 6;
    2 1.66667 2.4 3.28571 4.22222 5.18182 6.15385;
    4 2.33333 2.8 3.57143 4.44444 5.36364 6.30769;
    6 3 3.2 3.85714 4.66667 5.54545 6.46154]'
icmp = [0 -1 -2 -3 -4 -5 -6;
    0 0 0 0 0 0 0;
    0 1 2 3 4 5 6;
    0 2 4 6 8 10 12]'
scmp = ["Parting", "is such", "sweet", "sorrow."]
vicmp = Array{Int32}[[3,2,1],[1,1,2,3,5,8,13,21,34,55,89,144]]
opq = Array{Uint8}[[0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x30],
                   [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x31],
                   [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x32],
                   [0x4f, 0x50, 0x41, 0x51, 0x55, 0x45, 0x33]]

const savedir = joinpath(tempdir(), "h5")
if !isdir(savedir)
    mkdir(savedir)
end
function getfile(name)
    file = joinpath(savedir, name)
    if !isfile(file)
        file = download(urlbase*name, file)
    end
    file
end

file = getfile("h5ex_t_floatatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert norm(a - fcmp) < 1e-5
close(fid)

file = getfile("h5ex_t_float.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert norm(d - fcmp) < 1e-5
close(fid)

file = getfile("h5ex_t_intatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert a == icmp
close(fid)

file = getfile("h5ex_t_int.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert d == icmp
close(fid)

file = getfile("h5ex_t_objrefatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
g = fid[a[1]]
@assert isa(g, HDF5Group)
ds2 = fid[a[2]]
ds2v = read(ds2)
@assert isa(ds2v, Array{Int32})
@assert isempty(ds2v)
close(fid)

file = getfile("h5ex_t_objref.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
g = fid[d[1]]
@assert isa(g, HDF5Group)
ds2 = fid[d[2]]
ds2v = read(ds2)
@assert isa(ds2v, Array{Int32})
@assert isempty(ds2v)
close(fid)

file = getfile("h5ex_t_stringatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert a == scmp
close(fid)

file = getfile("h5ex_t_string.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert d == scmp
close(fid)

file = getfile("h5ex_t_vlenatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert a == vicmp
close(fid)

file = getfile("h5ex_t_vlen.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert d == vicmp
close(fid)

file = getfile("h5ex_t_vlstringatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert a == scmp
close(fid)

file = getfile("h5ex_t_vlstring.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert d == scmp
close(fid)

file = getfile("h5ex_t_opaqueatt.h5")
fid = h5open(file, "r")
dset = fid["DS1"]
a = a_read(dset, "A1")
@assert a.tag == "Character array"
@assert a.data == opq
close(fid)

file = getfile("h5ex_t_opaque.h5")
fid = h5open(file, "r")
d = read(fid, "DS1")
@assert d.tag == "Character array"
@assert d.data == opq
close(fid)


text:
"""

Plot the traces and densities of the chains.

$(SIGNATURES)

## Positional arguments
```julia
* `df::DataFrame` # DataFrame holding the draws
* `pars::Vector{Symbol}` # Vector of symbols to be included
```

## Keyword arguments
```julia
* `no_of_chains = 4` # Number of chains
* 'no_of_draws = 1000' # Number of draws (after warmup) in each chain
```

"""
function plot_chains(df::DataFrame, pars::Vector{Symbol};
    no_of_chains=4, no_of_draws=1000)

    df[!, :chain] = repeat(collect(1:no_of_chains); inner=no_of_draws)
    df[!, :chain] = categorical(df.chain)

    fig = Figure()
    for i in 1:length(pars)
        let
            plt = data(df) * visual(Lines) * mapping(pars[i]; color=:chain)
            axis = (; ylabel="$(pars[i])", xlabel="Iteration", title="Traces")
            draw!(fig[i, 1], plt; axis)
        end
        let
            plt = data(df) * mapping(pars[i]; color=:chain) * AlgebraOfGraphics.density()
            axis = (; title="Density $(pars[i])")
            draw!(fig[i, 2], plt; axis)
        end
    end
    return fig
end

export
    plot_chains


text:
#types used by kRPC and type conversion utilities

#part of runtime support
#indicate class and enum membership
module kRPCTypes
import ..krpc
abstract type Class end
abstract type Enum end
struct kStream
    id::UInt64
end
struct kEvent end
struct Service
    status::krpc.schema.Status
end
struct Services
    services::krpc.schema.Services
end
end

#type helper functions

function getWireValue(arg::Union{Int32,Int64,Float32,Float64,UInt32,UInt64})
    opb = PipeBuffer()
    ProtoBuf.write_fixed(opb, arg)
    return readavailable(opb)
end

function getWireValue(arg::Bool)
    opb = PipeBuffer()
    ProtoBuf.write_bool(opb, arg)
    return readavailable(opb)
end

function getWireValue(arg::String)
    opb = PipeBuffer()
    ProtoBuf.write_string(opb, arg)
    return readavailable(opb)
end

function getWireValue(arg::Array{UInt8,1})
    opb = PipeBuffer()
    ProtoBuf.write_bytes(opb, arg)
    return readavailable(opb)
end

function getWireValue(arg::Array{T, 1}) where T
    opb = PipeBuffer()
    ProtoBuf.writeproto(opb, krpc.schema.List(items=map(getWireValue, arg)))
    return readavailable(opb)
end

function getWireValue(arg::Dict{K,V}) where {K, V}
    opb = PipeBuffer()
    op = Array{Array{UInt8,1},1}()
    for kv in arg
        k,v = kv
        push!(op, krpc.schema.DictionaryEntry(key=getWireValue(k), value=getWireValue(v)))
    end
    ProtoBuf.writeproto(opb, krpc.schema.Dictionary(entries=op))
    return readavailable(opb)
end

function getWireValue(arg::Set{T}) where T
    opb = PipeBuffer()
    ProtoBuf.writeproto(opb, krpc.schema.Set(items=map(getWireValue, arg)))
    return readavailable(opb)
end

function getWireValue(arg::Tuple)
    opb = PipeBuffer()
    ProtoBuf.writeproto(opb, krpc.schema.Tuple_(items=convert(Array{Array{UInt8,1},1}, collect(map(getWireValue, arg)))))
    return readavailable(opb)
end

function getWireValue(arg::kRPCTypes.Class)
    res = getWireValue(arg.id)
    return res
end

function getWireValue(arg::kRPCTypes.Enum)
    return getWireValue(Remote.getEnumVal(arg))
end

function getWireValue(arg::Request{S, P, R}) where {S, P, R}
    opb = PipeBuffer()
    ProtoBuf.writeproto(opb, make_kerbal_call(arg))
    return readavailable(opb)
end

function getWireValue(arg::kRPCTypes.kStream)
    opb = PipeBuffer()
    ProtoBuf.writeproto(opb, krpc.schema.Stream(id=arg.id))
    return readavailable(opb)
end

function getJuliaValue(conn, value::krpc.schema.Response, rtype::Type)
    return getJuliaValue(value.return_value, rtype)
end
function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{T}) where {T<:kRPCTypes.Class}
    return rtype(conn, ProtoBuf.read_varint(PipeBuffer(value), Int64))
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Nothing)
    return Nothing()
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{AbstractString})
    return ProtoBuf.read_string(PipeBuffer(value))
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{kRPCTypes.kStream})
    backing = readproto(PipeBuffer(value), krpc.schema.Stream())
    return kRPCTypes.kStream(backing.id)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Union{Type{UInt32},Type{Int32},Type{UInt64},Type{Int64},Type{Float32},Type{Float64}})
    return ProtoBuf.read_fixed(PipeBuffer(value), rtype)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Bool})
    return ProtoBuf.read_bool(PipeBuffer(value))
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Array{UInt8,1}})
    return ProtoBuf.read_bytes(PipeBuffer(value))
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Array{T,1}}) where T
    res = readproto(PipeBuffer(value), krpc.schema.List())
    return T[getJuliaValue(item,T) for item in res.items]
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{T}) where T <: (Tuple{Vararg{T,N} where T where N})
    res = readproto(PipeBuffer(value), krpc.schema.Tuple())
    return ([getJuliaValue(conn,param[1],param[2]) for param in zip(res.items, rtype.parameters)]...,)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Dict{K,V}}) where {K,V}
    res = readproto(PipeBuffer(value), krpc.schema.Dictionary())
    return Dict{K,V}(getJuliaValue(conn,v.key,K) => getJuliaValue(conn,v.value,V) for v in res.entries)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Set{T}}) where {T}
    res = readproto(PipeBuffer(value), krpc.schema.Set())
    return Set{T}(getJuliaValue(conn,item,T) for item in res.items)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{Union{T, Nothing}}) where {T<:kRPCTypes.Class}
    res = ProtoBuf.read_varint(PipeBuffer(value), Int64)
    if res == 0
        return Nothing()
    else
        return T(res)
    end
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{kRPCTypes.Service})
    res = readproto(PipeBuffer(value), krpc.schema.Status())
    return kRPCTypes.Service(res)
end

function getJuliaValue(conn, value::Array{UInt8, 1}, rtype::Type{kRPCTypes.Services})
    res = readproto(PipeBuffer(value), krpc.schema.Services())
    return kRPCTypes.Services(res)
end


text:
"""

```
center_to_origin!(x::Vector{T}, center::T) where T
```

Translates atoms of vectory in x array such that center is in the origin. (`x` is a vector of vectors).

"""
function center_to_origin!(x::Vector{T}, center::T) where T
  for i in eachindex(x)
    x[i] = x[i] - center
  end
  nothing
end


text:
include("SigGenInterface.jl")
include("ParseConfigFiles.jl")

NamedTuple(::Missing) = (object = "missing",)
NamedTuple(d::Dict) = (dict_json_string = json(d),)
Base.convert(T::Type{NamedTuple}, x::Dict) = T(x)
Dict(nt::NamedTuple) = JSON.parse(nt.dict_json_string)

"""
    ssd_write(filename::AbstractString, sim::Simulation)

Converts a [`Simulation`](@ref) to a `NamedTuple` and writes it to a HDF5 file
with a given `filename` using [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl).

## Arguments
* `filename::AbstractString`: Filename of the HDF5 file.
* `sim::Simulation`: [`Simulation`](@ref) that should be written to the HDF5 file.

## Example
```julia
using HDF5
using LegendHDF5IO
using SolidStateDetectors
sim = Simulation(SSD_examples[:InvertedCoax])
simulate!(sim)
ssd_write("example_sim.h5", sim)
```

!!! warn
    If a file with `filename` already exists, it will be overwritten by this method.

!!! note
    In order to use this method, the packages [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) and
    [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl) have to be
    loaded before loading SolidStateDetectors.jl.

See also [`ssd_read`](@ref).
"""
function ssd_write end
export ssd_write

"""
    ssd_read(filename::AbstractString, ::Type{Simulation})

Reads a [`Simulation`](@ref) from a HDF5 file with a given `filename`
using [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl).

## Arguments
* `filename::AbstractString`: Filename of the HDF5 file.

## Example
```julia
using HDF5
using LegendHDF5IO
using SolidStateDetectors
sim = ssd_read("example_sim.h5", Simulation)
```

!!! note
    In order to use this method, the packages [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) and
    [LegendHDF5IO.jl](https://github.com/legend-exp/LegendHDF5IO.jl) have to be
    loaded before loading SolidStateDetectors.jl.

See also [`ssd_write`](@ref).
"""
function ssd_read end
export ssd_read