text:
@testset "Hamiltonian splitting" begin
    import GEMPIC: set_x!, set_v!
    import GEMPIC: set_weights!, get_charge, add_charge!
    import GEMPIC: operatorHp1, operatorHp2, operatorHE, operatorHB

    # Tolerance for comparison of real numbers: set it here!
    eqv_tol = 1.0e-14

    # Set parameters
    n_particles = 2
    eta_min = 0.0
    eta_max = 4π
    num_cells = 10
    delta_t = 0.1
    degree_smoother = 3
    rnd_seed = 10

    mesh = OneDGrid(eta_min, eta_max, num_cells)

    pg = ParticleGroup{1,2}(n_particles; common_weight=1.0)

    particle_info_ref = reshape(
        [
            11.780972450961723,
            5.4977871437821380,
            -1.5341205443525459,
            0.15731068461017067,
            0.15731068461017067,
            -1.5341205443525459,
            6.8636759376074723,
            5.7026946767517002,
        ],
        n_particles,
        4,
    )

    # Initialize particles from particle_info_ref
    xi = zeros(2)
    for i_part in 1:n_particles
        xi[1] = particle_info_ref[i_part, 1]
        set_x!(pg, i_part, xi[1])
        xi[1:2] = particle_info_ref[i_part, 2:3]
        set_v!(pg, i_part, xi)
        xi[1] = particle_info_ref[i_part, 4]
        set_weights!(pg, i_part, xi[1])
    end

    # Initialize kernel smoothers
    kernel_smoother_1 = ParticleMeshCoupling1D(
        mesh, n_particles, degree_smoother - 1, :galerkin
    )

    kernel_smoother_0 = ParticleMeshCoupling1D(
        mesh, n_particles, degree_smoother, :galerkin
    )

    # Initialize Maxwell solver
    maxwell_solver = Maxwell1DFEM(mesh, degree_smoother)

    n_dofs = kernel_smoother_0.n_dofs

    ex = zeros(Float64, (n_dofs))
    efield_1 = ones(Float64, (n_dofs))
    efield_2 = ones(Float64, (n_dofs))
    bfield = ones(Float64, (n_dofs))
    efield_ref = zeros(Float64, (n_dofs, 2))
    bfield_ref = zeros(Float64, (n_dofs))
    rho = zeros(Float64, (n_dofs))

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

    # Solve Poisson problem
    compute_e_from_rho!(efield_1, maxwell_solver, rho)

    propagator = HamiltonianSplitting(
        maxwell_solver,
        kernel_smoother_0,
        kernel_smoother_1,
        pg,
        [efield_1, efield_2],
        bfield,
    )

    operatorHp1(propagator, delta_t)

    # Compare to reference
    # Particle information after operatorV application
    particle_info_ref = reshape(
        [
            11.627560396526469,
            5.5135182122431550,
            -1.5341205443525459,
            0.15731068461017067,
            0.31072273904542569,
            -1.5498516128135633,
            6.8636759376074723,
            5.7026946767517002,
        ],
        n_particles,
        4,
    )

    # Compare computed values to reference values
    for i_part in 1:n_particles
        xi = get_x(pg, i_part)[1]

        @test xi ≈ particle_info_ref[i_part, 1]

        vi = get_v(pg, i_part)

        @test vi[1] ≈ particle_info_ref[i_part, 2]
        @test vi[2] ≈ particle_info_ref[i_part, 3]

        wi = get_charge(pg, i_part)

        @test wi[1] ≈ particle_info_ref[i_part, 4]
    end

    operatorHp2(propagator, delta_t)

    # Compare to reference
    # Particle information after operatorV application
    particle_info_ref = reshape(
        [
            11.627560396526469,
            5.5135182122431550,
            -1.5030482704480033,
            2.3255233288143329e-003,
            0.31072273904542569,
            -1.5498516128135633,
            6.8636759376074723,
            5.7026946767517002,
        ],
        n_particles,
        4,
    )

    # Compare computed values to reference values
    for i_part in 1:n_particles
        xi = get_x(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 1]

        xi = get_v(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 2]
        @test xi[2] ≈ particle_info_ref[i_part, 3]

        xi[1] = get_charge(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 4]
    end

    operatorHE(propagator, delta_t)

    # Compare to reference
    # Particle information after operatorV application
    particle_info_ref = reshape(
        [
            11.627560396526469,
            5.5135182122431550,
            -1.4311461960364338,
            2.7611337138493959e-002,
            0.39511488788223620,
            -1.3903131848268702,
            6.8636759376074723,
            5.7026946767517002,
        ],
        n_particles,
        4,
    )

    # Compare computed values to reference values
    for i_part in 1:n_particles
        xi = get_x(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 1]

        xi = get_v(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 2]
        @test xi[2] ≈ particle_info_ref[i_part, 3]

        xi[1] = get_charge(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 4]
    end

    operatorHB(propagator, delta_t)

    # Compare to reference
    # Particle information after operatorV application
    particle_info_ref = reshape(
        [
            11.627560396526469,
            5.5135182122431550,
            -1.4311461960364338,
            2.7611337138493959e-002,
            0.39511488788223620,
            -1.3903131848268702,
            6.8636759376074723,
            5.7026946767517002,
        ],
        n_particles,
        4,
    )

    # Compare computed values to reference values
    for i_part in 1:n_particles
        xi = get_x(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 1]

        xi = get_v(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 2]
        @test xi[2] ≈ particle_info_ref[i_part, 3]

        xi = get_charge(pg, i_part)

        @test xi[1] ≈ particle_info_ref[i_part, 4]
    end

    bfield_ref = [
        0.96241425319343565,
        1.0704842857730226,
        0.87749470807080665,
        1.0585080035598995,
        1.0312067899733146,
        0.98524486372490983,
        0.99163866889590213,
        1.0359265946162064,
        0.96576229888372844,
        1.0213195333087741,
    ]

    efield_ref = reshape(
        [
            0.32668967300827889,
            0.21111816665325256,
            -3.2797353014133206,
            0.93182059869704992,
            2.4099520529926513,
            -0.54074935301675420,
            -6.6508143734355582e-002,
            -4.4526750013638647,
            2.7314012140039408,
            2.4952249586987651,
            1.2896538770270727,
            0.45428428470806997,
            1.9372031866314203,
            1.2380946765497305,
            0.83732795515227731,
            1.0244885213955153,
            1.1187987978241094,
            0.68789307296286761,
            1.0938576576751671,
            0.85201507883794003,
        ],
        num_cells,
        2,
    )

    @test maximum(abs.(bfield .- bfield_ref)) ≈ 0.0 atol = 1e-14
    @test maximum(abs.(efield_1 .- efield_ref[:, 1])) ≈ 0.0 atol = 1e-14
    @test maximum(abs.(efield_2 .- efield_ref[:, 2])) ≈ 0.0 atol = 1e-14
end


text:
using AutoTM, nGraph
using Test
using Flux

# TODO: Delayed until putting IO into either PMM/DRAM is handled correctly.
#include("profiler.jl")

include("utils.jl")
#include("verify.jl")
include("benchmarker.jl")

ext: jl

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module IJuliaTimeMachine

    import IJulia

    include("varchive.jl")
    include("can_copy.jl")
    include("the_past.jl")
    include("spawn.jl")

    export @past, vars

    #================================

    __init__

    =#

    function __init__()
        dontsave(VX)
        dontsave(VX.past)
        dontsave(VX.store)
        dontsave(IJulia.Out)

        IJulia.push_preexecute_hook(tm_cleanup)
        saving!()

        notify_terminal!()

        nothing
    end

    """
        unhook()

    Remove the pre and postexecute hooks created by Time Machine.
    Made for when the Time Machine is causing strange errors.
    These typically cause strange errors to happen whenever a cell is executed.
    """
    function unhook()
        try
            IJulia.pop_preexecute_hook(tm_cleanup)
        catch _
            nothing
        end

        try
            IJulia.pop_postexecute_hook(save_state)
        catch _
            nothing
        end

    end

    #===========================================

    Variables that change behavior of the time machine

    =#

    debug_mode = false
    function debug!(x::Bool = true)
        global debug_mode = x
        println("debug: $(x)")
    end

    saving = false

    """
        saving!(bool)

    Turn saving of state on or off.  Works by pushing or poping an IJulia postexecute_hook.
    True by default.
    """
    function saving!(x = true)
        global saving = x
        if saving
            IJulia.push_postexecute_hook(save_state)
        else
            IJulia.pop_postexecute_hook(save_state)
        end
    end

    notify_jupyter = false
    notify_terminal = false
    notify_gui = false

    """
        notify_jupyter!(bool)

    If true, Print notifications about finishing jobs to the current jupyter output cell.
    """
    function notify_jupyter!(x::Bool = true)
        global notify_jupyter = x
    end

    """
    notify_terminal!(bool)

    If true, Print notifications about finishing jobs to the terminal from which jupyter was started.
    """
    function notify_terminal!(x::Bool = true)
        global notify_terminal = x
    end

    #================================================

    Variables the user might want to consult

    spawned lists the cells that have been spawned off
    running is the subset that are still running
    finished is the ones that completed
    tasks is gives a reference to the task run in a given cell.
    =#
    const spawned = Set()
    const running = Set()
    const finished = Vector{Int}()
    const tasks = Dict{Int,Task}()

    #==================================

    Intended for internal use only

    =#

    """
    `vars` is a dictionary of the variables at the time.
    `ans` is the value of `ans`, which is stored separately.
    """
    #=
    struct IJulia_State
        vars::Dict
        ans::Any
    end

    const past = Dict{Int,IJulia_State}()
=#

    const VX = Varchive()

    # DontSave is a set of variables not to save
    const DontSave = Set()

    """
        dontsave(x)

    Do not save variable `x` in the history.
    """
    function dontsave(x)
        push!(DontSave, objectid(x))
        nothing
    end

    struct Queue_Pair
        n::Int
        out::Any
    end

    const out_queue = Vector{Queue_Pair}()
    const past_queue = Vector{Queue_Pair}()

    # If a process is spawned in a thread, it will need to save its results.
    # We do this in a two-phase process.  First, it puts its results into a queue.
    # Then, during the start of an execution, we run a process in the main thread that checks
    # the queue and puts its results into Out an past.
    # The part that is executed during a thread has a lock.
    tm_lock = Threads.SpinLock()

end

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


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 = "libtirpc"
version = v"1.3.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://downloads.sourceforge.net/project/libtirpc/libtirpc/$(version)/libtirpc-$(version).tar.bz2", "e24eb88b8ce7db3b7ca6eb80115dd1284abc5ec32a8deccfed2224fc2532b9fd"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libtirpc-*/

install_license COPYING

#musl does not have this header, so we manually carry it around, queue.h file contents taken from https://github.com/dbmail/dbmail/blob/master/compatibility/queue.h
if [[ ${target} == *-musl* ]]; then
    mkdir sys
    cp $WORKSPACE/srcdir/queue.h $WORKSPACE/srcdir/libtirpc-*/sys/
fi

./configure \
--prefix=${prefix} \
--build=${MACHTYPE} \
--host=${target} \
--enable-static=no \
--enable-shared=yes \
--disable-gssapi

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental = true)
#this is really only supported for linux platforms
filter!(Sys.islinux, platforms)

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

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("OpticsUtils.jl")
include("GasOptics.jl")
include("CloudOptics.jl")
include("GrayOpticsKernels.jl")

