text:
# DESCRIPTION
# Converts pilot stick percentages to swashplate inputs.
#
# ------------------------------------------------------------------------------

function control_mixing(u)
    # convert stick percentages to inches
    uin = perc2in*u
    # input to mixer
    umix = LNKGAIN*uin
    # servo inputs
    servos=MIXGAIN*umix
    # servo inputs to blade deflections (MR cyclic, MR collective, TR collective)
    control = zeros(4)
    control[1:3] = SWASHGAIN*servos[1:3]+SWASHBIAS
    # MR cyclic convection: β_1c = -A1, β_1s = -A2
    control[1:2] = -control[1:2]
    # add tail rotor collective
    control[4] = TRGAIN*servos[4]+TRBIAS
    #
    return control
end

ext: jl

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-SimpleHel

max_forks_repo_name: umbertosaetti/J-SimpleHel


text:
maxiter = 1000

using Stopping

stp = NLPStopping(nlp, NLPAtX(nlp.meta.x0)  )
stp.meta.optimality_check = unconstrained_check
stp.meta.max_iter = maxiter
stp.meta.rtol = 0

using OneDmin

function test_Stp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed stp = algo(nlp; kwargs...)
    end

    return stats, stp
end

using Test
@info log_header([:name, :time, :iter, :f, :dual], [String, Float64, Int, Float64, Float64],
                 hdr_override=Dict(:name => "Solver", :f=>"f(x)", :dual=>"‖∇f‖"))

println("\n CG   ")

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_HZ)

