text:
"""
Provides an API to programmatically construct a [RequireJS](https://requirejs.org/) script.
"""
module JSDependencies
using JSON

"""
    struct RemoteLibrary

Declares a remote JS dependency that should be declared in the RequireJS configuration shim.

# Fields

* `name`: a unique name for the dependency, used to refer to it in other dependencies and
  snippets
* `url`: full remote URL from where the dependency can be loaded from
* `deps`: a list of the library's dependencies (becomes the `deps` configuration in the
  RequireJS shim)
* `exports`: sets the `exports` config in the resulting RequireJS shim

# Constructors

```julia
RemoteLibrary(name::AbstractString, url::AbstractString; deps=String[], exports=nothing)
```
"""
struct RemoteLibrary
    name :: String
    url :: String
    # The following become part of the shim
    deps :: Vector{String}
    exports :: Union{Nothing, String}

    function RemoteLibrary(name::AbstractString, url::AbstractString; deps=String[], exports=nothing)
        new(name, url, deps, exports)
    end
end

"""
    struct Snippet

Declares a JS code snipped that should be loaded with RequireJS. This gets wrapped in
`require([deps...], function(args...) {script...})` in the output.

# Fields

* `deps`: names of the [`RemoteLibrary`](@ref) dependencies of the snippet
* `args`: the arguments of the callback function, corresponding to the library objects
  of the dependencies, in the order of `deps`
* `js`: the JS code of the function that gets used as the function body of the callback

# Constructors

```julia
Snippet(deps::AbstractVector, args::AbstractVector, js::AbstractString)
```
"""
struct Snippet
    deps :: Vector{String}
    args :: Vector{String}
    js :: String

    function Snippet(deps::AbstractVector, args::AbstractVector, js::AbstractString)
        new(deps, args, js)
    end
end

"""
    struct RequireJS

Declares a single RequireJS configuration/app file.

# Fields

* `libraries`: a dictionary of [`RemoteLibrary`](@ref) declarations (keys are the library
  names)
* `snippets`: a list of JS snippets ([`Snippet`](@ref))

# Constructors

```julia
RequireJS(libraries::AbstractVector{RemoteLibrary}, snippets::AbstractVector{Snippet} = Snippet[])
```

# API

* The `push!` function can be used to add additional libraries and snippets.
*
"""
struct RequireJS
    libraries :: Dict{String, RemoteLibrary}
    snippets :: Vector{Snippet}

    function RequireJS(libraries::AbstractVector, snippets::AbstractVector = Snippet[])
        all(x -> isa(x, RemoteLibrary), libraries) || throw(ArgumentError("Bad element types for `libraries`: $(typeof.(libraries))"))
        all(x -> isa(x, Snippet), snippets) || throw(ArgumentError("Bad element types for `snippets`: $(typeof.(snippets))"))
        r = new(Dict(), [])
        for library in libraries
            push!(r, library)
        end
        for snippet in snippets
            push!(r, snippet)
        end
        return r
    end
end

function Base.push!(r::RequireJS, lib::RemoteLibrary)
    if lib.name in keys(r.libraries)
        error("Library already added.")
    end
    r.libraries[lib.name] = lib
end

Base.push!(r::RequireJS, s::Snippet) = push!(r.snippets, s)

"""
    verify(r::RequireJS; verbose=false) -> Bool

Checks that none of the dependencies are missing (returns `false` if some are). If `verbose`
is set to `true`, it will also log an error with the missing dependency.
"""
function verify(r::RequireJS; verbose=false)
    isvalid = true
    for (name, lib) in r.libraries
        for dep in lib.deps
            if !(dep in keys(r.libraries))
                verbose && @error("$(dep) of $(name) missing from libraries")
                isvalid = false
            end
        end
    end
    for s in r.snippets
        for dep in s.deps
            if !(dep in keys(r.libraries))
                verbose && @error("$(dep) missing from libraries")
                isvalid = false
            end
        end
    end
    return isvalid
end

"""
    writejs(io::IO, r::RequireJS)
    writejs(filename::AbstractString, r::RequireJS)

Writes out the [`RequireJS`](@ref) object as a valid JS that can be loaded with a `<script>`
tag, either into a stream or a file. It will contain all the configuration and snippets.
"""
function writejs end

function writejs(filename::AbstractString, r::RequireJS)
    open(filename, "w") do io
        writejs(io, r)
    end
end

function writejs(io::IO, r::RequireJS)
    write(io, """
    // Generated by Documenter.jl
    requirejs.config({
      paths: {
    """)
    for (name, lib) in r.libraries
        url = endswith(lib.url, ".js") ? replace(lib.url, r"\.js$" => "") : lib.url
        write(io, """
            '$(jsescape(lib.name))': '$(jsescape(url))',
        """)
    end
    write(io, "  }")

    shim = shimdict(r)
    isempty(shim) ? write(io, '\n') : write(io, ",\n  shim: ", json_jsescape(shim, 2))
    write(io, "});\n")

    for s in r.snippets
        args = join(s.args, ", ") # Note: not string literals => no escaping
        deps = join(("\'$(jsescape(d))\'" for d in s.deps), ", ")
        write(io, """
        $("/"^80)
        require([$(deps)], function($(args)) {
        $(s.js)
        })
        """)
    end
end

function shimdict(r::RequireJS)
    shim = Dict{String,Any}()
    for (name, lib) in r.libraries
        @assert name == lib.name
        libshim = shimdict(lib)
        if libshim !== nothing
            shim[name] = libshim
        end
    end
    return shim
end

function shimdict(lib::RemoteLibrary)
    isempty(lib.deps) && (lib.exports === nothing) && return nothing
    shim = Dict{Symbol,Any}()
    if !isempty(lib.deps)
        shim[:deps] = lib.deps
    end
    if lib.exports !== nothing
        shim[:exports] = lib.exports
    end
    return shim
end

"""
    parse_snippet(filename::AbstractString) -> Snippet
    parse_snippet(io::IO) -> Snippet

Parses a JS snippet file into a [`Snippet`](@ref) object.

# Format

The first few lines are parsed to get the dependencies and argument variable names of the
snippet. They need to match `^//\\s*([a-z]+):` (i.e. start with `//`, optional whitespace, a
lowercase identifier, and a colon). Once the parser hits a line that does not match that
pattern, it will assume that it and all the following lines are the actual script.

Only lowercase letters are allowed in the identifiers. Currently only `libraries` and
`arguments` are actually parsed and lines with other syntactically valid identifiers are
ignored. For `libraries` and `arguments`, the value (after the colon) must be a comma
separated list.

A valid snippet file would look like the following. Note that the list of arguments can be
shorter than the list of dependencies.

```js
// libraries: jquery, highlight, highlight-julia, highlight-julia-repl
// arguments: \$, hljs

// Initialize the highlight.js highlighter
\$(document).ready(function() {
  hljs.initHighlighting();
})
```
"""
function parse_snippet end

parse_snippet(filename::AbstractString; kwargs...) = open(filename, "r") do io
    parse_snippet(io; kwargs...)
end

function parse_snippet(io::IO)
    libraries = String[]
    arguments = String[]
    lineno = 1
    while true
        pos = position(io)
        line = readline(io)
        m = match(r"^//\s*([a-z]+):(.*)$", line)
        if m === nothing
            seek(io, pos) # undo the last readline() call
            break
        end
        if m[1] == "libraries"
            libraries = strip.(split(m[2], ","))
            if any(s -> match(r"^[a-z-_]+$", s) === nothing, libraries)
                error("Unable to parse a library declaration '$(line)' on line $(lineno)")
            end
        elseif m[1] == "arguments"
            arguments = strip.(split(m[2], ","))
        end
        lineno += 1
    end
    snippet = String(read(io))
    Snippet(libraries, arguments, snippet)
end

"""
Replaces some of the characters in the string with escape sequences so that the strings
would be valid JS string literals, as per the
[ECMAScript® 2017 standard](https://www.ecma-international.org/ecma-262/8.0/index.html#sec-literals-string-literals).
Note that it always escapes both potential `"` and `'` closing quotes.
"""
function jsescape(s)
    b = IOBuffer()
    # From the ECMAScript® 2017 standard:
    #
    # > All code points may appear literally in a string literal except for the closing
    # > quote code points, U+005C (REVERSE SOLIDUS), U+000D (CARRIAGE RETURN), U+2028 (LINE
    # > SEPARATOR), U+2029 (PARAGRAPH SEPARATOR), and U+000A (LINE FEED).
    #
    # https://www.ecma-international.org/ecma-262/8.0/index.html#sec-literals-string-literals
    #
    # Note: in ECMAScript® 2019 (10th edition), U+2028 and U+2029 do not actually need to be
    # escaped anymore:
    #
    # > Updated syntax includes /--/ allowing U+2028 (LINE SEPARATOR) and U+2029 (PARAGRAPH
    # > SEPARATOR) in string literals to align with JSON.
    #
    # https://www.ecma-international.org/ecma-262/10.0/index.html#sec-intro
    #
    # But we'll  keep these escapes around for now, as not all JS parsers may be compatible
    # with the latest standard yet.
    for c in s
        if c === '\u000a'     # LINE FEED,       i.e. \n
            write(b, "\\n")
        elseif c === '\u000d' # CARRIAGE RETURN, i.e. \r
            write(b, "\\r")
        elseif c === '\u005c' # REVERSE SOLIDUS, i.e. \
            write(b, "\\\\")
        elseif c === '\u0022' # QUOTATION MARK,  i.e. "
            write(b, "\\\"")
        elseif c === '\u0027' # APOSTROPHE,      i.e. '
            write(b, "\\'")
        elseif c === '\u2028' # LINE SEPARATOR
            write(b, "\\u2028")
        elseif c === '\u2029' # PARAGRAPH SEPARATOR
            write(b, "\\u2029")
        else
            write(b, c)
        end
    end
    String(take!(b))
end

"""
    json_jsescape(args...)

Call `JSON.json(args...)` to generate a `String` of JSON, but then also escape two Unicode
characters to get valid JS (since [JSON is not a JS subset](http://timelessrepo.com/json-isnt-a-javascript-subset)).

!!! note
    Technically, starting with ECMAScript® 2019 (10th edition), this is no longer necessary.
    The JS standard was changed in a way that all valid JSON is also valid JS.
"""
function json_jsescape(args...)
    escapes = ('\u2028' => "\\u2028", '\u2029' => "\\u2029")
    reduce(replace, escapes, init=JSON.json(args...))
end

end


text:
/* Testing that main must return. */

/* All functions should return a value of their value type. This is not a valid Javalette program: */

int main() {
}

ext: jl

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
# We want to explore fast ways to subset a large, square N*N matrix in julia dependent on
#   1. N = number of dimensions of large matrix / data set
#   2. S = subset size = floor(log(N)) - floor(10*log(N))
#   3. Method used to extract subset of dimensions/features

function extract_with_getindex(matrix, subset)
  matrix[subset,:][:,subset]
end

function extract_with_getindex2(matrix, subset)
  matrix[subset,subset]
end

function extract_in_loop(matrix, subset)
  m = length(subset)
  res = zeros(m, m)
  for i in 1:m
    si = subset[i]
    for j in 1:m
      res[i, j] = matrix[si, subset[j]]
    end
  end
  res
end

# Do some runs to ensure everything is compiled
m = randn(10,10)
ss = sort(shuffle(collect(1:10))[1:4])
r1 = extract_in_loop(m, ss)
r2 = extract_with_getindex(m, ss)
r3 = extract_with_getindex2(m, ss)

num_features(n) = rand(int(floor(log(n))):int(floor(10*log(n))))

NumRepeats = 10
for n in [100, 1000, 10000]
  times1 = zeros(NumRepeats)
  times2 = zeros(NumRepeats)
  times3 = zeros(NumRepeats)
  for i in 1:NumRepeats
    matrix = randn(n, n)
    ss = sort(shuffle(collect(1:n))[1:num_features(n)])
    tic()
      res1 = extract_in_loop(matrix, ss)
    times1[i] = toq()
    tic()
      res2 = extract_with_getindex(matrix, ss)
    times2[i] = toq()
    tic()
      res3 = extract_with_getindex2(matrix, ss)
    times3[i] = toq()
    if res1 != res2 || res1 != res3 || res2 != res3
      println("Results differ!")
    end
  end
  println("n = $(n), extract_in_loop = $(mean(times1)), extract_with_getindex = $(mean(times2)), extract_with_getindex2 = $(mean(times3)), getindex2 factor faster = $(mean(times1)/mean(times3))")
end


text:
module Flamenco

import Sockets

include("Server.jl")
include("Writer.jl")

"""
    close_server()

Stop Flamenco server
"""
function close_server()
    Server.close()
end

"""
    start_server(
        f::Function,
        host::Union{Sockets.IPAddr, String}=Sockets.localhost,
        port::Int64=8000
    ) :: Server.Request

Start Flamenco server on specified host/port
`f` is function that parsed `request` is forwarded to
"""
function start_server(
    f::Function,
    host::Union{Sockets.IPAddr, String}=Sockets.localhost,
    port::Int64=8000
) :: Server.Request
    request::Server.Request = Server.Request("", "", Dict(), "")

    try
        request = Server.start(f, host, port)
    catch e
        if e isa InterruptException
            @warn "Server interrupted"
        else
            rethrow(e)
        end
    end

    return request
end