function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::AtmosphericState{FT},
    glaycol,
    sf::AbstractSourceLW{FT},
    igpt,
    lkp::LookUpLW{I, FT},
    lkp_cld::LookUpCld,
) where {I <: Int, FT <: AbstractFloat}
    ibnd = lkp.major_gpt2bnd[igpt]

    compute_optical_props_kernel!(op, as, glaycol, sf, igpt, lkp) # longwave gas optics
    add_cloud_optics_2stream(op, as, lkp, lkp_cld, glaycol, ibnd, igpt) # add cloud optics
    return nothing
end

function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::AtmosphericState{FT},
    glaycol,
    sf::AbstractSourceLW{FT},
    igpt,
    lkp::LookUpLW{I, FT},
) where {I <: Int, FT <: AbstractFloat}
    t_sfc = as.t_sfc[glaycol[2]]
    t_lev = as.t_lev
    compute_optical_props_kernel!(op, as, glaycol, igpt, lkp, sf, t_lev, t_sfc)
    return nothing
end

function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::AtmosphericState{FT},
    glaycol,
    igpt,
    lkp::LookUpSW{I, FT},
    lkp_cld::LookUpCld,
) where {I <: Int, FT <: AbstractFloat}
    ibnd = lkp.major_gpt2bnd[igpt]

    compute_optical_props_kernel!(op, as, glaycol, igpt, lkp) # shortwave gas optics
    add_cloud_optics_2stream(op, as, lkp, lkp_cld, glaycol, ibnd, igpt) # add cloud optics
    return nothing
end

function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::AtmosphericState{FT},
    glaycol::Tuple{I, I},
    igpt::I,
    lkp::LookUpLW{I, FT},
    src_args...,
) where {I <: Int, FT <: AbstractFloat}

    vmr = as.vmr
    col_dry = as.col_dry[glaycol...]
    p_lay = as.p_lay[glaycol...]
    t_lay = as.t_lay[glaycol...]
    ibnd = lkp.major_gpt2bnd[igpt]
    τ, ssa = compute_τ_ssa_lw_src!(lkp, vmr, col_dry, igpt, ibnd, p_lay, t_lay, glaycol, src_args...)

    @inbounds op.τ[glaycol...] = τ

    if op isa TwoStream
        @inbounds op.ssa[glaycol...] = ssa
        @inbounds op.g[glaycol...] = FT(0) # initializing asymmetry parameter
    end
    return nothing
end

function compute_optical_props_kernel!(
    op::AbstractOpticalProps{FT},
    as::AtmosphericState{FT},
    glaycol::Tuple{I, I},
    igpt::I,
    lkp::LookUpSW{I, FT},
) where {I <: Int, FT <: AbstractFloat}

    vmr = as.vmr
    col_dry = as.col_dry[glaycol...]
    p_lay = as.p_lay[glaycol...]
    t_lay = as.t_lay[glaycol...]
    ibnd = lkp.major_gpt2bnd[igpt]
    τ, ssa = compute_τ_ssa_lw_src!(lkp, vmr, col_dry, igpt, ibnd, p_lay, t_lay, glaycol)

    @inbounds op.τ[glaycol...] = τ

    if op isa TwoStream
        @inbounds op.ssa[glaycol...] = ssa
        @inbounds op.g[glaycol...] = FT(0) # initializing asymmetry parameter
    end
    return nothing
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
export Semiring, Probability_semiring, Tropical_semiring, Log_semiring

# A semiring type that you can create. Currently the onus is on the user to
# make sure the semiring makes sense.
immutable Semiring
	addition::Function
	multiplication::Function
	ring_zero::Float64
	ring_one::Float64
end

Probability_semiring = Semiring(+, *, 0, 1)
Tropical_semiring = Semiring(min, +, Inf, 0)
Log_semiring = Semiring((x,y) -> -log(e^-x + e^-y), +, Inf, 0)


text:
include("../examples/nonconvex_smooth_noncompact.jl")
function lyap_exp(w, η)
	N = 20000
	le = 0.0
	for n = 1:N
		le += log(abs(dnext(w, η)))/N
		w = next(w, η)
	end
	return le
end
function run(w, η, n, s)
	for i = 1:n
		w = next(w, η, s)
	end
	return w
end

max_stars_repo_name: nishaChandramoorthy/neuralOMET

max_issues_repo_name: nishaChandramoorthy/neuralOMET

max_forks_repo_name: nishaChandramoorthy/neuralOMET


text:
"""
    getDocAnnotation(omc:: OMCSession, name:: String)

Returns documentation string from model annotation and strips `<html>` tag if
necessary.
"""
function getDocAnnotation(omc:: OMCSession, model:: String)
    htmldoc = sendExpression(omc, "getDocumentationAnnotation($model)")[1]
    ishtml = match(r"\s*<\s*html[^>]*>(.*)\s*</\s*html\s*>\s*"s, htmldoc)
    if !isnothing(ishtml)
        htmldoc = ishtml.captures[1]
    end
    return htmldoc
end

"""
    getcode(omc:: OMCSession, model:: String)

Returns the code of the given model or class as a string.
"""
function getcode(omc:: OMCSession, model:: String)
    cwd = sendExpression(omc, "cd()")
    tmpfile = joinpath(cwd, "export_$model.mo")
    sendExpression(omc, "saveModel(\"$(moescape(tmpfile))\", $model)")
    err = getErrorString(omc)
    if !isempty(err)
        throw(MoSTError("Could not export model code", err))
    end
    return read(tmpfile, String)
end

"""
    getequations(omc:: OMCSession, model::String)

Returns all equations of the given model as a list of strings with
presentation MathML syntax.
"""
function getequations(omc:: OMCSession, model::String)
    res = sendExpression(omc, "dumpXMLDAE($model, addMathMLCode=true)")
    err = getErrorString(omc)
    if !isempty(err)
        throw(MoSTError("Could not save model as independent XML file", err))
    end
    res = py"extract_equations"(res[2], xslt_dir=joinpath(@__DIR__, "..", "res"))
    return res
end

"""
    getfunctions(omc:: OMCSession, model:: String)

Returns all functions of the given model as a twodimensional array where each
row contains the function name, the function header, and the modelica definition
of the function in this order.
"""
function getfunctions(omc:: OMCSession, model:: String)
    res = sendExpression(omc, "dumpXMLDAE($model, addMathMLCode=true)")
    err = getErrorString(omc)
    if !isempty(err)
        throw(MoSTError("Could not save model as independent XML file", err))
    end
    funcs = py"extract_functions"(res[2])
    for i in 1:size(funcs)[1]
        # removes function name from end of header (i.e. from return type)
        if endswith(funcs[i,2], funcs[i,1])
            funcs[i,2] = funcs[i,2][1:end-length(funcs[i,1])]
        end
        # removes trailing dot from function type (present in versions < 1.17.0)
        if endswith(funcs[i,2], ".")
            funcs[i,2] = funcs[i,2][1:end-1]
        end
        funcs[i, 3] = replace(funcs[i, 3], "\"Inline if necessary\"" => "")
    end
    if isempty(funcs)
        return Array{String}(undef, 0, 2)
    else
        return funcs
    end
end

function algorithm(code:: AbstractString)
    funcname = match(
        r"^\s*function ([\w.$]+)",
        split(code, r"\r?\n")[1]
    ).captures[1]
    return match(
        Regex("algorithm\\n((?:.|\\n)+)end \\Q$funcname\\E;"),
        code
    ).captures[1]
end

function replacefuncsincode(code:: AbstractString, aliases:: Dict)
    for a in keys(aliases)
        code = replace(
            code,
            Regex("\\Q$a\\E\\(") => "$(aliases[a])("
        )
    end
    return code
end

function uniquefunctions(funcs:: Array)
    alg2func = Dict{String, Array{String, 2}}()
    for i in 1:size(funcs)[1]
        a = algorithm(funcs[i, 3])
        alg2func[a] = get(alg2func, a, Array{String}(undef, 0, 3))
        alg2func[a] = [alg2func[a]; permutedims(funcs[i, :])]
    end
    aliases = Dict()
    res = Array{String}(undef, 0, 3)
    if size(funcs)[1] == 0
        return aliases, res
    end
    for a in keys(alg2func)
        group = alg2func[a]
        prototype = group[1, :]
        res = [res; permutedims(prototype)]
        for i in 1:size(group)[1]
            aliases[group[i, 1]] = prototype[1,1]
        end
    end
    # update uses of functions in other functions
    for i in 1:size(res)[1]
        res[i, 3] = replacefuncsincode(res[i, 3], aliases)
    end
    if maximum([size(x)[1] for x in values(alg2func)]) > 1
        # there were some simplifications => repeat
        (na, res) = uniquefunctions(res)
        # update aliases
        for a in keys(aliases)
            if haskey(na, aliases[a])
                aliases[a] = na[aliases[a]]
            end
        end
    end
    return (aliases, res)
end

"""
    uniquehierarchy(names:: Array{<: AbstractString})

Returns a dictionary that maps the given names to the minimal hierarchical
postfix that is required to make the name unique.
A hierarchical postfix is a postfix that begins after a dot.
For example, `ab.cd.ef` has the hierarchical postfixes `ef`, `cd.ef` and
`ab.cd.ef`.
"""
function uniquehierarchy(names:: Array{<: AbstractString})
    function postfix(n, i)
        join(split(n, ".")[end-i+1:end], ".")
    end
    restoredollar = Dict(replace(x, '$' => '.') => x for x in names)
    res = Dict()
    remaining = Set(keys(restoredollar))
    i = 1
    while !isempty(remaining)
        # count names ending with .postfix
        # == 0: Postfix is full name of one of the elements => must be unique
        # == 1: Postfix is unique
        #  > 1: Postfix is not unique => continue
        namesdone = filter(
            x -> count(a -> endswith(a, '.' * postfix(x, i)), remaining) <= 1,
            remaining
        )
        for n in namesdone
            res[restoredollar[n]] = postfix(n, i)
            delete!(remaining, n)
        end
        i += 1
    end
    return res
end

"""
    getvariables(omc:: OMCSession, model::String)

Returns an array of dictionaries that contain the following keys describing
the variables and parameters of the given model.

* "name": the name of the variable
* "variability": is the variable a parameter, a constant, or a variable?
* "type": the data type of the variable (usually "Real")
* "label": the string label attached to the variable
* "quantity": the kind of physical quantity modeled by the variable
* "unit": the unit of the variable in Modelica unit syntax
* "initial": the initial value of the variable, if existent
* "bindExpression": the expression used to fix the variable value, if existent
* "aliasof": the name of the other variable of which this variable is an alias

An empty string is used as value for keys which are not applicable for the given variable.
"""
function getvariables(omc:: OMCSession, model::String)
    res = sendExpression(omc, "dumpXMLDAE($model, addMathMLCode=true)")
    err = getErrorString(omc)
    if !isempty(err)
        throw(MoSTError("Could not save model as independent XML file", err))
    end
    res = py"extract_variables"(res[2], xslt_dir=joinpath(@__DIR__, "..", "res"))
    return res
end

"""
    mdescape(s:: String)

Escapes characters that have special meaning in Markdown with a backslash.
"""
function mdescape(s:: String)
    escape_chars = "\\`*_#+-.!{}[]()\$"
    return join([(x in escape_chars ? "\\$x" : x) for x in s])
end

