text:
# The Dual Arrhenius and Michealis-Menten (DAMM) model, Davidson et al. 2012

"""
    DAMM(x::VecOrMat{<: Real}, p::NTuple{7, Float64})
Calculate respiration as a function of soil temperature and moisture.
# Examples
```julia-repl
julia> Tₛ = [18.0, 22.0] # soil temperature [°C]
julia> θ = [0.35, 0.22] # soil moisture [m³ m⁻³]
julia> x = hcat(Tₛ, θ)
julia> p = (1e9, 64.0, 3.46e-8, 2.0e-3, 0.4, 0.0125, 1.0) # αₛₓ, Eaₛₓ, kMₛₓ, kMₒ₂, Sxₜₒₜ, Q10kM
julia> DAMM(x, p)
  1.6 # μmolCO₂ m⁻² s⁻¹
  2.8 # μmolCO₂ m⁻² s⁻¹
```
"""
function DAMM(x::VecOrMat{<: Real}, p::NTuple{7, Float64})
# Independent variables
  Tₛ = x[:, 1]°C # Soil temperature
  Tₛ = Tₛ .|> K # Tₛ in Kelvin
  θ = x[:, 2]m^3*m^-3 # Soil moisture, m³ m⁻³
# Parameters
  αₛₓ = p[1]mg*cm^-3*hr^-1 # Pre-exponential factor, mgC cm⁻³ h⁻¹
  Eaₛₓ = p[2]kJ*mol^-1 # Activation energy, kJ mol⁻¹
  kMₛₓ = p[3]g*cm^-3 # Michaelis constant, gC cm⁻³
  kMₒ₂ = p[4]L*L^-1 # Michaelis constant for O₂, L L⁻¹
  porosity = p[5] # 1 - soil buld density / soil particle density
  Sxₜₒₜ = p[6]g*cm^-3 # Total soil carbon, gC cm⁻³
  Q10Km = p[7]
  # DAMM model
  Tref = 293.15K
  Vmax = @. αₛₓ * exp(-Eaₛₓ/(R * Tₛ)) # Maximum potential rate of respiration
  Sₓ = @. pₛₓ * Sxₜₒₜ * Dₗᵢ * θ^3 # All soluble substrate, gC cm⁻³
  MMₛₓ = @. Sₓ / (kMₛₓ * Q10Km^(ustrip(Tₛ - Tref)/10.0) + Sₓ) # Availability of substrate factor, 0-1
  airfilled_porosity = porosity .- θ
  airfilled_porosity[airfilled_porosity .< 0.0] .= NaN
  O₂ = @. Dₒₐ * O₂ₐ * (airfilled_porosity^(4/3)) # Oxygen concentration
  MMₒ₂ = @. O₂ / (kMₒ₂ + O₂) # Oxygen limitation factor, 0-1
  Rₛₘ = @. Vmax * MMₛₓ * MMₒ₂ # Respiration, mg C cm⁻³ hr⁻¹
  Rₛₘₛ = @. Rₛₘ * Eₛ # Respiration, effective depth 10 cm, mg C cm⁻² hr⁻¹
  Rₛₛₘ = Rₛₘₛ .|> molC*m^-2*s^-1
  Rₛ = Rₛₛₘ .|> μmolCO₂*m^-2*s^-1
  Rₛ = Float64.(ustrip(Rₛ)) # no units
  return Rₛ
end


text:
# Include this startup file prior to running Julia code

# Add project module locations to path
push!(LOAD_PATH, abspath(joinpath("src","distributions")))
push!(LOAD_PATH, abspath(joinpath("src","samplers")))

ext: jl


text:
function genproperties(encoder::Dict{Symbol,Symbol})
    decoder = Dict(collect(values(encoder)) .=> collect(keys(encoder)))
    decoder_diac = [encoder[Symbol(s[1])] for s in split(AR_DIACS_REGEX.pattern, "|")]
    idx = map(x -> in(x[1], SP_REGEX_CHARS), string.(decoder_diac))
    decoder_diac[idx] = Symbol.(raw"\\" .* string.(decoder_diac[idx]))
    decoder_diac = Regex(join(string.(decoder_diac), "|"))

    return decoder, decoder_diac
end

function Transliterator(x::Bool) end

abstract type AbstractEncoder end

"""
    @transliterator(dict, name)

Create a custom transliterator using an input `dict` (`Dict` object) with its corresponding
`name` as `String` object. This will automatically update the transliterator inside all
functions like `arabic`, `verses`, and `encode`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> old_keys = collect(keys(BW_ENCODING))
julia> new_vals = reverse(collect(values(BW_ENCODING)))
julia> my_encoder = Dict(old_keys .=> new_vals)
julia> basmala = arabic(verses(crpsdata[1][1])[1])
julia> @transliterator my_encoder "MyEncoder"
julia> encode(basmala)
"\"S%gAS zppj[KS zp`j[&gA[r]S zp`j[&SkAS"
```
"""
macro transliterator(dict, name)
    T = Symbol(uppercasefirst(name))
    esc(quote
            prop = genproperties(eval($dict))
            struct $T <: AbstractEncoder
                encode::Dict{Symbol,Symbol}
                decode::Dict{Symbol,Symbol}
                rx_diacs::Regex
                rx_ardiacs::Regex
            end
            QuranTree.Transliterator() = $T($dict, prop[1], prop[2], $AR_DIACS_REGEX)
            function Base.show(io::IO, t::$T)
                println(io, $T, ":")
                println(io, " ├ encode: ", Transliterator().encode)
                println(io, " ├ decode: ", Transliterator().decode)
                println(io, " ├ rx_diacs: ", Transliterator().rx_diacs)
                println(io, " └ rx_ardiacs: ", Transliterator().rx_ardiacs)
            end
        end
    )
end

"""
    @transliterator(symbl)

Fallback to the default `Buckwalter` transliterator.
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps);
julia> tnzldata = table(tnzl);
julia> old_keys = collect(keys(BW_ENCODING))
julia> new_vals = reverse(collect(values(BW_ENCODING)))
julia> my_encoder = Dict(old_keys .=> new_vals)
julia> basmala = arabic(verses(crpsdata[1][1])[1])
julia> @transliterator my_encoder "MyEncoder"
julia> encode(basmala)
"\"S%gAS zppj[KS zp`j[&gA[r]S zp`j[&SkAS"
julia> @transliterator :default
julia> encode(basmala)
"bisomi {ll~ahi {lr~aHoma`ni {lr~aHiymi"
```
"""
macro transliterator(symbl)
    esc(quote
        if string($symbl) === "default"
            eval(:(@transliterator(BW_ENCODING, "Buckwalter")))
        else
            throw(DomainError("Expects :default, got " * string(symbl)))
        end
    end)
end

@transliterator BW_ENCODING "Buckwalter"


text:
"""
Structures which contain a handle (opaque pointer) as primary data. Those structures are usually defined as mutable so that a finalizer can be registered. Any Handle structure is automatically converted to its handle data on `ccall`s, through `unsafe_convert`.
"""
abstract type Handle end

Base.unsafe_convert(::Type{<: Ptr}, handle::Handle) = handle.h

"""
Connection to the X server.
"""
mutable struct Connection <: Handle
    """
    Opaque handle to the connection, used for API calls.
    """
    h
    function Connection(h)
        conn = new(h)
        Base.finalizer(xcb_disconnect, conn)
        conn
    end
end

"""
Connection setup handle and data.
"""
struct Setup <: Handle
    """
    Handle to the setup, used for API calls.
    """
    h
    """
    Setup value, obtained when dereferencing its handle.
    """
    value::xcb_setup_t
end

"""
Initialize a connection to the X server.
"""
function Connection(; display=nothing)
    display = isnothing(display) ? C_NULL : display
    connection = Connection(xcb_connect(display, C_NULL))
    check(connection)
end

@enum XCB_CONN_ERROR_CODE begin
    XCB_CONN_ERROR = 1
    XCB_CONN_CLOSED_EXT_NOTSUPPORTED = 2
    XCB_CONN_CLOSED_MEM_INSUFFICIENT = 3
    XCB_CONN_CLOSED_REQ_LEN_EXCEED = 4
    XCB_CONN_CLOSED_PARSE_ERR = 5
    XCB_CONN_CLOSED_INVALID_SCREEN = 6
    XCB_CONN_CLOSED_FDPASSING_FAILED = 7
end

"""
Check that the connection to the X server was successful. Throws a [`ConnectionError`](@ref) if the connection failed.
"""
function check(connection::Connection)
    code = xcb_connection_has_error(connection)
    code ≠ 0 && throw(ConnectionError("XCB connection not successful ($(XCB_CONN_ERROR_CODE(code)))", code))
    connection
end

function Setup(connection::Connection)
    stp = xcb_get_setup(connection)
    Setup(stp, unsafe_load(stp))
end

"""
Check that the flush was successful, throwing a [`FlushError`](@ref) if the code is negative.
"""
function check_flush(code)
    if code <= 0
        throw(FlushError(code))
    end
end

Base.flush(connection::Connection) = check_flush(xcb_flush(connection))

"""
Check that the request was successfully handled by the server, throwing a [`RequestError`](@ref) if the request failed.

The severity level can be set via `level` with one of the two following values:

    - `:warn`: logs a warning with `@warn` (default).
    - `:error`: raise a `RequestError`.
"""
function check_request(conn, request; level=:warn)
    @assert level ∈ (:warn, :error) "Unknown severity level $level"
    errcode_ptr = xcb_request_check(conn, request)
    if errcode_ptr ≠ C_NULL
        errcode = Base.unsafe_load(errcode_ptr).error_code
        msg = unsafe_string(xcb_event_get_error_label(errcode))
        error_msg = "Request unsuccessful: (error code $errcode: $msg)"
        if level == :error
            throw(RequestError(error_msg))
        elseif level == :warn
            @warn error_msg
        end
    end
end

"""
Flush a connection attached to a request `expr`.

The connection is taken to be the first argument of `expr`. `expr` can be a call to `XCB.@check`.

# Examples
```
julia> @macroexpand @flush xcb_unmap_window(win.conn, win.id)
quote
    xcb_unmap_window(win.conn, win.id)
    (flush)(win.conn)
end
```
"""
macro flush(expr)
    if expr.head == :macrocall && expr.args[1] == Symbol("@check")
        _expr = expr.args[3]
        expr = _expr isa QuoteNode && _expr.value ∈ [:error, :warn] ? expr.args[4] : _expr
    end
    conn = expr.args[2]
    quote
        $(esc(expr))
        $(esc(flush))($(esc(conn)))
    end
end