@info log_row(Any["CG_HZ", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_HS)

@info log_row(Any["CG_HS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = false, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_PR)

@info log_row(Any["CG_PR", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(CG_generic, nlp, stp=stp, scaling = true, strongWolfe = true, LS_algo=bracket_N, CG_formula=formula_FR)

@info log_row(Any["CG_FR", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

println("\n bfgs   ")

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B)

@info log_row(Any["bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(bfgs_Stop, nlp, stp=stp)

@info log_row(Any["bfgsS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)

function test_noStp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    iter=0
    f=0.0
    normg=0.0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed iter, f, normg, B = algo(nlp; kwargs...)
    end

    return stats, iter, f, normg
end

stats, iter, f, g = test_noStp(bfgs, nlp, scaling = true, maxiter = maxiter, Lp = Inf)

@info log_row(Any["bfgs", stats.time,  iter, f, g])
@test g < 1e-6

println("\n L-bfgs   ")

mem = 70
reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(L_bfgs_StopLS, nlp, stp=stp, LS_algo=bracket_B, mem = mem)

@info log_row(Any["L-bfgsSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(L_bfgs_Stop, nlp, stp=stp, mem = mem)

@info log_row(Any["L-bfgsS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)

function test_noStp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    iter=0
    f=0.0
    normg=0.0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed iter, f, normg, B = algo(nlp; kwargs...)
    end

    return stats, iter, f, normg
end

stats, iter, f, g = test_noStp(L_bfgs, nlp, scaling = true, maxiter = maxiter, Lp = Inf, mem = mem)

@info log_row(Any["L-bfgs", stats.time,  iter, f, g])
@test g < 1e-6

println("\n Newton   ")

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(Newton_StopLS, nlp, stp=stp)

@info log_row(Any["NwtSLS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)
reinit!(stp)

stats, stp = test_Stp(Newton_Stop, nlp, stp=stp)

@info log_row(Any["NwtS", stats.time,  stp.meta.nb_of_stop, stp.current_state.fx, norm(stp.current_state.gx, Inf)])
@test norm(stp.current_state.gx, Inf) < 1e-6

reset!(nlp)

function test_noStp(algo::Function, nlp; kwargs...) where T

    stats = @timed bidon=0
    iter=0
    f=0.0
    normg=0.0
    logger = Logging.ConsoleLogger(stderr,Logging.Warn)
    Logging.with_logger(logger) do
        stats = @timed iter, f, normg = algo(nlp; kwargs...)
    end

    return stats, iter, f, normg
end

stats, iter, f, g = test_noStp(Newton_Spectral, nlp,  maxiter = maxiter, Lp = Inf, ϵ = 1e-6)

@info log_row(Any["Nwt", stats.time,  iter, f, g])
@test g < 1e-6

max_stars_repo_name: tmigot/LSDescent

max_issues_repo_name: tmigot/LSDescent

max_forks_repo_name: tmigot/LSDescent


text:
@doc raw"""
Holds the tableau of a stochastic implicit split partitioned Runge-Kutta method.

`qdrift`, `pdrift1`, `pdrift2` hold the RK coefficients for the drift parts,
and `qdiff`, `pdiff1`, `pdiff2` hold the RK coefficients for the diffusion part of the SDE.

Order of the tableau is not included, because unlike in the deterministic
setting, it depends on the properties of the noise (e.g., the dimension of
the Wiener process and the commutativity properties of the diffusion matrix)

Orders stored in `qdrift` and `qdiff` are understood as the classical orders of these methods.
"""
struct TableauSISPRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift ::CoefficientsRK{T}
    qdiff  ::CoefficientsRK{T}
    pdrift1::CoefficientsRK{T}
    pdrift2::CoefficientsRK{T}
    pdiff1 ::CoefficientsRK{T}
    pdiff2 ::CoefficientsRK{T}

    function TableauSISPRK{T}(name, s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2) where {T}
        @assert s == qdrift.s == qdiff.s == pdrift1.s == pdrift2.s == pdiff1.s == pdiff2.s
        new(name, s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2)
    end
end

function TableauSISPRK(name::Symbol, qdrift::CoefficientsRK{T},
                                     qdiff::CoefficientsRK{T},
                                     pdrift1::CoefficientsRK{T}, pdrift2::CoefficientsRK{T},
                                     pdiff1::CoefficientsRK{T}, pdiff2::CoefficientsRK{T}) where {T}
    TableauSISPRK{T}(name, qdrift.s, qdrift, qdiff, pdrift1, pdrift2, pdiff1, pdiff2)
end

function TableauSISPRK(name::Symbol, qorder_drift::Int, qa_drift::Matrix{T}, qb_drift::Vector{T}, qc_drift::Vector{T},
                                      qorder_diff::Int, qa_diff::Matrix{T}, qb_diff::Vector{T}, qc_diff::Vector{T},
                                      porder1_drift::Int, pa1_drift::Matrix{T}, pb1_drift::Vector{T}, pc1_drift::Vector{T},
                                      porder2_drift::Int, pa2_drift::Matrix{T}, pb2_drift::Vector{T}, pc2_drift::Vector{T},
                                      porder1_diff::Int, pa1_diff::Matrix{T}, pb1_diff::Vector{T}, pc1_diff::Vector{T},
                                      porder2_diff::Int, pa2_diff::Matrix{T}, pb2_diff::Vector{T}, pc2_diff::Vector{T}) where {T}
    TableauSISPRK{T}(name, length(qc_drift), CoefficientsRK(name, qorder_drift, qa_drift, qb_drift, qc_drift),
                                              CoefficientsRK(name, qorder_diff, qa_diff, qb_diff, qc_diff),
                                              CoefficientsRK(name, porder1_drift, pa1_drift, pb1_drift, pc1_drift),
                                              CoefficientsRK(name, porder2_drift, pa2_drift, pb2_drift, pc2_drift),
                                              CoefficientsRK(name, porder1_diff, pa1_diff, pb1_diff, pc1_diff),
                                              CoefficientsRK(name, porder2_diff, pa2_diff, pb2_diff, pc2_diff))
end

# TODO function readTableauSFIRKFromFile(dir::AbstractString, name::AbstractString)

@doc raw"""
Parameters for right-hand side function of implicit Runge-Kutta methods.

`A` - if positive, the upper bound of the Wiener process increments; if `A=0.0`, no truncation.
"""
mutable struct ParametersSISPRK{DT, TT, D, M, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equ::ET
    tab::TableauSISPRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    A::DT

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersSISPRK{DT,D,M}(equ::ET, tab::TableauSISPRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}, A::DT) where {DT, TT, D, M, ET <: NamedTuple}
        @assert M == length(ΔW) == length(ΔZ)
        new{DT, TT, D, M, tab.s, ET}(equ, tab, Δt, ΔW, ΔZ, A, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Structure for holding the internal stages `Q`, the values of the drift vector
and the diffusion matrix evaluated at the internal stages `V=v(Q)`, `B=B(Q)`,
and the increments `Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW`.
"""
struct IntegratorCacheSISPRK{DT,D,M,S} <: SDEIntegratorCache{DT,D,M}
    Q ::Vector{Vector{DT}}
    P ::Vector{Vector{DT}}
    V ::Vector{Vector{DT}}
    F1::Vector{Vector{DT}}
    F2::Vector{Vector{DT}}
    B ::Vector{Matrix{DT}}
    G1::Vector{Matrix{DT}}
    G2::Vector{Matrix{DT}}
    Y ::Vector{Vector{DT}}
    Z ::Vector{Vector{DT}}

    v ::Vector{DT}
    f1::Vector{DT}
    f2::Vector{DT}
    b ::Matrix{DT}
    g1::Matrix{DT}
    g2::Matrix{DT}
    y ::Vector{DT}
    z ::Vector{DT}

    Δy::Vector{DT}
    Δz::Vector{DT}

    function IntegratorCacheSISPRK{DT,D,M,S}() where {DT,D,M,S}

        # create internal stage vectors
        Q  = create_internal_stage_vector(DT, D, S)
        P  = create_internal_stage_vector(DT, D, S)
        V  = create_internal_stage_vector(DT, D, S)
        F1 = create_internal_stage_vector(DT, D, S)
        F2 = create_internal_stage_vector(DT, D, S)
        B  = create_internal_stage_matrix(DT, D, M, S)
        G1 = create_internal_stage_matrix(DT, D, M, S)
        G2 = create_internal_stage_matrix(DT, D, M, S)
        Y  = create_internal_stage_vector(DT, D, S)
        Z  = create_internal_stage_vector(DT, D, S)

        # create velocity and update vector
        v  = zeros(DT,D)
        f1 = zeros(DT,D)
        f2 = zeros(DT,D)
        b  = zeros(DT,D,M)
        g1 = zeros(DT,D,M)
        g2 = zeros(DT,D,M)
        y  = zeros(DT,D)
        z  = zeros(DT,D)

        # create temporary arrays
        Δy  = zeros(DT,M)
        Δz  = zeros(DT,M)

        new(Q, P, V, F1, F2, B, G1, G2, Y, Z, v, f1, f2, b, g1, g2, y, z, Δy, Δz)
    end
end

function Integrators.IntegratorCache{ST}(params::ParametersSISPRK{DT,TT,D,M,S}; kwargs...) where {ST,DT,TT,D,M,S}
    IntegratorCacheSISPRK{ST,D,M,S}(; kwargs...)
end

@inline Integrators.CacheType(ST, params::ParametersSISPRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = IntegratorCacheSISPRK{ST,D,M,S}

"Stochastic implicit partitioned Runge-Kutta integrator."
struct IntegratorSISPRK{DT, TT, D, M, S,
                PT <: ParametersSISPRK{DT,TT},
                ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S}
    params::PT
    solver::ST
    caches::CacheDict{PT}

    function IntegratorSISPRK(params::ParametersSISPRK{DT,TT,D,M,S}, solver::ST, caches) where {DT,TT,D,M,S,ST}
        new{DT, TT, D, M, S, typeof(params), ST}(params, solver, caches)
    end

    function IntegratorSISPRK{DT,D,M}(equations::NamedTuple, tableau::TableauSISPRK{TT}, Δt::TT; K::Int=0) where {DT,TT,D,M}
        # get number of stages
        S = tableau.s

        # K - the integer in the bound A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; K=0 no truncation
        K==0 ? A = 0.0 : A = sqrt( 2*K*Δt*abs(log(Δt)) )

        # create params
        params = ParametersSISPRK{DT,D,M}(equations, tableau, Δt, zeros(DT,M), zeros(DT,M), A)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, 2*D*S, params, caches)

        # create integrator
        IntegratorSISPRK(params, solver, caches)
    end

    function IntegratorSISPRK(equation::SPSDE{DT,TT}, tableau::TableauSISPRK{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorSISPRK{DT, ndims(equation), equation.m}(get_function_tuple(equation), tableau, Δt; kwargs...)
    end
end

@doc raw"""
This function computes initial guesses for `Y`, `Z` and assigns them to `int.solver.x`.

For SISPRK we are NOT IMPLEMENTING an InitialGuess.

SIMPLE SOLUTION
The simplest initial guess for `Y`, `Z` is 0.
"""
function initial_guess!(int::IntegratorSISPRK{DT,TT}, sol::AtomicSolutionPSDE{DT,TT},
                        cache::IntegratorCacheSISPRK{DT}=int.caches[DT]) where {DT,TT}
    int.solver.x .= 0
end

@doc raw"""
Unpacks the data stored in `x = (Y[1][1], Y[1][2], ... Y[1][D], Y[2][1], ..., Z[1][1], Z[1][2], ... Z[1][D], Z[2][1], ...)`
into `Y`, `Z`, calculates the internal stages `Q`, `P`, the values of the RHS
of the SDE ( `vi(Q,P)`, `fi(Q,P)`, `Bi(Q,P)` and `Gi(Q,P)` ), and assigns them
to `V[i]`, `F[i]`, `B[i]` and `G[i]`.
Unlike for FIRK, here
`Y = Δt a_drift v(Q,P) + a_diff B(Q,P) ΔW`,
`Z = Δt â1_drift f1(Q,P) + Δt â2_drift f2(Q,P) + â1_diff G1(Q,P) ΔW + â2_diff G2(Q,P) ΔW`.
"""
function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSISPRK{ST},
                         params::ParametersSISPRK{DT,TT,D,M,S}) where {ST,DT,TT,D,M,S}

    local tqᵢ::TT       #times for the q internal stages
    local tpᵢ₁::TT      #times for the p internal stages
    local tpᵢ₂::TT      #times for the p internal stages

    Q  = cache.Q
    P  = cache.P
    V  = cache.V
    F1 = cache.F1
    F2 = cache.F2
    B  = cache.B
    G1 = cache.G1
    G2 = cache.G2
    Y  = cache.Y
    Z  = cache.Z

    # copy x to Y, Z and calculate Q, P
    for i in 1:S
        for k in 1:D
            Y[i][k] = x[D*(  i-1)+k]
            Z[i][k] = x[D*(S+i-1)+k]
        end
        Q[i] .= params.q .+ Y[i]
        P[i] .= params.p .+ Z[i]
    end

    # compute Vi = v(Qi,Pi), Fi = f(Qi,Pi), Bi=B(Qi,Pi) and Gi=G(Qi,Pi)
    for i in 1:S
        tqᵢ  = params.t + params.Δt * params.tab.qdrift.c[i]
        tpᵢ₁ = params.t + params.Δt * params.tab.pdrift1.c[i]
        tpᵢ₂ = params.t + params.Δt * params.tab.pdrift2.c[i]

        # calculates v(t,Q,P), f1(t,Q,P), f2(t,Q,P) and assigns to V, F1, F2
        params.equ[:v](tqᵢ, Q[i], P[i], V[i])
        params.equ[:f1](tpᵢ₁, Q[i], P[i], F1[i])
        params.equ[:f2](tpᵢ₂, Q[i], P[i], F2[i])

        # calculates B(t,Q,P), G1(t,Q,P), G2(t,Q,P) and assigns to the matrices B[i], G1[i], G2[i]
        params.equ[:B](tqᵢ, Q[i], P[i], B[i])
        params.equ[:G1](tpᵢ₁, Q[i], P[i], G1[i])
        params.equ[:G2](tpᵢ₂, Q[i], P[i], G2[i])
    end
end

"Compute stages of stochastic implicit split partitioned Runge-Kutta methods."
function Integrators.function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersSISPRK{DT,TT,D,M,S},
                caches::CacheDict) where {ST,DT,TT,D,M,S}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(x, cache, params)

    local y1::ST
    local y2::ST
    local z1::ST
    local z2::ST

    # compute b = - (Y-AV)
    for i in 1:S
        for k in 1:D
            y1 = 0
            y2 = 0
            z1 = 0
            z2 = 0
            for j in 1:S
                y1 += params.tab.qdrift.a[i,j]  * cache.V[j][k]  * params.Δt
                y2 += params.tab.qdrift.â[i,j]  * cache.V[j][k]  * params.Δt
                z1 += params.tab.pdrift1.a[i,j] * cache.F1[j][k] * params.Δt
                   +  params.tab.pdrift2.a[i,j] * cache.F2[j][k] * params.Δt
                z2 += params.tab.pdrift1.â[i,j] * cache.F1[j][k] * params.Δt
                   +  params.tab.pdrift2.â[i,j] * cache.F2[j][k] * params.Δt
                for l in 1:M
                    y1 += params.tab.qdiff.a[i,j]  * cache.B[j][k,l]  * params.ΔW[l]
                    y2 += params.tab.qdiff.â[i,j]  * cache.B[j][k,l]  * params.ΔW[l]
                    z1 += params.tab.pdiff1.a[i,j] * cache.G1[j][k,l] * params.ΔW[l]
                       +  params.tab.pdiff2.a[i,j] * cache.G2[j][k,l] * params.ΔW[l]
                    z2 += params.tab.pdiff1.â[i,j] * cache.G1[j][k,l] * params.ΔW[l]
                       +  params.tab.pdiff2.â[i,j] * cache.G2[j][k,l] * params.ΔW[l]
                end
            end
            b[D*(  i-1)+k] = - cache.Y[i][k] + (y1 + y2)
            b[D*(S+i-1)+k] = - cache.Z[i][k] + (z1 + z2)
        end
    end
end

function update_solution!(sol::AtomicSolutionPSDE{T}, tab::TableauSISPRK{T}, Δt::T, ΔW::Vector{T},
                          cache::IntegratorCacheSISPRK{T}) where {T}

    update_solution!(sol, cache.V, cache.F1, cache.F2, cache.B, cache.G1, cache.G2,
                     tab.qdrift.b, tab.qdiff.b, tab.pdrift1.b, tab.pdrift2.b, tab.pdiff1.b, tab.pdiff2.b,
                     Δt, ΔW, cache.Δy, cache.Δz)

    update_solution!(sol, cache.V, cache.F1, cache.F2, cache.B, cache.G1, cache.G2,
                     tab.qdrift.b̂, tab.qdiff.b̂, tab.pdrift1.b̂, tab.pdrift2.b̂, tab.pdiff1.b̂, tab.pdiff2.b̂,
                     Δt, ΔW, cache.Δy, cache.Δz)
end

"Integrate PSDE with a stochastic implicit partitioned Runge-Kutta integrator."
function Integrators.integrate_step!(int::IntegratorSISPRK{DT,TT}, sol::AtomicSolutionPSDE{DT,TT},
                                     cache::IntegratorCacheSISPRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess and assign to int.solver.x
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # compute final update
    update_solution!(sol, tableau(int), timestep(int), int.params.ΔW, cache)
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libchemfiles"], :libchemfiles),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/Chemfiles_jll.jl/releases/download/Chemfiles-v0.10.2+0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.aarch64-linux-gnu-cxx03.tar.gz", "8c3bf277846e8d2c6278388ea6ed754d0b3699669d4ed1ce267bed8d8f2ffc2e"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.aarch64-linux-gnu-cxx11.tar.gz", "cda398a69e0b2f77126833dce519b61eec4c379a203801fcc06506b287372fab"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.aarch64-linux-musl-cxx03.tar.gz", "b5daf195d1b17bcc47a961240af433c48e61288b1d708a6df3142187ae89a795"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.aarch64-linux-musl-cxx11.tar.gz", "f0dfe558cc109c8c7df682630c46e72241c3369678a477722bf06ed064ddd4f3"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.armv7l-linux-gnueabihf-cxx03.tar.gz", "5b8e3246d0d90b7ee3fc3370ef3c10e448842955812f16a55ed798e16f497c24"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.armv7l-linux-gnueabihf-cxx11.tar.gz", "59dc17ffd3f5bac1a7c719baab9db39949c5ced802ea646daa7679052a68d598"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.armv7l-linux-musleabihf-cxx03.tar.gz", "cd601578fe322d7f1e073a7d5cf399ac9b738e188478e5d0962ccb541df535dd"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.armv7l-linux-musleabihf-cxx11.tar.gz", "aa0346081873b1e51c605596c4bba415fe59324a87d9cd6f40de57768404e2e3"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-linux-gnu-cxx03.tar.gz", "e571164d9b9f29bd18d319ccdd30f52f9b82a913d6a4432d4bbd126e84da2f79"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-linux-gnu-cxx11.tar.gz", "cd276c499c83f73f396e49371974a9feb5681c8247bfa9b23f1bf4b3a0fde066"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-linux-musl-cxx03.tar.gz", "1adc817ab41f204d68f7366589d7c5b43943e47489a5c02d4339a95b9396da22"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-linux-musl-cxx11.tar.gz", "5cf1609b12ed33d6aa440e3bec248214cc8d0e0fe73acf3e5572f8fb1fcc1382"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-w64-mingw32-cxx03.tar.gz", "c037b41ec4feb247c14aed582f58225dac7702419c6370777bb15d28ced99434"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.i686-w64-mingw32-cxx11.tar.gz", "dddeffbed57bb14ff5ada1ce9d17dc2e3be1bc84aa9562e49bd2bbead7892e07"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.powerpc64le-linux-gnu-cxx03.tar.gz", "9119d88ab98582949eda926c77ef52294a2066cd93b7b9407584c5ca6e2c46b7"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.powerpc64le-linux-gnu-cxx11.tar.gz", "a58fcfc747aa3ae56841283a385d619356aaa599ecaf71ce93aadb165e720d6c"),
    MacOS(:x86_64) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-apple-darwin.tar.gz", "5cfbcde3c3979f65333ad4c0b894fe9a3922e7aed13e8c10c2fafa3165dff7da"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-linux-gnu-cxx03.tar.gz", "449f6ea2e64e753faca57ef11b1b7a73af357a396eacf51072b36e7a1d8aebd3"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-linux-gnu-cxx11.tar.gz", "cfa42814faad11e26b10c18bdc74f6538b03cfcab1b59188525fc57f013e2820"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-linux-musl-cxx03.tar.gz", "a9d33deea6a101bf39d1b63dccddeec7ccc062c0da6372946cb17f0a2fc2424f"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-linux-musl-cxx11.tar.gz", "77a333de0f4b2a230f1d9944be5c62069a9dc0657e46b2aad0b2e58fee7d72c2"),
    FreeBSD(:x86_64) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-unknown-freebsd.tar.gz", "7540ef7c208d2f883ef90fdcf3c69fc3854e2766cad162602594bb769c0aab9c"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-w64-mingw32-cxx03.tar.gz", "5dc4779eae1351552018cc38462d49e4bc85a8baf8bf1bd8eade6fdab72fb858"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Chemfiles.v0.10.2.x86_64-w64-mingw32-cxx11.tar.gz", "5542b93bc22bac6f2fec222b7f926ce50d9bfe6eb3019903d335208b08fc56ca"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
module SpringCollab2020FlagToggleWater

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/FlagToggleWater" FlagToggleWater(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight, hasBottom::Bool=false, flag::String="flag_toggle_water", inverted::Bool=false)

const placements = Ahorn.PlacementDict(
    "Flag Toggle Water (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagToggleWater,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::FlagToggleWater) = 8, 8
Ahorn.resizable(entity::FlagToggleWater) = true, true

Ahorn.selection(entity::FlagToggleWater) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleWater, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.0, 0.0, 1.0, 0.4), (0.0, 0.0, 1.0, 1.0))
end

end


text:
#
# Factors Indexing
#

Base.getindex(ϕ::Factor, ::Colon) = Factor(ϕ.dimensions, ϕ.potential[:])

Base.getindex(ϕ::Factor, Is::Pair{Symbol}...) = ϕ[Assignment(Is...)]

function Base.getindex(ϕ::Factor, a::Assignment)
    (I, dims) = _translate_index(ϕ, a)
    return Factor(dims, ϕ.potential[I...])
end

Base.setindex!(ϕ::Factor, v, ::Colon) = (ϕ.potential[:] = v)

Base.setindex!(ϕ::Factor, v, Is::Pair{Symbol}...) = setindex!(ϕ, v, Assignment(Is...))

function Base.setindex!(ϕ::Factor, v, a::Assignment)
    (I, _) = _translate_index(ϕ, a)
    return ϕ.potential[I...] = v
end

# assignment to index (as array)
function _translate_index(ϕ::Factor, a::Assignment)
    ind = Array{Any}(ndims(ϕ))
    # all dimensions are accessed by default
    ind[:] = Colon()

    I = Array{Any}(ndims(ϕ))
    I[:] = Colon()
    dims = copy(scope(ϕ))
    for (i, d) in enumerate(scope(ϕ))
        n = name(d)
        if haskey(a, n)
            (ind, new_d) = update(d, a[n])
            I[i] = ind
            dims[i] = new_d
        end
    end

    dims = filter(d -> ~isempty(d), dims)
    return (I, dims)
end

"""
    ind2sub(ϕ::Factor, i::Integer)

Return the subsript into `potential`
"""
Base.ind2sub(ϕ::Factor, i::Integer) = ind2sub(size(ϕ), i)

"""
    sub2ind(ϕ::Factor, i::Integer...)

Return the index into `potential`
"""
Base.sub2ind(ϕ::Factor, i::Integer...) = sub2ind(size(ϕ), i...)

"""
    at2sub(ϕ::Factor, at::Vararg{Any, K}) -> subscript

Convert an assignment tuple to a subscript.
"""
at2sub{N, K}(ϕ::Factor{N}, vs::Vararg{Any, K}) =
    ntuple(Val{min(N, K)}) do k
        d = ϕ.dimensions[k]
        v = vs[k]

        return indexin(v, d)
    end

"""
    sub2at(ϕ::Factor, i::Integer...) -> Tuple

Convert a subscript into an assignment tuple.
"""
sub2at{N, K}(ϕ::Factor{N}, I::Vararg{Integer, K}) =
    ntuple(k -> ϕ.dimensions[k][I[k]], Val{min(N, K)})

"""
    at2a(ϕ::Factor, at::Vararg{Any, K}) -> Assignment

Convert an assignment tuple to an Assignment.
"""
at2a(ϕ::Factor, at::Vararg{Any}) = Dict( s => v for (s, v) in zip(names(ϕ), at) )

"""
    a2at(ϕ::Factor, a::Assignment) -> Tuple
    a2at(ϕ::Factor, a::Pair{Symbol}...) -> Tuple

Convert an assignment to an assignment tuple.
"""
a2at(ϕ::Factor, pairs::Pair{Symbol}...) = a2at(ϕ, Assignment(pairs...))
a2at{N}(ϕ::Factor{N}, a::Assignment) =
    ntuple(Val{N}) do k
        d = ϕ.dimensions[k]
        n = name(d)

        @boundscheck haskey(a, n) || assignment_missing_error(n)
        return a[n]
    end

"""
    a2sub(ϕ::Factor, a::Assignment) -> subscript
    a2sub(ϕ::Factor, a::Pair{Symbol}...) -> subscript

Convert an assignment to a subscript. All variables in scope specified as a single value.
"""
a2sub(ϕ::Factor, pairs::Pair{Symbol}...) = a2sub(ϕ, Assignment(pairs...))
a2sub(ϕ::Factor, a::Assignment) = at2sub(ϕ, a2at(ϕ, a)...)

"""
    sub2a(ϕ::Factor, i::Integer...) -> Assignment

Convert a subscript into an assignment.
"""
sub2a(ϕ::Factor, I::Integer...) = at2a(ϕ, sub2at(ϕ, I...)...)


text:
"""
precision(ArbFloat)           # show the current default precision
setprecision(ArbFloat, 120)   # change the current default precision
setprecision(ArbFloat, 53+7)  # akin to setprecision(BigFloat, 53)

ArbFloat(12)       # use the default precision, at run time
ArbFloat{200}(12)  # use specified precision, at run time
ArbFloat(200,"12") # use the specified precision, at run time
@ArbFloat(12)      # use the default precision, at compile time
@ArbFloat(200,12)  # use specified precision, at compile time

@ArbFloat(1.2345) == ArbFloat("1.2345")

          remember to do this        and           to avoid this

    goodValue = @ArbFloat(1.2345)         wrongValue = ArbFloat(1.2345);
        1.234500000000000000                   1.2344999999999999307

     ArbFloat(12345)/ArbFloat(1000)        ArbFloat(12.345)/ArbFloat(10)
        1.234500000000000000                   1.234500000000000064

```
setprecision(ArbFloat, 80)

exp1 = exp(ArbFloat(1));
stringsmall(exp1),stringcompact(exp1),string(exp1),stringall(exp1)
> ("2.7182818","2.71828182845905","2.71828182845904523536029","2.71828182845904523536029")
showall_pm(exp1)
> 2.718281828459045235360286±3.3216471534462276e-24
bounds(exp1)
> ( 2.71828182845904523536028,  2.718281828459045235360293 )

setprecision(ArbFloat, 116); # the initial default precision
fuzzed_e = tan(atanh(tanh(atan(exp(one(ArbFloat))))))
> 2.718281828459045235360287
showall(fuzzed_e)
> 2.7182818284590452353602874713527
bounds(fuzzed_e)
> ( 2.718281828459045235360287,
    2.718281828459045235360287 )
> they are not really the same ...
lo, hi = bounds(fuzzed_e); showall(lo,hi)
> ( 2.7182818284590452353602874713526543,
    2.7182818284590452353602874713526701 )

smartstring(fuzzed_e)
> "2.7182818284590452353602874713527-"
```
"""
module ArbFloats

export ArbFloat,      # co-matched decimal rounding, n | round(hi,n,10) == round(lo,n,10)
       @ArbFloat,     # converts string form of argument, precision is optional first arg
       simeq, nsime, prec, preceq, succ, succeq, # non-strict total ordering comparisons
       (≃), (≄), (≺), (⪯), (≻), (⪰),           #    matched binary operators
       upperbound, lowerbound, bounds,
       midpoint, radius, midpoint_radius,
       bounding_midpoint, bounding_radius, bounding_midpoint_radius,
       stringsmall, stringcompact, stringmedium, stringlarge, stringall,
       stringsmall_pm, stringcompact_pm, string_pm,
       stringlarge_pm, stringall_pm, string_exact,
       showsmall, showcompact, showlarge, showall, showpretty,
       showsmall_pm, showcompact_pm, show_pm,
       showlarge_pm, showall_pm,
       stringpretty, smartvalue, smartstring, showsmart,
       readable, show_readable, ReadableNumStyle,
       two, three, four, copymidpoint, copyradius, deepcopyradius,
       get_emax, get_emin, bounded, boundedrange,
       fmod, decompose,
       integerpart, decimalpart, fractionalpart, smartmodf,
       isexact, notexact,
       isposinf, isneginf,
       notnan, notinf, notposinf, notneginf, notfinite,
       iszero, notzero, nonzero, isone, notone, notinteger,
       ispositive, notpositive, isnegative, notnegative,
       includes_integer, excludes_integer, includes_zero, excludes_zero,
       includes_positive, excludes_positive, includes_negative, excludes_negative,
       includes_nonpositive,  includes_nonnegative,
       areequal, notequal, approxeq, (≊),
       narrow, overlap, donotoverlap,
       contains, iscontainedby, doesnotcontain, isnotcontainedby,
       absz, absz2, invsqrt, pow, root,
       tanpi, cotpi, logbase, sincos, sincospi, sinhcosh,
       doublefactorial, risingfactorial, rgamma, agm, polylog,
       relative_error, relative_accuracy, midpoint_precision, trimmed,
       PI,SQRTPI,LOG2,LOG10,EXP1,EULER,CATALAN,KHINCHIN,GLAISHER,APERY, # constants
       get_midpoint_digits_shown, get_radius_digits_shown,  # some interface control
       set_midpoint_digits_shown, set_radius_digits_shown,
       isolate_nonnegative_content, isolate_positive_content, # for interval algorithms
       force_nonnegative_content, force_positive_content,
       sort_intervals                                         # uses weak total ordering over intervals

import Base: stdout,
    hash, convert, promote_rule, isa,
    string, show, parse,
    finalizer, decompose, precision, setprecision,
    typemin, typemax, floatmin, floatmax,
    copy, deepcopy,
    size, length,
    zero, one, isinteger,
    ldexp, frexp, modf, eps,
    isequal, isless, (==),(!=),(<),(<=),(>=),(>),
    min, max, minmax,
    typemax, typemin, floatmax, floatmin,
    float, nextfloat, prevfloat,
    isnan, isinf, isfinite, issubnormal,
    signbit, sign, flipsign, copysign, abs, abs2,
    (+),(-),(*),(/),(\),(%),(^), inv, sqrt, hypot,
    trunc, round, ceil, floor,
    fld, cld, div, mod, rem, divrem, fldmod,
    muladd, fma,
    exp, expm1, log, log1p, log2, log10,
    sin, cos, tan, csc, sec, cot, asin, acos, atan,
    sinh, cosh, tanh, csch, sech, coth, asinh, acosh, atanh,
    sinc, factorial,
    in, union, intersect,
    rand, randn, sort,
    BigInt, BigFloat, Rational

import Base.Rounding: rounding_raw, setrounding_raw, rounding, setrounding

using Serialization
using SpecialFunctions

import SpecialFunctions: gamma, lgamma, digamma, zeta

if isdefined(Base, :iszero)
  import Base:iszero
else
  export iszero
end

NotImplemented(info::AbstractString="") = error(string("this is not implemented\n\t",info,"\n"))

include("support/libarb.jl")                 # for precompiled libraries
#include("support/NemoLibs.jl")                 # for precompiled libraries

using ReadableNumbers

include("type/ArbCstructs.jl")
include("type/MagFloat.jl")
include("type/ArfFloat.jl")
include("type/ArbFloat.jl")
include("type/ArbInterval.jl")

include("basics/primitive.jl")
include("basics/IEEEfp.jl")

include("basics/predicates.jl")
include("basics/convert.jl")
include("basics/compare.jl")
include("basics/sort.jl")

include("basics/string.jl")
include("basics/smartstring.jl")
include("basics/show.jl")
include("basics/serialize.jl")

include("math/arith.jl")
include("math/round.jl")
include("math/elementary.jl")
include("math/constants.jl")
include("math/special.jl")

include("math/rounding.jl")
include("math/rand.jl")
include("math/arrayops.jl")

#=
# precision is significand precision, significand_bits(FloatNN) + 1, for the hidden bit
=#
const ArbFloat16 = ArbFloat{ 11}  # read   2 ? 3 or fewer decimal digits to write the same digits ( 16bit Float)
const ArbFloat32 = ArbFloat{ 24}  # read   6 ? 7 or fewer decimal digits to write the same digits ( 32bit Float)
const ArbFloat64  = ArbFloat{ 53}  # read  15 ?15 or fewer decimal digits to write the same digits ( 64bit Float)
const ArbFloat128 = ArbFloat{113}  # read  33 ?34 or fewer decimal digits to write the same digits (128bit Float)
const ArbFloat256 = ArbFloat{237}  # read  71 ?71 or fewer decimal digits to write the same digits (256bit Float)
const ArbFloat512 = ArbFloat{496}  # read 148?149 or fewer decimal digits to write the same digits (512bit Float)

end # ArbFloats


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

name = "NetCDF"
version = v"4.7.4"

# Collection of sources required to build NetCDF
sources = [
    ArchiveSource("https://github.com/Unidata/netcdf-c/archive/v$(version).zip",
                  "170c9c9020f8909811b06e1034d5ea9288b3d5bd90793e3dd27490191faa7566")
]

# HDF5.h in /workspace/artifacts/805ccba77cd286c1afc127d1e45aae324b507973/include
# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/netcdf-c-*

export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"
export LDFLAGS_MAKE="${LDFLAGS}"

if [[ ${target} == *-mingw* ]]; then
    export LIBS="-lhdf5-0 -lhdf5_hl-0 -lcurl-4 -lz"
    # linking fails with: "libtool:   error: can't build x86_64-w64-mingw32 shared library unless -no-undefined is specified"
    # unless -no-undefined is added to LDFLAGS
    LDFLAGS_MAKE="${LDFLAGS} ${LIBS} -no-undefined"

elif [[ "${target}" == *-apple-* ]]; then
    # this file is referenced by hdf.h by not installed
    touch ${includedir}/features.h
fi

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-utilities \
    --enable-shared \
    --disable-static
make LDFLAGS="${LDFLAGS_MAKE}" -j${nproc}
make install
nc-config --all
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# Set equal to the supported platforms in HDF5
platforms = [
    Platform("x86_64", "linux"),
    # HDF5_jll on armv7l should use the same glibc as the root filesystem
    # before it can be used
    # https://github.com/JuliaPackaging/Yggdrasil/pull/1090#discussion_r432683488
    # Platform("armv7l", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="HDF5_jll", version="1.12.0")),
    Dependency("Zlib_jll"),
    Dependency("LibCURL_jll", v"7.71.1"),
    # The following libraries are dependencies of LibCURL_jll which is now a
    # stdlib, but the stdlib doesn't explicitly list its dependencies
    Dependency("LibSSH2_jll", v"1.9.0"),
    Dependency("MbedTLS_jll", v"2.16.8"),
    Dependency("nghttp2_jll", v"1.40.0"),
]

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

max_stars_repo_name: Affie/Yggdrasil

max_issues_repo_name: Affie/Yggdrasil

max_forks_repo_name: Affie/Yggdrasil


text:
import Test

Test.@testset "Test basic functionality" begin
    include("TestVariables.jl")
    include("TestFirstOrder.jl")
    include("TestFirstOrder2.jl")
    include("TestSource.jl")
    include("TestStateSelection.jl")
    include("TestFilterCircuit.jl")
    include("TestFilterCircuit2.jl")
end

Test.@testset "Test units, uncertainties" begin
    include("TestUnits.jl")
    include("TestUncertainties.jl")
    include("TestUnitsAndUncertainties.jl")

    include("TestTwoInertiasAndIdealGear.jl")
    include("TestTwoInertiasAndIdealGearWithUnits.jl")
    include("TestTwoInertiasAndIdealGearWithUnitsAndUncertainties.jl")
    include("TestTwoInertiasAndIdealGearWithMonteCarlo.jl")
    Test.@test_skip include("TestTwoInertiasAndIdealGearWithUnitsAndMonteCarlo.jl")  # MonteCarlo and Unitful do not yet work together
end

Test.@testset "Test model components" begin
    include("TestSingularLRRL.jl")
    include("TestStateSpace.jl")
    include("TestParameter.jl")
    include("TestHeatTransfer.jl")
end

Test.@testset "Test events, etc." begin
    include("TestSimpleStateEvents.jl")
    include("TestSynchronous.jl")
    include("TestInputOutput.jl")
    include("TestPathPlanning.jl")
    include("TestExtraSimulateKeywordArguments.jl")
    Test.@test_skip include("TestBouncingBall.jl")
end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
include("test_vehicles.jl")
include("test_trajdatas.jl")
include("test_scenes.jl")
include("test_scene_records.jl")


text:
using RangeArrays
using Test

# Regular RangeMatrix tests

R = RangeMatrix(UnitRange{Int}[1:10, 11:20, 21:30, 31:40])

@test vec(collect(R)) == collect(1:40)
@test R[1,1] == 1

for i=1:40
    @test R[i] == i
end
for (i,j) in zip(eachindex(R), 1:length(R))
    @test R[i] == j
end

@test_throws BoundsError R[0]
@test_throws BoundsError R[41]

@test_throws BoundsError R[:, 0]
@test R[:, 1] == 1:10
@test R[:, 2] == 11:20
@test R[:, 3] == 21:30
@test R[:, 4] == 31:40
@test_throws BoundsError R[:, 5]

@test R[:, 1:2] == RangeMatrix(1:10, 11:20)
@test R[:, 1:4] == R[:,:] == RangeMatrix(1:10, 11:20, 21:30, 31:40)
@test R[:, 3:4] == RangeMatrix(21:30, 31:40)

# RepeatedRangeMatrix tests
R = RepeatedRangeMatrix(1:10, 0:10:30)

@test vec(collect(R)) == collect(1:40)
@test R[1,1] == 1

for i=1:40
    @test R[i] == i
end
for (i,j) in zip(eachindex(R), 1:length(R))
    @test R[i] == j
end

@test_throws BoundsError R[0]
@test_throws BoundsError R[41]

@test_throws BoundsError R[:, 0]
@test R[:, 1] == 1:10
@test R[:, 2] == 11:20
@test R[:, 3] == 21:30
@test R[:, 4] == 31:40
@test_throws BoundsError R[:, 5]

@test R[:, 1:2] == RangeMatrix(1:10, 11:20)
@test R[:, 1:4] == R[:,:] == RangeMatrix(1:10, 11:20, 21:30, 31:40)
@test R[:, 3:4] == RangeMatrix(21:30, 31:40)

A = 1:100
@test_throws BoundsError A[RangeMatrix(0:9, 20:29)]
@test_throws BoundsError A[RangeMatrix(20:29, 0:9)]
@test_throws BoundsError A[RangeMatrix(20:29, 92:101)]
@test_throws BoundsError A[RangeMatrix(92:101, 20:29)]
@test_throws BoundsError A[RangeMatrix(-100:100)]
@test A[RangeMatrix(1:10, 91:100)] == [1:10 91:100]

@test_throws BoundsError A[RepeatedRangeMatrix(1:10, [-1,33])]
@test_throws BoundsError A[RepeatedRangeMatrix(1:10, [33,-1])]
@test_throws BoundsError A[RepeatedRangeMatrix(1:10, [91,33])]
@test_throws BoundsError A[RepeatedRangeMatrix(1:10, [33,91])]
@test_throws BoundsError A[RepeatedRangeMatrix(-100:100, [0])]
@test A[RepeatedRangeMatrix(1:10, [0,90])] == [1:10 91:100]


text:
push!(LOAD_PATH, "/Users/ortner/gits/Atoms.jl")

using JuLIP, TightBinding
import AtJuLIP
TB = TightBinding

nkpoints = (0,0,0)
beta = 10.0

# compute a configuration
at = bulk("Si", cubic = true) * 2
rattle!(at, 0.02)

# JuLIP TB model
tbj = TB.NRLTB.NRLTBModel(:Si, FermiDiracSmearing(beta), orbitals = :sp, cutoff=:energyshift)

# JulIP matrices
Hj, Mj = hamiltonian(tbj, at)
Hj = full(Hj)
Mj = full(Mj)

# Atoms.jl TB Model
tba = AtJuLIP.JuLIPTB(:Si, nkpoints = nkpoints)

# Atoms Matrices
Ha, Ma = hamiltonian(tba, at)
Ha = TB.NRLTB.half_eV * full(Ha)
Ma = full(Ma)

@show vecnorm(Hj - Ha, Inf)
@test vecnorm(Hj - Ha, Inf) < 1e-12
@show vecnorm(Mj - Ma, Inf)
@test vecnorm(Mj - Ma, Inf) < 1e-12
ae = full(Ha) |> eigvals |> real |> sort
je = full(Hj) |> eigvals |> real |> sort
@show vecnorm(ae - je, Inf)


text:
module DataTools

export averaging,
    firstitem,
    firstitems,
    inc1,
    lastitem,
    lastitems,
    meanvar,
    modifying,
    nitems,
    oncol,
    rightif

using Base: HasLength, HasShape, IteratorSize
using InitialValues: InitialValues
using Accessors: @optic, PropertyLens, modify, set
using StaticNumbers: static
using Statistics: Statistics, mean, std, var
using Tables: Tables
using Transducers:
    Composition,
    Count,
    IdentityTransducer,
    Map,
    MapSplat,
    Scan,
    Take,
    TakeLast,
    Transducers,
    combine,
    complete,
    extract_transducer,
    next,
    opcompose,
    reducingfunction,
    right,
    start

include("utils.jl")
include("oncol.jl")
include("modifying.jl")
include("reductions.jl")
include("reducers.jl")

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end DataTools

end


text:
using Revise
using Bem2d

"""
    constvsquadelements()

Compare farfield displacements and stresses from constant displacement
and quadratic displacement elements.  The point of this is to highlight
the fact the quadratic elements can be used to create linear slip elements
that can accurately model a linear displacement gradient with no spurious
singularities.
"""
function constvsquadelements()
    mu = 3e10
    nu = 0.25

    # Create a flat fault
    els = Elements(Int(1e5))
    nels = 20
    x1, y1, x2, y2 = discretizedline(-10e3, 0, 10e3, 0, nels)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Convenience dictionary for element names
    idx = getidxdict(els)

    # Set fault slip
    xcenters = els.xcenter[1:els.endidx]
    ycenters = els.ycenter[1:els.endidx]
    xnodes = sort(els.xnodes[1:els.endidx, :][:])
    ynodes = sort(els.ynodes[1:els.endidx, :][:])

    # # Constant x-slip only
    # constxslip = ones(nels)
    # constyslip = zeros(nels)
    # quadxslip = repeat(constxslip, 1, 3)
    # quadyslip = repeat(constyslip, 1, 3)

    # Linear x-slip only
    slope = 0.001
    constxslip = slope .* xcenters
    constyslip = zeros(nels)
    quadxslip = transpose(reshape(slope .* xnodes, 3, nels))
    quadyslip = zeros(size(quadxslip))

    # Observation coordinates for far-field calculation
    npts = 100
    obswidth = 3e3
    xobs, yobs = obsgrid(-obswidth, -obswidth, obswidth, obswidth, npts)

    # Displacements and stresses
    dispconst, stressconst = constdispstress(slip2dispstress, xobs, yobs, els, idx["fault"], constxslip, constyslip, mu, nu)
    dispquad, stressquad = quaddispstress(slip2dispstress, xobs, yobs, els, idx["fault"], quadxslip, quadyslip, mu, nu)

    # Plot
    close("all")
    plotfields(els, reshape(xobs, npts, npts), reshape(yobs, npts, npts),
               dispconst, stressconst, "constant slip - constant slip element")
    plotfields(els, reshape(xobs, npts, npts), reshape(yobs, npts, npts),
               dispquad, stressquad, "constant slip - quadratic slip element")
    plotfields(els, reshape(xobs, npts, npts), reshape(yobs, npts, npts),
               dispconst - dispquad, stressconst - stressquad, "residuals")
    show()

    return nothing
end
constvsquadelements()


text:
function gpu_energy!(
    pos::AbstractArray,
    forces::AbstractArray,
    N::Integer,
    L::Real,
    rc::Real,
    a::Real,
    b::Real,
    λ::Integer,
    temp::Real,
    full_ener::AbstractArray,
    vir::AbstractArray
)
    total_energy = 0.0f0
    virial = 0.0f0
    force = 0.0f0
    ener = 0.0f0

    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    @inbounds for i = index:stride:N
        # Reset variables
        total_energy = 0.0f0
        virial = 0.0f0

        for j = 1:N

            if i == j
                continue
            end

            xij = pos[1, i] - pos[1, j]
            yij = pos[2, i] - pos[2, j]
            zij = pos[3, i] - pos[3, j]

            # Periodic boundaries
            xij -= L * round(xij / L)
            yij -= L * round(yij / L)
            zij -= L * round(zij / L)

            # Compute distance
            Δpos = xij * xij + yij * yij + zij * zij
            Δpos = CUDA.sqrt(Δpos)

            if Δpos < rc
                if Δpos < b
                    # * Energy computation
                    ener = CUDA.pow(1.0f0 / Δpos, λ) - CUDA.pow(1.0f0 / Δpos, λ - 1.0f0)
                    ener *= a / temp
                    ener += 1.0f0 / temp

                    # * Force computation
                    force = λ * CUDA.pow(1.0f0 / Δpos, λ + 1.0f0)
                    force -= (λ - 1.0f0) * CUDA.pow(1.0f0 / Δpos, λ)
                    force *= a / temp
                else
                    force = 0.0f0
                    ener = 0.0f0
                end

                # * Update the energy
                total_energy += ener

                # * Compute forces
                fx = force * xij / Δpos
                fy = force * yij / Δpos
                fz = force * zij / Δpos

                # * Update the forces
                forces[1, i] += fx
                forces[2, i] += fy
                forces[3, i] += fz

                forces[1, j] -= fx
                forces[2, j] -= fy
                forces[3, j] -= fz

                # * Compute the virial, avoid double counting
                virial += 0.5f0 * (fx * xij + fy * yij + fz * zij)
            end
        end

        # * Save the values in their respective arrays
        full_ener[i] = total_energy
        vir[i] = virial
    end

    return nothing
end

function gpu_ermak!(
    positions::AbstractArray,
    forces::AbstractArray,
    N::Integer,
    L::Real,
    τ::Real,
    rnd_matrix::AbstractArray,
    pbc::Bool
)

    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x

    for j = index:stride:N
        positions[1, j] += (forces[1, j] * τ) + rnd_matrix[1, j]
        positions[2, j] += (forces[2, j] * τ) + rnd_matrix[2, j]
        positions[3, j] += (forces[3, j] * τ) + rnd_matrix[3, j]

        if pbc
            positions[1, j] -= L * round(positions[1, j] / L)
            positions[2, j] -= L * round(positions[2, j] / L)
            positions[3, j] -= L * round(positions[3, j] / L)
        end
    end

    return nothing
end


text:
module test_45_Save_svg_figure

import ModiaResult
ModiaResult.@usingModiaPlot

# Test the Figure operations
println("\n... test test_45_Save_svg_figure.jl:\n")

t = range(0.0, stop=10.0, length=100)
result = Dict{String,Any}("time" => t, "phi" => sin.(t))
ModiaResult.printResultInfo(result)

plot(result, "phi", figure=2)
saveFigure(2, "test_saveFigure.svg")

end


text:
function varargout = sharperatio(varargin)
# sharperatio MATLAB code for sharperatio.fig
#      sharperatio, by itself, creates a new sharperatio or raises the existing
#      singleton*.
#
#      H = sharperatio returns the handle to a new sharperatio or the handle to
#      the existing singleton*.
#
#      sharperatio('CALLBACK',hObject,eventData,handles,...) calls the local
#      function named CALLBACK in sharperatio.M with the given input arguments.
#
#      sharperatio('Property','Value',...) creates a new sharperatio or raises the
#      existing singleton*.  Starting from the left, property value pairs are
#      applied to the GUI before sharperatio_OpeningFcn gets called.  An
#      unrecognized property name or invalid value makes property application
#      stop.  All inputs are passed to sharperatio_OpeningFcn via varargin.
#
#      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
#      instance to run (singleton)".
#
# See also: GUIDE, GUIDATA, GUIHANDLES

# Edit the above text to modify the response to help sharperatio

# Last Modified by GUIDE v2.5 13-Nov-2013 10:44:25

# Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @sharperatio_OpeningFcn, ...
                   'gui_OutputFcn',  @sharperatio_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
# End initialization code - DO NOT EDIT

# --- Executes just before sharperatio is made visible.
function sharperatio_OpeningFcn(hObject, ~, handles, varargin)
# This function has no output args, see OutputFcn.
# hObject    handle to figure
# eventdata  reserved - to be defined in a future version of MATLAB
# handles    structure with handles and user data (see GUIDATA)
# varargin   command line arguments to sharperatio (see VARARGIN)

# Choose default command line output for sharperatio
handles.output = hObject;

# Update handles structure
guidata(hObject, handles);

# UIWAIT makes sharperatio wait for user response (see UIRESUME)
# uiwait(handles.figure1);

#File: sharperatio.jl
#Author:
#Date: 29-June-2020

#import XL data
[names,~,returns] = getXLData('hfma_matlab_data.xlsx','10 Hedge Funds 2');
#optimisation parameters
N = size(returns,2); #size of portfolio (N)
R = mean(returns)'*252; #annualised mean (R)
W0 = repmat(1/N,N,1); #initial weights (W0)
rf = 2.0; #risk-free rate

covar = cov(returns); #covariance matrix
#annualise variance for monthly returns
for i=1:N
    for j=1:N
        if i == j
            covar(i,j)=covar(i,j)*252;
        end
    end
end
#VCV matrix
VCV = covar;

#portfolio statistics
pR = W0'*R; #return
pVar  = W0'*VCV*W0; #variance
pSharpe = (pR - rf)./sqrt(pVar);
pStats = [pR;pVar;pSharpe];

#linear contstraints and bounds
A = [ones(1,N);-ones(1,N)];
b = [1;-1];
lb = zeros(N,1);
ub = ones(N,1)*0.5;

#maximise sharperatio ratio optimisation
options = optimset('LargeScale', 'off', 'Algorithm','sqp');
[W,fval,ExitFlag] = fmincon(@(W) maxSharpe(W,R,rf,VCV),W0,A,b,[],[],lb,ub,[],options);
disp('Maximise Sharpe ratio');
    if (ExitFlag ~= 1)
        disp('Optimization did not converge!');
    else
        disp('Optimization converged');
    end

#optimised results
oStats = -fval;

#OUTPUT
#format data
R = num2cell(R);
    for i = 1:numel(R)
        R{i} = sprintf('#12.3f',R{i});
    end
#set table data
set(handles.table_mean_returns,'ColumnName','Monthly|Returns|(#)','RowName',names,'Data',R);
#format data
W0 = num2cell(W0);
    for i = 1:numel(W0)
        W0{i} = sprintf('#12.2f',W0{i});
    end
set(handles.table_initial_weights,'ColumnName','Initial|Weights|(W0)','RowName',names,'Data',W0);

#format data
pStats = num2cell(pStats);
    for i = 1:numel(pStats)
        pStats{i} = sprintf('#12.3f',pStats{i});
    end
#set table data
set(handles.table_portfolio_statistics,'ColumnName','Annual|(#)','RowName',{'Port. Mean','Port. Variance','Port. Sharpe'},'Data',pStats);
#format data
VCV = num2cell(VCV);
    for i = 1:numel(VCV)
        VCV{i} = sprintf('#8.3f',VCV{i});
    end
#set size column widths
cwidth = num2cell(repmat(50,1,10));
#set table output
set(handles.table_VCV_matrix,'ColumnName',names,'ColumnWidth',cwidth,'RowName',names,'Data',VCV);

#format data
W = num2cell(W);
    for i = 1:numel(W)
        W{i} = sprintf('#12.3f',W{i});
    end
set(handles.table_optimised_weights,'ColumnName','Optimised|Weights|(W)','RowName',names,'Data',W);

#format data
oStats = num2cell(oStats);
    for i = 1:numel(oStats)
        oStats{i} = sprintf('#12.3f',oStats{i});
    end
#set table data
set(handles.table_optimised_results,'ColumnName','Maximised|(#)','RowName','Port. Sharpe','Data',oStats);

# --- Outputs from this function are returned to the command line.
function varargout = sharperatio_OutputFcn(~, ~, handles)
# varargout  cell array for returning output args (see VARARGOUT);
# hObject    handle to figure
# eventdata  reserved - to be defined in a future version of MATLAB
# handles    structure with handles and user data (see GUIDATA)

# Get default command line output from handles structure
varargout{1} = handles.output;


text:
# Jets/WaveFD bridge

# 2D source illumination
function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, nx::Int, ntrec::Int, nthreads::Int) where {T}
    srcillum!(zeros(T, nz, nx), filename, compressor, ginsu, ntrec, nthreads)
end

function srcillum!(γ::AbstractArray{T,2}, filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T}
    io = open(filename)
    field_ginsu = zeros(T, size(ginsu, interior=interior))
    for it = 1:ntrec
        WaveFD.compressedread!(io, compressor, it, field_ginsu)
        WaveFD.srcillum_helper(field_ginsu, nthreads)
        super!(γ, ginsu, field_ginsu, accumulate=true, interior=interior)
    end
    close(io)
    γ
end

# 3D source illumination
function srcillum(filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ::Type{T}, ginsu::Ginsu, nz::Int, ny::Int, nx::Int, ntrec::Int, nthreads::Int) where {T}
    srcillum!(zeros(T, nz, ny, nx), filename, compressor, interior, ginsu, ntrec, nthreads)
end

function srcillum!(γ::AbstractArray{T,3}, filename::AbstractString, compressor::WaveFD.Compressor, interior::Bool, ginsu::Ginsu, ntrec::Int, nthreads::Int) where {T}
    io = open(filename)
    field_ginsu = zeros(T, size(ginsu, interior=interior))
    for it = 1:ntrec
        WaveFD.compressedread!(io, compressor, it, field_ginsu)
        WaveFD.srcillum_helper(field_ginsu, nthreads)
        super!(γ, ginsu, field_ginsu, accumulate=true, interior=interior)
    end
    close(io)
    γ
end

@doc """
    srcillum(J)
    srcillum(J, m)
    srcillum!(y, J, m)

Compute and return the source illumination for `Jets` operator `J`. The source illumination
is defined as the sum of squares of the wavefield amplitudes everywhere in the model over
the time of the finite difference evolution.

If `J::Jop` is a `Jets` block operator or distributed block operator, the source
illuminations from all blocks will be accumulated with a simple sum.

`srcillum(J)` creates and returns an array with the source illumination from `J::JopLn`,
using the current location defined in the `Jet`.

`srcillum(J, m)` creates and returns an array with the source illumination from `J::Jop`
for the location `m`.

`srcillum!(y, J, m)` zeros the passed array `y` and then accumulates to the source
illumination from `J::Jop` at the location `m` into `y`.
"""
function srcillum(J::JopLn)
    s = zeros(eltype(J), size(domain(J))[1:end-1])
    srcillum!(s, J)
end

function srcillum(J::Jop, m::AbstractArray)
    s = zeros(eltype(J), size(domain(J))[1:end-1])
    srcillum!(s, J, m)
end

@doc (@doc srcillum(J))
srcillum!(γ::AbstractArray, J::Jop) = srcillum!(γ, J, jet(J).mₒ)

# composite operators, TODO: the try-catch here is a bit of a kludge
function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(Jets.JetComposite_f!)},T<:Jop{J}}
    for Aᵢ in state(A).ops
        try
            srcillum!(γ, Aᵢ, m)
            break
        catch e
            if length(m) == 0
                rethrow(e)
            end
        end
    end
    γ
end

# block operators
function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(Jets.JetBlock_f!)}, T<:Jop{J}}
    γ .= 0
    for Aᵢ in state(A).ops
        if length(m) == 0
            srcillum!(γ, Aᵢ, jet(Aᵢ).mₒ)
        else
            srcillum!(γ, Aᵢ, m)
        end
    end
    γ
end

# distributed block operators
function srcillum!(γ::AbstractArray, A::T, m::AbstractArray) where {D,R,J<:Jet{D,R,typeof(DistributedJets.JetDBlock_f!)},T<:Jop{J}}
    γ .= 0
    pids = procs(A)
    γ_partials = ArrayFutures(γ, DistributedJets.addmasterpid(pids))

    _srcillum!(γ_partials, A, m) = begin srcillum!(localpart(γ_partials), localpart(A), m); nothing end
    @sync for pid in pids
        @async remotecall_fetch(_srcillum!, pid, γ_partials, A, m)
    end
    reduce!(γ_partials)
end

export srcillum, srcillum!

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
include("solver_preamble.jl")
import CSDP
factory = optimizer_with_attributes(CSDP.Optimizer, "printlevel" => 0)
config = MOI.Test.Config(atol=1e-4, rtol=1e-4)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
function calculation()
    result = 0.0
    n = 1000_000
    for i in 1:n
        result += sin(rand())
    end
    return result
end

show(calculation())

@time include("1_5_func_better.jl")
459310.8032323998  0.026413 seconds (2.48 k allocations: 124.761 KB)

max_stars_repo_path: Section 1/1_5_func_better.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_5_func_better.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_5_func_better.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
abstract type AbstractPolytopalMesh end
abstract type AbstractCell{dim,V,F} end

"""
A `FaceIndex` wraps an (Int, Int) and defines a face by pointing to a (cell, face).
"""
struct FaceIndex
    cellidx::Int
    idx::Int
end

"""
A `EdgeIndex` wraps an (Int, Int) and defines an edge by pointing to a (cell, edge).
"""
struct EdgeIndex
    cellidx::Int
    idx::Int
end

# Vertices
struct Node{dim,T}
    x::Tensors.Vec{dim, T}
end

@inline get_coordinates(node::Node) = node.x

#--------------- Cells
struct Cell{dim, N, M, L}
    nodes::NTuple{N,Int}
end

#Common cell types
const LineCell = Cell{1,2,1,0}
@inline get_cell_name(::LineCell) = "Line"
@inline reference_edge_nodes(::Type{LineCell}) = ((1,2),)

const TriangleCell = Cell{2,3,3,1}
@inline get_cell_name(::TriangleCell) = "Triangle"
@inline reference_edge_nodes(::Type{TriangleCell}) = ((2,3),(3,1),(1,2))

const RectangleCell = Cell{2,4,4,1}
@inline get_cell_name(::RectangleCell) = "Rectangle"
@inline reference_edge_nodes(::Type{RectangleCell}) = ((1,2),(2,3),(3,4),(4,1))

# API
@inline getnnodes(cell::Cell{dim,N}) where {dim,N} = N
@inline getnedges(cell::Cell{dim,N,M}) where {dim,N,M} = M
@inline getnfaces(cell::Cell{dim,N,M,P}) where {dim,N,M,P} = P

# ----------------- Mesh
struct PolytopalMesh{dim,T,C} <: AbstractPolytopalMesh
    cells::Vector{C}
    nodes::Vector{Node{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{FaceIndex}}
    edgesets::Dict{String,Set{EdgeIndex}}
    nodesets::Dict{String,Set{Int}}
end

function PolytopalMesh(cells,
              nodes;
              cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
              facesets::Dict{String,Set{FaceIndex}}=Dict{String,Set{FaceIndex}}(),
              edgesets::Dict{String,Set{EdgeIndex}}=Dict{String,Set{EdgeIndex}}(),
              nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}()) where {dim}
    return PolytopalMesh(cells, nodes, cellsets, facesets, edgesets, nodesets)
end

# API
getdimension(mesh::PolytopalMesh{dim}) where {dim} = dim
@inline getncells(mesh::PolytopalMesh) = length(mesh.cells)
@inline getnnodes(mesh::PolytopalMesh) = length(mesh.nodes)
@inline getnodesidx(mesh, cell_idx) = mesh.cells[cell_idx].nodes
@inline getnodeset(mesh::PolytopalMesh, set::String) = mesh.nodesets[set]

"""
function getcoords(mesh, node_idx::Int)
Return a Tensor.Vec with the coordinates of node with index `node_idx`
"""
@inline getnodecoords(mesh::PolytopalMesh, node_idx::Int) = mesh.nodes[node_idx].x
@inline getnodecoords(mesh::PolytopalMesh{1}, node_idx::Int) = mesh.nodes[node_idx].x[]

"""
    getnodescoords(mesh::PolytopalMesh, cell_idx)
Return a vector with the coordinates of the nodes of cell number `cell`.
"""
@inline function getnodescoords(mesh::PolytopalMesh{dim,T}, cell_idx::Int) where {dim,T}
    N = getnnodes(mesh.cells[cell_idx])
    coords = Vector{Tensors.Vec{dim,T}}(undef, N)
    for (i,j) in enumerate(mesh.cells[cell_idx].nodes)
        coords[i] = mesh.nodes[j].x
    end
    return coords
end

function get_nodes_matrix(mesh::PolytopalMesh{dim,T,C}) where {dim,T,C}
    nodes_m = Matrix{T}(undef,length(mesh.nodes),dim)
    for (k,node) in enumerate(mesh.nodes)
        nodes_m[k,:] = node.x
    end
    nodes_m
end
function get_conectivity_list(mesh::PolytopalMesh{dim,T,C}) where {dim,T,C}
    cells_m = Vector()
    for k = 1:getncells(mesh)
        push!(cells_m,mesh.cells[k].nodes)
    end
    cells_m
end

function cell_volume(mesh::PolytopalMesh{2}, cell_idx::Int)
    N = getnnodes(mesh.cells[cell_idx])
    verts = getnodescoords(mesh,cell_idx)
    return 0.5*abs(sum(verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2] for j ∈ 1:N))
end

function cell_volume(mesh::PolytopalMesh{1}, cell_idx::Int)
    verts = getnodescoords(mesh,cell_idx)
    return abs(verts[2][1] - verts[1][1])
end

function cell_centroid(mesh::PolytopalMesh{2}, cell_idx::Int)
    verts = getnodescoords(mesh,cell_idx)
    Ve = cell_volume(mesh, cell_idx)
    N = getnnodes(mesh.cells[cell_idx])
    xc = 1/(6*Ve)*sum((verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2])*(verts[j][1]+verts[mod1(j+1,N)][1]) for j ∈ 1:N)
    yc = 1/(6*Ve)*sum((verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2])*(verts[j][2]+verts[mod1(j+1,N)][2]) for j ∈ 1:N)
    return Tensors.Vec{2}((xc,yc))
end

function cell_diameter(mesh::PolytopalMesh{dim,T}, cell_idx::Int) where {dim,T}
    K = mesh.cells[cell_idx]
    verts = getnodescoords(mesh,cell_idx)
    h = zero(T)
     for k in reference_edge_nodes(typeof(K))
        mσ = norm(verts[k[2]] - verts[k[1]])
        h = max(h, mσ)
    end
    h
end

function cell_centers(mesh::PolytopalMesh{1})
    xn = get_nodes_matrix(mesh)
    return [xn[i] + 0.5*(xn[i+1]-xn[i]) for i in 1:getncells(mesh)]
end


text:
@testset "Fundamentals" begin
    symbol = "IBM"

    @testset "Overview" begin
        data = company_overview(symbol)
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 59
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Income Statement" begin
        data = income_statement(symbol)
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Balance Sheet" begin
        data = balance_sheet(symbol)
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Cash Flow" begin
        data = cash_flow(symbol)
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end
end


text:
module SimplePlanarGraphs
using SimpleGraphs, LinearAlgebra

export quick_planar_check

"""
`quick_planar_check(G::SimpleGraph)` performs the most basic check for a
graph to be planar (or not). Specifically, we return true if the number of
vertices `n` is less than 5 or if the number of edges is at most `3n-6`.
"""
quick_planar_check(G::SimpleGraph)::Bool = NV(G)<5 || NE(G)<=3*NV(G)-6

include("count_cross.jl")
include("face_finder.jl")
include("random_triangulation.jl")

end


text:
export FileInfo

const FileInfo = Dict{String,String}

function pause()
  a = readline()
  (length(a) == 0) && return
  (a[1]=='q') && exit(0)
end

function get_arg(args,
                 key,
                 default=nothing)
  if haskey(args,key)
    return args[key]
  end
  if isnothing(default)
    throw(ArgumentError("key $key not found in argument dictionary"))
  end
  return default
end


text:
using Plots
using CSV

function general_setup()
    gr()
    fntsm = Plots.font(pointsize = 12)
    fntlg = Plots.font(pointsize = 18)
    default(titlefont = fntlg, guidefont = fntlg, tickfont = fntsm, legendfont = fntsm)
end

function plot_agent_losses(csv_path, plot_file; lower_bound = 1.0e-15, upper_bound = 1.0)
    general_setup()
    losses = DataFrame(CSV.File(csv_path))
    plot()
    for col in names(losses)
        plot!(
            losses[!, col],
            yaxis = (:log10, [lower_bound, upper_bound]),
            label = col,
            line = (2, :solid),
            legend = :bottomleft,
        )
    end
    xaxis!("Iteration count")
    yaxis!("Objective gap")
    savefig(plot_file)
end

function plot_step_sizes(csv_path, plot_file; lower_bound = 1.0, upper_bound = 1.0e15)
    general_setup()
    losses = DataFrame(CSV.File(csv_path))
    # lower_bound = minimum([s for s in losses[!, "step_size"] if s > 0.0]) / 2
    plot(
        losses[!, "step_size"],
        line = (3, :solid),
        yaxis = (:log10, [lower_bound, upper_bound]),
        label = false,
    )
    xaxis!("Iteration count")
    yaxis!("Ideal step size")
    savefig(plot_file)
end

function plot_objective(csv_path, plot_file; lower_bound = 1.0e-15, upper_bound = 1.0)
    general_setup()
    losses = DataFrame(CSV.File(csv_path))
    # lower_bound = minimum([s for s in losses[!, "objective"] if s > 0.0]) / 2
    plot(
        losses[!, "objective"],
        line = (3, :solid),
        yaxis = (:log10, [lower_bound, lower_bound]),
        color = 2,
        label = false,
    )
    xaxis!("Iteration")
    yaxis!("Objective gap")
    savefig(plot_file)
end


text:
module PoissonTests

using Test
using Gridap
import Gridap: ∇
#using LinearAlgebra

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)
order = 2

const h = (domain[2]-domain[1]) / partition[1]
const γ = 10

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[7,8])
add_tag_from_tags!(labels,"nitsche",6)

Ω = Triangulation(model)
Γn = BoundaryTriangulation(model,labels,tags="neumann")
Γd = BoundaryTriangulation(model,labels,tags="nitsche")

degree = order
dΩ = Measure(Ω,degree)
dΓn = Measure(Γn,degree)
dΓd = Measure(Γd,degree)

nn = get_normal_vector(Γn)
nd = get_normal_vector(Γd)

# Using automatic differentiation
u_scal(x) = x[1]^2 + x[2]
f_scal(x) = - Δ(u_scal)(x)

#u_scal(x) = x[1]^2 + x[2]
#∇u_scal(x) = VectorValue( 2*x[1], one(x[2]) )
#Δu_scal(x) = 2
#f_scal(x) = - Δu_scal(x)
#∇(::typeof(u_scal)) = ∇u_scal

scalar_data = Dict{Symbol,Any}()
scalar_data[:valuetype] = Float64
scalar_data[:u] = u_scal
scalar_data[:f] = f_scal

# Using automatic differentiation
u_vec(x) = VectorValue( x[1]^2 + x[2], 4*x[1] - x[2]^2 )
f_vec(x) = - Δ(u_vec)(x)

#u_vec(x) = VectorValue( x[1]^2 + x[2], 4*x[1] - x[2]^2 )
#∇u_vec(x) = TensorValue( 2*x[1], one(x[2]), 4*one(x[1]), - 2*x[2] )
#Δu_vec(x) = VectorValue( 2, -2 )
#f_vec(x) = - Δu_vec(x)
#∇(::typeof(u_vec)) = ∇u_vec

vector_data = Dict{Symbol,Any}()
vector_data[:valuetype] = VectorValue{2,Float64}
vector_data[:u] = u_vec
vector_data[:f] = f_vec

for data in [ vector_data, scalar_data ]

  T = data[:valuetype]
  u = data[:u]
  f = data[:f]

  for domain_style in (ReferenceDomain(),PhysicalDomain())

    cell_fe = FiniteElements(domain_style,model,lagrangian,T,order)
    V = TestFESpace(Ω,cell_fe,dirichlet_tags="dirichlet",labels=labels)
    U = TrialFESpace(V,u)

    uh = interpolate(u, U)

    a(u,v) =
      ∫( ∇(v)⊙∇(u) )*dΩ +
      ∫( (γ/h)*v⊙u  - v⊙(nd⋅∇(u)) - (nd⋅∇(v))⊙u )*dΓd

    l(v) =
      ∫( v⊙f )*dΩ +
      ∫( v⊙(nn⋅∇(uh)) )*dΓn +
      ∫( (γ/h)*v⊙uh - (nd⋅∇(v))⊙u )*dΓd

    op = AffineFEOperator(a,l,U,V)
    uh = solve(op)

    e = u - uh

    l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
    h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

    el2 = l2(e)
    eh1 = h1(e)
    ul2 = l2(uh)
    uh1 = h1(uh)

    @test el2/ul2 < 1.e-8
    @test eh1/uh1 < 1.e-7

  end
end

end # module


text:
# Here are tests which test what interface the solvers require.

################################################################################
# This is to test a scalar-like state variable
# (due to @acroy: https://gist.github.com/acroy/28be4f2384d01f38e577)

import Base: +, -, *, /, .+, .-, .*, ./

const delta0 = 0.
const V0 = 1.
const g0 = 0.

# define custom type ...
immutable CompSol
  rho::Matrix{Complex128}
  x::Float64
  p::Float64

  CompSol(r,x,p) = new(copy(r),x,p)
end

# ... which has to support the following operations
# to work with odeX
Base.norm(y::CompSol, p::Float64) = maximum([Base.norm(y.rho, p) abs(y.x) abs(y.p)])
Base.norm(y::CompSol) = norm(y::CompSol, 2.0)

+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s)
*(s::Real, y1::CompSol) = y1*s
/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)

### new for PR #68
Base.abs(y::CompSol) = norm(y, 2.) # TODO not needed anymore once https://github.com/JuliaLang/julia/pull/11043 is in current stable julia
Base.zero(::Type{CompSol}) = CompSol(complex(zeros(2,2)), 0., 0.)
ODE.isoutofdomain(y::CompSol) = any(isnan, vcat(y.rho[:], y.x, y.p))

# Because the new RK solvers wrap scalars in an array and because of
# https://github.com/JuliaLang/julia/issues/11053 these are also needed:
.+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
.-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
.*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s)
.*(s::Real, y1::CompSol) = y1*s
./(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)

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

# define RHSs of differential equations
# delta, V and g are parameters
function rhs(t, y, delta, V, g)
  H = [[-delta/2 V]; [V delta/2]]

  rho_dot = -im*H*y.rho + im*y.rho*H
  x_dot = y.p
  p_dot = -y.x

  return CompSol( rho_dot, x_dot, p_dot)
end

# inital conditons
rho0 = zeros(2,2);
rho0[1,1]=1.;
y0 = CompSol(complex(rho0), 2., 1.)

# solve ODEs
endt = 2*pi;

t,y1 = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., endt]) # used as reference
print("Testing interface for scalar-like state... ")
for solver in solvers
    # these only work with some Array-like interface defined:
    if solver in [ODE.ode23s, ODE.ode4s_s, ODE.ode4s_kr]
        continue
    end
    t,y2 = solver((t,y)->rhs(t, y, delta0, V0, g0), y0, linspace(0., endt, 500))
    @test norm(y1[end]-y2[end])<0.1
end
println("ok.")

################################################################################
# TODO: test a vector-like state variable, i.e. one which can be indexed.


text:
"""

```
mddf_frame_self!(iframe::Int, framedata::FrameData, options::Options, RNG, R::Result)
```

Computes the MDDF for a single frame, modifies the data in the `R` (`Result`) structure.

"""
function mddf_frame_self!(
    iframe::Int,
    framedata::FrameData,
    options::Options,
    RNG,
    R::Result,
)

    # Simplify code by assigning some shortened names
    trajectory = framedata.trajectory
    volume_frame = framedata.volume_frame
    rdf_count_random_frame = framedata.rdf_count_random_frame
    md_count_random_frame = framedata.md_count_random_frame
    dc = framedata.dc
    dmin_mol = framedata.dmin_mol
    dref_mol = framedata.dref_mol
    x_solvent_random = framedata.x_solvent_random
    lc_solvent = framedata.lc_solvent
    solute = trajectory.solute
    solvent = trajectory.solvent
    x_solute = trajectory.x_solute
    x_solvent = trajectory.x_solvent

    # Reset counters for this frame
    reset!(volume_frame)
    @. rdf_count_random_frame = 0.0
    @. md_count_random_frame = 0.0

    # get pbc sides in this frame
    sides = getsides(trajectory, iframe)

    volume_frame.total = sides[1] * sides[2] * sides[3]
    R.volume.total = R.volume.total + volume_frame.total

    R.density.solute = R.density.solute + (solute.nmols / volume_frame.total)
    R.density.solvent = R.density.solvent + (solvent.nmols / volume_frame.total)

    # Add the box side information to the box structure, in this frame
    box = Box(options.lcell, sides, R.cutoff)

    # Will wrap everything relative to the reference atom of the first molecule
    # and move everything such that that center is in the origin. This is important
    # to simplify the computation of cell indexes, as the minimum coordinates are
    # automatically -side/2 at each direction
    solute_center = x_solute[R.irefatom]
    wrap!(x_solvent, sides, solute_center)
    center_to_origin!(x_solvent, solute_center)

    # Initialize linked cells
    initcells!(x_solvent, box, lc_solvent)

    # Check if the cutoff is not too large considering the periodic cell size
    if R.cutoff > sides[1] / 2.0 || R.cutoff > sides[2] / 2.0 || R.cutoff > sides[3] / 2.0
        error(
            """
        in MDDF: cutoff or dbulk > periodic_dimension/2 in frame: $iframe
                 max(cutoff,dbulk) = $(R.cutoff)
                 sides read from file = $(sides)
                 If sides are zero it is likely that the PBC information is not available
                 in the trajectory file.
      """,
        )
    end

    n_solvent_in_bulk = 0.0
    local n_dmin_in_bulk
    for isolvent = 1:solvent.nmols-1
        # We need to do this one solute molecule at a time to avoid exploding the memory
        # requirements
        x_this_solute = viewmol(isolvent, x_solvent, solvent)

        # Compute all distances between solute and solvent atoms which are smaller than the
        # cutoff (this is the most computationally expensive part), the distances are returned
        # in the dc structure
        cutoffdistances_self!(
            R.cutoff,
            x_this_solute,
            x_solvent,
            lc_solvent,
            box,
            dc,
            solvent,
            isolvent,
        )

        # For each solute molecule, update the counters (this is highly suboptimal, because
        # within updatecounters there are loops over solvent molecules, in such a way that
        # this will loop with cost nsolute*nsolvent. However, I cannot see an easy solution
        # at this point with acceptable memory requirements
        n_dmin_in_bulk, n_dref_in_bulk =
            updatecounters!(R, solute, solvent, dc, dmin_mol, dref_mol)
        n_solvent_in_bulk += n_dref_in_bulk
    end
    # The normalization below is tricky. The number that comes out from updatecounters is the
    # sum, for every solvent molecule (minus one) of the distances that were not found to be in
    # the solute domain. The total number of distances is n^2, because the sum inside updatecounter
    # is made for all molecules (the same function is used for cross-distribution), but we called updatecounters
    # only (n-1) times, so the actual sum of the distances considered is n(n-1). From this set
    # the number of distances returned in n_dref_in_bulk is r=n(n-1)-nd/2, where nd is the
    # number of distances in the domain. Thus, the number of distances in the domain, considering
    # symmetric terms, is nd=2n(n-1)-2r. The average number of distances in the domain, per
    # molecule, is thus nd/n=2(n-1)-2r/n. Finally, the number of solvent molecules in
    # bulk, for each molecule, is the total number of other molecules, (n-1), minus the
    # number of molecules in the domain, that is (n-1)-nd/n=(n-1)-2(n-1)+2r/n, which
    # finally simplifies to 2r/n-(n-1), which is the equation below.
    n_solvent_in_bulk = 2 * n_solvent_in_bulk / solvent.nmols - (solvent.nmols - 1)

    #
    # Computing the random-solvent distribution to compute the random minimum-distance count
    #
    for i = 1:options.n_random_samples
        # generate random solvent box, and store it in x_solvent_random
        for j = 1:solvent.nmols
            # Choose randomly one molecule from the bulk, if there are actually bulk molecules
            if n_dmin_in_bulk != 0
                jmol =
                    dmin_mol[random(RNG, solvent.nmols-n_dmin_in_bulk+1:solvent.nmols)].jmol
            else
                jmol = random(RNG, 1:solvent.nmols)
            end
            # Indexes of this molecule in the x_solvent array
            x_ref = viewmol(jmol, x_solvent, solvent)
            # Indexes of the random molecule in random array
            x_rnd = viewmol(j, x_solvent_random, solvent)
            # Generate new random coordinates (translation and rotation) for this molecule
            random_move!(x_ref, R.irefatom, sides, x_rnd, RNG)
        end

        # wrap random solvent coordinates to origin
        wrap!(x_solvent_random, sides)

        # Initialize linked cells
        initcells!(x_solvent_random, box, lc_solvent)

        # Choose randomly one solute molecule to be the solute in this sample
        i_rand_mol = random(RNG, 1:solvent.nmols)
        x_this_solute = viewmol(i_rand_mol, x_solvent, solvent)

        # Compute all distances between solute and solvent atoms which are smaller than the
        # cutoff (this is the most computationally expensive part), the distances are returned
        # in the dc structure (here we do not use the self function)
        cutoffdistances!(R.cutoff, x_this_solute, x_solvent_random, lc_solvent, box, dc)

        # Update the counters and get the number of solvent molecules in bulk
        updatecounters!(
            R,
            rdf_count_random_frame,
            md_count_random_frame,
            solvent,
            dc,
            dmin_mol,
            dref_mol,
        )

    end # random solvent sampling

    # Update global counters with the data of this frame
    update_counters_frame!(
        R,
        rdf_count_random_frame,
        md_count_random_frame,
        volume_frame,
        solute,
        solvent,
        n_solvent_in_bulk,
    )

    nothing
end


text:
using ZXCalculus
using Graphs, Multigraphs
using GraphPlot: gplot
using Colors
using ZXCalculus: qubit_loc

function Multigraph2Graph(mg::Multigraph)
    g = SimpleGraph(nv(mg))
    vs = sort!(vertices(mg))
    for me in edges(mg)
        add_edge!(g, searchsortedfirst(vs, src(me)), searchsortedfirst(vs, dst(me)))
    end
    multiplicities = ["×$(mul(mg, vs[src(e)], vs[dst(e)]))" for e in edges(g)]
    for i = 1:length(multiplicities)
        if multiplicities[i] == "×1"
            multiplicities[i] = ""
        end
    end
    return g, multiplicities
end

ZX2Graph(zxd::ZXDiagram) = Multigraph2Graph(zxd.mg)
function ZX2Graph(zxg::ZXGraph)
    g = SimpleGraph(nv(zxg.mg))
    vs = sort!(vertices(zxg.mg))
    for me in edges(zxg.mg)
        add_edge!(g, searchsortedfirst(vs, src(me)), searchsortedfirst(vs, dst(me)))
    end
    # multiplicities = ["$(mul(mg, src(e), dst(e)))" for e in edges(g)]
    multiplicities = [ZXCalculus.is_hadamard(zxg, vs[src(e)], vs[dst(e)]) ? "×2" : "" for e in edges(g)]
    return g, multiplicities
end

function et2color(et::String)
    et == "" && return colorant"black"
    return colorant"blue"
end

function st2color(S::SpiderType.SType)
    S == SpiderType.Z && return colorant"green"
    S == SpiderType.X && return colorant"red"
    S == SpiderType.H && return colorant"yellow"
    S == SpiderType.In && return colorant"lightblue"
    S == SpiderType.Out && return colorant"gray"
end

ZX2nodefillc(zxd) = [st2color(zxd.st[v]) for v in sort!(vertices(zxd.mg))]

function ZX2nodelabel(zxd)
    nodelabel = String[]
    for v in sort!(vertices(zxd.mg))
        zxd.st[v] == SpiderType.Z && push!(nodelabel, "[$(v)]\n$(print_phase(zxd.ps[v]))")
        zxd.st[v] == SpiderType.X && push!(nodelabel, "[$(v)]\n$(print_phase(zxd.ps[v]))")
        zxd.st[v] == SpiderType.H && push!(nodelabel, "[$(v)]")
        zxd.st[v] == SpiderType.In && push!(nodelabel, "[$(v)]")
        zxd.st[v] == SpiderType.Out && push!(nodelabel, "[$(v)]")
    end
    return nodelabel
end

function print_phase(p)
    if typeof(p) <: Rational
        return "$(p.num)π/$(p.den)"
    else
        return "$p π"
    end
end

function layout2locs(zxd::ZXDiagram{T,P}) where {T,P}
    lo = zxd.layout
    spider_seq = ZXCalculus.spider_sequence(zxd)
    vs = sort!(spiders(zxd))
    locs = Dict()
    nqubit = lo.nbits
    frontier_v = ones(T, nqubit)
    frontier_locs = ones(nqubit)

    while sum([frontier_v[i] <= length(spider_seq[i]) for i = 1:nqubit]) > 0
        for q = 1:nqubit
            if frontier_v[q] <= length(spider_seq[q])
                v = spider_seq[q][frontier_v[q]]
                nb = neighbors(zxd, v)
                if length(nb) <= 2
                    locs[v] = (Float64(frontier_locs[q]), Float64(q))
                    frontier_locs[q] += 1
                    frontier_v[q] += 1
                else
                    v1 = nb[[qubit_loc(zxd, u) != q for u in nb]][1]
                    if spider_type(zxd, v1) == SpiderType.H
                        v1 = setdiff(neighbors(zxd, v1), [v])[1]
                    end
                    if sum([findfirst(isequal(u), spider_seq[qubit_loc(zxd, u)]) != frontier_v[qubit_loc(zxd, u)] for u in [v, v1]]) == 0
                        x = maximum(frontier_locs[min(qubit_loc(zxd, v), qubit_loc(zxd, v1)):max(qubit_loc(zxd, v), qubit_loc(zxd, v1))])
                        for u in [v, v1]
                            locs[u] = (Float64(x), Float64(qubit_loc(zxd, u)))
                            frontier_v[qubit_loc(zxd, u)] += 1
                        end
                        for q in min(qubit_loc(zxd, v), qubit_loc(zxd, v1)):max(qubit_loc(zxd, v), qubit_loc(zxd, v1))
                            frontier_locs[q] = x + 1
                        end
                    end
                end
            end
        end
    end
    for v in vs
        if !haskey(locs, v)
            v1, v2 = neighbors(zxd, v)
            x1, y1 = locs[v1]
            x2, y2 = locs[v2]
            locs[v] = ((x1+x2)/2, (y1+y2)/2)
        end
    end
    locs_x = [locs[v][1] for v in vs]
    locs_y = [locs[v][2] for v in vs]
    return locs_x, locs_y
end

function layout2locs(zxg::ZXGraph{T,P}) where {T,P}
    lo = zxg.layout
    spider_seq = ZXCalculus.spider_sequence(zxg)
    vs = sort!(spiders(zxg))
    locs = Dict()
    nqubit = lo.nbits
    frontier_v = ones(T, nqubit)
    frontier_locs = ones(nqubit)
    phase_gadget_loc = 1.0

    for v in vs
        if qubit_loc(zxg, v) !== nothing
            y = qubit_loc(zxg, v)
            x = findfirst(isequal(v), spider_seq[y])
            locs[v] = (Float64(x), Float64(y))
        else
            locs[v] = nothing
        end
    end
    for v in vs
        if locs[v] === nothing
            nb = neighbors(zxg, v)
            if length(nb) == 1
                gads = [v]
                u = v
                w = setdiff(neighbors(zxg, u), gads)[1]
                while locs[w] === nothing
                    push!(gads, w)
                    u = w
                    w = setdiff(neighbors(zxg, u), gads)[1]
                end
                push!(gads, w)
                for j = 1:(length(gads) - 1)
                    locs[gads[length(gads)-j]] = (phase_gadget_loc, Float64(nqubit + j))
                end
                phase_gadget_loc += 1
            end
        end
    end
    for v in vs
        if locs[v] === nothing
            # println(v)
            locs[v] = (phase_gadget_loc, Float64(nqubit + 1))
            phase_gadget_loc += 1
        end
    end
    locs_x = [locs[v][1] for v in vs]
    locs_y = [locs[v][2] for v in vs]
    return locs_x, locs_y
end

function plot(zxd::ZXDiagram; size_x=nothing, size_y=nothing, kwargs...)
    g, edgelabel = ZX2Graph(zxd)
    nodelabel = ZX2nodelabel(zxd)
    nodefillc = ZX2nodefillc(zxd)
    edgelabelc = colorant"black"
    if zxd.layout.nbits > 0
        locs_x, locs_y = layout2locs(zxd)
        if size_x === nothing
            size_x = maximum(locs_x) - minimum(locs_x)
        end
        if size_y === nothing
            size_y = maximum(locs_y) - minimum(locs_y)
        end
        set_default_graphic_size(3size_x*cm, 3size_y*cm)
        composition = gplot(g,
            locs_x, locs_y;
            nodelabel = nodelabel, edgelabel = edgelabel, edgelabelc = edgelabelc, nodefillc = nodefillc,
            NODESIZE = 1/(2size_x),
            kwargs...
            # EDGELINEWIDTH = 8.0 / sqrt(nv(g))
        )
        # draw(SVG("test.svg", size_x*cm, size_y*cm), composition)
    else
        gplot(g;
            nodelabel = nodelabel, edgelabel = edgelabel, edgelabelc = edgelabelc, nodefillc = nodefillc,
            kwargs...
            # NODESIZE = 0.35 / sqrt(nv(g)), EDGELINEWIDTH = 8.0 / sqrt(nv(g))
        )
    end
end
function plot(zxd::ZXGraph; size_x=nothing, size_y=nothing, kwargs...)
    g, edge_types = ZX2Graph(zxd)

    nodelabel = ZX2nodelabel(zxd)
    nodefillc = ZX2nodefillc(zxd)
    edgestrokec = et2color.(edge_types)
    if zxd.layout.nbits > 0
        locs_x, locs_y = layout2locs(zxd)
        if size_x === nothing
            size_x = maximum(locs_x) - minimum(locs_x)
        end
        if size_y === nothing
            size_y = maximum(locs_y) - minimum(locs_y)
        end
        set_default_graphic_size(3size_x*cm, 3size_y*cm)
        gplot(g,
            locs_x, locs_y;
            nodelabel = nodelabel,
            edgestrokec = edgestrokec,
            nodefillc = nodefillc,
            NODESIZE = 1/(2size_x),
            kwargs...
            )
    else
        gplot(g;
            nodelabel = nodelabel,
            edgestrokec = edgestrokec,
            nodefillc = nodefillc,
            kwargs...
            )
    end
end


text:
module Exports

###############################################################################
using ..MeshCore: AbsAttrib, VecAttrib, FunAttrib
using ..MeshCore: datavaluetype

export AbsAttrib, VecAttrib, FunAttrib
export datavaluetype

###############################################################################
using ..MeshCore: AbsShapeDesc, NoSuchShapeDesc, P1ShapeDesc, L2ShapeDesc, Q4ShapeDesc, H8ShapeDesc, T3ShapeDesc, T4ShapeDesc, L3ShapeDesc, T6ShapeDesc, Q8ShapeDesc
using ..MeshCore: manifdim, nvertices, nfacets, nridges, n1storderv, nshifts
using ..MeshCore: nfeatofdim
using ..MeshCore: SHAPE_DESC, NoSuchShape, P1, L2, T3, T4, Q4, H8, L3, T6, Q8

export AbsShapeDesc, NoSuchShapeDesc, P1ShapeDesc, L2ShapeDesc, Q4ShapeDesc, H8ShapeDesc, T3ShapeDesc, T4ShapeDesc, L3ShapeDesc, T6ShapeDesc, Q8ShapeDesc
export manifdim, n1storderv, nfacets, nridges, nshifts, nvertices
export nfeatofdim
export SHAPE_DESC, NoSuchShape, P1, L2, T3, T4, Q4, H8, L3, T6, Q8

###############################################################################
using ..MeshCore: ShapeColl
using ..MeshCore: shapedesc, nshapes, manifdim, nvertices, facetdesc, nfacets, facetconnectivity, ridgedesc, nridges, ridgeconnectivity, attribute

export ShapeColl
export shapedesc, nshapes, manifdim, nvertices, facetdesc, nfacets, facetconnectivity, ridgedesc, nridges, ridgeconnectivity, attribute

###############################################################################
using ..MeshCore: IncRel
using ..MeshCore: indextype, nrelations, nentities, ir_code
using ..MeshCore: ir_boundary, ir_subset
using ..MeshCore: ir_skeleton, ir_transpose, ir_bbyfacets, ir_bbyfacets, ir_identity

export IncRel
export indextype, nrelations, nentities, ir_code
export ir_boundary, ir_subset
export ir_skeleton, ir_transpose, ir_bbyfacets, ir_bbyfacets, ir_identity

end # module