"""
    variabletable(vars:: Array{Dict{Any, Any},1})

Creates a Markdown table from an array of variable descriptions as returned
by [`getvariables(omc:: OMCSession, model::String)`](@ref).
"""
function variabletable(vars:: Array{Dict{Any, Any},1})
    header = """
    |name|unit|label|value|
    |----|----|-----|-----|
    """
    lines = []
    for v in vars
        if length(v["aliasof"]) > 0
            continue # exclude aliases from table
        end
        value = if length(v["initial"]) == 0
            v["bindExpression"]
        else
            v["initial"]
        end
        vals = [v["name"], v["unit"], v["label"], value]
        vals = [mdescape(x) for x in vals]
        push!(lines, "|$(join(vals, "|"))|")
    end
    table = "$(header)$(join(lines, "\n"))"
    return Markdown.parse(table)
end

function replacefuncnames(equation:: AbstractString, replacements:: Dict)
    res = equation
    for n in keys(replacements)
        res = replace(
            res,
            Regex("<mi>\\s*\\Q$n\\E\\s*</mi>")
            => SubstitutionString("<mi>$(replacements[n])</mi>")
        )
    end
    return res
end

const ellipse = "_"

function equationlist(equations:: Array{<: AbstractString}, vars:: Array{Dict{Any, Any}, 1}, funcs:: Array)
    function hierarchify(equations:: Array)
        res = Dict()
        for (pref, eq) in equations
            levels = split(pref, ".")
            target = res
            for l in levels
                if l == "" continue end
                if !haskey(target, l) target[l] = Dict() end
                target = target[l]
            end
            if !haskey(target, "") target[""] = [] end
            push!(target[""], eq)
        end
        return res
    end
    function htmlify(hierarchy:: Dict)
        entries = []
        for k in sort(collect(keys(hierarchy)))
            v = hierarchy[k]
            if k == ""
                for e in v
                    push!(entries, "
$e")
                end
            else
                push!(entries, "
Within group $k (prefix $ellipse indicates shortened variable name)\n$(htmlify(v))")
            end
        end
        """

        $(join(entries, "\n"))

        """
    end
    (funcdict, funcs) = uniquefunctions(funcs)
    equations = [replacefuncnames(e, funcdict) for e in equations]
    funcdict = uniquehierarchy(funcs[1:end, 1])
    equations = [replacefuncnames(e, funcdict) for e in equations]
    equations = collect(map(explicify, equations))
    aliases = aliasdict(vars)
    prefixes = [commonhierarchy(e, aliases) for e in equations]
    pruneprefixes!(prefixes)
    equations = [(p, deprefix(e, p)) for (e, p) in zip(equations, prefixes)]
    return htmlify(hierarchify(equations))
end

function functionlist(funcs:: Array)
    (funcdict, funcs) = uniquefunctions(funcs)
    replacements = uniquehierarchy(funcs[1:end, 1])
    res = ["Functions:"]
    for i in 1:size(funcs)[1]
        fun = funcs[i, 1]
        rep = replacements[fun]
        code = replacefuncsincode(funcs[i, 3], replacements)
        cleaned = strip(replace(code, fun => rep))
        push!(res, "```modelica\n$cleaned\n```")
    end
    return Markdown.parse(join(res, "\n\n"))
end

# extend Documenter with new code block type @modelica
abstract type ModelicaBlocks <: Documenter.Expanders.ExpanderPipeline end
Documenter.Selectors.order(::Type{ModelicaBlocks}) = 5.0
Documenter.Selectors.matcher(::Type{ModelicaBlocks}, node, page, doc) = Documenter.Expanders.iscode(node, "@modelica")
function Documenter.Selectors.runner(::Type{ModelicaBlocks}, x, page, doc)
    cd(page.workdir) do
        result = ""
        modelnames = []
        result = []
        nullary(x) = Regex("%\\s*$x\\s*")
        unary(x) = Regex("%\\s*$x\\s*=\\s*(.*)")
        magics = Dict(
            "modeldir" => unary("modeldir"),
            "outdir" => unary("outdir"),
            "nocode" => nullary("nocode"),
            "noequations" => nullary("noequations"),
            "noinfo" => nullary("noinfo"),
            "omcargs" => unary("omcargs"),
            "headlevel" => unary("headlevel"),
            "libs" => unary("libs")
        )
        magicvalues = Dict()
        # get list of models and magic values
        for (line) in split(x.code, '\n')
            if startswith(line, '%')
                try
                    matches = [(k, match(v, line)) for (k, v) in magics]
                    matches = filter(p -> !isnothing(p[2]), matches)
                    if isempty(matches)
                        throw(MoSTError("magic line type in '$line' not recognized", ""))
                    end
                    name, rmatch = first(matches)
                    magicvalues[name] = if isempty(rmatch.captures)
                        true
                    else
                        string(rmatch.captures[1])
                    end
                catch err
                    push!(doc.internal.errors, :eval_block)
                    @warn("""
                        invalid magic line starting with '%' in `@modelica` block in $(Documenter.Utilities.locrepr(page.source))
                        ```$(x.language)
                        $(x.code)
                        ```
                        """, exception = err)
                end
            else
                push!(modelnames, String(line))
            end
        end
        # communicate with OMC to obtain documentation and equations
        try
            modeldir = get(magicvalues, "modeldir", "../..")
            outdir = get(magicvalues, "outdir", joinpath(modeldir, "../out"))
            withOMC(outdir, modeldir) do omc
                omcargs = get(magicvalues, "omcargs", "")
                sendExpression(omc, "setCommandLineOptions(\"$(moescape(omcargs))\")")
                for libstr in split(get(magicvalues, "libs", ""), " ")
                    if isempty(libstr)
                        continue
                    end
                    if occursin("@", libstr)
                        lib, version = split(libstr, "@")
                    else
                        lib = libstr
                        version = "latest"
                    end
                    installAndLoad(omc, lib; version=version)
                end
                for (model) in modelnames
                    # TODO automatically decide what to do based on class type
                    # load model without all extra checks
                    loadModel(omc, model; check=false)
                    # add header to result
                    hn = get(magicvalues, "headlevel", 3)
                    header = Documenter.Utilities.mdparse("$("#"^hn) $model")
                    push!(result, header)
                    # get documentation as HTML string
                    if !get(magicvalues, "noinfo", false)
                        htmldoc = getDocAnnotation(omc, model)
                        push!(result, Documenter.Documents.RawHTML(htmldoc))
                    end
                    # get model code
                    if !get(magicvalues, "nocode", false)
                        rawcode = getcode(omc, model)
                        push!(result, Documenter.Utilities.mdparse("```modelica\n$rawcode\n```\n"))
                    end
                    # TODO alternative way to get equations through getNthEquationItem()
                    # get model equations
                    if !get(magicvalues, "noequations", false)
                        equations = getequations(omc, model)
                        vars = getvariables(omc, model)
                        funcs = getfunctions(omc, model)
                        htmleqs = equationlist(equations, vars, funcs)
                        push!(result, Documenter.Documents.RawHTML(htmleqs))
                        funclist = functionlist(funcs)
                        push!(result, funclist)
                        vartab = variabletable(vars)
                        push!(result, vartab)
                    end
                end
            end
        catch err
            if err isa MoSTError
                push!(doc.internal.errors, :eval_block)
                @warn("""
                    failed to evaluate `@modelica` block in $(Documenter.Utilities.locrepr(page.source))
                    ```$(x.language)
                    $(x.code)
                    ```
                    """, exception = err)
            end
            rethrow(err)
        end
        page.mapping[x] = Documenter.Documents.MultiOutput(result)
    end
end

function commonprefix(str:: AbstractString...)
    # based on https://rosettacode.org/wiki/Longest_common_prefix#Julia
    if isempty(str) return "" end
    i = 1
    ref = str[1]
    while all(i ≤ length(s) && s[i] == ref[i] for s in str)
        i += 1
    end
    return ref[1:i-1]
end

function commonprefix(aliasgroups:: Set{<:AbstractString}...)
    if isempty(aliasgroups) return "" end
    res = []
    # if all aliasgroups have >= 2 elements, iterate over smallest group
    smallest = reduce((x, y) -> if length(x) < length(y) x else y end, aliasgroups)
    for ref in smallest
        push!(res, commonprefix(collect(aliasgroups), ref))
    end
    largestprefix = reduce((x, y) -> if length(x) > length(y) x else y end, res)
    return largestprefix
end

function commonprefix(aliasgroups:: Array{<:Set{<:AbstractString},1}, ref:: AbstractString)
    if isempty(aliasgroups) return "" end
    aliasgroups = deepcopy(aliasgroups)
    i = 1
    while all(length(g) > 0 for g in aliasgroups) && length(ref) ≥ i
        for g in aliasgroups
            filter!(x -> length(x) ≥ i && x[i] == ref[i], g)
        end
        i += 1
    end
    return ref[1:i-1]
end

function commonhierarchy(str:: Union{AbstractString, Set{<:AbstractString}}...)
    pref = commonprefix(str...)
    i = findlast('.', pref)
    if isnothing(i) return "" end
    return pref[1:i-1]
end

function findvarnames(str:: AbstractString)
    idnames = findidentifiers(str)
    fnames = findfuncnames(str)
    varnames = setdiff(idnames, fnames)
    return varnames
end

function findfuncnames(str:: AbstractString)
    pattern = r"<mi>\s*([\w.$]+)\s*<\/mi>\s*</mrow>\s*<mo></mo>"
    funcnames = [x.captures[1] for x in eachmatch(pattern, str)]
    return funcnames
end

function findidentifiers(str:: AbstractString)
    mi = r"<mi>\s*([\w.$]+)\s*<\/mi>"
    identifiers = [x.captures[1] for x in eachmatch(mi, str)]
    return identifiers
end

const commonconsts = Set(["e", "π", "c", "h", "G", "F", "R"])

function commonhierarchy(str:: AbstractString, aliases:: Dict{<:AbstractString, <:Set{<:AbstractString}})
    varnames = findvarnames(str)
    varnames = setdiff(varnames, commonconsts)
    aliasgroups = [get(aliases, n, Set{String}()) ∪ Set([n]) for n in varnames]
    pref = commonhierarchy(aliasgroups...)
    return pref
end

function deprefix(str:: AbstractString, pref:: AbstractString)
    de = replace(
        str,
        Regex("<mi>\\s*$pref\\.([\\w.]+)\\s*</mi>")
        => SubstitutionString("<mi>$ellipse\\1</mi>")
    )
    return de
end

function explicify(mml:: AbstractString)
    # use visible "dot operator" instead of "invisible times"
    # justification: our variables have multi-character names
    replace(mml, "" => "⋅")
end

function pruneprefixes!(prefs:: Array{<: AbstractString})
    counts = Dict(k => 0 for k in prefs)
    for p in prefs
        counts[p] += 1
    end
    while !all(([counts[p] for p in prefs] .> 1) .| (prefs .== ""))
        # index of longest single entry
        pi = argmax(map(x -> if counts[x] == 1 length(x) else 0 end, prefs))
        delete!(counts, prefs[pi])
        # remove last hierarchical level from prefix
        prefs[pi] = join(split(prefs[pi], ".")[1:end-1], ".")
        counts[prefs[pi]] = get(counts, prefs[pi], 0) + 1
    end
end

function aliasdict(vars:: Array{Dict{Any, Any},1})
    aliases = Dict{String, Set{String}}()
    for v in vars
        alias = v["name"]
        original = v["aliasof"]
        if length(original) == 0 continue end
        aliases[original] = get(aliases, original, Set{String}())
        push!(aliases[original], alias)
    end
    return aliases
end

function __init__doc()
    # import only used to install lxml automatically
    pyimport_conda("lxml.etree", "lxml")
    py"""
    import lxml.etree as et
    import lxml.objectify as lo
    import io
    import os

    def load_ctop(dirname):
        xslt = et.parse(os.path.join(dirname, "ctop.xsl"))
        tf = et.XSLT(xslt)
        return tf

    def cleanup_mathml(dom):
        # reference: https://stackoverflow.com/a/18160164
        for e in dom.getiterator():
            if not hasattr(e.tag, "find"):
                continue
            if e.tag.endswith("math"):
                continue
            i = e.tag.find('}')
            if i >= 0:
                e.tag = e.tag[i+1:]
        lo.deannotate(dom.getroot(), cleanup_namespaces=True)

    def fix_function_applications(dom, class_name, ns={}):
        functions = [str(x) for x in dom.xpath("/dae/functions/function/@name")]
        applies = dom.xpath("//mml:apply/*[1]", namespaces=ns)
        for app in applies:
            tag_name = et.QName(app).localname.replace("_dollar_", "$")
            if tag_name in functions:
                app.tag = et.QName(ns["mml"], "ci")
                app.text = tag_name

    def extract_variables(fname, xslt_dir="."):
        dom = et.parse(fname)
        vars = dom.xpath("//variable")
        result = []
        for v in vars:
            isalias = str(v.getparent().getparent().tag) == "aliasVariables"
            vdict = {
                "name": v.get("name"),
                "variability": v.get("variability"),
                "type": v.get("type"),
                "label": v.get("comment"),
                "quantity": v.xpath("string(attributesValues/quantity/@string)"),
                "unit": v.xpath("string(attributesValues/unit/@string)"),
                "initial": v.xpath("string(attributesValues/initialValue/@string)"),
                "bindExpression": v.xpath("string(bindExpression/@string)")
            }
            vdict["aliasof"] = vdict["bindExpression"] if isalias else None
            for k in vdict:
                if vdict[k] is None:
                    vdict[k] = ""
            # if len(vdict["bindExpression"]) > 0:
            #     vdict["bindExpression"] = c2p(vdict["bindExpression"], xslt_dir=xslt_dir)
            result.append(vdict)
        return result

    def c2p(content, xslt_dir="."):
        single = False
        if not isinstance(content, list):
            content = [content]
            single = True
        content_to_pres = load_ctop(xslt_dir)
        presentation = [content_to_pres(c) for c in content]
        for p in presentation:
            cleanup_mathml(p)
        if single:
            return presentation[0]
        else:
            return presentation

    def extract_equations(fname, xslt_dir="."):
        dom = et.parse(fname)
        ns = {"mml": "http://www.w3.org/1998/Math/MathML"}
        fix_function_applications(dom, os.path.splitext(os.path.basename(fname))[0], ns=ns)
        mathdoms = dom.xpath("/dae/equations/equation/MathML/mml:math", namespaces=ns)
        newdoms = c2p(mathdoms, xslt_dir=xslt_dir)
        return [et.tostring(x) for x in newdoms]

    def extract_functions(fname):
        dom = et.parse(fname)
        functions = dom.xpath("//function")
        res = []
        for f in functions:
            name = f.get("name")
            impl = f[0].text
            headi = impl.find("function " + name)
            head = impl[:headi]
            impl = impl[headi:]
            res.append([name, head, impl])
        return res
    """
end

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
include("experiments/rarefaction_wave.jl")
include("experiments/sod_shock_tube.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# Macro for checking arguments
macro check_args(K, param, cond, desc=string(cond))
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(K)), ": ", $(string(param)), " = ", $(esc(param)), " does not ",
                "satisfy the constraint ", $(string(desc)), ".")))
        end
    end