"""
Check the value returned by the function call `request` with [`check_request`](@ref).

The severity level (`:error` or `:warn`) can be supplied as first argument.
By default, the severity level is `:warn`.

Wraps `request` with [`check_request`](@ref). The [`Connection`](@ref) argument is taken as the first argument of the function call expression `request`. The request is transformed to be checkable, through the functions xcb_*_checked (or xcb_* if there exists a xcb_*_unchecked version). If no checkable substitute is found, an `ArgumentError` is raised.

TODO: `@macroexpand` example
"""
macro check(level, request)
    conn = request.args[2]
    request_fun = string(request.args[1])
    module_prefix = string(@__MODULE__) * "."
    has_module_prefix = startswith(request_fun, module_prefix)
    has_module_prefix && (request_fun = request_fun[5:end])

    # get checked version of the request function
    if endswith(request_fun, "_unchecked")
        request_fun_checked = replace(request_fun, "_unchecked" => "")
    elseif !endswith(request_fun, "_checked")
        if isdefined(@__MODULE__, Symbol(request_fun * "_unchecked"))
            request_fun_checked = request_fun
        else
            request_fun_checked = request_fun * "_checked"
        end
    else
        request_fun_checked = request_fun
    end
    # restore module prefix and add it to `check_request` as well if relevant
    if has_module_prefix
        request_fun_checked = module_prefix * request_fun_checked
        check_request_fun = Meta.parse(module_prefix * "check_request")
    else
        check_request_fun = :check_request
    end
    # change the request function in the original expr
    if isdefined(@__MODULE__, Symbol(replace(request_fun_checked, module_prefix => "")))
        request.args[1] = Meta.parse(request_fun_checked)
    else
        throw(ArgumentError("Function $request_fun does not have a checked version available."))
    end
    :($(esc(check_request_fun))($(esc(conn)), $(esc(request)), level=$level))
end

macro check(request) esc(:(@check(:warn, $request))) end

Base.showerror(io::IO, error::FlushError) = print("FlushError: server returned code $(error.code)")

function Base.show(io::IO, setup::xcb_setup_t)
    desc = ["Version"]
    vals = [VersionNumber(setup.protocol_major_version, setup.protocol_minor_version, setup.release_number)]
    println(io, "Setup")
    prettyprint(io, desc, vals)
end

function Base.show(io::IO, screen::xcb_screen_t)
    desc = ["White pixel value", "Black pixel value", "Width", "Height"]
    vals = [format_bignumber(screen.white_pixel), format_bignumber(screen.black_pixel), "$(screen.width_in_pixels) pixels, $(screen.width_in_millimeters / 10) cm", "$(screen.height_in_pixels) pixels, $(screen.height_in_millimeters / 10) cm"]
    println(io, "Screen details")
    prettyprint(io, desc, vals)
end

function format_bignumber(number)
    str = string(number)
    seq = map(zip(length(str):-1:1, str)) do (i, x)
        if i % 3 == 0 && i ≠ length(str)
            sep = ','
        else
            sep = ""
        end
        sep * x
    end
    join(seq)
end

prettyprint(io, desc, vals) = (description = join.(zip(rpad.(desc, 20), vals), ": "); print(io, join("    " .* description, "\n")))


text:
# # make file only for local make of the document.
# this result in errors in Travis
#
# ## To build the documentation locally
#
# ### Running inside Julia REPL:
#  cd to docs folder `cd("path to docs\\MatLang\\docs")` and run the following command:
# ```
# include("makeLocal.jl")
# ```
#
# ### Running inside OS Terminal:
# cd to docs folder using OS terminal and run the following command (julia path should be added to OS path):
# ```
# julia --color=yes makeLocal.jl
# ```

using Pkg
pkg"activate .."
push!(LOAD_PATH,"../src/")

#

using Documenter, MatLang