"""
    write_response(status_code::Int16, headers::Dict{String, String}; body::String) :: String

Return HTTP response String
"""
function write_response(status_code::Int16, headers::Dict{String, String}; body::String) :: String
    response::String = ""

    if length(body) > 0
        response = Writer.write(status_code, headers; body=body)
    else
        response = Writer.write(status_code, headers)
    end

    return response
end

end


text:
"Calculate the Hamming Distance between two DNA strands"
function distance(a, b)
    if length(a)!=length(b)
         throw(ArgumentError("diffrend strands length"))
    end
    isempty(a) || isempty(b) ? 0 : sum([a[i]!=b[i] for i in [1:length(a)...] ])
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
using Luxor, QRCode

qrc = qrcode("https://github.com/yakir12/Scallops.jl")

@pdf begin
    tiles = Tiler(250, 250, size(qrc)..., margin=0)
    setline(0.5)
    squares = first.(tiles)
    @layer begin
        for (n, pt) in enumerate(box(O, 140, 140, vertices=true))
            sethue([Luxor.julia_purple, Luxor.julia_blue,  Luxor.julia_red,  Luxor.julia_green][n])
            setopacity(0.75)
            circle(pt, 50, :fill)
        end
    end
    for x in eachindex(qrc)
        if qrc[x]
            sethue("black")
            box(squares[x], tiles.tilewidth, tiles.tileheight, :fill)
        end
    end
end 250 250 "/home/yakir/documents/formalia/Presentations/ICIV2019/betterposter-latex-template/img/qrcode"


text:
abstract type AbstractCamera end

# placeholder if no camera is present
struct EmptyCamera <: AbstractCamera end

@enum RaymarchAlgorithm begin
    IsoValue # 0
    Absorption # 1
    MaximumIntensityProjection # 2
    AbsorptionRGBA # 3
    AdditiveRGBA # 4
    IndexedAbsorptionRGBA # 5
end

include("interaction/iodevices.jl")

"""
This struct provides accessible `PriorityObservable`s to monitor the events
associated with a Scene.

Functions that act on a `PriorityObservable` must return `Consume()` if the function
consumes an event. When an event is consumed it does
not trigger other observer functions. The order in which functions are exectued
can be controlled via the `priority` keyword (default 0) in `on`.

Example:
```
on(events(scene).mousebutton, priority = Int8(20)) do event
    if is_correct_event(event)
        do_something()
        return Consume()
    end
    return
end
```

## Fields
$(TYPEDFIELDS)
"""
struct Events
    """
    The area of the window in pixels, as a `Rect2D`.
    """
    window_area::PriorityObservable{IRect2D}
    """
    The DPI resolution of the window, as a `Float64`.
    """
    window_dpi::PriorityObservable{Float64}
    """
    The state of the window (open => true, closed => false).
    """
    window_open::PriorityObservable{Bool}

    """
    Most recently triggered `MouseButtonEvent`. Contains the relevant
    `event.button` and `event.action` (press/release)

    See also [`ispressed`](@ref).
    """
    mousebutton::PriorityObservable{MouseButtonEvent}
    """
    A Set of all currently pressed mousebuttons.
    """
    mousebuttonstate::Set{Mouse.Button}
    """
    The position of the mouse as a `NTuple{2, Float64}`.
    Updates once per event poll/frame.
    """
    mouseposition::PriorityObservable{NTuple{2, Float64}} # why no Vec2?
    """
    The direction of scroll
    """
    scroll::PriorityObservable{NTuple{2, Float64}} # why no Vec2?

    """
    Most recently triggered `KeyEvent`. Contains the relevant `event.key` and
    `event.action` (press/repeat/release)

    See also [`ispressed`](@ref).
    """
    keyboardbutton::PriorityObservable{KeyEvent}
    """
    Contains all currently pressed keys.
    """
    keyboardstate::Set{Keyboard.Button}

    """
    Contains the last typed character.
    """
    unicode_input::PriorityObservable{Char}
    """
    Contains a list of filepaths to files dragged into the scene.
    """
    dropped_files::PriorityObservable{Vector{String}}
    """
    Whether the Scene window is in focus or not.
    """
    hasfocus::PriorityObservable{Bool}
    """
    Whether the mouse is inside the window or not.
    """
    entered_window::PriorityObservable{Bool}
end

function Events()
    mousebutton = PriorityObservable(MouseButtonEvent(Mouse.none, Mouse.release))
    mousebuttonstate = Set{Mouse.Button}()
    on(mousebutton, priority = typemax(Int8)) do event
        set = mousebuttonstate
        if event.action == Mouse.press
            push!(set, event.button)
        elseif event.action == Mouse.release
            delete!(set, event.button)
        else
            error("Unrecognized Keyboard action $(event.action)")
        end
        # This never consumes because it just keeps track of the state
        return Consume(false)
    end

    keyboardbutton = PriorityObservable(KeyEvent(Keyboard.unknown, Keyboard.release))
    keyboardstate = Set{Keyboard.Button}()
    on(keyboardbutton, priority = typemax(Int8)) do event
        set = keyboardstate
        if event.key != Keyboard.unknown
            if event.action == Keyboard.press
                push!(set, event.key)
            elseif event.action == Keyboard.release
                delete!(set, event.key)
            elseif event.action == Keyboard.repeat
                # set should already have the key
            else
                error("Unrecognized Keyboard action $(event.action)")
            end
        end
        # This never consumes because it just keeps track of the state
        return Consume(false)
    end

    return Events(
        PriorityObservable(IRect(0, 0, 0, 0)),
        PriorityObservable(100.0),
        PriorityObservable(false),

        mousebutton, mousebuttonstate,
        PriorityObservable((0.0, 0.0)),
        PriorityObservable((0.0, 0.0)),

        keyboardbutton, keyboardstate,

        PriorityObservable('\0'),
        PriorityObservable(String[]),
        PriorityObservable(false),
        PriorityObservable(false),
    )
end

# Compat only
function Base.getproperty(e::Events, field::Symbol)
    if field == :mousebuttons
        try
            error()
        catch ex
            bt = catch_backtrace()
            @warn(
                "`events.mousebuttons` is deprecated. Use `events.mousebutton` to " *
                "react to `MouseButtonEvent`s instead and ``."
            )
            Base.show_backtrace(stderr, bt)
            println(stderr)
        end
        mousebuttons = Node(Set{Mouse.Button}())
        on(getfield(e, :mousebutton), priority=typemax(Int8)-1) do event
            mousebuttons[] = getfield(e, :mousebuttonstate)
            return Consume(false)
        end
        return mousebuttons
    elseif field == :keyboardbuttons
        try
            error()
        catch ex
            bt = catch_backtrace()
            @warn(
                "`events.keyboardbuttons` is deprecated. Use " *
                "`events.keyboardbutton` to react to `KeyEvent`s instead."
            )
            Base.show_backtrace(stderr, bt)
            println(stderr)
        end
        keyboardbuttons = Node(Set{Keyboard.Button}())
        on(getfield(e, :keyboardbutton), priority=typemax(Int8)-1) do event
            keyboardbuttons[] = getfield(e, :keyboardstate)
            return Consume(false)
        end
        return keyboardbuttons
    elseif field == :mousedrag
        try
            error()
        catch ex
            bt = catch_backtrace()
            @warn(
                "`events.mousedrag` is deprecated. Use `events.mousebutton` or a " *
                "mouse state machine (`addmouseevents!`) instead."
            )
            Base.show_backtrace(stderr, bt)
            println(stderr)
        end
        mousedrag = Node(Mouse.notpressed)
        on(getfield(e, :mousebutton), priority=typemax(Int8)-1) do event
            if (event.action == Mouse.press) && (length(e.mousebuttonstate) == 1)
                mousedrag[] = Mouse.down
            elseif mousedrag[] in (Mouse.down, Mouse.pressed)
                mousedrag[] = Mouse.up
            end
            return Consume(false)
        end
        on(getfield(e, :mouseposition), priority=typemax(Int8)-1) do pos
            if mousedrag[] in (Mouse.down, Mouse.pressed)
                mousedrag[] = Mouse.pressed
            elseif mousedrag[] == Mouse.up
                mousedrag[] = Mouse.notpressed
            end
            return Consume(false)
        end
        return mousedrag
    else
        getfield(e, field)
    end
end

mutable struct Camera
    pixel_space::Node{Mat4f0}
    view::Node{Mat4f0}
    projection::Node{Mat4f0}
    projectionview::Node{Mat4f0}
    resolution::Node{Vec2f0}
    eyeposition::Node{Vec3f0}
    steering_nodes::Vector{ObserverFunction}
end

"""
Holds the transformations for Scenes.
## Fields
$(TYPEDFIELDS)
"""
struct Transformation <: Transformable
    parent::RefValue{Transformable}
    translation::Node{Vec3f0}
    scale::Node{Vec3f0}
    rotation::Node{Quaternionf0}
    model::Node{Mat4f0}
    flip::Node{NTuple{3, Bool}}
    align::Node{Vec2f0}
    # data conversion node, for e.g. log / log10 etc
    transform_func::Node{Any}
    function Transformation(translation, scale, rotation, model, flip, align, transform_func)
        return new(
            RefValue{Transformable}(),
            translation, scale, rotation, model, flip, align, transform_func
        )
    end
end

"""
`PlotSpec{P<:AbstractPlot}(args...; kwargs...)`

Object encoding positional arguments (`args`), a `NamedTuple` of attributes (`kwargs`)
as well as plot type `P` of a basic plot.
"""
struct PlotSpec{P<:AbstractPlot}
    args::Tuple
    kwargs::NamedTuple
    PlotSpec{P}(args...; kwargs...) where {P<:AbstractPlot} = new{P}(args, values(kwargs))
end

PlotSpec(args...; kwargs...) = PlotSpec{Combined{Any}}(args...; kwargs...)

Base.getindex(p::PlotSpec, i::Int) = getindex(p.args, i)
Base.getindex(p::PlotSpec, i::Symbol) = getproperty(p.kwargs, i)

to_plotspec(::Type{P}, args; kwargs...) where {P} =
    PlotSpec{P}(args...; kwargs...)

to_plotspec(::Type{P}, p::PlotSpec{S}; kwargs...) where {P, S} =
    PlotSpec{plottype(P, S)}(p.args...; p.kwargs..., kwargs...)

plottype(::PlotSpec{P}) where {P} = P

struct ScalarOrVector{T}
    sv::Union{T, Vector{T}}
end

Base.convert(::Type{<:ScalarOrVector}, v::AbstractVector{T}) where T = ScalarOrVector{T}(collect(v))
Base.convert(::Type{<:ScalarOrVector}, x::T) where T = ScalarOrVector{T}(x)
Base.convert(::Type{<:ScalarOrVector{T}}, x::ScalarOrVector{T}) where T = x

function collect_vector(sv::ScalarOrVector, n::Int)
    if sv.sv isa Vector
        if length(sv.sv) != n
            error("Requested collected vector with $n elements, contained vector had $(length(sv.sv)) elements.")
        end
        sv.sv
    else
        fill(sv.sv, n)
    end
end

"""
    GlyphCollection

Stores information about the glyphs in a string that had a layout calculated for them.
"""
struct GlyphCollection
    glyphs::Vector{Char}
    fonts::Vector{FTFont}
    origins::Vector{Point3f0}
    extents::Vector{FreeTypeAbstraction.FontExtent{Float32}}
    scales::ScalarOrVector{Vec2f0}
    rotations::ScalarOrVector{Quaternionf0}
    colors::ScalarOrVector{RGBAf0}
    strokecolors::ScalarOrVector{RGBAf0}
    strokewidths::ScalarOrVector{Float32}

    function GlyphCollection(glyphs, fonts, origins, extents, scales, rotations,
            colors, strokecolors, strokewidths)

        n = length(glyphs)
        @assert length(fonts)  == n
        @assert length(origins)  == n
        @assert length(extents)  == n
        @assert attr_broadcast_length(scales) in (n, 1)
        @assert attr_broadcast_length(rotations)  in (n, 1)
        @assert attr_broadcast_length(colors) in (n, 1)

        rotations = convert_attribute(rotations, key"rotation"())
        fonts = [convert_attribute(f, key"font"()) for f in fonts]
        colors = convert_attribute(colors, key"color"())
        strokecolors = convert_attribute(strokecolors, key"color"())
        strokewidths = Float32.(strokewidths)
        new(glyphs, fonts, origins, extents, scales, rotations, colors, strokecolors, strokewidths)
    end
end


text:
function fileio_save(f, x; kwargs...)
    get(kwargs, :plain, false) || error("Pass plain = true to select LigthBSON.jl over BSON.jl")
    bson_write(f.filename, x)
    nothing
end

function fileio_load(f)
    buf = read(f.filename)
    reader = BSONReader(buf)
    hastag = false
    hastype = false
    foreach(reader) do field
        hastag |= field.first == "tag"
        hastype |= field.first == "type"
        nothing
    end
    hastag && hastype && error("BSON.jl document detected, aborting LightBSON.jl load")
    reader[Any]
end


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

name = "acl"
version = v"2.2.53"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://download.savannah.nongnu.org/releases/acl/acl-2.2.53.tar.gz",
                  "06be9865c6f418d851ff4494e12406568353b891ffe1f596b34693c387af26c7")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/acl-*/