end

# Take highest Float among possibilities
# function promote_float(Tₖ::DataType...)
#     if length(Tₖ) == 0
#         return Float64
#     end
#     T = promote_type(Tₖ...)
#     return T <: Real ? T : Float64
# end

check_dims(K,X,Y,featdim,obsdim) = check_dims(X,Y,featdim,obsdim) && (size(K) == (size(X,obsdim),size(Y,obsdim)))

check_dims(X,Y,featdim,obsdim) = size(X,featdim) == size(Y,featdim)

feature_dim(obsdim::Int) = obsdim == 1 ? 2 : 1


text:
using Documenter
using ReinforcementLearning, ReinforcementLearningEnvironments

"filter concrete subtypes"
is_concrete_type_of(t) = x ->begin
    println(x)
    x isa Type && x <: t && !isabstracttype(x)
end

makedocs(
    modules = [ReinforcementLearning],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-149861753-1",
        canonical = "https://juliareinforcementlearning.github.io/ReinforcementLearning.jl/latest/",
        assets = [
            "assets/favicon.ico",
            "assets/custom.css"
            ],
        ),
    sitename = "ReinforcementLearning.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "A Quick Example" => "a_quick_example.md",
        "Overview" => "overview.md",
        "Manual" => [
            "Components" => [
                "Agents" => "components/agents.md",
                "Environments" => "components/environments.md",
                "Buffers" => "components/buffers.md",
                "Policies" => "components/policies.md",
                "Learners" => "components/learners.md",
                "Approximators" => "components/approximators.md",
                "Action Selectors" => "components/action_selectors.md",
                "Environment Models" => "components/environment_models.md"
            ],
            "Core" => "core.md",
            "Utils" => "utils.md"
        ],
        "Tips for Developers" => "tips_for_developers.md",
        "Experiments" => [
            "Play Atari Games with DQN" => "experiments/atari_dqn.md"
        ]],
)