makedocs(;
    modules=[MatLang],
	format = Documenter.HTML(
        prettyurls = prettyurls = get(ENV, "CI", nothing) == "true",
    ),
	pages=[
        "Home" => "index.md",
		"Functions" =>
			["Functions List" => "functions.md",
			"Functions Explanations and Examples" =>
				["Entering Commands" => "functions/Entering_Commands.md",
				"Matrices and Arrays" => "functions/Matrices_and_Arrays.md",
				"Numeric Types" => "functions/Numeric_Types.md"],],
        "Development - Contribution" => "development.md",
        "Native Julia noteworthy differences from MATLAB" => "juliavsmatlab.md"
    ],
    repo="https://github.com/juliamatlab/MatLang/blob/{commit}{path}#L{line}",
    sitename="MatLang",
    authors="Amin Yahyaabadi - JuliaMatlab",
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
# deploydocs(;
#     repo="github.com/juliamatlab/MatLang.git",
# )


text:
export AbstractParticle, Particle, MagneticParticle
####################################################
## Particles
####################################################
"""
    AbstractParticle
Particle supertype.
"""
abstract type AbstractParticle{T<:AbstractFloat} end
eltype(p::AbstractParticle{T}) where {T} = T

"""
```julia
Particle(ic::Vector{T}) #where ic = [x0, y0, φ0]
Particle(x0, y0, φ0)
Particle(pos::SVector, vel::SVector)
```
Create a particle with initial conditions `x0, y0, φ0`. It propagates as
a straight line.

The field `current_cell` shows at which cell of a periodic billiard is the particle
currently located.
"""
mutable struct Particle{T<:AbstractFloat} <: AbstractParticle{T}
    pos::SVector{2,T}
    vel::SVector{2,T}
    current_cell::SVector{2,T}
    function Particle(
        pos::SVector{2,T}, vel::SVector{2,T}, cc::SVector{2,T}) where{T<:AbstractFloat}
        new{T}(pos, normalize(vel), cc)
    end
end

Base.copy(p::Particle) =
Particle(p.pos, p.vel, p.current_cell)

function Particle(ic::AbstractVector{S}) where {S<:Real}
    T = S<:Integer ? Float64 : S
    φ0 = ic[3]
    pos = SVector{2,T}(ic[1:2]); vel = SVector{2,T}(cossin(φ0)...)
    return Particle(pos, vel, SVector{2,T}(0,0))
end
Particle(x::Real, y::Real, φ::Real) = Particle(collect(promote(x,y,φ)))
Particle() = Particle(rand(), rand(), rand()*2π)
function Particle(pos::SV{T}, vel::SV{T}) where {T}
    S = T<:Integer ? Float64 : T
    return Particle(pos, vel, SVector{2,S}(0.0, 0.0))
end
show(io::IO, p::Particle{T}) where {T} =
print(io, "Particle{$T}\n",
"position: $(p.pos+p.current_cell)\nvelocity: $(p.vel)")

"""
```julia
MagneticParticle(ic::AbstractVector{T}, ω::Real) # where ic = [x0, y0, φ0]
MagneticParticle(x0, y0, φ0, ω)
MagneticParticle(pos::SVector, vel::SVector, ω)
MagneticParticle(p::AbstractParticle, ω)
```
Create a *magnetic* particle with initial conditions `x0, y0, φ0` and angular
velocity `ω`. It propagates as a circle instead of a line,
with radius `1/abs(ω)`.

The field `current_cell` shows at which cell of a periodic billiard is the particle
currently located.
"""
mutable struct MagneticParticle{T<:AbstractFloat} <: AbstractParticle{T}
    pos::SVector{2,T}
    vel::SVector{2,T}
    current_cell::SVector{2,T}
    omega::T
    r::T
    center::SVector{2, T}
    function MagneticParticle(pos::SVector{2,T}, vel::SVector{2,T},
        current_cell::SVector{2,T}, ω::T) where {T<:AbstractFloat}
        if ω==0
            throw(ArgumentError("Angular velocity of magnetic particle cannot be 0."))
        end
        r = 1/ω
        c::SV{T} = pos - r*SV{T}(vel[2], -vel[1])
        new{T}(pos, normalize(vel), current_cell, ω, abs(1/ω), c)
    end
end

function Base.getproperty(p::MagneticParticle, s::Symbol)
    if s == :ω
        return Base.getfield(p, :omega)
    else
        return Base.getfield(p, s)
    end
end

Base.copy(p::MagneticParticle{T}) where {T} =
MagneticParticle(p.pos, p.vel, p.current_cell, p.omega)

function MagneticParticle(ic::AbstractVector{T}, ω::Real) where {T<:Real}
    φ0 = ic[3]
    S = T<:Integer ? Float64 : T
    pos = SVector{2,S}(ic[1:2]); vel = SVector{2,S}(cossin(φ0)...)
    return MagneticParticle(pos, vel, SVector{2,S}(0,0), convert(S,ω))
end
function MagneticParticle(x0::Real, y0::Real, φ0::Real, ω::Real)
    a = collect(promote(x0, y0, φ0, ω))
    MagneticParticle(a[1:3], a[4])
end
MagneticParticle() = MagneticParticle([rand(), rand(), rand()*2π], 1.0)
function MagneticParticle(pos::SV{T}, vel::SV{T}, ω::T) where {T}
    S = T<:Integer ? Float64 : T
    return MagneticParticle(pos, vel, SVector{2,S}(0,0), convert(S,ω))
end

MagneticParticle(p::AbstractParticle, ω) = MagneticParticle(p.pos, p.vel, p.current_cell, ω)

show(io::IO, p::MagneticParticle{T}) where {T} =
print(io, "MagneticParticle{$T}\n",
"position: $(p.pos+p.current_cell)\nvelocity: $(p.vel)\nang. velocity: $(p.omega)")

"""
    find_cyclotron(p::MagneticParticle)
Return the center of cyclotron motion of the particle.
"""
@inline find_cyclotron(p::MagneticParticle{T}) where {T} =
    c::SV{T} = p.pos - (1/p.omega)*SV{T}(p.vel[2], -p.vel[1])

@inline cyclotron(p) = p.center, p.r
function cyclotron(p, use_cell)
    pos = use_cell ? p.pos + p.current_cell : p.pos
    return pos, p.r
end


text:
module ParametronTest

include("util.jl")
include("parameter.jl")
include("functions.jl")
include("lazyexpression.jl")
include("model.jl")
include("debug.jl")

end # module


text:
"Type stored in `Compound.tuples`."
const ElementTuple = Tuple{String, Int}

"Regex to match `{...}` charge string."
const CHARGEREGEX = r"{(.*)}"

"""
Stores chemical compound's elements and charge in a structured way.

!!! info
    Electron is stored as `"e"`.
"""
struct Compound
    tuples::Vector{ElementTuple}
    charge::Int
end

"""
Constructs a compound from `str`.

An element begins with an uppercase unicode letter and
ends with a lowercase unicode letter or a unicode symbol.

!!! info
    An element can also begin with a symbol if
    the symbol is the first character (e.g. `"⬡H"`).

Parsing is insensitive to whitespace and underscores (`_`),
but also to state symbols (`(s)`, `(l)`, `(g)`, `(aq)`).
Special parsing is implemented for:
- parens (e.g. `"(CH3COO)2Mg"`)
- compounds with `"*"` (e.g. `"CuSO4 * 5H2O"`)
- electrons (`"e"`)
Charge is in the form `"{±n}"` or `"{n±}"`.
It is automatically deduced for electron (`"e"`).

# Examples
```jldoctest
julia> Compound("H2O(s)")
H2O

julia> Compound("H3O{+}")
H3O{+}

julia> Compound("(CH3COO)2Mg")
C4H6O4Mg

julia> Compound("CuSO4 * 5H2O")
CuSO9H10

julia> Compound("⬡Cl")
⬡Cl
```
"""
function Compound(str::AbstractString)
    str = replace(str, [' ', '_'] => "")
    str = replace(str, r"\((s|l|g|aq)\)$" => "") # remove state symbols
    Compound(_elementtuples(str), _charge(str))
end

"Extracts element tuples from compound's string."
function _elementtuples(str::AbstractString)
    str = replace(str, CHARGEREGEX => "")
    if str ∈ ("", "e")
        return [("e", 1)]
    end

    tuples = ElementTuple[]

    # Add 1 to elements and parens without a coefficient
    str = replace(str,  r"(?<x>\p{L}|\p{S}|\))(?=(\p{Lu}|\(|\)|\*|$))" => s"\g<x>1")

    # Expand parens
    for substr ∈ eachmatch(r"\(((\p{L}|\d)+)\)(\d+)", str)
        capture = replace(
            substr.captures[1],
            isdigit => x -> string(parse(Int, substr.captures[3]) * parse(Int, x))
        )
        str = replace(str, substr.match => capture)
    end

    # Expand compounds with '*', e.g. CuSO4*5H2O
    str = split(str, '*')
    for (i, substr) ∈ enumerate(str)
        if isdigit(substr[1])
            k, substr = match(r"(^\d+)(.+)", substr).captures
            substr = replace(substr, isdigit => x -> string(parse(Int, k) * parse(Int, x)))
            str[i] = substr
        end
    end
    str = join(str)

    str = split(str, r"(?=\p{Lu})")
    for element ∈ str
        element, k = split(element, r"(?=\d)(?<!\d)")
        index = findfirst(x -> x[1] == element, tuples)

        if isnothing(index)
            push!(tuples, (element, parse(Int, k)))
        else
            tuples[index] = (element, tuples[index][2] + parse(Int, k))
        end
    end

    return tuples
end

"Extracts charge from compound's string into a number of specified type."
function _charge(str::AbstractString)
    if str == "e"
        return -1
    end

    strmatch = match(CHARGEREGEX, str)
    if isnothing(strmatch)
        return 0
    else
        str = strmatch.captures[1]
        if str ∈ ("-", "+")
            str *= "1"
        elseif str[end] ∈ ('-', '+')
            str = str[end] * str[1:end-1]
        end
        return parse(Int, str)
    end
end

"""
Constructs a compound with `cc"str"` syntax, instead of `Compound(str)`.

# Examples
```jldoctest
julia> cc"H3O{+1}"
H3O{+}
```
"""
macro cc_str(str) Compound(str) end

"""
Checks whether two compounds are chemically equal.

# Examples
```jldoctest
julia> cc"MgOHOH" == cc"Mg(OH)2"
true
```
"""
function Base.:(==)(compound_1::Compound, compound_2::Compound)
    return sort(compound_1.tuples) == sort(compound_2.tuples) &&
        compound_1.charge == compound_2.charge
end

"""
Creates a string to represent the compound.

All elements are displayed only once (e.g. `"H2O"` and not `"HOH"`),
in the order in which they were originally given (e.g. `"MgO2H2"` from `cc"Mg(OH)2"`),
with coefficients equal to 1 not displayed (e.g. `"H"` and not `"H1"`).

# Examples
```jldoctest
julia> string(cc"CuSO4 * 5 H2O")
"CuSO9H10"
```
"""
function Base.string(compound::Compound)
    if compound == cc"e"
        return "e"
    end

    str = ""
    for (element, k) ∈ compound.tuples
        str *= element
        if k ≠ 1
            str *= string(k)
        end
    end
    if hascharge(compound)
        str *= "{"
        if compound.charge > 0
            str *= "+"
        end
        if compound.charge == -1
            str *= "-"
        elseif compound.charge ≠ 1
            str *= string(compound.charge)
        end
        str *= "}"
    end
    return str
end

"Displays the compound using [`Base.string(::Compound)`](@ref)."
function Base.show(io::IO, compound::Compound)
    print(io, string(compound))
end

"""
Returns compound's elements as strings.

# Examples
```jldoctest
julia> elements(cc"CH3COOH")
3-element Array{String,1}:
 "C"
 "H"
 "O"
```
"""
function elements(compound::Compound)
    return first.(compound.tuples)
end

"True if the compound's charge is nonzero."
function hascharge(compound::Compound)
    return compound.charge ≠ 0
end


text:
################################################################################
#
#          EllCrv/Finite.jl : Elliptic curves over finite fields
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# (C) 2016 Tommy Hofmann
# (C) 2016 Robin Ammon
# (C) 2016 Sofia Brenner
#
################################################################################

export hasse_interval, order, order_via_bsgs, order_via_legendre,
       order_via_schoof, rand, elem_order_bsgs

################################################################################
#
#  Random point
#
################################################################################

doc"""
***
    rand(E::EllCrv) -> EllCrvPt
> Returns a random point on the elliptic curve $E$ defined over a finite field.
> It is assumed that $E$ is given in short form.
"""
# only works for short form
function rand(E::EllCrv)
  R = base_field(E)

  if E.short == false
    error("does not work for long form")
  end

  while true
  # choose random x-coordinate and check if it is a square in F_q
  # if not, choose new x-coordinate
    x = rand(R)
    square = x^3 + E.coeff[1]*x + E.coeff[2]

    a = issquare(square)
    if a[1] == true # square is a square in F_q, so have found point on the curve
      y = a[2]
      P = E([x, y])
      return P
    end
  end
end

################################################################################
#
# Order via Legendre symbol
#
################################################################################

doc"""
***
    order_via_legendre(E::EllCrv{Generic.Res{fmpz}) -> fmpz

Calculates the number of points on an elliptic curve $E$ over a finite field
$\mathbf Z/p\mathbf Z$ using the Legendre symbol. It is assumed that $p$ is
prime.
"""
# Th. 4.14
function order_via_legendre(E::EllCrv{Generic.Res{fmpz}})
  R = base_field(E)
  p = characteristic(R)
  q = order(R)
  grouporder = FlintZZ(0)

  p == 0 && error("Base field must be finite")

  if p != q
    error("Finite field must have degree 1")
  end

  if E.short == false
    E = short_weierstrass_model(E)[1]
  end

  x = FlintZZ(0)

  while x < p
    C = x^3 + (E.coeff[1])*x + (E.coeff[2])
    Cnew = C.data # convert to fmpz
    a = jacobi(Cnew, p) # can be used to compute (C/F_p) since p prime
    grouporder = grouporder + a
    x = x + 1
  end

  grouporder = grouporder + p + 1
  return grouporder
end

################################################################################
#
#  Hasse inverval
#
################################################################################

doc"""
***
    hasse_interval(E::EllCrv) -> Array{fmpz, 1}

> Given an elliptic curve $E$ over a finite field $\mathbf F$, returns an array
> `[l, b]` > of integers, such that $l \leq \#E(\mathbf F) \leq b$ using
> Hasse's theorem.
"""
function hasse_interval(E::EllCrv)
  R = base_field(E)
  characteristic(R) == 0 && error("Base field must be finite")
  q = order(R)
  s = isqrtrem(q)[1] # sqrt(q) rounded down

  l = q + 1 - 2*(s + 1)
  b = q + 1 + 2*(s + 1)

  return [l, b]
end

doc"""
***
    elem_order_bsgs(P::EllCrvPt) -> fmpz

> Calculates the order of a point P on an elliptic curve given over a finite
> field using BSGS.
"""
# section 4.3.4
function elem_order_bsgs(P::EllCrvPt)
  R = base_field(P.parent)
  p = characteristic(R)
  p == 0 && error("Base field must be finite")
  q = order(R) # R = F_q

  # step 1
  Q = Int(q + 1) * P

  # step 2
  m = Int( ceil(Int(q)^(1//4)) )

  list_points = []
  for j = 0:m
    S = j*P
    push!(list_points, S)
  end

  # step 3
  k = -m
  H = (2*m)*P
  M = FlintZZ(0) # initialize M, so that it is known after the while loop

  while k < m + 1
    Snew = Q + (k*H)

    i = 1
    while i < m + 2 # is there a match?
      if Snew == list_points[i]
        M = q + 1 + (2*m*k) - (i-1)

        if M != 0
          i = m + 2 # leave both while loops
          k = m + 1
        else
          i = i + 1 # search on if M == 0
        end

      elseif Snew == -(list_points[i])
        M = q + 1 + (2*m*k) + (i-1) # step 4

        if M != 0
          i = m + 2 # leave both while loops
          k = m + 1
        else
          i = i + 1 # search on if M == 0
        end
      else
        i = i + 1
      end
    end

    k = k + 1
  end

  # step 5
  gotostep5 = true
  while gotostep5
    fac = factor(M)
    primefactors = []
    for i in fac
      push!(primefactors, i[1])
    end
    r = length(primefactors)

    # step 6
    i = 1
    while i < (r + 1)
      T = Int(divexact(M, primefactors[i]))*P
      if T.isinfinite == true
        M = divexact(M, primefactors[i])
        i = r + 2  # leave while-loop
      else
        i = i + 1
      end
    end

    if i == r + 1  # then (M/p_i)P != oo for all i
      gotostep5 = false
    end
  end

  return M
end

doc"""
***
    order(P::EllCrvPt) -> fmpz

Given a point on an elliptic curve over a finite field, returns the order
of this point.
"""
order(P::EllCrvPt) = elem_order_bsgs(P)

################################################################################
#
#  Order via BSGS
#
################################################################################

doc"""
***
    order_via_bsgs(E::EllCrv) -> Array{fmpz, 1}

Calculates candidates for the number of points on an elliptic curve $E$ given
over a finite field $\mathbf F_q$, using the baby step giant step method. If
$q$ prime, $q > 229$, then the order is determined uniquely by this algorithm.
It is assumed that the characteristic is not 2.
"""
function order_via_bsgs(E::EllCrv)
  R = base_field(E)
  p = characteristic(R)
  p == 0 && error("Base field must be finite")

  q = order(R)

  if (p == 2)
    error("Characteristic must not be 2")
  end

  if E.short == false
    E = short_weierstrass_model(E)[1]
  end

  Nposs = FlintZZ(1)
  h = hasse_interval(E)
  l = h[1]
  b = h[2]

  counter = 0
  runwhile = true

  # if Nposs is greater than interval length, there is only one multiple of
  # Nposs in the interval, so stop
  # Also, if lcm does not change in 10(?) consecutive loops, leave while loop
  while (Nposs <= (b-l)) && (runwhile == true)

    P = rand(E)
    ord = elem_order_bsgs(P)

    if lcm(ord, Nposs) == Nposs
      counter = counter + 1
    else
      Nposs = lcm(ord, Nposs)
      counter = 0
    end

    if counter > 10 # maybe take another criterion
      runwhile = false
    end

  end

  if runwhile == false # could not determine group order uniquely
    candidates = fmpz[]
    Ncand = divrem(l, Nposs)[1]*Nposs
    if Ncand != 0
      push!(candidates, Ncand)
    end

    Ncand = Ncand + Nposs

    while Ncand < b # compute all possible group orders using Hasse's theorem
      push!(candidates, Ncand)
      Ncand = Ncand + Nposs
    end

    output = candidates

  else # group order is determined
    N = (divrem(l+1, Nposs)[1] + 1) * Nposs
    output = [N]
  end

  if (p == q) && (p > 229) && (length(output) != 1)
  # group order is uniquely determined, but by quadratic twist

    d = R(0)
    boolie = true
    while boolie # get a quadratic non-residue mod p
      d = rand(R)
      if issquare(d)[1] == false
        boolie = false
      end
    end

    Eprime = EllipticCurve([E.coeff[1]*d^2, E.coeff[2]*d^3]) # quadratic twist
    bb = order_via_bsgs(Eprime)[1]
    output = [2*p + 2 - bb]
  end

  return output
end

################################################################################
#
#  Schoof's algorithm
#
################################################################################

function fn_from_schoof(E::EllCrv, n::Int, x)

  R = base_field(E)
  S, y = PolynomialRing(parent(x),"y")

  f = psi_poly_field(E, n, x, y)

 # println("f: $f, $(degree(f))")

  g = x^3 + E.coeff[1]*x + E.coeff[2]

  if isodd(n)
    return replace_all_squares(f, g)
  else
    return replace_all_squares(divexact(f, y), g)
  end
end

doc"""
***
    order_via_schoof(E::EllCrv) -> fmpz

Given an elliptic curve $E$ elliptic curve given a finite field $\mathbf F$,
this function computes the order of $E(\mathbf F)$ using Schoof's algorithm
The characteristic must not be $2$ or $3$.
"""
function order_via_schoof(E::EllCrv)
  R = base_field(E)
  q = order(R)
  p = characteristic(R)

  if (p == 2) || (p == 3) || (p == 0)
    error("Characteristic must not be 2 or 3 or 0")
  end

  if E.short == false
    E = short_weierstrass_model(E)[1]
  end

  # step 1: get suitable set S of prime numbers
  sqrt_q = isqrtrem(q)[1]
  S = prime_set(4*(sqrt_q + 1), p)

  L = length(S)
  product = 1

  # step 2: compute t (mod l) for every l in S
  t_mod_l = []
  for i = 1:L
    a = S[i]
    push!(t_mod_l, t_mod_prime(S[i], E))
    product = product * S[i]
  end

  # step 3: use chinese remainder theorem to get value of t
  t = 0
  for i = 1:L
    n_i = div(product, S[i])
    B = ResidueRing(FlintZZ, S[i])
    M_i = inv(B(n_i))
    M_i = M_i.data
    t = t + (M_i * n_i * t_mod_l[i])
  end

  t = mod(t, product)
  if t > product//2
    t = t - product
  end

  return (q + 1 - t)::fmpz
end

#prime_set(M::Nemo.fmpz, char::Nemo.fmpz) -> Array{Nemo.fmpz}
#  returns a set S of primes with:
# 1) char not contained in S
# 2) product of elements in S is greater than M
function prime_set(M, char)
  S = Nemo.fmpz[]

  p = 1
  product = 1

  while product < M
    p = next_prime(p)

    if p != char
      push!(S,p)
      product = product * p
    end
  end

  return S
end

# t_mod_prime(l::Nemo.fmpz, E::EllCrv) -> Nemo.fmpz
# determines the value of t modulo some prime l (used in Schoof's algorithm)
function t_mod_prime(l, E)
  R = base_field(E)
  q = order(R)
  q_int = Int(q)
  l = Int(l)

  S, x = PolynomialRing(R, "x")
  T, y = PolynomialRing(S, "y")
  Z = ResidueRing(FlintZZ, l)

  f = x^3 + E.coeff[1]*x + E.coeff[2]
  fl = fn_from_schoof(E, l, x)
  U = ResidueRing(S, fl)

  PsiPoly = [] # list of psi-polynomials
  for i = -1:(l + 1)
    push!(PsiPoly, psi_poly_field(E,i,x,y)) # Psi[n] is now found in PsiPoly[n+2]
  end

  Fnschoof = [] # list of the fn- functions # Psi[n] is now found in PsiPoly[n+2]
  for i = -1:(l + 1)
    push!(Fnschoof, fn_from_schoof(E,i,x))
  end

  # case where l == 2. value of t mod l determined by some gcd, see p. 124
  if l == 2
    x_q = powmod(x, q_int, f)
    ggt = gcd(f, x_q - x)
    if ggt == 1
      t = FlintZZ(1)
    else
      t = FlintZZ(0)
    end

    return t
  end

  # case where l != 2
  k = mod(q, l) # reduce q mod l
  k_mod = Z(k)

  fk = Fnschoof[k+2]
  fkme = Fnschoof[k+1]
  fkpe = Fnschoof[k+3]
  fkpz = Fnschoof[k+4]

  # is there a nonzero P = (x,y) in E[l] with phi^2(P) == +-kP ?
  if mod(k,2) == 0
    g = U( (powmod(x, q_int^2, fl) - x) * fk^2 * f + fkme * fkpe).data
    ggT = gcd(g, fl)
  else
    g = U( (powmod(x, q_int^2, fl) - x) * fk^2 + fkme * fkpe * f).data
    ggT = gcd(g, fl)
  end

  if ggT != 1 # case 1
    if jacobi(FlintZZ(k), FlintZZ(l)) == -1
      return FlintZZ(0)
    else
      # need square root of q (mod l)
      w = issquare(k_mod)[2]
      if w.data < 0
        w = w + l
      end
      w_int = Int(w.data)

      fw = Fnschoof[w_int+2]
      fwme = Fnschoof[w_int+1]
      fwpe = Fnschoof[w_int+3]

      if mod(w_int, 2) == 0
        g = U((powmod(x,q_int,fl) - x) * fw^2*f + fwme*fwpe).data # reduce mod fl
        ggT = gcd(g, fl)
      else
        g = U((powmod(x,q_int,fl) - x) * fw^2 + fwme*fwpe*f).data
        ggT = gcd(g, fl)
      end

      if ggT == 1
        return FlintZZ(0)
      else
        fwmz = Fnschoof[w_int]
        fwpz = Fnschoof[w_int+4]

        if mod(w_int, 2) == 0
          g = U(4 * powmod(f,div(q_int + 3, 2),fl)*fw^3 - (fwpz * fwme^2) + (fwmz*fwpe^2)).data
          ggT2 = gcd(g, fl)
        else
          g = U(4 * powmod(f,div(q_int - 1, 2),fl)*fw^3 - (fwpz * fwme^2) + (fwmz*fwpe^2)).data
          ggT2 = gcd(g, fl)
        end

        if ggT2 == 1
          return -2*(w.data)
        else
          return 2*(w.data)
        end
      end
    end
  else # case 2
    Fkmz = PsiPoly[k]
    Fkme = PsiPoly[k+1]
    Fk = PsiPoly[k+2]
    Fkpe = PsiPoly[k+3]
    Fkpz = PsiPoly[k+4]

    alpha = Fkpz*psi_power_mod_poly(k-1,E,x,y,2,fl) - Fkmz*psi_power_mod_poly(k+1,E,x,y,2,fl) - 4*powmod(f,div(q_int^2+1,2),fl)*psi_power_mod_poly(k,E,x,y,3,fl)
    beta = ((x - powmod(x, (q_int^2), fl))*psi_power_mod_poly(k,E,x,y,2,fl)- Fkme*Fkpe)*4*y*Fk

    tau = 1
    while tau < l

      ftaumz = PsiPoly[tau]
      ftaume = PsiPoly[tau+1]
      ftau = PsiPoly[tau+2]
      ftaupe = PsiPoly[tau+3]
      ftaupz = PsiPoly[tau+4]

      fntaumz = Fnschoof[tau]
      fntaume = Fnschoof[tau+1]
      fntaupe = Fnschoof[tau+3]
      fntaupz = Fnschoof[tau+4]

      gammahelp = powmod(fntaupz*fntaume^2- fntaumz * fntaupe^2,q_int,fl)

      if mod(tau, 2) == 0
        gamma = y * powmod(f,div(q_int-1,2),fl) * gammahelp
      else
        gamma = powmod(f,q_int,fl) * gammahelp
      end

      monster1 = ((Fkme*Fkpe - psi_power_mod_poly(k,E,x,y,2,fl)*(powmod(x, q_int^2, fl) + powmod(x, q_int, fl) + x)) * beta^2 + psi_power_mod_poly(k,E,x,y,2,fl)*alpha^2) * psi_power_mod_poly(tau,E,x,y,2*q_int,fl) + psi_power_mod_poly(tau-1, E, x,y,q_int,fl)*psi_power_mod_poly(tau+1, E, x,y,q_int, fl)*beta^2*psi_power_mod_poly(k,E,x,y,2,fl)

      if divrem(degree(monster1), 2)[2] == 1
        monster1 = divexact(monster1, y)
      end

      monster1_1 = replace_all_squares_modulo(monster1, f, fl)
      monster1_2 = U(monster1_1).data # monster1_1 reduced

      if monster1_2 != 0
        tau = tau + 1
      else
        monster2 = 4*y*powmod(f,div(q_int-1,2),fl)*psi_power_mod_poly(tau,E,x,y,3*q_int,fl) * (alpha * (((2*powmod(x, q_int^2, fl) + x) * psi_power_mod_poly(k,E,x,y,2,fl) - Fkme*Fkpe )*beta^2 - alpha^2*psi_power_mod_poly(k,E,x,y,2,fl)) - y*powmod(f,div(q_int^2-1,2),fl)*beta^3 * Fk^2) - beta^3*psi_power_mod_poly(k,E,x,y,2,fl)*gamma

        if divrem(degree(monster2), 2)[2] == 1
          monster2 = divexact(monster2, y)
        end

        monster2_1 = replace_all_squares_modulo(monster2, f,fl)
        monster2_2 = U(monster2_1).data # monster2_1 mod fl

        if monster2_2 != 0
          tau = tau + 1
        else
          return tau
        end
      end
    end
  end
end

# computes psi_n^power mod g
function psi_power_mod_poly(n, E, x, y, power, g)

    fn = fn_from_schoof(E, n, x)
    f = x^3 + E.coeff[1]*x + E.coeff[2]
    p = powmod(fn,power,g)

    if mod(n, 2) == 0
        if mod(power, 2) == 0
            p1 = powmod(f, div(power, 2), g)
        else
            p1 = powmod(f, div(power - 1, 2), g) * y
        end
    else p1 = 1
    end

    return p * p1
end

#special_order2(E::EllCrv{NemoResidue}) -> Nemo.fmpz
#> counts points on an elliptic curve E given over F_2
function _special_order2(E)
  R = base_field(E) # should be Z/2Z
  ord = FlintZZ(1)

  for i = 0:1
    for j = 0:1
      if ison_curve(E, [R(i), R(j)])
        ord = ord + 1
      end
    end
  end

  return ord
end

#special_order3(E::EllCrv{NemoResidue}) -> Nemo.fmpz
# counts points on an elliptic curve E given over F_3
function _special_order3(E)
  R = base_field(E) # should be Z/3Z
  ord = FlintZZ(1)

  for i = 0:2
    for j = 0:2
      if ison_curve(E, [R(i), R(j)])
        ord = ord + 1
      end
    end
  end

  return ord
end

################################################################################
#
#  Point counting
#
################################################################################

doc"""
***
    order(E::EllCrv{NemoResidue}) -> Nemo.fmpz

Given an elliptic curve $E$ over a finite field $\mathbf F$, computes
$\#E(\mathbf F)$.
"""
function order(E::EllCrv)
  R = base_field(E)
  p = characteristic(R)
  q = order(R)

  p == 0 && error("Characteristic must be nonzero")

  # char 2
  if p == 2
    if q > 2
      error("Don't have algorithm for char = 2 and not F_2") # legendre is the only algorithm that can deal with char = 2, but q must be equal to p
    else
      return _special_order2(E)
    end
  end

  # char 3
  if p == 3
    if q > 3
      error("Don't have algorithm for char = 3 and not F_3")
    else
      return _special_order3(E)
    end
  end

  A = order_via_bsgs(E)
  if length(A) == 1
    return A[1]
  end
  return order_via_schoof(E) # bsgs may only return candidate list
end


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

# An exception type used internally to signal that an unsatisfiable
# constraint was detected
struct UnsatError <: Exception
    p0::Int
end

# Some parameters to drive the decimation process
mutable struct MaxSumParams
    dec_interval::Int # number of iterations between decimations
    dec_fraction::Float64 # fraction of nodes to decimate at every decimation
                          # step

    function MaxSumParams()
        accuracy = parse(Int, get(ENV, "JULIA_PKG_RESOLVE_ACCURACY",
                                  # Allow for `JULIA_PKGRESOLVE_ACCURACY` for backward
                                  # compatibility with Julia v1.7-
                                  get(ENV, "JULIA_PKGRESOLVE_ACCURACY", "1")))
        accuracy > 0 || error("JULIA_PKG_RESOLVE_ACCURACY must be > 0")
        dec_interval = accuracy * 5
        dec_fraction = 0.05 / accuracy
        return new(dec_interval, dec_fraction)
    end
end

# Messages has the cavity messages and the total fields, and
# gets updated iteratively (and occasionally decimated) until
# convergence
mutable struct Messages
    # cavity incoming messages: for each package p0,
    #                           for each neighbor p1 of p0,
    #                           msg[p0][p1] is a vector of length spp[p0]
    #                           messages are normalized (i.e. the max is always 0)
    msg::Vector{Vector{Field}}

    # overall fields: for each package p0,
    #                 fld[p0] is a vector of length spp[p0]
    #                 fields are not normalized
    fld::Vector{Field}

    # backup of the initial value of fld, to be used when resetting
    initial_fld::Vector{Field}

    function Messages(graph::Graph)
        np = graph.np
        spp = graph.spp
        gadj = graph.gadj
        gconstr = graph.gconstr
        req_inds = graph.req_inds
        ignored = graph.ignored
        pvers = graph.data.pvers
        pdict = graph.data.pdict

        ## generate wveights (v0 == spp[p0] is the "uninstalled" state)
        vweight = [[VersionWeight(v0 < spp[p0] ? pvers[p0][v0] : v"0") for v0 = 1:spp[p0]] for p0 = 1:np]

        # external fields: favor newest versions over older, and no-version over all;
        #                  explicit requirements use level l1 instead of l2
        fv(p0, v0) = p0 ∈ req_inds ?
            FieldValue(0, vweight[p0][v0],     zero(VersionWeight), (v0==spp[p0])) :
            FieldValue(0, zero(VersionWeight), vweight[p0][v0],     (v0==spp[p0]))
        fld = [[fv(p0, v0) for v0 = 1:spp[p0]] for p0 = 1:np]

        initial_fld = [copy(f0) for f0 in fld]

        # allocate cavity messages
        msg = [[Field(undef, spp[p0]) for j1 = 1:length(gadj[p0])] for p0 = 1:np]

        msgs = new(msg, fld, initial_fld)

        reset_messages!(msgs, graph)

        return msgs
    end
end

# enforce constraints and normalize fields;
# zero out all cavity messages
function reset_messages!(msgs::Messages, graph::Graph)
    msg = msgs.msg
    fld = msgs.fld
    initial_fld = msgs.initial_fld
    np = graph.np
    spp = graph.spp
    gconstr = graph.gconstr
    ignored = graph.ignored
    for p0 = 1:np
        ignored[p0] && continue
        map(m->fill!(m, zero(FieldValue)), msg[p0])
        copyto!(fld[p0], initial_fld[p0])
        gconstr0 = gconstr[p0]
        for v0 = 1:spp[p0]
            gconstr0[v0] || (fld[p0][v0] = FieldValue(-1))
        end
        fld[p0] .-= maximum(fld[p0])
    end
    return msgs
end

mutable struct SolutionTrace
    # the solution is built progressively by a decimation process
    solution::Vector{Int}
    num_nondecimated::Int

    best::Vector{Int}
    staged::Union{Tuple{Int,Int},Nothing}

    function SolutionTrace(graph::Graph)
        np = graph.np
        solution = zeros(Int, np)
        num_nondecimated = np

        best = zeros(Int, np)
        staged = nothing

        return new(solution, num_nondecimated, best, staged)
    end
end

function update_solution!(strace::SolutionTrace, graph::Graph)
    np = graph.np
    ignored = graph.ignored
    gconstr = graph.gconstr
    solution = strace.solution
    best = strace.best
    nnd = np
    fill!(solution, 0)
    for p0 in findall(ignored)
        s0 = findfirst(gconstr[p0])
        solution[p0] = s0
        nnd -= 1
    end
    strace.num_nondecimated = nnd
    if nnd ≤ sum(best .== 0)
        copyto!(best, solution)
        strace.staged = nothing
        return true
    else
        return false
    end
end

# This is the core of the max-sum solver:
# for a given node p0 (i.e. a package) updates all
# input cavity messages and fields of its neighbors
function update!(p0::Int, graph::Graph, msgs::Messages)
    gadj = graph.gadj
    gmsk = graph.gmsk
    adjdict = graph.adjdict
    ignored = graph.ignored
    spp = graph.spp
    np = graph.np
    cavfld = graph.cavfld
    newmsg = graph.newmsg
    diff = graph.diff

    msg = msgs.msg
    fld = msgs.fld

    maxdiff = zero(FieldValue)

    gadj0 = gadj[p0]
    msg0 = msg[p0]
    fld0 = fld[p0]
    spp0 = spp[p0]
    adjdict0 = adjdict[p0]

    # iterate over all neighbors of p0
    for j0 in 1:length(gadj0)
        p1 = gadj0[j0]
        ignored[p1] && continue
        j1 = adjdict0[p1]
        #@assert j0 == adjdict[p1][p0]
        bm1 = gmsk[p1][j1]
        spp1 = spp[p1]
        msg1 = msg[p1]

        # compute the output cavity field p0->p1
        resize!(cavfld, length(fld0))
        cavfld .= fld0 .- msg0[j0]

        # keep the old input cavity message p0->p1
        oldmsg = msg1[j1]

        # init the new message to minus infinity
        resize!(newmsg, spp1)
        fill!(newmsg, FieldValue(-1))

        # compute the new message by passing cavfld
        # through the constraint encoded in the bitmask
        # (equivalent to:
        #    newmsg = [maximum(cavfld[bm1[:,v1]]) for v1 = 1:spp1]
        # )
        # This is hot code for the resolver
        @inbounds for v1 = 1:spp1, v0 = 1:spp0
            bm1[v0, v1] || continue
            newmsg[v1] = max(newmsg[v1], cavfld[v0])
        end
        m = maximum(newmsg)
        validmax(m) || return Unsat(p0) # No state available without violating some
                                        # hard constraint

        # normalize the new message
        @inbounds for i in 1:length(newmsg)
            newmsg[i] -= m
        end

        resize!(diff, spp1)
        diff .= newmsg .- oldmsg
        maxdiff = max(maxdiff, maximum(abs, diff))

        # update the field of p1
        fld[p1] .+= diff

        # put the newly computed message in place
        copy!(msg1[j1], newmsg)
    end
    return maxdiff
end

# A simple shuffling machinery for the update order in iterate!()
# (wouldn't pass any random quality test but it's arguably enough)
mutable struct NodePerm
    p::Vector{Int}
    step::Int64
    NodePerm(np::Integer) = new(collect(1:np), 1)
end

function Random.shuffle!(perm::NodePerm)
    p = perm.p
    for j = length(p):-1:2
        k = perm.step % j + 1
        p[j], p[k] = p[k], p[j]
        perm.step += isodd(j) ? 1 : k
    end
    #@assert isperm(p)
end

Base.iterate(perm::NodePerm, state...) = iterate(perm.p, state...)

# Call update for all nodes (i.e. packages) in
# random order
function iterate!(graph::Graph, msgs::Messages, perm::NodePerm)
    maxdiff = zero(FieldValue)
    shuffle!(perm)
    for p0 in perm
        graph.ignored[p0] && continue
        maxdiff0 = update!(p0, graph, msgs)
        if maxdiff0 isa Unsat
            return maxdiff0
        end
        maxdiff = max(maxdiff, maxdiff0)
    end
    return maxdiff
end

function decimate1!(p0::Int, graph::Graph, strace::SolutionTrace, msgs::Messages)
    solution = strace.solution
    fld = msgs.fld
    adjdict = graph.adjdict
    gmsk = graph.gmsk
    gconstr = graph.gconstr

    @assert solution[p0] == 0
    @assert !graph.ignored[p0]
    fld0 = fld[p0]
    s0 = argmax(fld0)
    # only do the decimation if it is consistent with
    # the constraints...
    gconstr[p0][s0] || return 0
    # ...and with the previously decimated nodes
    for p1 in findall(solution .> 0)
        haskey(adjdict[p0], p1) || continue
        s1 = solution[p1]
        j1 = adjdict[p0][p1]
        gmsk[p1][j1][s0,s1] || return 0
    end
    solution[p0] = s0
    strace.num_nondecimated -= 1
    return s0
end

function decimate!(graph::Graph, strace::SolutionTrace, msgs::Messages, n::Integer)
    np = graph.np
    gconstr = graph.gconstr
    ignored = graph.ignored
    fld = msgs.fld

    @assert n ≥ 1
    dtrace = Tuple{Int,Int}[]
    dec = 0

    fldorder = sort(findall(.!(ignored)), by=p0->secondmax(fld[p0], gconstr[p0]))
    for p0 in fldorder
        s0 = decimate1!(p0, graph, strace, msgs)
        s0 == 0 && continue
        push!(dtrace, (p0,s0))
        dec += 1
        dec == n && break
    end

    return dtrace
end

function clean_forbidden!(graph::Graph, msgs::Messages)
    np = graph.np
    gconstr = graph.gconstr
    ignored = graph.ignored
    fld = msgs.fld
    affected = Tuple{Int,Int}[]
    removed = 0

    for p0 = 1:np
        ignored[p0] && continue
        fld0 = fld[p0]
        gconstr0 = gconstr[p0]
        for v0 in findall(gconstr0)
            validmax(fld0[v0]) && continue
            push!(affected, (p0,v0))
            removed += 1
        end
    end
    return affected
end

# Iterative solver: run iterate!() until convergence
# (occasionally calling decimate())
function maxsum(graph::Graph)
    params = MaxSumParams()

    perm = NodePerm(graph.np)
    strace = SolutionTrace(graph)
    msgs = Messages(graph)

    push_snapshot!(graph)
    # gconstr_sav = graph.gconstr
    # ignored_sav = graph.ignored
    ok = converge!(graph, msgs, strace, perm, params)
    # @assert graph.gconstr ≡ gconstr_sav
    # @assert graph.ignored ≡ ignored_sav
    pop_snapshot!(graph)
    if ok
        @assert strace.best == strace.solution
        @assert strace.num_nondecimated == 0
        @assert all(strace.solution .> 0)
        @assert strace.staged ≡ nothing
    else
        @assert strace.staged ≢ nothing
    end
    return ok, strace.best, strace.staged
end

function try_simplify_graph_soft!(graph, sources)
    try
        simplify_graph_soft!(graph, sources, log_events = false)
    catch err
        err isa ResolverError || rethrow()
        return false
    end
    return true
end

struct Unsat
    p0::Int
end

function converge!(graph::Graph, msgs::Messages, strace::SolutionTrace, perm::NodePerm, params::MaxSumParams)
    is_best_sofar = update_solution!(strace, graph)

    # this is the base of the recursion: the case when
    # the solver has succeeded in decimating everything
    strace.num_nondecimated == 0 && return true

    reset_messages!(msgs, graph)

    # perform some maxsum iterations, then decimate one node.
    # If failure happens during this process, we bail (return false)
    it = 0
    for it = 1:params.dec_interval
        maxdiff = iterate!(graph, msgs, perm)
        if maxdiff isa Unsat
            if is_best_sofar
                p0 = maxdiff.p0
                s0 = findlast(graph.gconstr[p0])
                strace.staged = (p0, s0)
            end
            return false
        end
        maxdiff == zero(FieldValue) && break
    end

    # if maxsum has found some forbidden states, remove
    # them and propagate the effect
    affected = clean_forbidden!(graph, msgs)

    isempty(affected) && @goto decimate

    sources = Set{Int}()
    for (p0,v0) in affected
        graph.gconstr[p0][v0] = false
        push!(sources, p0)
    end

    if !try_simplify_graph_soft!(graph, sources)
        # found an implicit contradiction
        is_best_sofar && (strace.staged = first(affected))
        return false
    end
    return converge!(graph, msgs, strace, perm, params)

    @label decimate

    ndec = max(1, round(Int, params.dec_fraction * strace.num_nondecimated))
    dtrace = decimate!(graph, strace, msgs, ndec)
    if isempty(dtrace)
        # decimation has failed, all candidate states are forbidden
        # (which shouldn't really happen, this is a failsafe)
        if is_best_sofar
            # pick the first decimation candidate
            smx(p1) = secondmax(msgs.fld[p1], graph.gconstr[p1])
            p0 = reduce((p1,p2)->(smx(p1)≤smx(p2) ? p1 : p2), findall(.!(graph.ignored)))
            s0 = argmax(fld[p0])
            strace.staged = dec_firstcandidate(graph, msgs)
        end
        return false
    end

    while true
        # decimation has succeeded, at least nominally.
        # We need to propagate its effects and check for contradictions

        lentr = length(dtrace)
        if lentr == 1 && is_best_sofar
            strace.staged = dtrace[1]
        end

        push_snapshot!(graph)
        # info("setting dtrace=$dtrace")
        for (p0,s0) in dtrace
            @assert !graph.ignored[p0]
            @assert graph.gconstr[p0][s0]
            fill!(graph.gconstr[p0], false)
            graph.gconstr[p0][s0] = true
            graph.ignored[p0] = true
        end

        # if decimation has produced an implicit contradiction, backtrack
        try_simplify_graph_soft!(graph, Set{Int}(first.(dtrace))) || @goto backtrack

        # otherwise, keep going...
        converge!(graph, msgs, strace, perm, params) && (pop_snapshot!(graph); return true)

        @label backtrack

        # warn("reverting dtrace=$dtrace")

        # if we're here, the last decimation step has been proven to lead
        # to an unsat situation at some point, we need to roll it back

        # revert the state of the constraints
        pop_snapshot!(graph)

        lentr == 1 && break
        # halve the dtrace
        deleteat!(dtrace, ((lentr÷2)+1):lentr)
    end

    @assert length(dtrace) == 1
    p0, s0 = dtrace[1]

    is_best_sofar && @assert strace.staged ≢ nothing

    # forbid the state used in the last attempt
    # (note that we're working on the "entry" snapshot here!)
    graph.gconstr[p0][s0] = false
    # check if we have finished all available possibilities
    any(graph.gconstr[p0]) || return false
    # if neither the s0 state nor its negation are valid, give up
    try_simplify_graph_soft!(graph, Set{Int}([p0])) || return false

    # keep going, with one possible state less...
    return converge!(graph, msgs, strace, perm, params)
end


text:
using Printf

export Grid, State

struct State
    x::Int # x-position
    y::Int # y-position
    function State(x::Int, y::Int)::State
        return new(x, y);
    end
end

@inline Base.:(==)(s1::State, s2::State)::Bool = (s1.x == s2.x) && (s1.y == s2.y);

struct Grid
    dims::Tuple{Int,Int}     # size of the grid
    null_state::State        # terminal state off the grid
    p_transition::Real       # probability of transitioning correctly
    rewards::Matrix{Float64} # reward matrix
    function Grid(dims::Tuple{Int,Int} = (10, 10);
                  null_state::State = State(-1, -1),
                  p_transition::Real = 0.70)::Grid
        # initialize empty rewards matrix
        rewards = zeros(dims)
        return new(dims, null_state, p_transition, rewards);
    end
end

function Base.getindex(g::Grid, i::Int, j::Int)::Float64
    if (1 ≤ i ≤ g.dims[1]) && (1 ≤ j ≤ g.dims[2])
        return g.rewards[i,j];
    end
    throw(BoundsError(g));
end

function Base.setindex!(g::Grid, v::Real, i::Int, j::Int)
    value = Float64(v)
    if (1 ≤ i ≤ g.dims[1]) && (1 ≤ j ≤ g.dims[2])
        g.rewards[i,j] = value;
    else
        throw(BoundsError(g));
    end
end

function Base.show(io::IO, g::Grid)::Nothing
    println(io, "\033[2J") # clear repl
    println(io, "\033[$(displaysize(stdout)[1])A")
    rows, cols = g.dims
    border = "\t" * " "^7 * "-"^(cols * 11)
    println(io, border)
    for i in rows:-1:1
        print("\t")
        for j in 1:cols
            if j == 1
                @printf "%-5i  | %5.1f   |" i g.rewards[j,i]
            else
                @printf "  %5.1f   |" g.rewards[j,i]
            end
        end
        println(io, "\n" * border)
    end
    print(io, "\t" * " "^4)
    for k in 1:cols
        @printf "    %5.0f  " k
    end
end

@inline statespace(g::Grid)::Vector{State} = [[State(x,y) for x=1:g.dims[1], y=1:g.dims[2]]..., g.null_state];
@inline stateindex(g::Grid, s::State)::Int = findall(x->x == s, statespace(g))[1];
@inline inbounds(g::Grid, s::State)::Bool = (1 ≤ s.x ≤ g.dims[1]) && (1 ≤ s.y ≤ g.dims[2]);

max_stars_repo_name: Joshua-Robison/GridWorld

max_issues_repo_name: Joshua-Robison/GridWorld

max_forks_repo_name: Joshua-Robison/GridWorld


text:
module PkgTemplates

using Dates
using InteractiveUtils
using LibGit2
using Mustache
using Pkg
using REPL.TerminalMenus
using URIParser

export
    # Template/package generation.
    Template,
    generate,
    interactive_template,
    generate_interactive,
    # Licenses.
    show_license,
    available_licenses,
    # Plugins.
    GitHubPages,
    GitLabPages,
    AppVeyor,
    TravisCI,
    GitLabCI,
    CirrusCI,
    Codecov,
    Coveralls,
    Citation

"""
A plugin to be added to a [`Template`](@ref), which adds some functionality or integration.
New plugins should almost always extend [`GenericPlugin`](@ref) or [`CustomPlugin`](@ref).
"""
abstract type Plugin end

include("licenses.jl")
include("template.jl")
include("generate.jl")
include("plugin.jl")
include(joinpath("plugins", "documenter.jl"))
include(joinpath("plugins", "coveralls.jl"))
include(joinpath("plugins", "appveyor.jl"))
include(joinpath("plugins", "codecov.jl"))
include(joinpath("plugins", "travisci.jl"))
include(joinpath("plugins", "gitlabci.jl"))
include(joinpath("plugins", "cirrusci.jl"))
include(joinpath("plugins", "githubpages.jl"))
include(joinpath("plugins", "gitlabpages.jl"))
include(joinpath("plugins", "citation.jl"))

const DEFAULTS_DIR = normpath(joinpath(@__DIR__, "..", "defaults"))
const BADGE_ORDER = [GitHubPages, GitLabPages, TravisCI, AppVeyor, GitLabCI, Codecov, Coveralls]

end


text:
const BUILD_DIR = dirname(@__FILE__)
const DEPS_FILE = joinpath(BUILD_DIR, "deps.jl")
const LIBCSHELL_DIR = joinpath(BUILD_DIR, "libcshell")

function clean_libcshell()
    run(`cargo clean`)
end

function build_libcshell()
    try
        output = readstring(`cargo build --verbose --release`)
        libc_rlib = match(r"libc=(.*liblibc.*\.rlib)", output).captures[1]
    catch e
        println(STDERR, "Failed to build the rust libc rlib due to the following error:")
        rethrow(e)
    end
end

function write_depsfile(text)
    open(DEPS_FILE, "w") do fp
        println(fp, text)
    end
end

function generate_deps()
    cd(LIBCSHELL_DIR)
    clean_libcshell()
    libc_rlib = build_libcshell()
    write_depsfile("const LIBC_RLIB = \"$libc_rlib\"")
end

generate_deps()


text:
@testset "236.lowest-common-ancestor-of-a-binary-tree.jl" begin
    root = TreeNode{Int}([3, 5, 1, 6, 2, 0, 8, nothing, nothing, 7, 4])
    @test lowest_common_ancestor_236(root, root.left, root.left.right.right) == root.left
    @test lowest_common_ancestor_236(root, root.left, root.right) == root
    @test lowest_common_ancestor_236(
        TreeNode{Int}([1, 2]), TreeNode{Int}([1, 2]), TreeNode(2)
    ) == TreeNode{Int}([1, 2])
end


text:
# simple arithmetic expression generator

# note: no recursion so as to avoid issue of infinite strings
@generator EXSimpleExprGen begin
    start() = expression()
    expression() = operand() * " " * operator() * " " * operand()
    operand() = (choose(Bool) ? "-" : "") * join(plus(digit))
    digit() = string(choose(Int,0,9))
    operator() = "+"
    operator() = "-"
    operator() = "/"
    operator() = "*"
end

@testset "simple arithmetic expression generator" begin

	gn = EXSimpleExprGen()

	@mtestset "emits a valid simple expression as a string" reps=Main.REPS alpha=Main.ALPHA begin
	    td = choose(gn)
	    @test typeof(td) <: AbstractString
	    @test ismatch(r"^-?[0-9]+ [+\-/*] -?[0-9]+$", td)
	    @mtest_values_vary td
	end

end


text:
int_rules_7_2_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.2.1*(a+b*arccosh(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ~x * (~a + ~b * acosh(~c * ~x)) ^ ~n - ~b * ~c * ~n * Antiderivative((~x * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (~c / (~b * (~n + 1))) * Antiderivative((~x * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => (1 / (~b * ~c)) * Subst(Antiderivative((~x) ^ ~n * sinh(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.2 Inverse hyperbolic cosine/7.2.1 (a+b arccosh(c x))^n.jl


text:
module strip2

using Random
using FFTW
using ConjGrad
using Plots
using LinearAlgebra
using ConjGrad

# gensky methods
export make_CMB_T_map
export make_CMB_pol_maps
export SZ_source
export PS_exp
export PS_poisson
export Plot_Sky

# Beam methods
export convolve_beam

# Noise sources
export white_noise
export atmospheric_noise
export one_over_f
export one_over_f2
export filter_noise
export observe_sky
export get_tod
export tod_noise
export binned_map
export destriper
export conjgrad
export denoise

# Analysis tools
export windowing!
export get_power_spectrum

abstract type TOD end

struct TOD_fake_pointing <: TOD
    time_order_data::AbstractArray{<:Number}
    pointing::AbstractArray{<:Number}
    noise::AbstractArray{<:Number}
end

include("gensky.jl")
include("beam.jl")
include("noise_map.jl")
include("filter.jl")
include("genpoints.jl")
include("mapmaker.jl")
include("analysis.jl")

end  # module strip2


text:
module JungleHelperEscapeSinkingPlatform

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Platform (Sinking, Escape) (JungleHelper)" => Ahorn.EntityPlacement(
        Maple.SinkingPlatform,
        "rectangle",
        Dict{String, Any}(
          "texture" => "JungleHelper/escape"
        )

    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
# Equivalent to PyCall.PyArray but
# does not require the array to be writable, so that we can read back jax data
function PyCall.PyArray(o::PyObject, readonly::Bool)
    if !readonly
        return PyCall.PyArray(o)
    else
        return PyArray_ReadOnly(o)
    end
end

function PyArray_ReadOnly(o::PyObject)
    info = PyArray_Info_ReadOnly(o)
    return PyCall.PyArray{eltype(info),length(info.sz)}(o, info)
end

function PyArray_Info_ReadOnly(o::PyObject)
    # n.b. the pydecref(::PyBuffer) finalizer handles releasing the PyBuffer
    pybuf = PyCall.PyBuffer(o, PyCall.PyBUF_ND_STRIDED & ~PyCall.PyBUF_WRITABLE)
    T, native_byteorder = PyCall.array_format(pybuf)
    sz = size(pybuf)
    strd = strides(pybuf)
    length(strd) == 0 && (sz = ())
    N = length(sz)
    isreadonly = pybuf.buf.readonly == 1
    return PyCall.PyArray_Info{T,N}(
        native_byteorder,
        sz,
        strd,
        pybuf.buf.buf,
        isreadonly,
        pybuf,
    )
end

PyCall.getindex(o::PyObject, i::PyObject) = PyCall.get(o, i)
PyCall.getindex(o::PyObject, ::Colon) = PyCall.get(o, :)

PyCall.PyObject(::Colon) = _pycolon

##
#=
const _pynamedtuples = Dict{Any, Any}()

struct PyStructSequenceField2
    name::Cstring
    doc::Cstring
end

struct PyStructSequenceFieldJl
    name::String
    doc::String
end
Base.cconvert(::Type{PyStructSequenceField2}, v::PyStructSequenceFieldJl) =
    PyStructSequenceField2(Base.cconvert(Cstring, v.name), Base.cconvert(Cstring, v.doc))

Base.unsafe_convert(::Type{PyStructSequenceField2}, v::Tuple) = PyStructSequenceField2(Base.unsafe_convert(Cstring, v[1]), Base.unsafe_convert(Cstring, v[2]))

struct PyStructSequenceDesc
    name::Cstring
    doc::Cstring
    fields::Vector{Ptr{PyStructSequenceField}}
    n_in_sequence::Cint
end

function PyStructSequenceDesc(T::Type{<:NamedTuple})
    name = "namedtuple_$(length(_pynamedtuples))"
    nm = Base.unsafe_convert(Cstring, name)
    names = [string(nm) for nm=T.names]
    Cnames = [Base.unsafe_convert(Cstring, nm) for nm=names]
    fields = [PyStructSequenceField(nm, C_NULL) for nm=Cnames]
    flds = [convert(Ptr{PyStructSequenceField},Base.pointer_from_objref(f)) for f in fields]
    push!(flds, C_NULL)
    return PyStructSequenceDesc(nm,
                                C_NULL,
                                flds,
                                length(fields))
end

function PyCall.prepycheck(o::PyObject)

    return false
end

function PyCall.PyObject(t::NamedTuple)
    T = typeof(t)

    if T ∈ keys(_pynamedtuples)
        pyT = _pynamedtuples[T]
    else
        pyT = PyStructSequenceDesc(T)

        oT = PyObject(PyCall.@pycheckn ccall((PyCall.@pysym :PyStructSequence_NewType), PyPtr, (Ptr{PyStructSequenceDesc},), Ref(pyT)))
        _pynamedtuples[T] = oT
    end

    #=len = lastindex(t) # lastindex, not length, because of julia#14924
    o = PyObject(PyCall.@pycheckn ccall((PyCall.@pysym :PyTuple_New), PyPtr, (Int,), len))
    for i = 1:len
        oi = PyObject(t[i])
        PyCall.@pycheckz ccall((PyCall.@pysym :PyTuple_SetItem), Cint, (PyPtr,Int,PyPtr),
                         o, i-1, oi)
        PyCall.pyincref(oi) # PyTuple_SetItem steals the reference
    end
    return o=#
end
=#


text:
using Fed.Config: BaseConfig, VanillaConfig, QuantizedConfig,
    QDiffConfig, QDiffStaticConfig, GDConfig, GDStaticConfig
using SHA
using CRC32c

"""
    hash_crc32(data)

Hash data with CRC32 and returns the hash as a `Vector{UInt8}`.
"""
function hash_crc32(data)
    hash = Vector{UInt8}(undef, 4)
    data = crc32c(data)
    hash[1] = (data & 0x000000ff)
    hash[2] = (data & 0x0000ff00) >> 8
    hash[3] = (data & 0x00ff0000) >> 16
    hash[4] = (data & 0xff000000) >> 24
    return hash
end

# config
const SERVERURL = "http://127.0.0.1:8080"
const NUM_COMM_ROUNDS = 100
const FRACTION_CLIENTS = 0.1f0
const NUM_TOTAL_CLIENTS = 100

const base_config = BaseConfig(
    SERVERURL,
    NUM_COMM_ROUNDS,
    FRACTION_CLIENTS,
    NUM_TOTAL_CLIENTS
)

newVanillaConfig() = VanillaConfig(base_config)
newQuantizedConfig(::Type{T}) where T <: Unsigned = QuantizedConfig{T}(base_config)
newQDiffConfig(::Type{T}, is_client::Bool) where T <: Unsigned =
    QDiffConfig{T}(base_config, 256, is_client)

newQDiffStaticConfig(::Type{T}, is_client::Bool) where T <: Unsigned =
    QDiffStaticConfig{T}(base_config, 256, is_client)

newGDConfig(::Type{T}, port::Int, is_client::Bool) where T <: Unsigned =
    GDConfig{T}(
        base_config,
        256,
        sha1,
        # hash_crc32,
        round(T, 0.125 * sizeof(T) * 8), # % of each weight which goes in the basis
        "127.0.0.1",
        port,
        is_client
    )

newGDStaticConfig(::Type{T}, port::Int, is_client::Bool) where T <: Unsigned =
    GDStaticConfig{T}(
        base_config,
        256,
        # sha1,
        hash_crc32,
        round(T, 0.6 * sizeof(T) * 8), # % of each weight which goes in the basis
        "127.0.0.1",
        port,
        is_client
    )

"""
    newconfig()

Returns a new current configuration.
"""
# newconfig(_, _) = newVanillaConfig()
# newconfig(_, _) = newQuantizedConfig(UInt8)
# newconfig(_, is_client) = newQDiffConfig(UInt8, is_client)
# newconfig(_, is_client) = newQDiffStaticConfig(UInt16, is_client)
newconfig(port, is_client) = newGDConfig(UInt8, port, is_client)
# newconfig(port, is_client) = newGDStaticConfig(UInt8, port, is_client)

"""
    newmodel()

Returns a new current model.
"""
newmodel() = initialize_mlp(28*28, 200, 10)


text:
"""
    LagSelection(data, maxlag; constant = true, trend = false) -> LagSelection

Calculate AIC, SIC, and HQC lag selection criteria for an unrestricted VAR.

# Arguments
- `data` : `DataFrame` or `TimeArray` of observations on endogenous variables
- `maxlag::Int` : the maximum number of lags to estimate
- `constant::Bool = true` : include an intercept term
- `trend::Bool = false` : include a linear trend

# Fields
- `maxlag::Int`
- `table::DataFrame`
- `selection::Dict`
"""
struct LagSelection
    maxlag::Int
    table::DataFrame
    selection::Dict
end

function LagSelection(data::DataFrame, maxlag::Int; constant::Bool = true, trend::Bool = false)
    if maxlag < 1
        error("maxlag must be ≥ 1")
    end

    table, selection = lagselect(data, maxlag, constant, trend)
    LagSelection(maxlag, table, selection)
end

function LagSelection(data::TimeArray, maxlag::Int; constant::Bool = true, trend::Bool = false)
    data = DataFrame(data)[:, 2:end]
    LagSelection(data, maxlag; constant = constant, trend = trend)
end

function lagselect(data, maxlag, constant, trend)
    critname = ["AIC", "HQIC", "SIC"]
    AIC = zeros(maxlag)
    HQC = zeros(maxlag)
    SIC = zeros(maxlag)
    for lag = 1:maxlag
        dmat = data[(maxlag-lag+1):end, :]
        v = VAR(dmat, lag; constant = constant, trend = trend)
        AIC[lag] = aic(v)
        HQC[lag] = hqc(v)
        SIC[lag] = sic(v)
    end
    table = DataFrame(lag = collect(1:maxlag), AIC = AIC, HQC = HQC, SIC = SIC)
    selection = Dict(cn => argmin(table[!, cn]) for cn in names(table)[2:end])
    return table, selection
end

function show(io::IO, ls::LagSelection)
    println(io, typeof(ls))
    println(io, "Maximum lags: ", ls.maxlag)
    for k in keys(ls.selection)
        print(io, k, ": ", rpad(ls.selection[k], 3))
    end
    print(io, '\n')
end


text:
type SumKernel <: Kernel
    kerns::Vector{Kernel}
    function SumKernel(args...)
        kerns = Array(Kernel, length(args))
        for (i,k) in enumerate(args)
            isa(k, Kernel) || throw(ArgumentError("All arguments of SumKernel must be Kernel objects"))
            kerns[i] = k
        end
        return new(kerns)
    end
end

function show(io::IO, sumkern::SumKernel, depth::Int = 0)
    pad = repeat(" ", 2 * depth)
    println(io, "$(pad)Type: $(typeof(sumkern))")
    for k in sumkern.kerns
        show(io, k, depth+1)
    end
end

function cov(sumkern::SumKernel, x::Vector{Float64}, y::Vector{Float64})
    s = 0.0
    for k in sumkern.kerns
        s += cov(k, x, y)
    end
    return s
end

# This slows down cov...

## function cov(X::Matrix{Float64}, sumkern::SumKernel)
##     d, nobsv = size(X)
##     s = zeros(nobsv, nobsv)
##     for k in sumkern.kerns
##         BLAS.axpy!(nobsv*nobsv, 1.0, cov(X,k), 1, s, 1)
##         #s += cov(X, k)
##     end
##     return s
## end

## function add_matrices!(X::AbstractMatrix, Y::AbstractMatrix)
##     m,n = size(X)
##     for i in 1:m, j in 1:n
##         X[i,j] += Y[i,j]
##     end
## end

function cov(sumkern::SumKernel, X::Matrix{Float64})
    d, nobsv = size(X)
    s = zeros(nobsv, nobsv)
    for k in sumkern.kerns
        #add_matrices!(s, cov(X,k))
        s[:,:] = s + cov(k, X)
    end
    return s
end

function get_params(sumkern::SumKernel)
    p = Array(Float64, 0)
    for k in sumkern.kerns
        append!(p, get_params(k))
    end
    p
end

get_param_names(sumkern::SumKernel) = composite_param_names(sumkern.kerns, :sk)

function num_params(sumkern::SumKernel)
    n = 0
    for k in sumkern.kerns
        n += num_params(k)
    end
    n
end

function set_params!(sumkern::SumKernel, hyp::Vector{Float64})
    i, n = 1, num_params(sumkern)
    length(hyp) == num_params(sumkern) || throw(ArgumentError("SumKernel object requires $(n) hyperparameters"))
    for k in sumkern.kerns
        np = num_params(k)
        set_params!(k, hyp[i:(i+np-1)])
        i += np
   p end
end

function grad_kern(sumkern::SumKernel, x::Vector{Float64}, y::Vector{Float64})
     dk = Array(Float64, 0)
      for k in sumkern.kerns
        append!(dk,grad_kern(k, x, y))
      end
    dk
end

function grad_stack!(stack::AbstractArray, X::Matrix{Float64}, sumkern::SumKernel)
    s = 1
    for kern in sumkern.kerns
        np = num_params(kern)
        grad_stack!(view(stack,:, :, s:(s+np-1)), X, kern)
        s += np
    end
    return stack
end

# Addition operators
function +(k1::SumKernel, k2::Kernel)
    kerns = [k1.kerns; k2]
    SumKernel(kerns...)
end
function +(k1::SumKernel, k2::SumKernel)
    kerns = [k1.kerns; k2.kerns]
    SumKernel(kerns...)
end
+(k1::Kernel, k2::Kernel) = SumKernel(k1,k2)
+(k1::Kernel, k2::SumKernel) = +(k2,k1)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module Repeat

export basismatrix

using Base: offset_if_vec
using StaticArrays:SVector,MVector,SMatrix,MMatrix
using DataFrames:DataFrame
using Colors
import LazySets
using LinearAlgebra:norm
import ..OpticSim #only LensletAssembly uses OpticSim. This doesn't seem like a great idea. Probably should move LensletAssembly somewhere else or at least remove the dependency.
import OpticSim: surfaceintersection
using OpticSim: virtualpoint,SphericalPolygon,processintersection,point
import Unitful

include("Lattice.jl")
include("HexagonalLattice.jl")
include("RectangularLattice.jl")
include("Array.jl")
include("Cluster.jl")
include("Multilens/Multilens.jl")
include("LensletAssembly.jl")

end #module
export Repeat


text:
# ---
# title: 10. Regular Expression Matching
# id: problem10
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming, Backtracking
# link: <https://leetcode.com/problems/regular-expression-matching/description/>
# hidden: true
# ---
#
# Given an input string (`s`) and a pattern (`p`), implement regular expression
# matching with support for `'.'` and `'*'` where:` `
#
#   * `'.'` Matches any single character.
#   * `'*'` Matches zero or more of the preceding element.
#
# The matching should cover the **entire** input string (not partial).
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aa", p = "a"
#     Output: false
#     Explanation: "a" does not match the entire string "aa".
#
#
# **Example 2:**
#
#
#
#     Input: s = "aa", p = "a*"
#     Output: true
#     Explanation:  '*' means zero or more of the preceding element, 'a'. Therefore, by repeating 'a' once, it becomes "aa".
#
#
# **Example 3:**
#
#
#
#     Input: s = "ab", p = ".*"
#     Output: true
#     Explanation:  ".*" means "zero or more (*) of any character (.)".
#
#
# **Example 4:**
#
#
#
#     Input: s = "aab", p = "c*a*b"
#     Output: true
#     Explanation:  c can be repeated 0 times, a can be repeated 1 time. Therefore, it matches "aab".
#
#
# **Example 5:**
#
#
#
#     Input: s = "mississippi", p = "mis*is*p*."
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 20`
#   * `0 <= p.length <= 30`
#   * `s` contains only lowercase English letters.
#   * `p` contains only lowercase English letters, `'.'`, and `'*'`.
#   * It is guaranteed for each appearance of the character `'*'`, there will be a previous valid character to match.
#
#
## @lc code=start
using LeetCode

function is_match(s::AbstractString, p::AbstractString)::Bool
    if isempty(p)
        isempty(s)
    elseif length(p) >= 2 && p[2] == '*'
        if !isempty(s) && (s[1] == p[1] || p[1] == '.')
            ## case 1: * means 1 or more of previous char, then the first char of `s` and `p` must match.
            is_match(SubString(s, 2), p)
        else
            ## case 2: * means zero of previous char, ignore it and match the rest
            is_match(s, SubString(p, 3))
        end
    else
        if isempty(s)
            false  ## length mismatch
        elseif (s[1] == p[1] || p[1] == '.')
            is_match(SubString(s, 2), SubString(p, 2))
        else
            false  ## first char mismatch
        end
    end
end
## @lc code=end


text:
module BenchMoments

using BenchmarkTools
using LatticeBoltzmann
using StaticArrays

import LatticeBoltzmann:
    density, velocity!, temperature, pressure, momentum_flux, deviatoric_tensor

suite = BenchmarkGroup()

for q in LatticeBoltzmann.Quadratures
    suite[string(q)] = BenchmarkGroup([string(q), "moments"])

    suite[string(q)]["density"] =
        @benchmarkable(density($q, f), setup = (f = copy($q.weights)))

    suite[string(q)]["velocity"] = @benchmarkable(
        velocity!($q, f, ρ, u),
        setup = (f = copy($q.weights);
        ρ = density($q, f);
        u = zeros(eltype(f), dimension($q)))
    )

    suite[string(q)]["pressure"] = @benchmarkable(
        pressure($q, f, ρ, u),
        setup = (f = copy($q.weights);
        ρ = density($q, f);
        u = zeros(eltype(f), dimension($q));
        velocity!($q, f, ρ, u))
    )

    suite[string(q)]["temperature"] = @benchmarkable(
        pressure($q, f, ρ, u),
        setup = (f = copy($q.weights);
        ρ = density($q, f);
        u = zeros(eltype(f), dimension($q));
        velocity!($q, f, ρ, u))
    )

    suite[string(q)]["momentum_flux"] = @benchmarkable(
        momentum_flux($q, f, ρ, u),
        setup = (f = copy($q.weights);
        ρ = density($q, f);
        u = zeros(eltype(f), dimension($q));
        velocity!($q, f, ρ, u))
    )

    suite[string(q)]["deviatoric_tensor"] = @benchmarkable(
        deviatoric_tensor($q, τ, f, ρ, u),
        setup = (τ = 1.0;
        f = copy($q.weights);
        ρ = density($q, f);
        u = zeros(eltype(f), dimension($q));
        velocity!($q, f, ρ, u))
    )

    suite[string(q)]["density", "static"] = @benchmarkable(
        density(static_q, f),
        setup = (
            f = $q.weights |> MVector{length($q.weights)};
            static_q = D2Q9(
                $q.abscissae |> SMatrix{dimension($q), length($q.weights)},
                tuple($q.weights...) |> SVector{length($q.weights)},
                $q.speed_of_sound_squared,
            )
        )
    )
    suite[string(q)]["velocity", "static"] = @benchmarkable(
        velocity!(static_q, f, ρ, u),
        setup = (
            f = $q.weights |> MVector{length($q.weights)};
            ρ = density($q, f);
            u = @MVector zeros(eltype(f), dimension($q));
            static_q = D2Q9(
                $q.abscissae |> SMatrix{dimension($q), length($q.weights)},
                tuple($q.weights...) |> SVector{length($q.weights)},
                $q.speed_of_sound_squared,
            )
        )
    )
end

end  # module
BenchMoments.suite


text:
pluralize(io, n) = n == 1 ? print(io, n) : print(io, n, "s")

function Base.summary{T,S,N}(sig::AbstractSignal{T,S,N})
    b = IOBuffer()
    N == 1 && print(b, length(sig), "-", sig.dims[1] == symbol("") ? "element" : sig.dims[1])
    N > 1 && print(b, join(size(sig), "x"))

    print(b, " ", typeof(sig).name)
    N > 1 && sig.dims[1] != symbol("") && (print(b, " with "); pluralize(b, size(sig, 1)); print(b, " ", sig.dims[1]))
    N > 2 && sig.dims[2] != symbol("") && (print(b, " across "); pluralize(b, size(sig, 2)); print(b, " ", sig.dims[2]))
    # N > 3 && ... # TODO: how do I show other dimension names? Do I at all?
    takebuf_string(b)
end

function Base.writemime(io::IO, m::MIME"text/plain", sig::Signal)
    print(io, summary(sig))
    length(sig) == 0 && return
    println(io, ":")
    print(io, "  Each ", !isempty(sig.dims) && sig.dims[1] != symbol("") ? sig.dims[1] : "element",
              " has ", length(sig.time), " datapoints from ", sig.time[1], " to ", sig.time[end])
    isa(sig, RegularSignal) && print(io, ", at ", samplingfreq(sig))

    show_signal(io, sig)
end

# If the signal contains nested signals display them nicely with sparklines
# I'm not sure I like this since it displays the shape of the data transposed
function show_signal{T,S<:Signal}(io::IO, sig::Signal{T, S}, limit_output::Bool=true)
    length(sig) > 1 && return # TODO: display multiple signals?

    # Determine screen size
    rows, cols = limit_output ? Base.tty_size() : (typemax(Int), typemax(Int))
    rows = min(rows-5, length(sig.time))

    # Gather the times
    ts = Array(ByteString, rows)
    for i=1:min(length(ts), rows)
        ts[i] = string(sig.time[i], ": ")
    end
    twidth = maximum([strwidth(t) for t in ts]) + 2

    # And display the signal for each time
    s = sig[1]
    for i=1:length(ts)
        println(io)
        print(io, lpad(ts[i], twidth))
        if length(s[i]) > cols - twidth
            spark(io, s[i][1:(cols - twidth - 2)])
            print(io, " …")
        else
            spark(io, s[i])
        end
    end

    # And display a continuation mark if we have more data
    if rows < length(sig.time)
        println(io)
        print(io, "  ⋮", lpad(": ", twidth-3), "⋮")
    end
end

function show_signal(io::IO, sig::Signal, limit_output::Bool=true)
    # TODO: How should generic signals of many channels display?
    println(io)
    Base.showarray(io, sig.data; header=false)
end


text:
module AutoPreallocation
using Cassette
using LinearAlgebra: LinearAlgebra

export avoid_allocations, record_allocations, freeze, preallocate, reinitialize!, @no_prealloc

include("record_types.jl")
include("recording.jl")
include("replaying.jl")
include("inference_fixes.jl")
include("no_prealloc.jl")
include("preallocate.jl")

end # module


text:
# Make sure this is precompiled on only one process first
using MPIArrays

julia_exe = joinpath(Sys.BINDIR, Base.julia_exename())
mpiexec = "mpiexec"

testdir = dirname(@__FILE__)

mpifiles = ["basics.jl", "matmul_mpi.jl"]
juliafiles = ["matmul_darray.jl"]
nprocs = clamp(Sys.CPU_THREADS ÷ 2, 2, 4)
ENV["OPENBLAS_NUM_THREADS"] = 1

for f in mpifiles
    cmd = `$mpiexec -n $nprocs $julia_exe --startup-file=no $(joinpath(testdir, f))`
    println("Executing $cmd")
    run(cmd)
end

for f in juliafiles
    cmd = `$julia_exe -p $(nprocs-1) --startup-file=no $(joinpath(testdir, f))`
    println("Executing $cmd")
    run(cmd)
end