export CPPFLAGS="-I${includedir} -fPIC -DPIC"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
install_license doc/COPYING*
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if Sys.islinux(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libacl", :libacl),
    ExecutableProduct("getfacl", :getfacl),
    ExecutableProduct("setfacl", :setfacl),
    ExecutableProduct("chacl", :chacl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Attr_jll", uuid="1fd713ca-387f-5abc-8002-d8b8b1623b73"))
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenSpecFun_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenSpecFun")
JLLWrappers.@generate_main_file("OpenSpecFun", UUID("efe28fd5-8261-553b-a9e1-b2916fc3738e"))
end  # module OpenSpecFun_jll


text:
###############################################################################
#
# Canopy4RT structure
#
###############################################################################
"""
    mutable struct Canopy4RT

A canopy struct for the radiation transfer module

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct Canopy4RT{FT<:AbstractFloat}
    # local storage of dimension information
    "Number of canopy layers"
    nLayer::Int = 20

    # canopy information
    "Leaf Area Index"
    LAI       ::FT = FT(3.0 )
    "Clumping factor"
    Ω         ::FT = FT(1.0 )
    "Structure factor a"
    clump_a   ::FT = FT(1.0 )
    "Structure factor b"
    clump_b   ::FT = FT(0.0 )
    "Leaf width"
    leaf_width::FT = FT(0.1 )
    "Vegetation height"
    hc        ::FT = FT(2.0 )
    "Leaf Inclination"
    LIDFa     ::FT = FT(0.0 )
    "Variation in leaf inclination"
    LIDFb     ::FT = FT(0.0 )
    "HotSpot parameter (still need to check!)"
    hot       ::FT = FT(0.05)

    # tree/canopy/leaf traits
    "Canopy height `[m]`"
    height::FT = FT(20.0  )
    "Canopy roughness `[m]`"
    z0m   ::FT = FT(1.0   )
    "Tree roughtnes `[m]`"
    z0h   ::FT = FT(-999.0)
    "Canopy displacement height `[m]`"
    d     ::FT = FT(-999.0)
    "m/sqrt(s) turbulent transfer coefficient"
    Cd    ::FT = FT(0.01  )

    # Some more derived parameters:
    "List of mean inclination angles `[°]`"
    litab    ::Array{FT,1} = collect(FT,5:10:85)
    "List of inclination angle boundaries `[°]`"
    litab_bnd::Array{FT,2} = [collect(0:10:80) collect(FT,10:10:90)]
    "List of mean azimuth angles `[°]`"
    lazitab  ::Array{FT,1} = collect(FT,5:10:355)

    # variables used for canopy_geometry!
    "Cosine of lazitab"
    cos_ttlo ::Array{FT,1} = cosd.(lazitab)
    "Cosine of lazitab - raa (relative azimuth angle), update with time"
    cos_philo::Array{FT,1} = cosd.(lazitab .- 0)
    "Cosine of litab"
    cos_ttli ::Array{FT,1} = cosd.(litab)
    "Sine of litab"
    sin_ttli ::Array{FT,1} = sind.(litab)
    "Cache for volome scatter function"
    vol_scatt::Array{FT,1} = ones(FT, 4)

    # This is changed afterwards, ignore here.
    "Inclination angles weight distribution"
    lidf::Array{FT,1} = dladgen(LIDFa, LIDFb, litab_bnd)
    "List of level location (level = layer + 1)"
    xl  ::Array{FT,1} = collect(FT, 0:-1.0/nLayer:-1)
    "1/nLayer"
    dx  ::FT = FT(1)/nLayer

    # local storage of dimension information
    "Number of azimuth angles"
    nAzi ::Int = length(lazitab)
    "Number of inclination angles"
    nIncl::Int = length(litab)
end

max_stars_repo_name: climate-machine/Land


text:
using Test, StorageColorTypes
using ColorTypes

const ENDIAN = StorageColorTypes.BE ? "big" : "little"
@info "This script is running on a $ENDIAN-endian machine."

@test isempty(detect_ambiguities(StorageColorTypes, ColorTypes, Base, Core))

@testset "types" begin
    include("types.jl")
end

@testset "conversions" begin
    include("conversions.jl")
end


text:
# For comparison with ffr-LFDFT

push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_p( NN, AA, BB, verbose=true )

  #print("\n")

  print("Using Kronecker product directly: nabla2_x")
  @time Lx_v2 = LF.LFx.D2jl ⊗ speye(Ny) ⊗ speye(Nz)

  print("Using Kronecker product directly: nabla2_y")
  @time Ly_v2 = speye(Nx) ⊗ LF.LFy.D2jl ⊗ speye(Nz)

  @printf("Using Kronecker product directly: nabla2_z")
  @time Lz_v2 = speye(Nx) ⊗ speye(Ny) ⊗ LF.LFz.D2jl

  print("Forming Laplacian:")
  @time ∇2 = Lx_v2 + Ly_v2

  #println(Lx_v2)
  #println(Ly_v2)
  #println(Lz_v2)
  println( Lx_v2 + Ly_v2 + Lz_v2 )

end

test_main( [3, 3, 5])


text:
# free_energy.jl

"""
Implementation of Feynman's original variational technique applied to the Polaron model.

See Feynman 1955:
http://dx.doi.org/10.1103/PhysRev.97.660
"""

# Equation 31: The <|X(t) - X(s)|^{-1}> * exp(-|t-w|) effective action.
A_integrand(x, v, w) = (w^2 * x + (v^2 - w^2) / v * (1 - exp(-v * x)))^(-0.5) * exp(-x)

A(v, w, α) = π^(-0.5) * α * v * QuadGK.quadgk(x -> A_integrand(x, v, w), 0, Inf)[1]

# Equation 33: Lowest Free energy E = -B - A where B = -3/(4v)*(v-w)^2.
free_energy(v, w, α) = (3 / (4 * v)) * (v - w)^2 - A(v, w, α)

"""
    Implementation of Osaka's extended treatment of Feynman's variational technique applied to the Polaron model; generalising it from the case at O^{∘}K to the case at finite temperature.

    See Osaka 1959:
    Progress of Theoretical Physics, Vol. 22, No.3, September 1959

    And Hellwarth 1999:
    https://doi.org/10.1103/PhysRevB.60.299

"""

# Equation 62d in Hellwarth.
Y(x, v, β) = 1 / (1 - exp(-v * β)) * (1 + exp(-v * β) - exp(-v * x) - exp(v * (x - β)))

# Integrand of Equation 62c in Hellwarth.
A_integrand(x, v, w, β) =  (exp(β - x) + exp(x)) / sqrt(w^2 * x * (1 - x / β) + Y(x, v, β) * (v^2 - w^2) / v)

# Equation 62c in Hellwarth.
err = false
A(v, w, α, β) = try
    α * v / (sqrt(π) * (exp(β) - 1)) * QuadGK.quadgk(x -> A_integrand(x, v, w, β), 0, β / 2)[1]
catch
    err = true
end

if err
    A(v, w, α, β) = α * v / (sqrt(π) * (exp(BigFloat(β)) - 1)) * QuadGK.quadgk(x -> A_integrand(x, v, w, BigFloat(β)), 0, BigFloat(β / 2))[1]
end

# Equation 62b in Hellwarth. Equation 20 in Osaka.
B(v, w, β) = 3 / β * (log(v / w) - 1 / 2 * log(2 * π * β) - log(sinh(v * β / 2) / sinh(w * β / 2)))

# Equation 62e in Hellwarth. Equation 17 in Osaka.
C(v, w, β) = 3 / 4 * (v^2 - w^2) / v * (coth(v * β / 2) - 2 / (v * β))

# Equation 62a in Hellwarth. In paragraph below Equation 22 in Osaka; has extra 1/β due to different definition of A, B & C.
free_energy(v, w, α, β) = -(A(v, w, α, β) + B(v, w, β) + C(v, w, β))

"""
Multiple Phonon Branches
"""


text:
module CS1_2
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Test
using JLD2

## CS1_2  - Case-study example: Monitoring flight actuator faults
#           Local measurements of control surface angles are used
println("Case study CS1_2")

## Part 1 - Model setup
# load matrices of the aircraft multiple-model SYSACM, actuator model ACT,
# output-feedback gain K, percentages of mass variations massi

cd(joinpath(pkgdir(FaultDetectionTools), "Examples"))
@load "cs1data.jld2"

N = size(SYSACM,1)
p, m = size(SYSACM[1].D)

# build minimal realizations of AC-models with actuator faults
# set dimensions
nom = 7;             # nominal system
# set primary actuator indices
#           [ ailerons,  elevators, stabilizer, ruder ]
act_prim = [ [1,2,15,16]; [17,19];     18;       20   ];
mu = size(ACT,1); md = m-mu; mf = length(act_prim);

# form augmented aircraft model with extended measurement set
ee = eye(m); p1 = p+mf;

sysact = similar(Vector{DescriptorStateSpace},N)
for i = 1:N
    sysact[i] = gir([ee[act_prim,:]; SYSACM[i]] * append(ACT,dss(eye(md))); atol = 1.e-6)
end

sysactf = fdimodset(sysact; c = Vector(1:mu), d = mu .+ Vector(1:md), f = act_prim)

# determine closed-loop stability margin
sdegcl = -5
for i = 1:N
    global sdegcl = max(sdegcl, maximum(real(eigvals(feedback(sysact[i],K,1:mu, mf.+(1:p); negative = false).A))))
end
@test sdegcl < 0
println("sdegcl = $sdegcl")

## Part 2 - Setup of the synthesis specifications

# compute the achievable strong specifications for constant faults
S_strong = fdigenspec(sysactf[nom], atol = 1.e-6, atol3 = 1.e-6, FDtol = 1.e-5, FDGainTol = 0.001, FDfreq = 0, sdeg = -0.05);

@test size(S_strong,1) == 2^mf-1

# define SFDI, the signatures for isolation of single faults
SFDI = eye(mf) .> 0;

## Part 3 - Synthesis using Procedure EFDI
# set options for least order synthesis with EFDISYN
Q, Rf, info = efdisyn( sysactf[nom], SFDI;  atol = 1.e-6, sdeg = -5, smarg = -1, simple = false, FDfreq = 0, FDGainTol = 0.001, rdim = 1);

##  Part 4 - Assesment of synthesis results
# form extended system [ Gu Gd Gf; I 0 0] with outputs extended with u
syse = similar(Vector{DescriptorStateSpace},N)
for i = 1:N
    syse[i] = [sysactf[i].sys;eye(mu,mu+md+mf)]
end

# build overall detector and fault detection system for open-loop system
Qtot = vcat(Q.sys...)
Rftilde = vcat(Rf.sys...)
# check synthesis conditions
# open-loop checks: Q*[Gu Gd;I 0] = 0 and Q*[Gf; 0] = Rftilde
# check of achieved structure matrix
@test iszero(Qtot*syse[nom][:,1:mu+md]; atol=1.e-5) &&
      iszero(Qtot*syse[nom][:,mu+md+1:end]-Rftilde; atol=1.e-5) &&
      isequal(fditspec_(Rftilde; FDfreq = 0, atol = 1.e-5)[:,:,1],SFDI)

# form extended closed-loop system with output feedback u = K*y+v
# and evaluate [Ru Rd Rf] for the closed-loop setup

sysefb = similar(Vector{DescriptorStateSpace},N)
Rtot = similar(Vector{DescriptorStateSpace},N)

for i = 1:N
    sysefb[i] = feedback(syse[i], K, 1:mu, mf.+(1:p); negative = false)
    Rtot[i] = Qtot*sysefb[i]
end

# check robustness by computing the H-inf norms
NormRu = similar(Vector{Float64},N)
NormRd = similar(Vector{Float64},N)
NormRfmRfnom = similar(Vector{Float64},N)
for i = 1:N
    NormRu[i] = glinfnorm(Rtot[i][:,sysactf[i].controls])[1]
    NormRd[i] = glinfnorm(Rtot[i][:,sysactf[i].disturbances])[1]
    NormRfmRfnom[i] = glinfnorm(Rtot[i][:,sysactf[i].faults]-Rftilde)[1]
end

norms = [NormRu NormRd NormRfmRfnom]
println(" Table 8.3  Robustness analysis results for the nominal synthesis with position measurements")
println(" NormRu      NormRd       NormRfmRfnom")
display(norms)
@test norm(NormRu,Inf) < 1.e-4 &&  norm(NormRd,Inf) < 1.e-4 &&   norm(NormRfmRfnom,Inf) < 1.e-4
println("orders = $(order.(Q.sys))")

end  # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using QuadraticAssignmentProblem

using ProgressMeter
using DataFrames
using CSV
using Dates
using Random

using Gurobi
qap_set_optimizer!(Gurobi.Optimizer)

Ns = collect(10:10:100)
iterations = 4

problem = (
	objective = :max,
	matrix_type = :known,
	density = missing,
)

algorithm = (name = "faq", method = faq)
# (name = "opt",    method = qap_exact)
# (name = "lin",    method = qap_linearization)
# (name = "lpr",    method = qap_lprounding)
# (name = "faq",    method = faq)
# (name = "rand",   method = qap_random)
# (name = "mms",    method = qap_mms)
# (name = "ns",     method = qap_ns)
# (name = "star",   method = qap_starpacking)
# (name = "dense",  method = qap_densemapping)

function evaluate(Ns, problem, algorithm, iterations)
	evaluate_step(5, problem, algorithm)

	results = []
	@showprogress for N in Ns, i in 1:iterations
		push!(results, evaluate_step(N, problem, algorithm))
	end

	save_results(results)

	return results
end

function evaluate_step(N, problem, algorithm)
	A, B, opt = generate_example(N, problem)
	uid = rand(UInt)

	alg = algorithm.method

	time = @elapsed P, _ = alg(A, B, problem.objective)
	obj = qap_objective(A, B, P)

	result = merge(problem, (;N, uid, method=algorithm.name, time, obj, ratio=obj/opt))
	return result
end

function generate_example(N, problem)
	opt = missing
	if problem.matrix_type == :uniform
		A = uniform_matrix(N)
		B = uniform_matrix(N)
	elseif problem.matrix_type == :metric
		A = metric_matrix(N)
		B = metric_matrix(N)
	elseif problem.matrix_type == :undirected_graph
		A = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=true)
		B = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=true)
	elseif problem.matrix_type == :directed_graph
		A = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=false)
		B = zeroone_matrix(N, p=problem.density, selfedges=false, symmetric=false)
	elseif problem.matrix_type == :known
		A, B, perm = generate_qap(N)
		opt = qap_objective(A, B, perm)
	else
		error("Invalid matrix type: $(problem.matrix_type)")
	end

	return A, B, opt
end

function save_results(results)

	d = Dates.format(now(), "yyyymmddHHMM")
	d = isdir("results/$d/") ? d * "_" * randstring(8) : d

	dir = joinpath("./results/", d)
	mkpath(dir)

	metric_cols = [:time, :obj, :ratio]
	problem_cols = setdiff(keys(results[1]), metric_cols, [:uid])

	results_df_long = DataFrame(results)
	results_df_long |> CSV.write(joinpath(dir, "results_long.csv"))

	mean(xs) = sum(xs) / length(xs)
	results_mean_df = combine(groupby(results_df_long, problem_cols), [m => mean => m for m in metric_cols]...)
	results_mean_df |> CSV.write(joinpath(dir, "results_mean.csv"))

	return
end

if abspath(PROGRAM_FILE) == @__FILE__
	evaluate(Ns, problem, algorithm, iterations)
end


text:
function foo_reorder(args...; node::Int64 = Int[], n::Int64 = Int[], m::Int64 = Int[],
    e1::Array{Int64, 1} = Int[], e2::Array{Int64, 1} = Int[],
    neworder::Array{Int64, 1} = Int[], L::Array{Int64, 1} = Int[], pos::Array{Int64, 1} = Int[])
    i  = node - n
    for j in 1:pos[i]
        k = L[i + m * j]
        neworder[iii] = k + 1
        if e2[k] > n
            foo_reorder(node = e2[k], n = n, m = m, e1 = e1, e2 = e2, neworder = neworder, L = L, pos = pos)
        end
        iii = iii + 1
    end
end


text:
module SpectatorV4 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

struct FeaturedGames <: AbstractDTO
    clientRefreshInterval
    gameList
    FeaturedGames(clientRefreshInterval, gameList) = new(clientRefreshInterval, gameList)
end

struct CurrentGameInfo <: AbstractDTO
    gameId
    gameStartTime
    platformId
    gameMode
    mapId
    gameType
    bannedChampions
    observers
    participants
    gameLength
    gameQueueConfigId
    CurrentGameInfo(gameId, gameStartTime, platformId, gameMode, mapId, gameType, bannedChampions, observers, participants, gameLength, gameQueueConfigId) = new(gameId, gameStartTime, platformId, gameMode, mapId, gameType, bannedChampions, observers, participants, gameLength, gameQueueConfigId)
end

"""
         featured_games(api_key::String,
                        platform::String ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::FeaturedGames
"""
function featured_games(api_key::String,
                        platform::String ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::FeaturedGames
    call_api(FeaturedGames, api_key, action, endpoint, "/lol/spectator/v4/featured-games", event, featured_games)
end

"""
         active_games(api_key::String,
                      platform::String,
                      encryptedSummonerId::String ;
                      endpoint::HTTP.URI = lol_api_server(platform),
                      action::Function = http_action,
                      event::Function = nothing_in_event)::CurrentGameInfo
"""
function active_games(api_key::String,
                      platform::String,
                      encryptedSummonerId::String ;
                      endpoint::HTTP.URI = lol_api_server(platform),
                      action::Function = http_action,
                      event::Function = nothing_in_event)::CurrentGameInfo
    call_api(CurrentGameInfo, api_key, action, endpoint, "/lol/spectator/v4/active-games/by-summoner/$encryptedSummonerId", event, active_games)
end

end # module LOLTools.SpectatorV4


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module ForwardNets

using LightGraphs

export
    ForwardNet,
    Node,
    Layer,
    Variable,
    Activation,
    Affine,
    LSTM,
    GRU,
    Conv2d,
    BatchNorm1,
    BatchNorm3,
    Concatenator,
    Reshaper,
    ReLU,
    ELU,
    TanH,
    SoftPlus,

    ForwardPass,

    read_binary_vec,
    infer_shape,
    convert_to_column_major_array,
    convert_to_column_major_array!,

    forward!,
    add_node!,
    restore!,
    indexof,
    lastindex,
    name,
    output,

    sigmoid,
    relu,
    elu,
    softplus,

    zero!,  # to zero out the hidden state of an LSTM

    calc_forwardpass,
    print_forward_pass

pkgdir = dirname(@__FILE__)

include(joinpath(pkgdir, "read_binary_vec.jl"))
include(joinpath(pkgdir, "utils.jl"))
include(joinpath(pkgdir, "node.jl"))
include(joinpath(pkgdir, "forwardnet.jl"))
include(joinpath(pkgdir, "variable.jl"))
include(joinpath(pkgdir, "layers.jl"))
include(joinpath(pkgdir, "recurrent.jl"))
include(joinpath(pkgdir, "conv.jl"))
include(joinpath(pkgdir, "batchnorm.jl"))
include(joinpath(pkgdir, "activations.jl"))
include(joinpath(pkgdir, "forwardpass.jl"))

include(joinpath(pkgdir, "deprecated.jl"))

end # module


text:
using MacroTools: MacroTools, isexpr, isline, @capture, @q
using Setfield: @set

struct VarInfo{S<:Union{Symbol,Nothing}}
    system::Symbol
    name::Symbol
    alias::Union{Symbol,Nothing}
    args::Vector
    kwargs::Vector
    body#::Union{Expr,Symbol,Nothing}
    state::S
    type::Union{Symbol,Expr,Nothing}
    tags::Dict{Symbol,Any}
    line::Expr
    linenumber::LineNumberNode
    docstring::String
end

Base.show(io::IO, v::VarInfo) = begin
    println(io, "system: $(v.system)")
    println(io, "name: $(v.name)")
    println(io, "alias: $(v.alias)")
    println(io, "func ($(repr(v.args)); $(repr(v.kwargs))) = $(repr(v.body))")
    println(io, "state: $(repr(v.state))")
    println(io, "type: $(repr(v.type))")
    for (a, b) in v.tags
        println(io, "tag $a = $(repr(b))")
    end
    println(io, "line: $(v.line)")
    println(io, "linenumber: $(v.linenumber)")
    println(io, "docstring: $(v.docstring)")
end

VarInfo(system::Symbol, line::Expr, linenumber::LineNumberNode, docstring::String, scope::Module, substs::Dict) = begin
    # name[(args..; kwargs..)][: alias] [=> body] [~ [state][::stype|<:dtype][(tags..)]]
    @capture(bindscope(line, scope), (decl_ ~ deco_) | decl_)
    @capture(deco,
        (state_::stype_(tags__)) | (::stype_(tags__)) | (state_::stype_) | (::stype_) |
        (state_<:dtype_(tags__)) | (<:dtype_(tags__)) | (state_<:dtype_) | (<:dtype_) |
        state_(tags__) | state_
    )
    @capture(decl, (def1_ => body_) | def1_)
    @capture(def1, (def2_: alias_) | def2_)
    @capture(def2, name_(args__; kwargs__) | name_(; kwargs__) | name_(args__) | name_)
    name = parsename(name, system)
    alias = parsealias(alias, system)
    args = parseargs(args, system)
    kwargs = parsekwargs(kwargs)
    body = parsebody(body)
    state = parsestate(state)
    type = parsetype(stype, dtype, state, scope, substs)
    tags = parsetags(tags; name, alias, args, kwargs, system, state, type)
    try
        VarInfo{typeof(state)}(system, name, alias, args, kwargs, body, state, type, tags, line, linenumber, docstring)
    catch
        error("unrecognized variable declaration: $line")
    end
end

#HACK: experimental support for scope placeholder `:$`
bindscope(l, s::Module) =  MacroTools.postwalk(x -> @capture(x, :$) ? nameof(s) : x, l)

parsename(name, system) = canonicalname(name, system)
parsealias(alias, system) = parsename(alias, system)
canonicalname(n::Symbol, s) = Symbol(isprivatename(n) ? "__$(s)_$(n)" : n)
canonicalname(n::Symbol, ::Union{S,Type{S}}) where {S<:System} = canonicalname(n, namefor(S))
canonicalname(n, _) = n
splitcanonicalname(n) = begin
    m = match(r"__(.+)__(.+)", string(n))
    isnothing(m) ? (nothing, n) : Symbol.((m[1], '_' * m[2]))
end
uncanonicalname(n) = splitcanonicalname(n)[2]
uncanonicalname(n, s) = begin
    cs, cn = splitcanonicalname(n)
    (s == cs) ? cn : Symbol(n)
end
canonicalnamelabel(n) = begin
    cs, cn = splitcanonicalname(n)
    isnothing(cs) ? string(cn) : "($cs) $cn"
end
isprivatename(n) = begin
    s = string(n)
    #HACK: support private variable name with single prefix `_` (i.e. _a => __S__a, __b => __b)
    startswith(s, "_") && !startswith(s, "__")
end
privatename(n) = Symbol(isprivatename(n) ? string(n)[2:end] : n)
bindname(ex, s) = MacroTools.postwalk(x -> x isa Symbol ? canonicalname(x, s) : x, ex)

#TODO: prefixscope to args/kwargs type specifier?
parseargs(args, system) = parsearg.(args, system)
#HACK: support private variable name in the dependency list
parsearg(a::Symbol, system) = Expr(:kw, privatename(a), canonicalname(a, system))
parsearg(a::Expr, system) = @capture(a, k_=v_) ? Expr(:kw, k, canonicalname(v, system)) : a
parseargs(::Nothing, _) = []

parsekwargs(kwargs) = kwargs
parsekwargs(::Nothing) = []

parsebody(body) = begin
    #HACK: disable return checking for now, too aggressive for local scope return
    #TODO: translate `return` to a local safe statement
    #MacroTools.postwalk(x -> @capture(x, return(_)) ? error("`return` is not allowed: $body") : x, body)
    body
end
parsebody(::Nothing) = nothing

parsestate(state) = typestate(Val(state))
typestate(::Val{S}) where {S} = Symbol(uppercasefirst(string(S)))
typestate(::Val{nothing}) = nothing

parsetype(::Nothing, ::Nothing, state, _, _) = typetag(Val(state))
parsetype(stype, ::Nothing, _, args...) = parsetype(Val(:static), stype, args...)
parsetype(::Nothing, dtype, _, args...) = parsetype(Val(:dynamic), dtype, args...)
parsetype(trait, type, scope, substs) = begin
    if haskey(substs, type)
        s = substs[type]
        type = s.type
        scope = s.scope
        #HACK: force static type when substituion takes place
        trait = Val(:static)
    end
    parsetype(type, scope, trait)
end
parsetype(type, scope, trait) = begin
    T = if @capture(type, elemtype_[])
        :(Vector{$(gentype(elemtype, scope, trait))})
    else
        gentype(type, scope, trait)
    end
end
gentype(type, scope, trait) = genactualtype(genscopedtype(type, scope), trait)
genscopedtype(type, scope) = begin
    l = Symbol[]
    add(t::Symbol) = push!(l, t)
    add(t) = nothing
    isexpr(type, :braces) && MacroTools.postwalk(type) do ex
        @capture(ex, $:|(T__)) && add.(T)
        ex
    end
    conv(t) = prefixscope(parsetypealias(t), scope)
    isempty(l) ? conv(type) : :(Union{$(conv.(l)...)})
end
genactualtype(type, ::Val{:static}) = :($C.typefor($type))
genactualtype(type, ::Val{:dynamic}) = type

parsetypealias(type::Symbol) = parsetypealias(Val(type), type)
parsetypealias(type) = type
parsetypealias(::Val{:int}, _) = :Int64
parsetypealias(::Val{:uint}, _) = :UInt64
parsetypealias(::Val{:float}, _) = :Float64
parsetypealias(::Val{:bool}, _) = :Bool
parsetypealias(::Val{:sym}, _) = :Symbol
parsetypealias(::Val{:str}, _) = :String
parsetypealias(::Val{:∅}, _) = :Nothing
parsetypealias(::Val{:_}, _) = :Missing
parsetypealias(::Val{:date}, _) = :($C.Dates.Date)
parsetypealias(::Val{:datetime}, _) = :($C.TimeZones.ZonedDateTime)
parsetypealias(::Val, type) = type

using Dates: Date
using TimeZones: ZonedDateTime

extractscope(x) = begin
    l = []
    if @capture(x, a_{c__})
        x = a
    end
    while true
        if @capture(x, a_.b_)
            push!(l, b)
            if isexpr(a)
                x = a
            else
                push!(l, a)
                break
            end
        else
            push!(l, x)
            break
        end
    end
    (; l=reverse(l), c)
end
genscope(lc) = genscope(lc.l, lc.c)
genscope(l, ::Nothing) = reduce((a, b) -> Expr(:., a, QuoteNode(b)), l)
genscope(l, c) = :($(genscope(l, nothing)){$(c...)})
prefixscope(x, p::Module) = prefixscope(x, nameof(p))
prefixscope(x, p::Symbol) = let (l, c) = extractscope(x)
    genscope([p, l...], c)
end

getsystem(m::Module, i) = reduce((a, b) -> getfield(a, b), split(String(i), ".") .|> Symbol, init=m)

parsetags(::Nothing; a...) = parsetags([]; a...)
parsetags(tags::Vector; system, state, type, a...) = begin
    s = Val(state)
    d = Dict{Symbol,Any}()
    for t in tags
        if @capture(t, k_=v_)
            d[k] = bindname(v, system)
        elseif @capture(t, @u_str(v_))
            d[:unit] = @q @u_str($v)
        else
            d[t] = true
        end
    end
    !haskey(d, :unit) && (d[:unit] = nothing)
    d[:_type] = esc(type)
    updatetags!(d, s; a...)
    d
end

typetag(::Val) = :Float64
typetag(::Val{:Advance}) = :Float64 #HACK: avoid unexpected promotion (i.e. Rational) when using :Int
typetag(::Val{:Flag}) = :Bool
typetag(::Val{:Produce}) = :(Vector{System})
typetag(::Val{:Provide}) = :(Cropbox.DataFrame)
typetag(::Val{nothing}) = nothing

supportedtags(::Val) = nothing

filterconstructortags(v::VarInfo) = begin
    stags = constructortags(Val(v.state))
    filter(v.tags) do (k, v)
        isnothing(stags) ? true : k in stags ||
        startswith(String(k), "_")
    end
end

updatetags!(d, ::Val; _...) = nothing

istag(v::VarInfo, t) = get(v.tags, t, false)
istag(v::VarInfo, t::Tuple) = any(istag.(Ref(v), t))
istag(v::VarInfo, t...) = istag(v, t)
istag(t) = Base.Fix2(istag, t)
hastag(v::VarInfo, t) = haskey(v.tags, t)
gettag(v::VarInfo, t, d=nothing) = get(v.tags, t, d)

Base.names(v::VarInfo) = let n = v.name, a = v.alias
    isnothing(a) ? [n] : [n, a]
end
linenumber(v::VarInfo, prefix="", postfix="") = begin
    n = v.linenumber
    @set n.file = Symbol(n.file, ":$prefix|", v.name, "|$postfix")
end

####

abstract type VarStep end
struct PreStep <: VarStep end
struct MainStep <: VarStep end
struct PostStep <: VarStep end

Base.print(io::IO, ::PreStep) = print(io, "∘")
Base.print(io::IO, ::MainStep) = print(io, "")
Base.print(io::IO, ::PostStep) = print(io, "⋆")

struct Node{I,S}
    info::I
    step::S
end

const VarNode = Node{VarInfo,VarStep}

prev(n::VarNode) = begin
    if n.step == MainStep()
        VarNode(n.info, PreStep())
    elseif n.step == PostStep()
        VarNode(n.info, MainStep())
    elseif n.step == PreStep()
        error("Pre-step node can't have a previous node: $n")
    end
end

####

const C = :Cropbox
const ϵ = @q begin end

genvartype(v::VarInfo{Nothing}) = gettag(v, :_type)
genvartype(v::VarInfo{Symbol}) = begin
    N = gettag(v, :_type)
    U = gettag(v, :unit)
    V = @q $C.valuetype($N, $U)
    T = genvartype(v, Val(v.state); N, U, V)
    istag(v, :ref) ? @q(StateRef{$V,$T}) : T
end

genfield(v::VarInfo) = begin
    type = genvartype(v)
    name = symname(v)
    docstring = isempty(v.docstring) ? ϵ : v.docstring
    alias = v.alias
    @q begin
        $docstring
        $name::$type
        $(isnothing(alias) ? ϵ : docstring)
        $(isnothing(alias) ? ϵ : :($alias::$type))
    end
end
genfields(infos) = [genfield(v) for v in infos]

genpredecl(name) = @q _names = $C.names.($C.mixincollect($(esc(name)))) |> reverse |> Iterators.flatten |> collect
gennewargs(infos) = names.(infos) |> Iterators.flatten |> collect

genoverride(v::VarInfo) = begin
    !isnothing(v.body) && error("`override` can't have funtion body: $(v.body)")
    gengetkwargs(v, nothing)
end

genextern(v::VarInfo, default) = gengetkwargs(v, default)

gengetkwargs(v::VarInfo, default) = begin
    K = [Meta.quot(n) for n in names(v)]
    K = names(v)
    @q $C.getbynames(__kwargs__, $K, $default)
end

getbynames(d, K, default=nothing) = begin
    for k in K
        haskey(d, k) && return d[k]
    end
    default
end

using DataStructures: OrderedSet
gendecl(N::Vector{VarNode}) = gendecl.(OrderedSet([n.info for n in N]))
gendecl(v::VarInfo) = begin
    name = Meta.quot(v.name)
    alias = Meta.quot(v.alias)
    decl = if istag(v, :override)
        genoverride(v)
    else
        value = istag(v, :extern) ? genextern(v, gendefault(v)) : gendefault(v)
        stargs = [:($(esc(k))=$v) for (k, v) in filterconstructortags(v)]
        @q $C.$(v.state)(; _name=$name, _alias=$alias, _value=$value, $(stargs...))
    end
    decl = istag(v, :ref) ? @q(StateRef($decl)) : decl
    gendecl(v, decl)
end
gendecl(v::VarInfo{Nothing}) = begin
    args = emitfuncargpair.(v.args; value=false)
    if istag(v, :option)
        push!(args, @q $(esc(:option)) = __kwargs__)
    end
    decl = if istag(v, :override)
        genoverride(v)
    elseif isnothing(v.body)
        @q $(esc(v.type))(; $(args...))
    else
        @q let $(args...); $(esc(v.body)) end
    end
    # implicit :expose
    decl = :($(esc(v.name)) = $decl)
    gendecl(v, decl)
end
gendecl(v::VarInfo, decl) = @q begin
    $(linenumber(v, "gendecl"))
    $(v.name) = $decl
    $(isnothing(v.alias) ? ϵ : :($(v.alias) = $(v.name)))
end

#HACK: @capture doesn't seem to support GlobalRef
const DOCREF = GlobalRef(Core, Symbol("@doc"))
isdoc(ex) = isexpr(ex, :macrocall) && ex.args[1] == DOCREF
gensource(v::VarInfo, name) = begin
    v.system != name && return ϵ
    if isempty(v.docstring)
        @q begin $(v.linenumber); $(v.line) end
    else
        Expr(:macrocall, DOCREF, v.linenumber, v.docstring, v.line)
    end
end
gensource(infos, name) = MacroTools.flatten(@q begin $(gensource.(infos, name)...) end)

updateconsts(consts, infos) = begin
    #HACK: keep type of Context in case needed for field construction (i.e. timeunit for Accumulate)
    context = filter(v -> v.name == :context, infos) |> only
    merge(consts, Dict(:__Context__ => context.type))
end
genconstpatches(consts, scope, incl) = begin
    @gensym CS
    consts0 = merge(constsof.(_mixincollect(mixinsof(scope, incl)))...)
    K = [keys(consts0)..., keys(consts)...]
    [
        :($CS = $(genconstbase(:(mixinsof($scope, $incl)), consts))),
        (:($k = $CS[$(Meta.quot(k))]) for k in K)...
    ]
end
genconstbase(M, consts) = :(merge(constsof.(_mixincollect($M))..., Dict($((@q($(Meta.quot(p[1])) => $(esc(p[2]))) for p in consts)...))))

genfieldnamesunique(infos) = Tuple(v.name for v in infos)
genfieldnamesalias(infos) = Tuple((v.name, v.alias) for v in infos)

genstruct(name, type, infos, consts, substs, incl, scope) = begin
    #FIXME: simple `gensym(name)` leads to very strange invalid redefinition of constant error on Julia 1.6.0
    _S = esc(gensym(Symbol(:_, name)))
    S = esc(name)
    T = esc(type)
    N = Meta.quot(name)
    nodes = sortednodes(infos)
    #HACK: field declarations inside block doesn't work as expected
    #fields = genfields(infos)
    fields = MacroTools.flatten(@q begin $(genfields(infos)...) end).args
    predecl = genpredecl(name)
    decls = gendecl(nodes)
    args = gennewargs(infos)
    source = gensource(infos, name)
    consts = updateconsts(consts, infos)
    constpatches = genconstpatches(consts, scope, incl)

    system = quote
        Core.@__doc__ abstract type $S <: $T end
        $S(; kw...) = $_S(; kw...)
        $C.typefor(::Type{<:$S}) = $_S
        let $(constpatches...)
            Core.@__doc__ mutable struct $_S <: $S
                $(fields...)
                function $_S(; __kwargs__...)
                    $predecl
                    $(decls...)
                    new($(args...))
                end
            end
        end
        $C.namefor(::Type{$_S}) = $C.namefor($S)
        $C.typefor(::Type{$_S}) = $_S
        $C.source(::Type{$_S}) = $(Meta.quot(source))
        $C.constsof(::Type{$_S}) = $(genconstbase(:(mixinsof($S)), consts))
        $C.substsof(::Type{$_S}) = $substs
        $C.mixinsof(::Type{$_S}) = $(mixinsof(scope, incl))
        $C.fieldnamesunique(::Type{$_S}) = $(genfieldnamesunique(infos))
        $C.fieldnamesalias(::Type{$_S}) = $(genfieldnamesalias(infos))
        $C.scopeof(::Type{$_S}) = $scope
        $C._update!($(esc(:self))::$_S, ::$C.MainStage) = $(genupdate(nodes, MainStage(); scope))
        $C._update!($(esc(:self))::$_S, ::$C.PreStage) = $(genupdate(infos, PreStage(); scope))
        $C._update!($(esc(:self))::$_S, ::$C.PostStage) = $(genupdate(infos, PostStage(); scope))
        Base.Docs.getdoc(::Type{$S}) = $C.getdoc($_S)
        Base.Docs.getdoc(::Type{$_S}) = $C.getdoc($_S)
        $S
    end
    system #|> MacroTools.flatten
end

source(s::S) where {S<:System} = source(S)
source(::Type{S}) where {S<:System} = source(typefor(S))
source(::Type{System}) = quote
    context ~ ::Context(override)
    config(context) => context.config ~ ::Config
end
source(::Type) = :()

constsof(s::S) where {S<:System} = constsof(S)
constsof(::Type{S}) where {S<:System} = constsof(typefor(S))
constsof(::Type{System}) = Dict()

substsof(s::S) where {S<:System} = substsof(S)
substsof(::Type{S}) where {S<:System} = substsof(typefor(S))
substsof(::Type{System}) = Dict()

mixinsof(s::S) where {S<:System} = mixinsof(S)
mixinsof(::Type{S}) where {S<:System} = mixinsof(typefor(S))
mixinsof(::Type{System}) = (System,)
mixinsof(::Type) = ()
mixinsof(scope::Module, incl) = Tuple(getsystem.(Ref(scope), incl))

using DataStructures: OrderedSet
mixincollect(s::S) where {S<:System} = mixincollect(S)
mixincollect(S::Type{<:System}, l=OrderedSet()) = begin
    S in l && return l
    push!(l, S)
    _mixincollect(mixinsof(S), l)
    #HACK: ensure mixins come before composite system
    #TODO: need testsets for mixins/mixincollect
    push!(delete!(l, S), S)
    l
end
_mixincollect(M::Tuple, l=OrderedSet()) = begin
    for m in M
        union!(l, mixincollect(m, l))
    end
    l
end
mixincollect(s) = ()

mixinof(s, SS::Type{<:System}...) = begin
    M = mixincollect(s)
    for S in SS
        for m in M
            m <: S && return S
        end
    end
    nothing
end

mixindispatch(s, SS::Type{<:System}...) = begin
    m = mixinof(s, SS...)
    n = isnothing(m) ? m : namefor(m)
    (s, Val(n))
end

typefor(s::Symbol, m::Module=Main) = getsystem(m, s) |> typefor
typefor(T) = T
vartype(::Type{S}, k) where {S<:System} = fieldtype(typefor(S), k) |> typefor

fieldnamesunique(::Type{System}) = ()
fieldnamesalias(::Type{System}) = ()
fieldnamesunique(::S) where {S<:System} = fieldnamesunique(S)
fieldnamesalias(::S) where {S<:System} = fieldnamesalias(S)
fieldnamesunique(::Type{S}) where {S<:System} = fieldnamesunique(typefor(S))
fieldnamesalias(::Type{S}) where {S<:System} = fieldnamesalias(typefor(S))

subsystemsof(::Type{System}) = ()
subsystemsof(::S) where {S<:System} = subsystemsof(S)
subsystemsof(::Type{S}) where {S<:System} = filter(n -> fieldtype(typefor(S), n) <: System, fieldnamesunique(S))

scopeof(::Type{System}) = @__MODULE__
scopeof(::S) where {S<:System} = scopeof(S)
scopeof(::Type{S}) where {S<:System} = scopeof(typefor(S))

abstract type UpdateStage end
struct PreStage <: UpdateStage end
struct MainStage <: UpdateStage end
struct PostStage <: UpdateStage end

Base.print(io::IO, ::PreStage) = print(io, "†")
Base.print(io::IO, ::MainStage) = print(io, "")
Base.print(io::IO, ::PostStage) = print(io, "‡")

_update!(S::Vector{<:System}, t) = begin
    for s in S
        update!(s, t)
    end
end
_update!(s, t) = s
update!(s, t=MainStage()) = _update!(s, t)

parsehead(head; scope) = begin
    # @system name[{patches..}][(mixins..)] [<: type]
    @capture(head, (decl_ <: type_) | decl_)
    @capture(decl, name_{patches__}(mixins__) | name_{patches__} | name_(mixins__) | name_)
    type = isnothing(type) ? :System : type
    patches = isnothing(patches) ? [] : patches
    mixins = isnothing(mixins) ? [] : mixins
    consts, substs = parsepatches(patches, scope)
    incl = [:System]
    for m in mixins
        push!(incl, Symbol(m))
    end
    (; name, consts, substs, incl, type)
end

parsepatches(patches::Vector, scope::Module) = begin
    consts = Dict()
    substs = Dict()
    for p in patches
        if @capture(p, o_ => n_)
            substs[o] = (; type=n, scope)
        elseif @capture(p, o_ = n_)
            consts[o] = n
        else
            error("unsupported patch format: $p")
        end
    end
    (consts, substs)
end

using DataStructures: OrderedDict, OrderedSet
gensystem(body; name, consts, substs, incl, type, scope, _...) = genstruct(name, type, geninfos(body; name, substs, incl, scope), consts, substs, incl, scope)
geninfos(body; name, substs, incl, scope, _...) = begin
    # update type substitution dict cascaded over all participating mixins
    substs = merge(substsof.(_mixincollect(mixinsof(scope, incl)))..., substs)
    con(b, s, sc) = begin
        d = OrderedDict{Symbol,VarInfo}()
        #HACK: default in case LineNumberNode is not attached
        ln = LineNumberNode(@__LINE__, @__FILE__)
        for l in b.args
            isline(l) && (ln = l; continue)
            if isdoc(l)
                lnn, ds, l = l.args[2:4]
                isline(lnn) && (ln = lnn)
            else
                ds = ""
            end
            v = VarInfo(s, l, ln, ds, sc, substs)
            n = v.name
            v0 = get(d, n, nothing)
            if !isnothing(v0)
                @warn "duplicate variable" system=s name=v.name alias=(v0.alias => v.alias) state=(v0.state => v.state)
            end
            d[n] = v
        end
        d
    end
    add!(d, b, s, sc) = begin
        for (n, v) in con(b, s, sc)
            if haskey(d, n)
                v0 = d[n]
                if v.state == :Hold
                    continue
                # support simple body replacement (i.e. `a => 1` without `~ ...` part)
                elseif isnothing(v.state) && isnothing(v.type)
                    v = @set v0.body = v.body
                elseif v0.alias != v.alias && v0.state != :Hold
                    @warn "variable replaced with inconsistent alias" name=v.name system=(v0.system => v.system) alias=(v0.alias => v.alias)
                end
            end
            d[n] = v
        end
    end
    combine() = begin
        d = OrderedDict{Symbol,VarInfo}()
        for i in incl
            S = getsystem(scope, i)
            for m in mixincollect(S)
                add!(d, source(m), namefor(m), scopeof(m))
            end
        end
        add!(d, body, name, scope)
        d
    end
    combine() |> values |> collect
end
geninfos(S::Type{<:System}) = geninfos(source(S); name=namefor(S), substs=substsof(S), incl=Symbol.(mixinsof(S)), scope=scopeof(S))

include("dependency.jl")
sortednodes(infos) = sort(dependency(infos))

"""
    @system name[{patches..}][(mixins..)] [<: type] [decl] -> Type{<:System}

Declare a new system called `name` with new variables declared in `decl` block using a custom syntax. The resultant system is subtype of `System` or a custom `type`. `mixins` allows reusing specification of existing systems to be pasted into the declaration of new system. `patches` may provide type substitution and/or constant definition needed for advanced use.

# Variable
    name[(args..; kwargs..)][: alias] [=> expr] [~ [state][::type|<:type][(tags..)]]

- `name`: variable name; usually short abbreviation.
- `args`: automatically bound depending variables
- `kwargs`: custom bound depending variables; used by `call` and `integrate`.
- `alias`: alternative name; long description.
- `expr`: state-specific code snippet; use begin-end block for multiple statements.
- `type`: internal data type; default is Float64 for many, but not all, variables.
- `tags`: state-specific options; `unit`, `min`/`max`, etc.

# States
- `hold`: marks a placeholder for variable shared between mixins.
- `wrap`: passes a state variable to other fucnction as is with no unwrapping its value.
- `advance`: manages a time-keeping variable; `time` and `tick` from `Clock`.
- `preserve`: keeps initially assigned value with no further updates; constants, parameters.
- `tabulate`: makes a two dimensional table with named keys; *i.e.* partitioning table.
- `interpolate`: makes a curve fuction interpolated with discrete values; *i.e.* soil characteristic curve.
- `track`: evaluates variable expression as is for each update.
- `remember`: keeps tracking variable until a certain condition is met; essentially `track` turning into `preserve`.
- `provide`: manages a table of time-series in DataFrame.
- `drive`: fetchs the current value from a time-series; maybe supplied by `provide`.
- `call`: defines a partial function bound with some variables.
- `integrate`: calculates integral using Gaussian method; not for time domain.
- `accumulate`: emulates integration of rate variable over time; essentially Euler method.
- `capture`: calculates difference between integration for each time step.
- `flag`: sets a boolean flag; essentially `track::Bool`.
- `produce`: attaches a new instance of system dynamically constructed; *i.e.* root structure growth.
- `bisect`: solves nonlinear equation using bisection method; *i.e.* gas-exchange model coupling.
- `solve`: solves polynomical equation symbolically; *i.e.* quadratic equations in photosynthesis model.

# Examples
```julia-repl
julia> @system S(Controller) begin
           a => 1 ~ preserve(parameter)
           b(a) ~ accumulate
       end
S
```
"""
macro system(head, body=:(begin end))
    scope = __module__
    gensystem(body; scope, parsehead(head; scope)...)
end

export @system, update!

gendefault(v::VarInfo) = gendefault(v, Val(v.state))
gendefault(v::VarInfo, ::Val) = gendefaultvalue(v)

gensample(v::VarInfo, x) = @q $C.sample($x)
genunitfy(v::VarInfo, x) = begin
    u = gettag(v, :unit)
    isnothing(u) ? x : @q $C.unitfy($x, $C.value($u))
end
genminmax(v::VarInfo, x) = begin
    l = gettag(v, :min)
    u = gettag(v, :max)
    #TODO: validate (min <= max)
    x = isnothing(l) ? x : @q max($(genunitfy(v, @q $C.value($l))), $x)
    x = isnothing(u) ? x : @q min($x, $(genunitfy(v, @q $C.value($u))))
    x
end
genround(v::VarInfo, x) = begin
    f = gettag(v, :round)
    isnothing(f) && return x
    f = if f isa Bool
        f ? :round : return x
    elseif f isa QuoteNode
        f.value
    else
        error("unsupported value for tag `round`: $f")
    end
    U = gettag(v, :roundunit)
    U = isnothing(U) ? gettag(v, :unit) : U
    N = gettag(v, :_type)
    if isnothing(U)
        #HACK: rounding functions with explicit type only supports Integer target
        # https://github.com/JuliaLang/julia/issues/37984
        @q convert($N, $f($x))
    else
        @q $f($C.valuetype($N, $U), $x)
    end
end
genwhen(v::VarInfo, x) = begin
    w = gettag(v, :when)
    isnothing(w) && return x
    N = gettag(v, :_type)
    U = gettag(v, :unit)
    i = gettag(v, :init)
    d = if isnothing(i)
        @q zero($C.valuetype($N, $U))
    else
        @q $C.unitfy($C.value($i), $C.value($U))
    end
    @q $C.value($w) ? $x : $d
end

genparameter(v::VarInfo) = begin
    @gensym o
    @q let $o = $C.option(config, _names, $(names(v)))
        ismissing($o) ? $(genbody(v)) : $o
    end
end
geninit(v::VarInfo) = @q $C.value($(gettag(v, :init)))

gendefaultvalue(v::VarInfo; parameter=false, init=false, sample=true, unitfy=true, minmax=true, round=true, when=true) = begin
    x0 = if parameter && istag(v, :parameter)
        genparameter(v)
    elseif init && hastag(v, :init)
        geninit(v)
    else
        genbody(v)
    end
    x = x0
    sample && (x = gensample(v, x))
    unitfy && (x = genunitfy(v, x))
    minmax && (x = genminmax(v, x))
    round && (x = genround(v, x))
    when && (x = genwhen(v, x))
    if istag(v, :optional)
        @q isnothing($x0) ? nothing : $x
    else
        x
    end
end

genupdate(nodes::Vector{VarNode}, ::MainStage; kw...) = @q begin
    $([genupdateinit(n.info; kw...) for n in nodes]...)
    $([genupdate(n; kw...) for n in nodes]...)
    self
end

genupdate(infos::Vector{VarInfo}, t::UpdateStage; kw...) = @q begin
    $([genupdateinit(v; kw...) for v in infos]...)
    $([genupdate(v, t; kw...) for v in infos]...)
    self
end

symname(v::VarInfo) = symname(v.system, v.name)
symname(s::Symbol, n::Symbol) = n #Symbol(:_, s, :__, n)
symstate(v::VarInfo) = symname(v) #Symbol(symname(v), :__state)
symsuffix(::T) where {T} = "__$T"
symlabel(v::VarInfo, t, s...) = Symbol(symname(v), symsuffix(t), s...)
symcall(v::VarInfo) = Symbol(v.name, :__call)

genupdateinit(v::VarInfo; kw...) = begin
    # implicit :expose
    @q begin
        $(v.name) = self.$(v.name)
        $(isnothing(v.alias) ? ϵ : :($(v.alias) = $(v.name)))
    end
end

genupdate(n::VarNode; kw...) = genupdate(n.info, n.step; kw...)
genupdate(v::VarInfo, t; kw...) = @q begin
    $(linenumber(v, "genupdate", t))
    @label $(symlabel(v, t))
    $(genupdate(v, Val(v.state), t; kw...))
end

genvalue(v::VarInfo) = @q $C.value($(symstate(v)))
genstore(v::VarInfo, val=nothing; unitfy=true, minmax=true, round=true, when=true) = begin
    isnothing(val) && (val = genbody(v))
    unitfy && (val = genunitfy(v, val))
    minmax && (val = genminmax(v, val))
    round && (val = genround(v, val))
    when && (val = genwhen(v, val))
    #TODO: remove redundant unitfy() in store!()
    @gensym s
    @q let $s = $(symstate(v))
        $C.store!($s, $val)
    end
end

genupdate(v::VarInfo, ::Val{nothing}, ::PreStep; kw...) = begin
    if istag(v, :context)
        @gensym c
        @q let $c = $(v.name)
            $C.update!(self, $C.PreStage())
            $c
        end
    end
end
genupdate(v::VarInfo, ::Val{nothing}, ::MainStep; kw...) = begin
    if istag(v, :override, :skip)
        nothing
    else
        @q $C.update!($(v.name))
    end
end
genupdate(v::VarInfo, ::Val{nothing}, ::PostStep; kw...) = begin
    if istag(v, :context)
        l = symlabel(v, PreStep())
        @gensym c cc
        @q let $c = $(v.name),
               $cc = $c.context
            $C.update!(self, $C.PostStage())
            if !isnothing($cc) && $C.value($c.clock.time) < $C.value($cc.clock.time)
                @goto $l
            end
            $c
        end
    end
end

genupdate(v::VarInfo, ::Val, ::PreStep; kw...) = nothing
genupdate(v::VarInfo, ::Val, ::MainStep; kw...) = istag(v, :override, :skip) ? nothing : genstore(v)
genupdate(v::VarInfo, ::Val, ::PostStep; kw...) = nothing

genupdate(v::VarInfo, ::Val, ::UpdateStage; kw...) = nothing
genupdate(v::VarInfo, ::Val{nothing}, ::PreStage; kw...) = @q $C.update!($(v.name), $C.PreStage())
genupdate(v::VarInfo, ::Val{nothing}, ::PostStage; kw...) = @q $C.update!($(v.name), $C.PostStage())

#TODO: merge extractfuncargdep() and extractfuncargkey()?
extractfuncargdep(v::Expr) = begin
    a = v.args
    if isexpr(v, :call)
        # detect boolean operators between state vars (i.e. `a`, `b` in `a && b`, `a || b`)
        if a[1] == :& || a[1] == :|
            extractfuncargdep.(a[2:3]) |> Iterators.flatten |> collect
        # detect variable inside wrapping function (i.e. `a` in `nounit(a.b, ..)`)
        else
            extractfuncargdep(a[2])
        end
    # detect shorthand syntax for calling value() (i.e. `a` in `a'` = `value(a)`)
    elseif isexpr(v, Symbol("'"))
        extractfuncargdep(a[1])
    # detect first callee of dot chaining (i.e. `a` in `a.b.c`)
    elseif isexpr(v, :., :ref)
        extractfuncargdep(a[1])
    else
        nothing
    end
end
extractfuncargdep(v::Symbol) = [v]
extractfuncargdep(v) = nothing

extractfuncargkey(v::Expr) = begin
    a = v.args
    if isexpr(v, :call)
        # detect boolean operators between state vars (i.e. `a`, `b` in `a && b`, `a || b`)
        if a[1] == :& || a[1] == :|
            error("missing function argument key: $v")
        # detect variable inside wrapping function (i.e. `b` in `nounit(a.b, ..)`)
        else
            extractfuncargkey(a[2])
        end
    # detect shorthand syntax for calling value() (i.e. `b` in `a.b'` = `value(a.b)`)
    elseif isexpr(v, Symbol("'"))
        extractfuncargkey(a[1])
    # detect last callee of dot chaining (i.e. `c` in `a.b.c`)
    elseif isexpr(v, :., :ref)
        extractfuncargkey(a[2])
    else
        error("unrecognized function argument key: $v")
    end
end
extractfuncargkey(v::QuoteNode) = extractfuncargkey(v.value)
extractfuncargkey(v::Symbol) = v

extractfuncargpair(a) = let k, v
    !@capture(a, k_=v_) && (k = a; v = a)
    extractfuncargkey(k) => v
end
emitfuncargpair(a; kw...) = emitfuncargpair(extractfuncargpair(a)...; kw...)
emitfuncargpair(k, v; value=true) = begin
    k = esc(k)
    v = value ? @q($C.value($v)) : v
    @q $k = $v
end

genbody(v::VarInfo, body=nothing) = begin
    if isnothing(v.body) && length(v.args) == 1
        a = v.args[1]
        #HACK: can't use or (|) syntax: https://github.com/MikeInnes/MacroTools.jl/issues/36
        if @capture(a, k_=v_)
            args = emitfuncargpair(k, v)
        elseif @capture(a, k_Symbol)
            args = emitfuncargpair(k, k)
        else
            @gensym k
            args = emitfuncargpair(k, a)
        end
        body = esc(k)
    else
        args = @q begin $(emitfuncargpair.(v.args)...) end
        isnothing(body) && (body = esc(v.body))
    end
    MacroTools.flatten(@q let $args; $body end)
end

extractfunckwargtuple(a) = let k, t, u
    @capture(a, k_::t_(u_) | k_::t_ | k_(u_))
    isnothing(k) && (k = a)
    (k, t, u)
end
emitfunckwargkey(a) = @q $(esc(extractfunckwargtuple(a)[1]))
emitfunckwargpair(a) = begin
    k, t, u = extractfunckwargtuple(a)
    v = esc(k)
    v = isnothing(u) ? @q($v) : @q($C.unitfy($v, $u))
    # Skip type assertion (maybe only needed for Call, not Integrate)
    #v = @q $v::$C.valuetype($(gencallargtype(t)), $u)
    @q $k = $v
end

genfunc(v::VarInfo; unitfy=true) = begin
    innerargs = @q begin $(emitfuncargpair.(v.args)...) end
    innerbody = MacroTools.flatten(@q let $innerargs; $(esc(v.body)) end)
    unitfy && (innerbody = @q $C.unitfy($innerbody, $C.value($(gettag(v, :unit)))))

    callargs = emitfunckwargkey.(v.kwargs)
    argsheader = emitfunckwargpair.(v.kwargs)

    @q function $(symcall(v))($(callargs...))
        let $(argsheader...)
            $innerbody
        end
    end
end


text:
using Printf
function MergeRestorefile( srcfile )
    ftgr = open( srcfile,"w" )
    IDX = 0
    while true
        IDX +=1
        srcfilename = @sprintf("%s.%03d",srcfile,IDX)
        if  ! ispath(srcfilename )
            break
        end
        try
            fsrc =open(srcfilename,"r")
            a= UInt8(0)
            while !eof( fsrc )
                a = read( fsrc,UInt8 )
                write( ftgr, a)
            end
            close(fsrc)
            println(IDX," pos:",position( ftgr ) )
        catch(e)
            println(e)
            break
        end
    end
    close(ftgr)
end

function sliptfile( srcfile )
    fobj= open( srcfile,"r" )
    # byteNum =0x00;
    maxBytesPerfile =100*1014*1024
    IDX = 0
    while !eof( fobj )
        IDX +=1
        tgrfilename = @sprintf("%s.%03d",srcfile,IDX)
        fout =open(tgrfilename,"w")
        segNum = 1;
        a= UInt8(0)
        while !eof( fobj )
            a = read( fobj,UInt8 )
            write( fout, a)
            segNum +=1
            if segNum > maxBytesPerfile
                break
            end
        end
        close(fout)
        println(IDX," pos:",position( fobj) )
    end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# read_rdb_hdr.jl

using Test: @test, @testset, @test_throws
using MIRTio: header_write
using MIRTio: rdb_hdr_26_002_def, read_rdb_hdr
using MIRTio: header_init, header_size, header_string, header_write

# random data self test

hd = rdb_hdr_26_002_def()
ht = header_init(hd) # random header values
@test ht isa NamedTuple

tname = tempname()

open(tname, "w") do fid
    header_write(fid, ht ; bytes = header_size(hd))
end

# intentionally fail due to random :rdbm_rev
@test_throws String read_rdb_hdr(tname)

# now set :rdbm_rev so that read test passes
ht = merge(ht, [:rdbm_rev => Float32(26.002)]) # instead of setindex

open(tname, "w") do fid
    header_write(fid, ht ; bytes = header_size(hd))
end

hr = read_rdb_hdr(tname)
@test hr isa NamedTuple
@test hr[:rdbm_rev] isa Float32
@test hr == header_string(ht)

# self test using a local pfile at UM
# UM-only test below here

file = "/n/ir71/d3/fessler/fmri-data-michelle-L+S/P97792.7"
isfile(file) && (@test read_rdb_hdr(file).dab[2] == 31)


text:
###
# transforms directionality
###
import NumericExtensions.negate!, NumericExtensions.negate!, NumericExtensions.multiply!, NumericExtensions.subtract!, NumericExtensions.add!, NumericExtensions.map1!, NumericExtensions.SinFun, NumericExtensions.ModFun

function alignJl!(Conc :: Matrix, Dir :: Matrix, Newdir :: Matrix, attraction :: Real, step :: Real) #concentration and direction
    d1,d2 = size(Conc)

    diff = zeros(Float64, 3,3)
    tmp = zeros(Float64, 3,3)
    potential = zeros(Float64, 3,3)

    for j in 1:d2
        for i in 1:d1

            west  = j == 1  ? d2 : j-1
            east  = j == d2 ? 1  : j+1
            north = i == d1 ? 1  : i+1
            south = i == 1  ? d1 : i-1

            # Getting the data and storing it in two eight vectors
            direction = Dir[i, j];

            s0 = Conc[south,west];
            s1 = Conc[south,j   ];
            s2 = Conc[south,east];
            s3 = Conc[i    ,east];
            s4 = Conc[north,east];
            s5 = Conc[north,j   ];
            s6 = Conc[north,west];
            s7 = Conc[i    ,west];

            conc = [s0, s1, s2, s3, s4, s5, s6, s7]

            s0 = Dir[south,west];
            s1 = Dir[south,j   ];
            s2 = Dir[south,east];
            s3 = Dir[i    ,east];
            s4 = Dir[north,east];
            s5 = Dir[north,j   ];
            s6 = Dir[north,west];
            s7 = Dir[i    ,west];

            dir = [s0, s1, s2, s3, s4, s5, s6, s7]

            # Calculate
            diff = negate!(dir)
            diff = add!(diff, direction);
            diff = map1!(ModFun(), diff, PI);
            diff = negate!(map1!(SinFun(), (multiply!(diff,2))));

            dtheta = sum(multiply!(diff, conc));

            dtheta = attraction * dtheta / 8;
            ndir = direction + dtheta * step;

            Newdir[i,j] = fmod(ndir,PI);
        end
    end
end


text:
using LinearAlgebra, ForwardDiff, Distributions
using Convex, SCS, ECOS
using IterativeSolvers

# problem setup
p = 60
N = 30
x_sol = max.(0.0, randn(p))
G = rand(N, p)
h = G * x_sol
ν_sol = randn(N)
λ_sol = rand(p)
λ_sol[x_sol .> 0.0] .= 0.0
c = -G' * ν_sol + λ_sol

"Convex.jl"
x = Variable(p)
prob = minimize(c' * x)
prob.constraints += h - G * x == 0.0
prob.constraints += x >= 0.0

@time solve!(prob, ECOS.Optimizer)

@show prob.status
@show prob.constraints
@show x.value
# @show prob.constraints[1].dual
# @show prob.constraints[2].dual
A = [-1.0 * Diagonal(ones(p)); G; -G]
b = [zeros(p); h; -h]
m = 2 * N + p
k = p + m + 1

Q = Array([zeros(p, p) A' c;
           -A zeros(m, m) b;
           -c' -b'      0.0])

function κ(z)
    # z_proj = zero(z)
    # z_proj[1:N] = z[1:N]
    # z_proj[N .+ (1:p)] = max.(0.0, z[N .+ (1:p)])

    z_proj = max.(0.0, z)
    return z_proj
    # return z
end

κ(randn(10))

function Pc(z)
    z_proj = zero(z)
    z_proj[1:p] = z[1:p]
    z_proj[p .+ (1:m)] = κ(z[p .+ (1:m)])
    z_proj[p + m + 1] = max(0.0, z[p + m + 1])

    return z_proj
end

function F(z)
    ũ = z[1:k]
    u = z[k .+ (1:k)]
    v = z[2 * k .+ (1:k)]

    [(I + Q) * ũ - (u + v);
     u - Pc(ũ - v);
     ũ - u]
end

function Ju(z)
    ũ = z[1:k]
    u = z[k .+ (1:k)]
    v = z[2k .+ (1:k)]

    JP = zeros(m, m)
    dif = ũ[p .+ (1:m)] - v[p .+ (1:m)]
    for i = 1:m
        if dif[i] >= 0.0
            JP[i, i] = 1.0
        else
            JP[i, i] = 0.0
        end
    end

    if z[k] - z[3k] >= 0.0
        ℓ = 1.0
    else
        ℓ = 0.0
    end

    [-I zeros(p, m) zeros(p) I zeros(p, m) zeros(p) I zeros(p, m) zeros(p);
     zeros(m, p) -JP zeros(m) zeros(m, p) I zeros(m) zeros(m, p) JP zeros(m);
     [zeros(1, p) zeros(1, m) -ℓ zeros(1, p) zeros(1, m) 1.0 zeros(1, p) zeros(1, m) ℓ]]
end

function J(z)
    [(I + Q) -I -I
     Ju(z);
     I -I zeros(k, k)]
end

ũ = zeros(k)
u = zeros(k)
v = zeros(k)
ũ[end] = 1.0
u[end] = 1.0
v[end] = 1.0

z = [ũ; u; v]
z = rand(3k)
F(z)
rank(Ju(z))
J(z)
rank(ForwardDiff.jacobian(F, z))
norm(J(z) - ForwardDiff.jacobian(F, z), Inf)
gmres(J(z), -F(z))
# (J(z)' * J(z)) \ (J(z)' * F(z))
# Δ = zero(z)
# Δ[k] = 1.0
# Δ[2k] = 1.0
# Δ[3k] = 1.0
# gmres!(Δ, J(z), -1.0 * F(z), abstol = 1.0)

function solve()
    ũ = zeros(k)
    u = zeros(k)
    v = zeros(k)
    ũ[end] = 1.0
    u[end] = 1.0
    v[end] = 1.0

    z = [ũ; u; v]
    # z = rand(3k)
    # Δ = zero(z)

    extra_iters = 0

    for i = 1:500
        _F = F(z)
        _J = J(z)
        Δ = gmres(_J, 1.0 * _F, abstol = 1.0e-12, maxiter = i + extra_iters)

        iter = 0
        α = 1.0
        while norm(F(z - α * Δ))^2.0 >= (1.0 - 0.001 * α) * norm(_F)^2.0 && α > 1.0e-4
            α = 0.5 * α
            # println("   α = $α")
            iter += 1
            if iter > 100
                @error "line search fail"
                # return z
                x = z[k .+ (1:p)]
                y = z[k + p .+ (1:m)]
                τ = z[k + k]
                κ = z[3k]
                println("τ = $τ")
                println("κ = $κ")
                return x ./ τ, z
            end
        end

        if α <= 1.0e-4
            extra_iters += 1
        end

        println("iter ($i) - norm: $(norm(F(z)))")

        z .-= α * Δ
    end

    # return z
    x = z[k .+ (1:p)]
    y = z[k + p .+ (1:m)]
    τ = z[k + k]
    κ = z[3k]
    println("τ = $τ")
    println("κ = $κ")
    return x ./ τ, z
end

x_sol, z_sol = solve()

x_sol
b - A * x_sol
norm(x_sol - x.value)

(Q - I) * Pc(z_sol[1:k]) + z_sol[1:k]


text:
export Radiation, Elements, Isotope

@enum DecayType α β γ n
# Write your package code here.
struct Radiation
    type::DecayType
    energy::Float32
    Radiation() = new(α,6.5)
end

struct Elements
    file
    Elements() = new(CSV.File("../Elements.csv"))
end
get_iso_row(z::Int,elems::Elements) = elems.file[elems.file.AtomicNumber .== z,:]
get_iso_row(z::Int) = get_iso_row(z,Elements())
get_iso_row(symbol::String,elems::Elements) = elems.file[elems.file.Symbol .== symbol,:]

struct Isotope
    Z
    A
    symbol::String
    name::String
    radioactive::Bool
    Isotope() = new(94,239,"Plutonium","Pu",true)
    function Isotope(z::Int,a::Int)
        isoRow = get_iso_row(z)
        result = first(isoRow)
        sym = result.Symbol
        nme = result.Element
        wkStr = result.Radioactive
        if wkStr !== nothing
            active = parsebool(String(wkStr))
        end
        a = result.AtomicMass
        new(z,a,sym,nme,active)
    end
    function Isotope(z::Int)
        isoRow = get_iso_row(z)
        result = first(isoRow)
        sym = result.Symbol
        nme = result.Element
        wkStr = result.Radioactive
        if wkStr !== nothing
            active = parsebool(String(wkStr))
        end
        a = result.AtomicMass
        new(z,a,sym,nme,active)
    end
end
function decay!(iso::Isotope,radType::DecayType)
    if radType == α
        newIso = Isotope(iso.Z-2,iso.A-4)
    elseif radType == β
        newIso = Isotope(iso.Z+1,iso.A)
    elseif radType == n
        newIso = Isotope(iso.Z,iso.A-1)
    else
        newIso = Isotope(iso.Z,iso.A)
    end
    return newIso
end


text:
function primal_inds(d)
    N = length(d)
    inds = zeros(Cint, N, 2)
    inds[1, 1] = 1
    inds[1, 2] = d[1]
    for n ∈ 2:N
        inds[n, 1] = inds[n-1, 2] + 1
        inds[n, 2] = inds[n, 1] + d[n] - 1
    end
    inds
end
@non_differentiable primal_inds(d)

function prob_prod(x, ind, primal_inds, n...)
    N = size(primal_inds, 1)
    length(n) == N && return 1.0
    prod(x[primal_inds[i, 1]+ind[i]-1] for i ∈ 1:N if i ∉ n)
end

function expected_cost(CT, x::Vector{T}, indices, primal_inds)::T where {T}
    val = sum(CT[ind] * prob_prod(x, ind, primal_inds) for ind ∈ indices)
end

function expected_cost(x, cost_tensor)
    N = length(x)
    d = size(cost_tensor)
    tensor_indices = CartesianIndices(cost_tensor)
    primal_indices = primal_inds(d)
    expected_cost(cost_tensor, reduce(vcat, x), tensor_indices, primal_indices)
end

function ChainRulesCore.rrule(::typeof(expected_cost), x, cost_tensor)
    N = length(x)
    d = size(cost_tensor)
    tensor_indices = CartesianIndices(cost_tensor)
    primal_indices = primal_inds(d)
    xx = reduce(vcat, x)
    value = expected_cost(cost_tensor, xx, tensor_indices, primal_indices)

    function expected_cost_pullback(∂value)
        ∂x = zeros(sum(d))
        for n ∈ 1:N
            grad!(∂x, cost_tensor, xx, n, tensor_indices, primal_indices)
        end
        ∂x .*= ∂value
        ∂x = [∂x[primal_indices[n, 1]:primal_indices[n, 2]] for n ∈ 1:N]

        ∂ct = zero(cost_tensor)
        for ind ∈ tensor_indices
            ∂ct[ind] = ∂value * prob_prod(xx, ind, primal_indices)
        end
        NoTangent(), ∂x, ∂ct
    end
    value, expected_cost_pullback
end

function grad!(f, CT, x, n, indices, primal_inds)
    f[primal_inds[n, 1]:primal_inds[n, 2]] .= 0.0
    for ind ∈ indices
        prob = prob_prod(x, ind, primal_inds, n)
        for i ∈ 1:primal_inds[n, 2]+1-primal_inds[n, 1]
            if ind[n] == i
                f[primal_inds[n, 1]+i-1] += CT[ind] * prob
            end
        end
    end
    Cint(0)
end

function jac(J, loc, T, x, n1, n2, i, j, indices, primal_inds)
    l1 = primal_inds[n1, 2] - primal_inds[n1, 1] + 1
    l2 = primal_inds[n2, 2] - primal_inds[n2, 1] + 1
    (n1 == n2) && return zeros(l1, l2)
    for ind ∈ indices
        if ind[n1] == i && ind[n2] == j
            J[loc] += T[ind] * prob_prod(x, ind, primal_inds, n1, n2)
        end
    end
    Cint(0)
end

struct Wrapper{L} <: Function
    tensors::Vector{Array{Float64,L}}
    tensor_indices::CartesianIndices{L,NTuple{L,Base.OneTo{Int64}}}
    primal_inds::Matrix{Cint}
    dual_inds::Vector{Cint}
    N::Cint
    m::NTuple{L,Cint}
    num_primals::Cint
end

function (T::Wrapper)(n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})
    ind = 0
    for (n, tensor) ∈ enumerate(T.tensors)
        f[ind+1:ind+T.m[n]] .= 0.0
        grad!(f, tensor, x, n, T.tensor_indices, T.primal_inds)
        for j ∈ ind+1:ind+T.m[n]
            f[j] -= x[T.dual_inds[n]]
        end
        ind += T.m[n]
    end
    for n ∈ 1:T.N
        f[ind+n] = -1.0
        for i ∈ T.primal_inds[n, 1]:T.primal_inds[n, 2]
            f[ind+n] += x[i]
        end
    end
    Cint(0)
end

function (T::Wrapper)(n::Cint,
    nnz::Cint,
    x::Vector{Cdouble},
    col::Vector{Cint},
    len::Vector{Cint},
    row::Vector{Cint},
    data::Vector{Cdouble})
    @assert n == length(x) == length(col) == length(len)
    @assert nnz == length(row) == length(data)

    col[1] = 1
    ind = 2
    for n in 1:T.N
        inc = 1
        for j in (1:T.N)
            if j ≠ n
                inc = inc + T.m[j]
            end
        end
        for _ in 1:T.m[n]
            col[ind] = col[ind-1] + inc
            len[ind-1] = inc
            ind += 1
        end
    end
    for n in 1:T.N
        if n ≠ T.N
            col[ind] = col[ind-1] + T.m[n]
        end
        len[ind-1] = T.m[n]
        ind += 1
    end

    ind = 0
    for n_col in 1:T.N
        for c in 1:T.m[n_col]
            for n_row in 1:T.N
                if n_row == n_col
                    continue
                end
                for p in 1:T.m[n_row]
                    row[ind+1] = Cint(T.primal_inds[n_row, 1] + p - 1)
                    data[ind+1] = 0.0
                    jac(data, ind + 1, T.tensors[n_row], x, n_row, n_col, p, c, T.tensor_indices, T.primal_inds)
                    ind += 1
                end
            end
            row[ind+1] = Cint(T.num_primals + n_col)
            data[ind+1] = 1.0
            ind += 1
        end
    end
    for n in 1:T.N
        row[ind+1:ind+T.m[n]] = T.primal_inds[n, 1]:T.primal_inds[n, 2]
        data[ind+1:ind+T.m[n]] .= -1.0
        ind += T.m[n]
    end
    return Cint(0)
end

function compute_equilibrium(cost_tensors::AbstractVector{<:AbstractArray{<:ForwardDiff.Dual{T}}}; kwargs...) where {T}
    # strip off the duals:
    cost_tensors_v = [ForwardDiff.value.(c) for c in cost_tensors]
    cost_tensors_p = [ForwardDiff.partials.(c) for c in cost_tensors]
    # forward pass
    res = compute_equilibrium(cost_tensors_v; kwargs...)
    # backward pass
    # 1. compute jacobian
    _back = _compute_equilibrium_pullback(res)
    # 2. project input-sensitivy through jacobian to yield output sensitivy
    x_p = let
        # output sensitivities stacked for all players
        x_p_stacked = sum(zip(cost_tensors_p, _back)) do (cost_tensor_p, ∂cost_tensor)
            inflated_cost_tensor_p = reshape(cost_tensor_p, 1, size(cost_tensor_p)...)
            dims = 2:ndims(inflated_cost_tensor_p)
            reshape(sum(∂cost_tensor .* inflated_cost_tensor_p; dims), size(∂cost_tensor, 1))
        end
        # unstacking the output sensitivities
        n_actions_per_player = size(first(cost_tensors))
        x_p_stacked_it = Iterators.Stateful(x_p_stacked)
        map(n_actions_per_player) do n_actions
            Iterators.take(x_p_stacked_it, n_actions) |> collect
        end
    end

    # 3. glue primal and dual results together into a ForwardDiff.Dual-valued result
    x_d = [ForwardDiff.Dual{T}.(xi_v, xi_p) for (xi_v, xi_p) in zip(res.x, x_p)]

    (; x = x_d, res.λ, res._deriv_info)
end

function compute_equilibrium(cost_tensors;
    initialization = nothing,
    ϵ = 0.0,
    silent = true,
    convergence_tolerance = 1e-6)

    N = Cint(length(cost_tensors))
    m = Cint.(size(cost_tensors[1]))
    @assert all(m == size(tensor) for tensor ∈ cost_tensors)
    primal_indices = primal_inds(m)

    num_primals::Cint = sum(m)
    dual_inds = Vector{Cint}(num_primals+1:num_primals+N)

    tensor_indices = CartesianIndices(cost_tensors[1])
    nnz = 2 * num_primals
    for n = 1:N
        nnz += m[n] * (num_primals - m[n])
    end

    wrapper! = Wrapper(cost_tensors, tensor_indices, primal_indices, dual_inds, N, m, num_primals)

    lb = [ϵ * ones(Cdouble, num_primals); -1e20 * ones(Cdouble, N)]
    ub = 1e20 * ones(Cdouble, num_primals + N)
    z = zeros(Cdouble, num_primals + N)
    if isnothing(initialization)
        for n ∈ 1:N
            start_ind = primal_indices[n, 1]
            end_ind = primal_indices[n, 2]
            z[start_ind:end_ind] .= 1.0 / m[n]
        end
    else
        z .= initialization
    end

    n = Cint(sum(m) + N)
    x = randn(Cdouble, n)
    f = zeros(Cdouble, n)

    status, vars::Vector{Cdouble}, info = PATHSolver.solve_mcp(
        wrapper!,
        wrapper!,
        lb,
        ub,
        z;
        nnz,
        convergence_tolerance,
        silent)

    x = [vars[primal_indices[n, 1]:primal_indices[n, 2]] for n ∈ 1:N]
    λ = vars[dual_inds]

    (; x, λ, _deriv_info = (; ϵ, wrapper!, nnz = Cint(nnz), tensor_indices, primal_indices))
end

function ChainRulesCore.rrule(::typeof(compute_equilibrium),
    cost_tensors;
    initialization = nothing,
    ϵ = 0.0,
    silent = true,
    convergence_tolerance = 1e-6)
    res = compute_equilibrium(cost_tensors; initialization, ϵ, silent, convergence_tolerance)

    _back = _compute_equilibrium_pullback(res)

    function compute_equilibrium_pullback(∂res)
        ∂self = NoTangent()

        ∂cost_tensors = let
            full_sensitivities = map(∂res.x, res.x) do r, x
                r isa ZeroTangent ? zeros(x) : r
            end
            derivs = reduce(vcat, full_sensitivities)

            map(_back) do ∂cost_tensor
                dropdims(sum(∂cost_tensor .* derivs; dims = 1); dims = 1)
            end
        end
        ∂self, ∂cost_tensors
    end

    res, compute_equilibrium_pullback
end

function _compute_equilibrium_pullback(res; bound_tolerance = 1e-6, singularity_tolerance = 1e-6)
    primals = vcat(res.x...)
    vars = [primals; res.λ]
    n = Cint(length(vars))
    N = length(res.λ)
    d = [length(xi) for xi ∈ res.x]
    starts = cumsum([0; d[1:end-1]])
    num_primals = sum(d)

    lb = [res._deriv_info.ϵ * ones(num_primals); -Inf * ones(N)]
    unbound_indices = (vars .> (lb .+ bound_tolerance))
    unbound_primals = unbound_indices[1:end-N]
    nup = sum(unbound_primals)
    ubmap = zeros(Int, num_primals)
    idx = 1
    for i ∈ 1:num_primals
        if unbound_primals[i]
            ubmap[i] = idx
            idx += 1
        end
    end
    col = zeros(Cint, n)
    len = zeros(Cint, n)
    row = zeros(Cint, res._deriv_info.nnz)
    data = zeros(Cdouble, res._deriv_info.nnz)
    res._deriv_info.wrapper!(n, res._deriv_info.nnz, vars, col, len, row, data)

    colptr = zeros(Cint, n + 1)
    colptr[1:n] .= col
    colptr[n+1] = col[n] + len[n]

    ∂cost_tensors = [zeros(sum(d), d...) for _ ∈ 1:N]

    nJ = Matrix{Cdouble}((SparseMatrixCSC{Cdouble,Cint}(n, n, colptr, row, data))[unbound_indices, unbound_indices])
    factorization = qr(nJ)
    if any(abs(r) ≤ singularity_tolerance for r ∈ diag(factorization.R))
        # Artificially returning zero-derivatives if solution is non-isolated.
        return ∂cost_tensors
    end
    nJi = factorization \ (-I)

    for ind ∈ res._deriv_info.tensor_indices
        for n ∈ 1:N
            if ubmap[starts[n]+ind[n]] > 0
                ∂cost_tensors[n][unbound_primals, ind] .= (nJi[1:nup, ubmap[starts[n]+ind[n]]] *
                                                           prob_prod(primals, ind, res._deriv_info.primal_indices, n))
            end
        end
    end
    ∂cost_tensors
end