deploydocs(
    repo = "github.com/JuliaReinforcementLearning/ReinforcementLearning.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module FDFD

using Peacock

include("solver_fdm.jl")
export Solver, solve, FDFDBasis

end


text:
@testset "system" begin
    include("system/core.jl")
    include("system/clock.jl")
    include("system/controller.jl")
    include("system/calendar.jl")
    include("system/store.jl")
end


text:
function starts_in(arc, vertex, arc_set)
    return(arc_set[arc, :start_node]==vertex)
end

function ends_in(arc, vertex, arc_set)
    return(arc_set[arc, :end_node]==vertex)
end

function solve_slave(graph, data, tour_pool, branch_pool, depot, dual_values, node_branch)
    """Construct a tour starting and ending at the same depot wich is
    continuous and respects the capacity of the vehicle."""

    ### Constants ###

    # Number of arcs/edges
    card_Ar = data[1,:req_arcs]
    card_Anr = data[1,:noreq_arcs]
    card_Er = data[1,:req_edges]
    card_Enr = data[1,:noreq_edges]

    # Number of nodes
    nb_nodes = data[1,:nodes]

    # From an edge, create two opposite arcs
    Er_swap = copy(graph["req_edges"])
    Er_swap[:start_node], Er_swap[:end_node] = Er_swap[:end_node], Er_swap[:start_node]
    Enr_swap = copy(graph["noreq_edges"])
    Enr_swap[:start_node], Enr_swap[:end_node] = Enr_swap[:end_node], Enr_swap[:start_node]

    # All arcs
    A = vcat(graph["req_arcs"], graph["noreq_arcs"], graph["req_edges"],
        graph["noreq_edges"], Er_swap, Enr_swap)

    # Required arcs
    R = vcat(graph["req_arcs"], graph["req_edges"], Er_swap)

    # Demand on required arcs
    demand = R[:demand]

    # Costs of edges
    serv_cost = R[:serv_cost]
    trav_cost = A[:trav_cost]

    ### Model definiton ###
    model = Model(with_optimizer(CPLEX.Optimizer, CPX_PARAM_SCRIND=0))

    ### Variables ###

    # Service variables
    @variable(model, X_Ar[1:card_Ar], Bin)
    @variable(model, X_Er[1:card_Er], Bin)
    @variable(model, X_Er_swap[1:card_Er], Bin)
    X_R = vcat(X_Ar, X_Er, X_Er_swap)
    X_tot = vcat(X_Ar, [0 for arc in 1:card_Anr],
        X_Er, [0 for arc in 1:card_Enr],
        X_Er_swap, [0 for arc in 1:card_Enr])

    # Passage without service variables
    @variable(model, Y_Ar[1:card_Ar], Int)
    @constraint(model, positivity_Ar, Y_Ar.>=0)
    @variable(model, Y_Anr[1:card_Anr], Int)
    @constraint(model, positivity_Anr, Y_Anr.>=0)
    @variable(model, Y_Er[1:card_Er], Int)
    @constraint(model, positivity_Er, Y_Er.>=0)
    @variable(model, Y_Enr[1:card_Enr], Int)
    @constraint(model, positivity_Enr, Y_Enr.>=0)
    @variable(model, Y_Er_swap[1:card_Er], Int)
    @constraint(model, positivity_Er_swap, Y_Er_swap.>=0)
    @variable(model, Y_Enr_swap[1:card_Enr], Int)
    @constraint(model, positivity_Enr_swap, Y_Enr_swap.>=0)
    Y_tot = vcat(Y_Ar, Y_Anr, Y_Er, Y_Enr, Y_Er_swap, Y_Enr_swap)

    # Flow variables
    @variable(model, F[1:(card_Ar+card_Anr+2*(card_Er+card_Enr))] >= 0)

    ### Objective ###
    sum_1 = try
        sum(dual_values["req_arcs"].*X_Ar)
    catch
        0
    end
    sum_2 = try
        sum(dual_values["req_edges"].*(X_Er+X_Er_swap))
    catch
        0
    end

    sum_3 = 0
    for (j,i) in enumerate(node_branch)
        if (branch_pool[i][2]=="req_arcs") & (branch_pool[i][1]==depot)
            sum_3 += try
                X_Ar[branch_pool[i][3]]*(-dual_values["branch_upper"][j]-dual_values["branch_lower"][j])
            catch
                0
            end
        elseif (branch_pool[i][2]=="req_edges") & (branch_pool[i][1]==depot)
            sum_3 += try
                X_Er[branch_pool[i][3]]*(-dual_values["branch_upper"][j]-dual_values["branch_lower"][j])
            catch
                0
            end
        elseif (branch_pool[i][2]=="req_edges_swap") & (branch_pool[i][1]==depot)
            sum_3 += try
                X_Er_swap[branch_pool[i][3]]*(-dual_values["branch_upper"][j]-dual_values["branch_lower"][j])
            catch
                0
            end
        end
    end

    @objective(model, Min,
        sum(X_R.*serv_cost) +
        sum(Y_tot.*trav_cost) +
        dual_values["vehicle"] -
        sum_1 -
        sum_2 +
        sum_3
        )

    ### Constraints ###

    # Continuity constraint
    @constraint(model, continuity[i=1:nb_nodes],
        sum(map(arc->starts_in(arc,i,A) ? Y_tot[arc] : 0, 1:nrow(A)))
        + sum(map(arc->starts_in(arc,i,R) ? X_R[arc] : 0, 1:nrow(R)))
        == sum(map(arc->ends_in(arc,i,A) ? Y_tot[arc] : 0, 1:nrow(A)))
        + sum(map(arc->ends_in(arc,i,R) ? X_R[arc] : 0, 1:nrow(R)))
        )

    # Dumping constraint
    @constraint(model, dumping,
        sum(map(arc->starts_in(arc,depot,A) ? Y_tot[arc] : 0, 1:nrow(A)))
        +  sum(map(arc->starts_in(arc,depot,R) ? X_R[arc] : 0, 1:nrow(R)))
        <= 1
        )

    # Flow constraint 1
    @constraint(model, flow_1[i=[node for node in 1:nb_nodes if node!=depot]],
        sum(map(arc->ends_in(arc,i,A) ? F[arc] : 0, 1:nrow(A)))
        - sum(map(arc->starts_in(arc,i,A) ? F[arc] : 0, 1:nrow(A)))
        == sum(map(arc->ends_in(arc,i,R) ? X_R[arc]*demand[arc] : 0, 1:nrow(R)))
        )

    # Flow constraint 2
    @constraint(model, flow_2,
        sum(map(arc->starts_in(arc,depot,A) ? F[arc] : 0, 1:nrow(A)))
        == sum(X_R.*demand)
        )

    # Flow constraint 3
    @constraint(model, flow_3,
        sum(map(arc->ends_in(arc,depot,A) ? F[arc] : 0, 1:nrow(A)))
        == sum(map(arc->ends_in(arc,depot,R) ? X_R[arc]*demand[arc] : 0, 1:nrow(R)))
        )

    # Capacity constraint
    @constraint(model, capacity[arc=1:nrow(A)],
        F[arc] <= data[1,:capa]*(Y_tot[arc]+X_tot[arc])
        )

    # Avoid servicing twice on edges
    @constraint(model, avoid_return[edge=1:data[1,:req_edges]],
        X_Er[edge]+X_Er_swap[edge] <= 1
        )

    ### Solve slave ###

    optimize!(model)

    ### Format output ###

    reduced_cost = objective_value(model)

    cost = sum(value.(X_R).*serv_cost) +sum(value.(Y_tot).*trav_cost)

    X = Dict(
        "req_arcs" => [value.(X_Ar[arc]) for arc in 1:length(X_Ar)],
        "req_edges" => [value.(X_Er[arc]) for arc in 1:length(X_Er)],
        "req_edges_swap" => [value.(X_Er_swap[arc]) for arc in 1:length(X_Er_swap)]
        )

    Y = Dict(
        "req_arcs" => [value.(Y_Ar[arc]) for arc in 1:length(Y_Ar)],
        "noreq_arcs" => [value.(Y_Anr[arc]) for arc in 1:length(Y_Anr)],
        "req_edges" => [value.(Y_Er[arc]) for arc in 1:length(Y_Er)],
        "noreq_edges" => [value.(Y_Enr[arc]) for arc in 1:length(Y_Enr)],
        "req_edges_swap" => [value.(Y_Er_swap[arc]) for arc in 1:length(Y_Er_swap)],
        "noreq_edges_swap" => [value.(Y_Enr_swap[arc]) for arc in 1:length(Y_Enr_swap)]
        )

    tour = Dict(
        "depot" => depot,
        "reduced_cost" => reduced_cost,
        "cost" => cost,
        "capacity" => sum(value.(X_R).*demand),
        "X" => X,
        "Y" => Y
        )

    tour_index = length(tour_pool)+1
    push!(tour_pool, tour)

    return(tour_index, tour_pool)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
module SessionTests

using Test
using GCP
using ..GCPTests

@testset "Session" begin

  @test begin
    sess = GCP.Session(GCPTests.CredFilename, ["cloud-platform"])
    sess != nothing
  end

end

end # module SessionTests

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
"""
```
init_subspec!(m::DSGEVAR)
```
initializes a model subspecification of a DSGEVAR. We utilize
subspecifications to distinguish between possible VARs, such as
specifying the observables and the number of lags.
"""
function init_subspec!(m::DSGEVAR)
    if subspec(m) == "ss0"
        return # initializes an empty DSGEVAR
    elseif subspec(m) == "ss1"
        ss1!(m)
    elseif subspec(m) == "ss2"
        ss2!(m)
    elseif subspec(m) == "ss3"
        ss3!(m)
    elseif subspec(m) == "ss10"
        ss10!(m)
    elseif subspec(m) == "ss11"
        ss11!(m)
    elseif subspec(m) == "ss12"
        ss12!(m)
    elseif subspec(m) == "ss13"
        ss13!(m)
    else
        error("DSGEVAR subspec $(subspec(m)) is not defined.")
    end

    return m
end

function ss1!(m::DSGEVAR)
    observables = [:obs_hours, :obs_gdpdeflator]
    lags        = 4
    λ           = .5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss2!(m::DSGEVAR)
    observables = [:obs_gdp, :obs_cpi]
    lags        = 4
    λ           = 0.5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss3!(m::DSGEVAR)
    observables = [:obs_gdp, :obs_cpi, :obs_nominalrate]
    lags        = 4
    λ           = 0.5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss10!(m::DSGEVAR)
    observables = [:obs_hours, :obs_gdpdeflator, :laborshare_t, :NominalWageGrowth]
    lags        = 4
    λ           = .5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss11!(m::DSGEVAR)
    observables = [:obs_hours, :π_t, :laborshare_t, :NominalWageGrowth]
    lags        = 4
    λ           = .5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss12!(m::DSGEVAR)
    observables = [:obs_hours, :π_t, :laborshare_t, :NominalWageGrowth, :Epi_t]
    lags        = 4
    λ           = .5
    update!(m; observables = observables, lags = lags, λ = λ)
end

function ss13!(m::DSGEVAR)
    observables = [:obs_spread, :obs_hours, :π_t, :laborshare_t, :NominalWageGrowth, :Epi_t]
    lags        = 4
    λ           = .5
    update!(m; observables = observables, lags = lags, λ = λ)
end


text:
# global constants used by more than one module

const DATA_DIR = expanduser("~/data/")
const CK_DATA_DIR_FACE_ALIGNED = DATA_DIR * "/CK/faces_aligned/"
const CK_DATA_DIR = expanduser("~/data/CK/")
const FACES_ALIGNED_SUBDIR = "faces_aligned"
const LABELES_SUBDIR = "labels"
const IMSIZE = (96, 96)
const MODEL_FILE = DATA_DIR * "facehunter/model.jld"
const NZS_FILE = DATA_DIR * "facehunter/nzs.jld"

max_stars_repo_name: dfdx/masque

max_issues_repo_name: dfdx/masque

max_forks_repo_name: dfdx/masque


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m7_rho3_59(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m7_rho3_59!(A_copy)
end

@inline function exp_ps_m7_rho3_59!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 B3 B4 B5 Ba8_6 Bb6 B6 Bb7 B7 Bb8 B8 T2k10
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing B3 with operation: mult
    mul!(memslots3,memslots2,memslots1)
    # Computing B4 with operation: mult
    mul!(memslots4,memslots3,memslots1)
    # Computing B5 with operation: mult
    mul!(memslots5,memslots4,memslots1)
    # Computing Ba8_6 = x*I+x*B5
    coeff1=0.40308943770549005
    coeff2=1.0
    memslots6 .= coeff2.*memslots5
    mul!(memslots6,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=7.647163731819816e-13
    coeff2=4.779477332387385e-14
    coeff3=2.8114572543455206e-15
    coeff4=1.5619206968586225e-16
    coeff5=8.22063524662433e-18
    coeff6=4.110317623312165e-19
    memslots7 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5
    mul!(memslots7,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots5,memslots7)
    # Deallocating Bb6 in slot 7
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B6
    coeff1=2.755731922398589e-7
    coeff2=2.505210838544172e-8
    coeff3=2.08767569878681e-9
    coeff4=1.6059043836821613e-10
    coeff5=1.1470745597729725e-11
    coeff6=1.0
    # Smart lincomb recycle B6
    memslots8 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots8
    mul!(memslots8,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots7,memslots5,memslots8)
    # Deallocating B5 in slot 5
    # Deallocating Bb7 in slot 8
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B7
    coeff1=0.008333333333333333
    coeff2=0.001388888888888889
    coeff3=0.0001984126984126984
    coeff4=2.48015873015873e-5
    coeff5=2.7557319223985893e-6
    coeff6=1.0
    # Smart lincomb recycle B7
    memslots7 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots7
    mul!(memslots7,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots5,memslots6,memslots7)
    # Deallocating Ba8_6 in slot 6
    # Deallocating Bb8 in slot 7
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B8
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=0.16666666666666666
    coeff5=0.041666666666666664
    coeff6=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B3 in slot 3
    # Deallocating B4 in slot 4
    # Deallocating B8 in slot 5
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
mutable struct AMGPreconditioner{T, TML<:AMG.MultiLevel}
    ml::TML
end
struct RugeStuben end
struct SmoothedAggregation end

for (t, f) in [(:RugeStuben, :ruge_stuben), (:SmoothedAggregation, :smoothed_aggregation)]
    @eval begin
        function AMGPreconditioner(::Type{$t}, A::AbstractMatrix)
            _A = get_data(A)
            ml = $f(_A)
            return AMGPreconditioner{$t, typeof(ml)}(ml)
        end

        function UpdatePreconditioner!(C::AMGPreconditioner{$t}, A::AbstractMatrix)
            _A = get_data(A)
            C.ml = $f(_A)
            return C
        end
    end
end

AMGPreconditioner(A::AbstractMatrix) = AMGPreconditioner(RugeStuben, A)
AMGPreconditioner{T}(A::AbstractMatrix) where T = AMGPreconditioner(T, A)

@inline function \(p::AMGPreconditioner, b)
    x = copy(b);
    return ldiv!(x, AMG.Preconditioner(p.ml), b)
end
@inline *(p::AMGPreconditioner, b) = AMG.Preconditioner(p.ml) * b
@inline ldiv!(p::AMGPreconditioner, b) = b .= p \ b
@inline function ldiv!(x, p::AMGPreconditioner, b)
    x .= b
    return ldiv!(x, AMG.Preconditioner(p.ml), b)
end
@inline mul!(b, p::AMGPreconditioner, x) = mul!(b, AMG.Preconditioner(p.ml), x)


text:
import jInv.ForwardShare.getSensMat
function getSensMat(m::Vector,pFor::SfSParam)
	if pFor.method == MATFree || pFor.method == MATBased
		error("getSensMat: Not supported or too expensive");
	end
	return pFor.Jacobian';
end

export getSensMatVec
function getSensMatVec(v::Vector,m::Vector,pFor::SfSParam)
if pFor.method == MATFree
	nshots = size(pFor.theta_phi_rad,1);
	traceLength = prod(pFor.ScreenMeshX2X3.n);
	n = pFor.ObjMesh.n;
	d = zeros(Float32,traceLength,nshots);
	v = reshape(v,tuple(n...));
	m = reshape(m,tuple(n...));
	vr = zeros(eltype(v),size(v));
	XT = zeros(0);
	XTT = zeros(0);
	for ii = 1:nshots
		(vr,XT,XTT) = rotateAndMove3D(v,pFor.theta_phi_rad[ii,:],(pFor.b[ii,:]./(pFor.ObjMesh.h)),false,vr,XT,XTT);
		# (mr,XT,XTT) = rotateAndMove3D(m,pFor.theta_phi_rad[ii,:],(pFor.b[ii,:]./(pFor.ObjMesh.h)),false,vr,XT,XTT);
		d[:,ii] = softMaxSensMatVec(pFor.SampleMatT,m[:],vr[:]);
	end
	# d = reshape(d,traceLength,nshots);
	return d;
else
	S = pFor.Jacobian;
	d = S'*v[:];
	nshots = size(pFor.theta_phi_rad,1);
	d = reshape(d,div(length(d),nshots),nshots);
end
end
export getSensTMatVec
function getSensTMatVec(v::Vector,m::Vector,pFor::SfSParam)
JtV = 0;
if pFor.method == MATFree
	error("TODO: not implemented.");
else
	S = pFor.Jacobian;
	JtV = S*v[:];
end
return JtV;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
include("Include.jl")

# Phase 1: build a default model dictionary -
path_model_file = "Coagulation.net"

# load the file -
model_buffer = read_model_file(path_model_file)

# build the default dictionary -
dd = build_default_model_dictionary(model_buffer)

# Phase 2: customize the dictionary -
SF = 1e9 # all concentrations are in nmol/L
sfa = dd["static_factors_array"]
# sfa[1] = (2.5e-9) # 1 TFPI
# sfa[2] = (3.4e-6) # 2 AT
# sfa[3] = (5e-12)  # 3 TF
# sfa[6] = (1e-9)   # 6 TRAUMA

sfa[1] = (2.5e-9) * SF        # 1 TFPI
sfa[2] = (3.4e-6) * SF        # 2 AT
sfa[3] = (5e-12) * SF         # 3 TF
sfa[6] = 0.005              # 6 TRAUMA
#sfa = SF*sfa

ℳ = dd["number_of_dynamic_states"]
xₒ = zeros(ℳ)
xₒ[1] = (1.4e-6)    # 1 FII
xₒ[2] = (1e-8)      # 2 FVII
xₒ[3] = (2e-8)      # 3 FV
xₒ[4] = (1.6e-7)    # 4 FX
xₒ[5] = (7e-10)     # 5 FVIII
xₒ[6] = (9e-8)      # 6 FIX
xₒ[7] = (1e-8)      # 7 FXI
xₒ[8] = (1e-8)      # 8 FXII
xₒ[9] = (1e-13)     # 9 FIIa
xₒ = SF * xₒ # convert to nmol

# update α -
α = dd["α"]
α[1] = 0.061
α[9] = 0.70

# update the G -
G = dd["G"]
idx = indexin(dd, "FVIIa")
G[idx, 4] = 0.1

# what is the index of TRAUMA?
idx = indexin(dd, "AT")
G[idx, 9] = 0.05

# what is the index of TFPI?
idx = indexin(dd, "TFPI")
G[idx, 1] = -0.65

# Phase 3: solve the model -
# setup the solver -
tspan = (0.0, 20.0)
prob = ODEProblem(balances, xₒ, tspan, dd; saveat = 0.01)
soln = solve(prob)

T = soln.t
X = hcat(soln.u...)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
function TidalDisplacement(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::PVCoordinatesProvider, arg4::PVCoordinatesProvider, arg5::IERSConventions, arg6::jboolean)
    return TidalDisplacement((jdouble, jdouble, jdouble, PVCoordinatesProvider, PVCoordinatesProvider, IERSConventions, jboolean), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function displacement(obj::TidalDisplacement, arg0::BodiesElements, arg1::Frame, arg2::Vector3D)
    return jcall(obj, "displacement", Vector3D, (BodiesElements, Frame, Vector3D), arg0, arg1, arg2)
end


text:
#actually this task is similar to #20
#every rectangle has only 1 North side

function great_painter(r::Robot)
    while (!isborder(r, Sud) || !isborder(r, West))
        moves!(r, West)
        moves!(r, Sud)
    end
    count = 0
    side = Ost
    while !isborder(r, Nord)
        move!(r, Nord)
        change = 1
        while change != 0 #checking the row
            change = move_full_detour(r, side)
            if isborder(r, Sud)
                count += 1
                while isborder(r, Sud) #"skipping" the barries
                    move!(r, side)
                end
            end
        end
        side = inverse(side)
    end
    println("количество прямоугольников: ", count)
end

function moves!(r::Robot,side::HorizonSide)
    num_steps = 0
    while !isborder(r, side)
        move!(r, side)
        num_steps += 1
    end
    return num_steps
end

function move_full_detour(r::Robot, side::HorizonSide)
    num_steps = 0
    while isborder(r, side) && !isborder(r, next(side))
        move!(r,next(side))
        num_steps += 1
    end
    change = 0
    if !isborder(r, side)
        move!(r, side)
        change += 1
    end
    if num_steps != 0
        while isborder(r, inverse(next(side)))
            move!(r, side)
            change += 1
        end
        for _ in 1:num_steps
            move!(r, inverse(next(side)))
        end
    end
    return change
end

next(side::HorizonSide)=(HorizonSide(mod(Int(side)+1,4)))

inverse(side::HorizonSide)=(HorizonSide(mod(Int(side)+2,4)))


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks

using ChainRulesCore
import ForwardDiff
using Interpolations: bounds

# helper to collect variable IdSet
function indiciesForValueReferences(fmu::FMU2,
                                    refs::Array{fmi2ValueReference})
    array = fmi2ValueReference[]
    for i in 1:length(fmu.modelDescription.valueReferences)
        if fmu.modelDescription.valueReferences[i] in refs
            push!(array, i)
        end
    end
    array
end

"""
Performs something similar to `fmiDoStep` for ME-FMUs (note, that fmiDoStep is for CS-FMUs only).
Event handling (state- and time-events) is supported. If you don't want events to be handled, you can disable event-handling for the NeuralFMU `nfmu` with the attribute `eventHandling = false`.

Optional, additional FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, additional FMU-values can be retrieved by keyword argument `getValueReferences`.

Function takes the current system state array ("x") and returns an array with state derivatives ("x dot") and optionally the FMU-values for `getValueReferences`.
Setting the FMU time via argument `t` is optional, if not set, the current time of the ODE solver around the NeuralFMU is used.
"""
function fmi2EvaluateME(fmu::FMU2,
        x::Array{<:Real},
        t,#::Real,
        setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
        setValues::Array{<:Real} = zeros(Real, 0),
        getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    _fmi2EvaluateME(fmu, x, t, setValueReferences, setValues, getValueReferences)
end
function fmi2EvaluateME(fmu::FMU2,
    x::Array{<:ForwardDiff.Dual{Tx, Vx, Nx}},
    t,#::Real,
    setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
    setValues::Array{<:Real} = zeros(Real, 0),
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tx, Vx, Nx}

    _fmi2EvaluateME_fd((Tx, Vx, Nx), (Tx, Vx, Nx), fmu, x, t, setValueReferences, setValues, getValueReferences)
end
function fmi2EvaluateME(fmu::FMU2,
    x::Array{<:ForwardDiff.Dual{Tx, Vx, Nx}},
    t,#::Real,
    setValueReferences::Array{fmi2ValueReference},
    setValues::Array{<:ForwardDiff.Dual{Tu, Vu, Nu}},
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tx, Vx, Nx, Tu, Vu, Nu}

    _fmi2EvaluateME_fd((Tx, Vx, Nx), (Tu, Vu, Nu), fmu, x, t, setValueReferences, setValues, getValueReferences)
end

# ForwardDiff backend using the existing ChainRulesCore.frule
# adapted from: https://discourse.julialang.org/t/chainrulescore-and-forwarddiff/61705/8
function _fmi2EvaluateME_fd(TVNx, TVNu, fmu, x, t, setValueReferences, setValues, getValueReferences)

    Tx, Vx, Nx = TVNx
    Tu, Vu, Nu = TVNu

    ȧrgs = [NoTangent(), NoTangent(), collect(ForwardDiff.partials(e) for e in x), ForwardDiff.partials(t), NoTangent(), collect(ForwardDiff.partials(e) for e in setValues), NoTangent()]
    args = [fmi2EvaluateME, fmu, collect(ForwardDiff.value(e) for e in x), ForwardDiff.value(t), setValueReferences, collect(ForwardDiff.value(e) for e in setValues), getValueReferences]

    # ToDo: Find a good fix!
    #ignore_derivatives() do @debug "From $(typeof(args[6]))"
    if typeof(args[6]) == Vector{Any}
        args[6] = convert(Vector{Float64}, args[6])
        #ignore_derivatives do @debug "To $(typeof(args[6]))"
    end

    ȧrgs = (ȧrgs...,)
    args = (args...,)

    @assert (typeof(args[3]) == Vector{Float64}) || (typeof(args[3]) == Vector{Float32}) "After conversion, `x` is still an invalid type `$(typeof(args[3]))`."

    y, _, dx, _, _, du, _ = ChainRulesCore.frule(ȧrgs, args...)

    if Vx != Float64
        Vx = Float64
    end

    if Vu != Float64
        Vu = Float64
    end

    # original function returns [dx_1, ..., dx_n, y_1, ..., y_m]
    # ToDo: Add sensitivities (partials) from du -> (dx, y)

    [collect( ForwardDiff.Dual{Tx, Vx, Nx}(y[i], dx[i]) for i in 1:length(dx) )...]
end

function _fmi2EvaluateME(fmu::FMU2,
                      x::Array{<:Real},
                      t::Real,
                      setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                      setValues::Array{<:Real} = zeros(Real, 0),
                      getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    @assert fmi2IsModelExchange(fmu) ["fmi2EvaluateME(...): As in the name, this function only supports ME-FMUs."]

    setter = (length(setValueReferences) > 0)
    getter = (length(getValueReferences) > 0)

    if setter
        @assert length(setValueReferences) == length(setValues) ["fmi2EvaluateME(...): `setValueReferences` and `setValues` need to be the same length!"]
    end

    comp = fmu.components[end]

    if setter
        fmi2SetReal(comp, setValueReferences, setValues)
    end

    fmi2SetContinuousStates(comp, x)

    # if t >= 0.0
    #     fmi2SetTime(comp, t)
    # end

    y = []
    if getter
        y = fmi2GetReal(comp, getValueReferences)
    end

    dx = fmi2GetDerivatives(comp)

    return [dx..., y...]
end

function evaluateJacobians(fmu::FMU2,
                            x::Array{<:Real},
                            t::Real = comp.t,
                            setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                            setValues::Array{<:Real} = zeros(Real, 0),
                            getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    setter = (length(setValueReferences) > 0)
    getter = (length(getValueReferences) > 0)

    comp = fmu.components[end]

    fmi2SetContinuousStates(comp, x)

    # if t >= 0.0
    #     fmi2SetTime(comp, t)
    # end

    rdx = vcat(fmu.modelDescription.derivativeValueReferences, getValueReferences)
    rx = fmu.modelDescription.stateValueReferences
    ru = setValueReferences

    if comp.jac_x != x || comp.jac_t != t || comp.jac_u != setValues

        # Jacobian ∂ ẋy / ∂ x
        if size(comp.jac_ẋy_x) != (length(rdx), length(rx))
            comp.jac_ẋy_x = zeros(length(rdx), length(rx))
        end
        comp.jacobianUpdate!(comp.jac_ẋy_x, comp, rdx, rx)

        # Jacobian ∂ ẋy / ∂ u
        if setter
            if size(comp.jac_ẋy_u) != (length(rdx), length(ru))
                comp.jac_ẋy_u = zeros(length(rdx), length(ru))
            end
            comp.jacobianUpdate!(comp.jac_ẋy_u, comp, rdx, ru)

        end

        comp.jac_u = setValues
        comp.jac_x = x
        comp.jac_t = t

    end

end

function ChainRulesCore.rrule(::typeof(fmi2EvaluateME),
                              fmu::FMU2,
                              x::Array{<:Real},
                              t::Real = comp.t,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2EvaluateME(fmu, x, t, setValueReferences, setValues, getValueReferences)

    function fmi2EvaluateME_pullback(ȳ)

        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.frule(fmi2EvaluateME, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        evaluateJacobians(fmu, x, t, setValueReferences, setValues, getValueReferences)
        comp = fmu.components[end]

        n_dx_x = @thunk(comp.jac_ẋy_x' * ȳ)
        n_dx_u = ZeroTangent()

        if setter
            n_dx_u = @thunk(comp.jac_ẋy_u' * ȳ)
        end

        f̄ = NoTangent()
        f̄mu = ZeroTangent()
        x̄ = n_dx_x
        t̄ = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return f̄, f̄mu, x̄, t̄, s̄etValueReferences, s̄etValues, ḡetValueReferences
    end
    return y, fmi2EvaluateME_pullback
end

function ChainRulesCore.frule((Δself, Δcomp, Δx, Δt, ΔsetValueReferences, ΔsetValues, ΔgetValueReferences),
                              ::typeof(fmi2EvaluateME),
                              fmu, #::FMU2Component,
                              x,#::Array{<:Real},
                              t,#::Real = comp.t,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    y = fmi2EvaluateME(fmu, x, t, setValueReferences, setValues, getValueReferences)

    function fmi2EvaluateME_pullforward(Δx, ΔsetValues)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.frule(fmi2EvaluateME, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        evaluateJacobians(fmu, x, t, setValueReferences, setValues, getValueReferences)
        comp = fmu.components[end]

        n_dx_x = comp.jac_ẋy_x * Δx
        n_dx_u = ZeroTangent()

        if setter
            n_dx_u = comp.jac_ẋy_u * ΔsetValues
        end

        f̄mu = ZeroTangent()
        x̄ = n_dx_x
        t̄ = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (f̄mu, x̄, t̄, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return (y, fmi2EvaluateME_pullforward(Δx, ΔsetValues)...)
end

"""
Performs a fmiDoStep for CS-FMUs (note, that fmiDoStep is for CS-FMUs only).

Optional, FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, FMU-values can be retrieved by keyword argument `getValueReferences`.

Function returns the FMU-values for the optional keyword argument `getValueReferences`.
The CS-FMU performs one macro step with step size `dt`. Dependent on the integrated numerical solver, the FMU may perform multiple (internal) micro steps if needed to meet solver requirements (stability/accuracy). These micro steps are hidden by FMI2.
"""
function fmi2DoStepCS(fmu::FMU2,
                      dt::Real,
                      setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                      setValues::Array{<:Real} = zeros(Real, 0),
                      getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))
    _fmi2DoStepCS(fmu, dt, setValueReferences, setValues, getValueReferences)
end
function fmi2DoStepCS(fmu::FMU2,
    dt,#::Real,
    setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
    setValues::Array{<:ForwardDiff.Dual{Tu, Vu, Nu}} = Array{ForwardDiff.Dual{Tu, Vu, Nu}}(),
    getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0) ) where {Tu, Vu, Nu}

    _fmi2DoStepCS_fd((Td, Vd, Nd), fmu, dt, setValueReferences, setValues, getValueReferences)
end

# Helper because keyword arguments are (currently) not supported by Zygote.
function _fmi2DoStepCS(fmu::FMU2,
                       dt::Real,
                       setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                       setValues::Array{<:Real} = zeros(Real, 0),
                       getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    comp = fmu.components[end]
    @assert fmi2IsCoSimulation(comp.fmu) ["fmi2DoStepCS(...): As in the name, this function only supports CS-FMUs."]
    @assert length(setValueReferences) == length(setValues) ["fmi2DoStepCS(...): `setValueReferences` ($(length(setValueReferences))) and `setValues` ($(length(setValues))) need to be the same length!"]

    if length(setValueReferences) > 0
        fmi2SetReal(comp, setValueReferences, setValues)
    end

    fmi2DoStep(fmu, dt)

    y = zeros(Float64, 0)

    if length(getValueReferences) > 0
        y = fmi2GetReal(comp, getValueReferences)
    end

    y
end

# ForwardDiff backend using the existing ChainRulesCore.frule
# adapted from: https://discourse.julialang.org/t/chainrulescore-and-forwarddiff/61705/8
function _fmi2DoStepCS_fd(TVNu,
                          fmu,
                          dt,
                          setValueReferences,
                          setValues,
                          getValueReferences)

    Tu, Vu, Nu = TVNu
    comp = fmu.components[end]

    ȧrgs = [NoTangent(), NoTangent(), ForwardDiff.partials(dt), NoTangent(), collect(ForwardDiff.partials(e) for e in setValues), NoTangent()]
    args = [fmi2DoStepCS, fmu, ForwardDiff.value(dt), setValueReferences, collect(ForwardDiff.value(e) for e in setValues), getValueReferences]

    # ToDo: Find a good fix!
    if typeof(args[5]) == Vector{Any}
        args[5] = convert(Vector{Float64}, args[5])
    end

    ȧrgs = (ȧrgs...,)
    args = (args...,)

    y, _, _, _, du, _ = ChainRulesCore.frule(ȧrgs, args...)

    if Vu != Float64
        Vu = Float64
    end

    # original function returns [y_1, ..., y_m]
    [collect( ForwardDiff.Dual{Tu, Vu, Nu}(y[i], du[i]) for i in 1:length(du) )...]
end

function ChainRulesCore.rrule(::typeof(fmi2DoStepCS),
                              fmu::FMU2,
                              dt::Real = comp.t,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    comp = fmu.components[end]
    y = fmi2DoStepCS(comp, dt, setValueReferences, setValues, getValueReferences)
    function fmi2DoStepCS_pullback(ȳ)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.rrule(fmi2DoStepCS, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        rdx = getValueReferences
        ru = setValueReferences

        n_dx_u = ZeroTangent()

        if getter
            if size(comp.jac_ẋy_u) != (length(rdx), length(ru))
                comp.jac_ẋy_u = zeros(length(rdx), length(ru))
            end
            comp.jacobianUpdate!(comp.jac_ẋy_u, comp, rdx, ru)

            n_dx_u = @thunk(comp.jac_ẋy_u' * ȳ)
        end

        f̄ = NoTangent()
        f̄mu = ZeroTangent()
        d̄t = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (f̄, f̄mu, d̄t, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return y, fmi2DoStepCS_pullback
end

function ChainRulesCore.frule((Δself, Δcomp, Δdt, ΔsetValueReferences, ΔsetValues, ΔgetValueReferences),
                              ::typeof(fmi2DoStepCS),
                              fmu, #::FMU2,
                              dt,#::Real = comp.t,
                              setValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0),
                              setValues::Array{<:Real} = zeros(Real, 0),
                              getValueReferences::Array{fmi2ValueReference} = zeros(fmi2ValueReference, 0))

    comp = fmu.components[end]
    y = fmi2DoStepCS(fmu, dt, setValueReferences, setValues, getValueReferences)
    function fmi2DoStepCS_pullforward(ΔsetValues)
        setter = (length(setValueReferences) > 0)
        getter = (length(getValueReferences) > 0)

        if setter
            @assert length(setValueReferences) == length(setValues) ["ChainRulesCore.frule(fmi2DoStepCS, ...): `setValueReferences` and `setValues` need to be the same length!"]
        end

        rdx = getValueReferences
        ru = setValueReferences

        n_dx_u = ZeroTangent()

        if getter
            if size(comp.jac_ẋy_u) != (length(rdx), length(ru))
                comp.jac_ẋy_u = zeros(length(rdx), length(ru))
            end
            comp.jacobianUpdate!(comp.jac_ẋy_u, comp, rdx, ru)

            n_dx_u = comp.jac_ẋy_u * ΔsetValues
        end

        f̄mu = ZeroTangent()
        d̄t = ZeroTangent()
        s̄etValueReferences = ZeroTangent()
        s̄etValues = n_dx_u
        ḡetValueReferences = ZeroTangent()

        return (f̄mu, d̄t, s̄etValueReferences, s̄etValues, ḡetValueReferences)
    end
    return (y, fmi2DoStepCS_pullforward(ΔsetValues)...)
end

"""
    fmi2InputDoStepCSOutput(comp::FMU2Component,
                            dt::Real,
                            u::Array{<:Real})

Sets all FMU inputs to `u`, performs a ´´´fmi2DoStep´´´ and returns all FMU outputs.
"""
function fmi2InputDoStepCSOutput(fmu::FMU2,
                                 dt::Real,
                                 u::Array{<:Real})

    @assert fmi2IsCoSimulation(fmu) ["fmi2InputDoStepCSOutput(...): As in the name, this function only supports CS-FMUs."]

    fmi2DoStepCS(fmu, dt,
                 fmu.modelDescription.inputValueReferences,
                 u,
                 fmu.modelDescription.outputValueReferences)
end

# FMU wrappers

function fmi2EvaluateME(comp::FMU2Component, args...; kwargs...)
    fmi2EvaluateME(comp.fmu, args...; kwargs...)
end

function fmi2DoStepCS(comp::FMU2Component, args...; kwargs...)
    fmi2DoStepCS(comp.fmu, args...; kwargs...)
end

function fmi2InputDoStepCSOutput(comp::FMU2Component, args...; kwargs...)
    fmi2InputDoStepCSOutput(comp.fmu, args...; kwargs...)
end


text:
import SExpressions.Parser: readsize, peek

@testset "readsize" begin
    @test readsize(stdin, 0x10) == 1
    @test readsize(stdin, 0x1000) == 2
    @test readsize(stdin, 0x10000000) == 4
    @test readsize(stdin, 0x1000000000000000) == 8
    @test readsize(stdin, 0x10000000000000000000000000000000) == 16
    @test readsize(stdin, Int8(1)) == 1
    @test readsize(stdin, Int16(1)) == 2
    @test readsize(stdin, Int32(1)) == 4
    @test readsize(stdin, Int64(1)) == 8
    @test readsize(stdin, Int128(1)) == 16
    @test readsize(stdin, Float16(0)) == 2
    @test readsize(stdin, 0.0f0) == 4
    @test readsize(stdin, 0.0) == 8
    @test readsize(stdin, 'x') == 1
    @test readsize(stdin, 'α') == 2
    @test readsize(stdin, '←') == 3
    @test readsize(stdin, '🍕') == 4
end

@testset "peek" begin
    buf = IOBuffer("Hello World")
    @test @inferred(peek(buf, UInt8)) == UInt8('H')
    @test @inferred(peek(buf, Char)) == 'H'
    @test @inferred(peek(buf, Char)) == 'H'  # test again
    @test @inferred(read(buf, Char)) == 'H'
    @test @inferred(peek(buf, Char)) == 'e'
    @test length(read(buf, String)) == 10
    @test_throws EOFError peek(buf, UInt8)
    @test_throws EOFError peek(buf, Char)
end


text:
struct Shell{M,N}
    exponents::SVector{N, Float64}
    angular_momentum::Int
    coefficients::SVector{M, SVector{N, Float64}}
end

struct Basis
    shells::Vector{Vector{Shell}}
    numbers::Vector{Int}
end

function Basis(numbers::AbstractVector{Int}, basis_dict::Dict)
    numbers = sort(unique(numbers))
    basis = Basis(Vector{Shell}[], Int[])
    for number in numbers
        shells = Shell[]
        for shell in basis_dict[string(number)]["electron_shells"]
            exponents = parse.(Float64, shell["exponents"])
            angular_momentum = shell["angular_momentum"]
            N = length(exponents)
            if length(angular_momentum) == length(shell["coefficients"])
                for i = 1:length(angular_momentum)
                    coefficients = [parse.(Float64, shell["coefficients"][i])]
                    push!(shells, Shell{1, N}(exponents, angular_momentum[i], coefficients))
                end
            elseif length(angular_momentum) == 1
                M = length(shell["coefficients"])
                coefficients = []
                for i = 1:M
                    push!(coefficients, parse.(Float64, shell["coefficients"][i]))
                end
                push!(shells, Shell{M, N}(exponents, angular_momentum[1], coefficients))
            end
        end
        push!(basis.shells, shells)
        push!(basis.numbers, number)
    end
    return basis
end

function Basis(numbers::AbstractVector{Int}, basis_str::String)
    file = joinpath(@__DIR__, "data/basis/$(basis_str).json")
    basis_dict = JSON.parsefile(file)["elements"]
    return Basis(numbers, basis_dict)
end

function Basis(nuclei::Nuclei, basis::Union{Dict,String})
    return Basis(nuclei.numbers, basis)
end


text:
function default_voigt_order(::Val{dim}) where {dim}
    dim == 1 && return @SVector[(1,1)]
    dim == 2 && return @SVector[(1,1), (2,2), (1,2), (2,1)]
    dim == 3 && return @SVector[(1,1), (2,2), (3,3), (2,3), (1,3), (1,2), (3,2), (3,1), (2,1)]
    throw(ArgumentError("no default Voigt order for `dim > 3`"))
end

"""
    tovoigt(A::Union{SecondOrderTensor, FourthOrderTensor}; [order])
    tovoigt(A::Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor}; [order, offdiagonal])

Convert a tensor to Voigt form.

Keyword arguments:
 - `offdiagscale`: Determines the scaling factor for the offdiagonal elements.
 - `order`: A vector of cartesian indices (`Tuple{Int, Int}`) determining the Voigt order.
   The default order is `[(1,1), (2,2), (3,3), (2,3), (1,3), (1,2), (3,2), (3,1), (2,1)]` for `dim=3`.

See also [`fromvoigt`](@ref).

# Examples
```jldoctest
julia> x = Mat{3,3}(1:9...)
3×3 Tensor{Tuple{3, 3}, Int64, 2, 9}:
 1  4  7
 2  5  8
 3  6  9

julia> tovoigt(x)
9-element StaticArrays.SVector{9, Int64} with indices SOneTo(9):
 1
 5
 9
 8
 7
 4
 6
 3
 2

julia> x = SymmetricSecondOrderTensor{3}(1:6...)
3×3 SymmetricSecondOrderTensor{3, Int64, 6}:
 1  2  3
 2  4  5
 3  5  6

julia> tovoigt(x; offdiagscale = 2,
                  order = [(1,1), (2,2), (3,3), (1,2), (1,3), (2,3)])
6-element StaticArrays.SVector{6, Int64} with indices SOneTo(6):
  1
  4
  6
  4
  6
 10
```
"""
function tovoigt end

@generated function tovoigt(x::SecondOrderTensor{dim};
                            order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim))) where {dim}
    L = ncomponents(x)
    exps = [:(x[order[$i]...]) for i in 1:L]
    quote
        @_propagate_inbounds_meta
        @assert length(order) == $L
        SVector{$L}($(exps...))
    end
end

@generated function tovoigt(x::FourthOrderTensor{dim};
                            order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim))) where {dim}
    L = Int(sqrt(ncomponents(x)))
    exps = [:(x[order[$i][1], order[$i][2], order[$j][1], order[$j][2]]) for i in 1:L, j in 1:L]
    quote
        @_propagate_inbounds_meta
        @assert length(order) == $L
        SMatrix{$L, $L}($(exps...))
    end
end

@generated function tovoigt(x::SymmetricSecondOrderTensor{dim, T};
                            order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim)),
                            offdiagscale::T = one(T)) where {dim, T}
    L = ncomponents(x)
    function tovoigt_element(I)
        ex = :(x[i,j])
        quote
            i, j = order[$I]
            (i == j ? $ex : $ex * offdiagscale)
        end
    end
    exps = map(tovoigt_element, 1:L)
    quote
        @_propagate_inbounds_meta
        @assert length(order) ≥ $L
        SVector{$L, T}($(exps...))
    end
end

@generated function tovoigt(x::SymmetricFourthOrderTensor{dim, T};
                            order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim)),
                            offdiagscale::T = one(T)) where {dim, T}
    L = Int(sqrt(ncomponents(x)))
    function tovoigt_element(I, J)
        ex = :(x[i,j,k,l])
        quote
            i, j = order[$I]
            k, l = order[$J]
            ((i == j && k == l) ? $ex :
             (i == j || k == l) ? $ex * offdiagscale :
                                  $ex * (offdiagscale * offdiagscale))
        end
    end
    exps = [tovoigt_element(i, j) for i in 1:L, j in 1:L]
    quote
        @_propagate_inbounds_meta
        @assert length(order) ≥ $L
        SMatrix{$L, $L, T}($(exps...))
    end
end

"""
    tomandel(A::Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor})

Convert a tensor to Mandel form which is equivalent to `tovoigt(A, offdiagscale = √2)`.

See also [`tovoigt`](@ref).
"""
@inline function tomandel(x::Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor})
    @_propagate_inbounds_meta
    tovoigt(x, offdiagscale = eltype(x)(√2))
end

"""
    fromvoigt(S::Type{<: Union{SecondOrderTensor, FourthOrderTensor}}, A::AbstractArray{T}; [order])
    fromvoigt(S::Type{<: Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor}}, A::AbstractArray{T}; [order, offdiagscale])

Converts an array `A` stored in Voigt format to a Tensor of type `S`.

Keyword arguments:
 - `offdiagscale`: Determines the scaling factor for the offdiagonal elements.
 - `order`: A vector of cartesian indices (`Tuple{Int, Int}`) determining the Voigt order.
   The default order is `[(1,1), (2,2), (3,3), (2,3), (1,3), (1,2), (3,2), (3,1), (2,1)]` for `dim=3`.

!!! note
    Since `offdiagscale` is the scaling factor for the offdiagonal elements in **Voigt form**,
    they are multiplied by `1/offdiagscale` in `fromvoigt` unlike [`tovoigt`](@ref).
    Thus `fromvoigt(tovoigt(x, offdiagscale = 2), offdiagscale = 2)` returns original `x`.

See also [`tovoigt`](@ref).

# Examples
```jldoctest
julia> fromvoigt(Mat{3,3}, 1.0:1.0:9.0)
3×3 Tensor{Tuple{3, 3}, Float64, 2, 9}:
 1.0  6.0  5.0
 9.0  2.0  4.0
 8.0  7.0  3.0

julia> fromvoigt(SymmetricSecondOrderTensor{3},
                 1.0:1.0:6.0,
                 offdiagscale = 2.0,
                 order = [(1,1), (2,2), (3,3), (1,2), (1,3), (2,3)])
3×3 SymmetricSecondOrderTensor{3, Float64, 6}:
 1.0  2.0  2.5
 2.0  2.0  3.0
 2.5  3.0  3.0
```
"""
function fromvoigt end

@generated function fromvoigt(::Type{TT},
                              v::AbstractVector;
                              order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim))) where {dim, TT <: Tensor{Tuple{dim, dim}}}
    S = Space(TT)
    L = ncomponents(S)
    inds = independent_indices(S)
    T = eltype(TT) == Any ? eltype(v) : eltype(TT)
    quote
        @_propagate_inbounds_meta
        @assert length(v) == length(order) == $L
        data = MVector{$L, $T}(undef)
        for I in 1:$L
            i, j = @inbounds order[I]
            index = $inds[i,j]
            @inbounds data[index] = v[I]
        end
        TT(Tuple(data))
    end
end

@generated function fromvoigt(::Type{TT},
                              v::AbstractMatrix;
                              order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim))) where {dim, TT <: Tensor{NTuple{4, dim}}}
    S = Space(TT)
    L = Int(sqrt(ncomponents(S)))
    inds = independent_indices(S)
    T = eltype(TT) == Any ? eltype(v) : eltype(TT)
    quote
        @_propagate_inbounds_meta
        @assert size(v) == ($L, $L)
        @assert length(order) == $L
        data = MVector{$(L*L), $T}(undef)
        for J in 1:$L, I in 1:$L
            @inbounds begin
                i, j = order[I]
                k, l = order[J]
            end
            index = $inds[i,j,k,l]
            @inbounds data[index] = v[I,J]
        end
        TT(Tuple(data))
    end
end

@generated function fromvoigt(::Type{TT},
                              v::AbstractVector{T};
                              order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim)),
                              offdiagscale::T = one(T)) where {dim, T, TT <: Tensor{Tuple{@Symmetry{dim, dim}}}}
    S = Space(TT)
    L = ncomponents(S)
    inds = independent_indices(S)
    T = eltype(TT) == Any ? eltype(v) : eltype(TT)
    quote
        @_propagate_inbounds_meta
        @assert length(v) == $L
        @assert length(order) ≥ $L
        data = MVector{$L, $T}(undef)
        for I in 1:$L
            @inbounds i, j = order[I]
            index = $inds[i,j]
            if i == j
                @inbounds data[index] = v[I]
            else
                @inbounds data[index] = v[I] / offdiagscale
            end
        end
        TT(Tuple(data))
    end
end

@generated function fromvoigt(::Type{TT},
                              v::AbstractMatrix{T};
                              order::AbstractVector{Tuple{Int, Int}} = default_voigt_order(Val(dim)),
                              offdiagscale::T = one(T)) where {dim, T, TT <: Tensor{NTuple{2, @Symmetry{dim, dim}}}}
    S = Space(TT)
    L = Int(sqrt(ncomponents(S)))
    inds = independent_indices(S)
    T = eltype(TT) == Any ? eltype(v) : eltype(TT)
    quote
        @_propagate_inbounds_meta
        @assert size(v) == ($L, $L)
        @assert length(order) ≥ $L
        data = MVector{$(L*L), $T}(undef)
        for J in 1:$L, I in 1:$L
            @inbounds begin
                i, j = order[I]
                k, l = order[J]
            end
            index = $inds[i,j,k,l]
            if i == j && k == l
                @inbounds data[index] = v[I,J]
            elseif i == j || k == l
                @inbounds data[index] = v[I,J] / offdiagscale
            else
                @inbounds data[index] = v[I,J] / (offdiagscale * offdiagscale)
            end
        end
        TT(Tuple(data))
    end
end

"""
    frommandel(S::Type{<: Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor}}, A::AbstractArray{T})

Create a tensor of type `S` from Mandel form.
This is equivalent to `fromvoigt(S, A, offdiagscale = √2)`.

See also [`fromvoigt`](@ref).
"""
@inline function frommandel(TT::Type{<: Union{SymmetricSecondOrderTensor, SymmetricFourthOrderTensor}}, v::AbstractArray{T}) where T
    @_propagate_inbounds_meta
    fromvoigt(TT, v, offdiagscale = T(√2))
end