text:
## Defines all HMAC functionality
## Pretty much everything in here learned from http://www.lysator.liu.se/~nisse/nettle/nettle.html#Keyed-hash-functions
## Also from reading the header files from libnettle source tree

import Base: show
export HMACState, update!, digest, digest!, hexdigest!, hexdigest

immutable HMACState
    hash_type::HashType
    outer::Array{UInt8,1}
    inner::Array{UInt8,1}
    state::Array{UInt8,1}
end

# Constructor for HMACState
function HMACState(name::AbstractString, key)
    hash_types = get_hash_types()
    name = uppercase(name)
    if !haskey(hash_types, name)
        throw(ArgumentError("Invalid hash type $name: call Nettle.get_hash_types() to see available list"))
    end

    # Construct HMACState object for this type and initialize using Nettle's init functions
    hash_type = hash_types[name]
    outer = Array(UInt8, hash_type.context_size)
    inner = Array(UInt8, hash_type.context_size)
    state = Array(UInt8, hash_type.context_size)
    ccall((:nettle_hmac_set_key,nettle), Void, (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Csize_t,Ptr{UInt8}),
        outer, inner, state, hash_type.ptr, sizeof(key), key)
    return HMACState(hash_type, outer, inner, state)
end

function update!(state::HMACState, data)
    ccall((:nettle_hmac_update,nettle), Void, (Ptr{Void},Ptr{Void},Csize_t,Ptr{UInt8}), state.state,
        state.hash_type.ptr, sizeof(data), data)
    return state
end

function digest!(state::HMACState)
    digest = Array(UInt8,state.hash_type.digest_size)
    ccall((:nettle_hmac_digest,nettle), Void, (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}, Csize_t,
        Ptr{UInt8}), state.outer, state.inner, state.state, state.hash_type.ptr, sizeof(digest), digest)
    return digest
end

# Take a digest, and convert it to a printable hex representation
hexdigest!(state::HMACState) = bytes2hex(digest!(state))

# The one-shot functions that makes this whole thing so easy
digest(hmac_name::AbstractString, key, data) = digest!(update!(HMACState(hmac_name, key), data))
hexdigest(hmac_name::AbstractString, key, data) = hexdigest!(update!(HMACState(hmac_name, key), data))

show(io::IO, x::HMACState) = write(io, "$(x.hash_type.name) HMAC state")


text:
getinputinds(fun::RD.AbstractFunction, n, m) = getinputinds(RD.functioninputs(fun), n, m)
getinputinds(::RD.StateControl, n, m) = 1:n+m
getinputinds(::RD.StateOnly, n, m) = 1:n
getinputinds(::RD.ControlOnly, n, m) = n .+ (1:m)

"""
    ProjectedNewtonSolver

Projected Newton Solver
Direct method developed by the REx Lab at Stanford University
Achieves machine-level constraint satisfaction by projecting onto the feasible subspace.
    It can also take a full Newton step by solving the KKT system.
This solver is to be used exlusively for solutions that are close to the optimal solution.
    It is intended to be used as a "solution polishing" method for augmented Lagrangian methods.
"""
struct ProjectedNewtonSolver{L<:DiscreteDynamics,O<:AbstractObjective,Nx,Nu,T,F} <: ConstrainedSolver{T}
    # Problem Info
    model::L
    obj::O
    x0::Vector{T}
    conset::PNConstraintSet{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    funsig::F

    # Data
    _data::Vector{T}  # storage for primal and dual variables
    Zdata::VectorView{T,Int}
    Z̄data::VectorView{T,Int}
    dY::Vector{T}                 # KKT solution vector
    Ydata::VectorView{T,Int}
    Atop::SparseMatrixCSC{T,Int}  # top Np rows of KKT matrix
    colptr::Vector{Int}
    rowval::Vector{Int}
    nzval::Vector{T}
    d::Vector{T}               # constraints
    b::Vector{T}               # rhs vector
    active::BitVector

    # Trajectories
    Z::SampledTrajectory{Nx,Nu,T,KnotPoint{Nx,Nu,VectorView{T,Int},T}}
    Z̄::SampledTrajectory{Nx,Nu,T,KnotPoint{Nx,Nu,VectorView{T,Int},T}}

    # Cost function
    hess::Vector{Matrix{T}}
    hessdiag::Vector{Diagonal{T,SubArray{T,1,Matrix{T}, Tuple{Vector{CartesianIndex{2}}}, false}}}
    grad::Vector{VectorView{T,Int}}

    # Dynamics
    ∇f::Matrix{Matrix{T}}          # dynamics Jacobians
    f::Vector{Vector{T}}           # dynamics values
    e::Vector{VectorView{T,Int}}   # dynamics residuals
    Iinit::Diagonal{T,Vector{T}}
    Ireg::Diagonal{T,Vector{T}}    # primal regularization

    # Indices
    ix::Vector{UnitRange{Int}}
    iu::Vector{UnitRange{Int}}
    iz::Vector{UnitRange{Int}}
    blocks::SparseBlocks
    hessblocks::Vector{SparseBlockIndex}
    ∇fblocks::Matrix{SparseBlockIndex}
    qdldl::Cqdldl.QDLDLSolver{T,Int}
end

function ProjectedNewtonSolver(prob::Problem{T}, opts::SolverOptions=SolverOptions{T}(),
                                stats::SolverStats=SolverStats(parent=solvername(ProjectedNewtonSolver));
                                use_static::Val{USE_STATIC}=Val(false),
                                kwargs...) where {T,USE_STATIC}
    nx,nu = RD.dims(prob)
    N = TO.horizonlength(prob)
    Nc = sum(num_constraints(prob.constraints))  # stage constraints
    @assert length(nx) == N
    @assert length(nu) == N
    Np = sum(nx) + sum(nu)
    Nd = sum(nx) + Nc       # number of duals (no constraints)

    set_options!(opts; kwargs...)
    if stats.parent == solvername(ProjectedNewtonSolver)
        reset!(stats, opts.iterations)
    end

    nn = [[0;0] reshape(cumsum(vec([nx'; nu'])), 2, :)]
    ix = map(2:N+1) do k
        nn[2,k-1] + 1 : nn[1,k]
    end
    iu = map(2:N) do k
        nn[1,k] + 1 : nn[2,k]
    end
    iz = map(2:N+1) do k
        nn[2,k-1] + 1 : nn[2,k]
    end

    _data = zeros(2Np + Nd + 2nu[1])  # leave extra room for terminal control
    Zdata = view(_data, 1:Np)
    Z̄data = view(_data, Np .+ (1:Np))
    Ydata = view(_data, 2Np + 2nu[1] .+ (1:Nd))
    Atop = spzeros(T, Np, Np + Nd)
    d = zeros(T, Nd)
    b = zeros(T, Np + Nd)  # allocate maximum size
    active = trues(Np + Nd)
    Ireg = Diagonal(fill(-one(T), Nd))

    Z = SampledTrajectory([KnotPoint{nx[k],nu[k]}(view(_data, iz[k]), prob.Z[k].t, prob.Z[k].dt) for k = 1:N])
    Z[end].dt = Inf
    Z̄ = SampledTrajectory([KnotPoint{nx[k],nu[k]}(view(_data, Np .+ iz[k]), prob.Z[k].t, prob.Z[k].dt) for k = 1:N])
    Z̄[end].dt = Inf
    dY = zeros(T, Np + Nd)

    # Create sparse blocks and initialize cost Hessian sparsity
    blocks = SparseBlocks(Np, Np + Nd)
    isdiag = prob.obj[1] isa TO.QuadraticCostFunction ? TO.is_diag(prob.obj[1]) : false
    for i in iz
        addblock!(blocks, i, i, isdiag)
    end
    id = Np .+ (1:Nd)  # dual indices

    # Create the constraint set
    # NOTE: this initializes the sparsity structure in Atop
    conset = PNConstraintSet(prob.constraints, Z̄, opts, Atop, d, view(active, id), blocks)

    # Initialize the storage for the full sparse KKT matrix
    #   These will get dynamically resized depending on the active set
    nnz_Atop = nnz(Atop)
    nnz_A = nnz_Atop + Nd  # add nnz for the dual regularization
    colptr = zeros(Int, Np + Nd + 1)
    rowval = zeros(Int, nnz_A)
    nzval = zeros(T, nnz_A)

    # Storage for the dynamics expansion
    ∇f = [j == 1 ? zeros(T, nx[k], nx[k]+nu[k]) : zeros(T, nx[k], nx[k]) for k = 1:N-1, j = 1:2]
    f = [zeros(T,nx[k]) for k = 1:N-1]
    e = [view(d, conset.cinds[end][k+1]) for k = 1:N-1]

    # Cache the sparse view blocks for the dynamics Jacobians
    izk(k,j) = j == 1 ? iz1(k) : iz2(k)
    iz1(k) = iz[k]    #(k-1) * (n + m) .+ (1:n+m)
    iz2(k) = ix[k+1]  #(k) * (n + m) .+ (1:n)
    ci(k) = conset.cinds[end][k+1] .+ Np
    ∇fblocks = [
        begin blocks[izk(k,j), ci(k)] end for k = 1:N-1, j = 1:2
    ]
    ∇fblocks = vcat(
        [blocks[ix[1], ci(0)] blocks[ix[1], ci(0)]],
        ∇fblocks
    )

    # Storage for the cost expansion
    hess = [zeros(T, length(i), length(i)) for i in iz]
    grad = [view(b,i) for i in iz]
    hessdiag = map(hess) do H
        inds = [CartesianIndex(i,i) for i = 1:size(H,1)]
        Diagonal(view(H, inds))
    end
    hessblocks = map(iz) do i  # cache for the Hessian block views
        blocks[i, i]
    end

    Iinit = -Diagonal(ones(nx[1]))
    Ireg = Diagonal(ones(nx[1]+nu[1]))  # TODO: allow varying sizes
    qdldl = QDLDLSolver{T}(Np + Nd, nnz_A, 2*nnz_A)

    # Check if in-place method works
    if !USE_STATIC
        try
            RD.discrete_dynamics!(prob.model[1], f[1], Z[1])
            f[1] .= 0
            funsig = RD.InPlace()
        catch
            @warn "In-place dynamics not found. Using StaticReturn. May allocate in ProjectedNewtonSolver."
            funsig = RD.StaticReturn()
        end
    else
        funsig = RD.StaticReturn()
    end

    ProjectedNewtonSolver(prob.model[1], prob.obj, Vector(prob.x0), conset, opts, stats, funsig,
        _data, Zdata, Z̄data, dY, Ydata, Atop, colptr, rowval, nzval, d, b, active,
        Z, Z̄, hess, hessdiag, grad, ∇f, f, e, Iinit, Ireg, ix, iu, iz, blocks, hessblocks,
        ∇fblocks, qdldl,
    )
end

# Getters
RD.dims(pn::ProjectedNewtonSolver{<:Any,<:Any,n,m}) where {n,m} = n,m,length(pn.ix)
TO.get_objective(pn::ProjectedNewtonSolver) = pn.obj
TO.get_model(pn::ProjectedNewtonSolver) = pn.model
TO.get_trajectory(pn::ProjectedNewtonSolver) = pn.Z
iterations(pn::ProjectedNewtonSolver) = pn.stats.iterations_pn
solvername(::Type{<:ProjectedNewtonSolver}) = :ProjectedNewton
num_primals(pn::ProjectedNewtonSolver) = length(pn.Zdata)
num_duals(pn::ProjectedNewtonSolver) = length(pn.Ydata)
TO.get_constraints(pn::ProjectedNewtonSolver) = pn.conset

function primalregularization!(pn::ProjectedNewtonSolver)
    ρ_primal = pn.opts.ρ_primal
    Ireg = pn.Ireg
    Ireg.diag .= ρ_primal
    for block in pn.hessblocks
        pn.Atop[block] .+= Ireg
    end
end

function cost_hessian!(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z)
    obj = pn.obj
    for k in eachindex(Z)
        RD.hessian!(obj.diffmethod[k], obj.cost[k], pn.hess[k], Z[k])

        # Copy to Sparse Array
        hessblock = pn.hessblocks[k]
        if hessblock.block.isdiag
            pn.Atop[hessblock] .= pn.hessdiag[k]
        else
            pn.Atop[hessblock] .= UpperTriangular(pn.hess[k])
        end
    end
end

function cost_gradient!(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z)
    obj = pn.obj
    for k in eachindex(Z)
        RD.gradient!(obj.diffmethod[k], obj.cost[k], pn.grad[k], Z[k])
    end
end

function dynamics_expansion!(pn::ProjectedNewtonSolver{<:Any,<:Any,Nx,Nu},
        Z::SampledTrajectory=pn.Z̄, ∇f=pn.∇f) where {Nx,Nu}
    diff = pn.opts.dynamics_diffmethod
    model = pn.model
    N = length(pn.ix)
    for k = 1:N - 1
        RD.jacobian!(pn.funsig, diff, model, ∇f[k,1], pn.f[k], Z[k])

        # Copy to Sparse Array
        ∇fview1 = pn.Atop[pn.∇fblocks[k+1,1]]'
        ∇fview1 .= ∇f[k,1]

        ∇fblock2 = pn.∇fblocks[k+1,2]
        pn.Atop[∇fblock2] .= pn.Iinit  # assumes explicit
    end
end

function dynamics_error!(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z̄)
    model = pn.model
    N = length(Z)
    for k = 1:N - 1
        RD.evaluate!(pn.funsig, model, pn.f[k], Z[k])
        RD.state_diff!(model, pn.e[k], pn.f[k], RD.state(Z[k+1]))
    end
end

function evaluate_constraints!(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z̄)
    ix = pn.ix
    n = RD.dims(pn)[1]

    # Initial condition
    pn.d[ix[1]] .= pn.x0 .- RD.state(Z[1])

    # Dynamics
    dynamics_error!(pn, Z)

    # Stage constraints
    evaluate_constraints!(pn.conset, Z)
end

function constraint_jacobians!(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z̄)
    ix = pn.ix
    # n = RD.state_dim(pn.model)

    # Initial condition
    # TODO: skip this since it's constant
    x0block = pn.∇fblocks[1]
    pn.Atop[x0block] .= pn.Iinit

    # Dynamics
    dynamics_expansion!(pn, Z)

    # Stage constraints
    constraint_jacobians!(pn.conset, Z)
end

function max_violation(pn::ProjectedNewtonSolver, Z::SampledTrajectory=pn.Z̄)
    evaluate_constraints!(pn, Z)
    update_active_set!(pn)
    max_violation(pn, nothing)
end

function norm_violation(pn::ProjectedNewtonSolver, p=1)
    return norm(pn.d[pn.active], p)
end

function max_violation(pn::ProjectedNewtonSolver, Z::Nothing)
    Np = num_primals(pn)
    c_max = zero(eltype(pn.d))
    for i in eachindex(pn.d)
        if pn.active[i + Np]
            c_max = max(c_max, abs(pn.d[i]))
        end
    end
    return c_max
end

function update_active_set!(pn::ProjectedNewtonSolver)
    update_active_set!(pn.conset)
    return nothing
end

function active_constraints(pn::ProjectedNewtonSolver)
    return sparse(pn.D)[pn.active, :], pn.d[pn.active]
end

ext: jl


text:
module RatFun
using Base, ApproxFun, RecipesBase

import ApproxFun: evaluate, dimension, domain, setdomain, PiecewiseSpace, DiracSpace,
                    PointSpace, ∞, FunTypes, components
import Base: +, -, *, /, getindex, broadcast, //, numerator, denominator

export RationalFun, inv

struct RationalFun{F1<:Fun,F2<:Fun}
    p::F1
    q::F2
end

//(a::FunTypes, b::FunTypes) = RationalFun(a,b)

domain(r::RationalFun) = domain(r.p)
setdomain(r::RationalFun,d) = RationalFun(setdomain(r.p,d),setdomain(r.q,d))

function evaluate(r::RationalFun,x)
    (r.p)(x)/(r.q)(x)
end

(r::RationalFun)(x) = evaluate(r,x)

numerator(r::RationalFun) = r.p
denominator(r::RationalFun) = r.q

*(r1::RationalFun,r2::RationalFun) = RationalFun(r1.p*r2.p,r1.q*r2.q)
*(r::RationalFun,a::Union{Number,Fun}) = RationalFun(r.p*a,r.q)
*(a::Union{Number,Fun},r::RationalFun) = RationalFun(a*r.p,r.q)

Base.inv(r::RationalFun) = RationalFun(r.q,r.p)

components(r::RationalFun) = RationalFun.(components(r.p),components(r.q))

broadcast(::typeof(/), r1::RationalFun, r2::RationalFun) = r1.*inv(r2)
broadcast(::typeof(/), a, r::RationalFun) = a.*inv(r)

(/)(r1::RationalFun,r2::RationalFun) = r1*inv(r2)
(/)(a,r::RationalFun) = a*inv(r)

(/)(r::RationalFun,a) = RationalFun(r.p,r.q*a)
broadcast(::typeof(/), r::RationalFun, a) = (1/a).*r

for op = (:+,:-)
  @eval begin
      $op(r1::RationalFun,r2::RationalFun) = RationalFun($op((r1.p*r2.q),(r2.p*r1.q)),r1.q*r2.q)
      broadcast(::typeof($op),r1::RationalFun,r2::RationalFun) = $op(r1,r2)
  end
end

Base.convert(::Type{Fun},r::RationalFun) = r.p/r.q

# The padding in this function can be improved
# No support for functions with poles within the domain
function plotptsvals(r::RationalFun)
    p = r.p
    q = r.q
    plen = ncoefficients(p)
    qlen = ncoefficients(q)
    if dimension(space(p)) == ∞ && dimension(space(q)) == ∞
        p=pad(p,10plen+10*qlen+1000)
        q=pad(q,10plen+10*qlen+1000)
        r = RationalFun(p,q)
    elseif dimension(space(p)) == ∞ && dimension(space(q)) < ∞
        p=pad(p,10plen+10*dimension(space(q))+1000)
        q=pad(q,10plen+10*dimension(space(q))+1000)
        r = RationalFun(p,q)
    elseif dimension(space(q)) == ∞ && dimension(space(p)) < ∞
        q=pad(q,dimension(space(p))+10qlen+1000)
        p=pad(p,dimension(space(p))+10qlen+1000)
        r = RationalFun(p,q)
    else # both finite dimensional.
        @assert dimension(space(p)) == dimension(space(q))
        p=pad(p,dimension(space(p)))
        q=pad(q,dimension(space(q)))
        r = RationalFun(p,q)
    end
    return points(r.p),values(r.p)./values(r.q)
end

@recipe function f(g::RationalFun{Fun{S,T,V1}}) where {S,T<:Real,V1}
    plotptsvals(g)
end

@recipe function f(x::AbstractVector{T},g::RationalFun{Fun{S,T,V1}}) where {S,T<:Real,V1}
    x,g(x)
end

@recipe function f(r::RationalFun{Fun{S1,T1,V1},
                                  Fun{S2,T2,V2}}) where {S1<:PiecewiseSpace,T1<:Real,S2<:PiecewiseSpace,T2<:Real,V1,V2}
    vp = components(r.p)
    vq = components(r.q)

    if !isempty(vp)
        @series begin
            primary --> true
            RationalFun(vp[1],vq[1])
        end
    end

    for k=2:length(vp)
        @series begin
            primary := false
            RationalFun(vp[k],vq[k])
        end
    end
end

# For dirac space, we draw a dotted line extending to infinity
@recipe function f(r::RationalFun{Fun{S1,T1,V1},
                                  Fun{S2,T2,V2}}) where {S1<:DiracSpace,T1<:Real,S2<:PointSpace,T2<:Real,V1,V2}
    p = r.p
    q = r.q
    pts=space(p).points
    n=length(pts)
    ws=pad(p.coefficients./q.coefficients,length(pts))
    @series begin
        primary --> true
        ones(2)*pts[1],[0,1]*ws[1]
    end

    if length(ws) > 1
        @series begin
            primary := false
            ones(2)*pts[2:end]',[0,1]*ws[2:end]'
        end
    end

    @series begin
        primary := false
        linestyle := :dot
        ones(2)*pts',[1,2]*ws'
    end
end

# for PointSpace, we draw just a line
@recipe function f(r::RationalFun{Fun{S1,T1,V1},
                                  Fun{S2,T2,V2}}) where {S1<:PointSpace,T1<:Real,S2<:PointSpace,T2<:Real,V1,V2}
    p = r.p
    q = r.q
    pts=space(p).points
    n=length(pts)
    ws=pad(p.coefficients./q.coefficients,length(pts))

    @series begin
        primary --> true
        ones(2)*pts[1],[0,1]*ws[1]
    end

    @series begin
        primary := false
        ones(2)*pts[2:end]',[0,1]*ws[2:end]'
    end
end

end #module


text:
struct Name{name}
end

"""
    Name(name)

Create a typed `Name`. Inverse of [`unname`](@ref)

```jlodctest
julia> using LightQuery

julia> Name(:a)
`a`
```
"""
@pure function Name(name)
    Name{name}()
end
export Name

"""
    unname(::Name{name}) where name

Inverse of [`Name`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred unname(Name(:a))
:a
```
"""
@inline function unname(::Name{name}) where {name}
    name
end
export unname

@inline function show(output::IO, ::Name{name}) where {name}
    print(output, '`', name, '`')
end

const Named{name} = Tuple{Name{name}, Any}

@inline function recursive_get(initial, name)
    throw(BoundsError(initial, (name,)))
end
@inline function recursive_get(initial, name, (check_name, value), rest...)
    if name === check_name
        value
    else
        recursive_get(initial, name, rest...)
    end
end

@inline function (::Name{name})(object) where {name}
    getproperty(object, name)
end
@inline function (name::Name)(data::Some{Named})
    recursive_get(data, name, data...)
end

@inline function get_pair(data, name)
    name, name(data)
end

@inline function (some_names::Some{Name})(data)
    partial_map(get_pair, data, some_names)
end
@inline function (::Tuple{})(::Some{Named})
    ()
end

@inline function isless(::Name{name1}, ::Name{name2}) where {name1, name2}
    isless(name1, name2)
end

# to override recusion limit on constant propagation
@pure function to_Names(some_names::Some{Symbol})
    map_unrolled(Name, some_names)
end
@inline function to_Names(them)
    map_unrolled(Name, Tuple(them))
end

@inline function NamedTuple(data::Some{Named})
    NamedTuple{map_unrolled(unname ∘ key, data)}(
        map_unrolled(value, data)
    )
end

function code_with_names(other)
    other
end
function code_with_names(symbol::QuoteNode)
    Name{symbol.value}()
end
function code_with_names(code::Expr)
    if @capture code data_.name_
        Expr(:call, Name{name}(), code_with_names(data))
    elseif @capture code name_Symbol = value_
        Expr(:tuple, Name{name}(), code_with_names(value))
    else
        Expr(code.head, map(code_with_names, code.args)...)
    end
end

"""
    macro name(code)

Switch to [`named_tuple`](@ref)s.

```jldoctest name
julia> using LightQuery

julia> using Test: @inferred

julia> data = @name (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
```

based on typed `Name`s.

```jldoctest name
julia> @name :a
`a`
```

`Name`s can be used as properties

```jldoctest name
julia> @name @inferred data.c
1

julia> @name data.g
ERROR: BoundsError: attempt to access ((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
[...]
```

and selector functions.

```jldoctest name
julia> @name @inferred (:c)(data)
1
```

Multiple names can be used as selector functions

```jldoctest name
julia> @name @inferred (:c, :f)(data)
((`c`, 1), (`f`, 1.0))
```

You can also convert back to `NamedTuple`s.

```jldoctest name
julia> @inferred NamedTuple(data)
(a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
```
"""
macro name(code)
    esc(code_with_names(code))
end
export @name

"""
    named_tuple(data)

Convert `data` to a `named_tuple` (see [`@name`](@ref)).

```jldoctest named_tuple
julia> using LightQuery

julia> using Test: @inferred

julia> @inferred named_tuple((a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0))
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
```

For stability working with arbitrary `struct`s, `propertynames` must constant propagate.

```jldoctest named_tuple
julia> struct MyType
            a::Int
            b::Float64
            c::Int
            d::Float64
            e::Int
            f::Float64
        end

julia> import Base: propertynames

julia> @inline propertynames(::MyType) = (:a, :b, :c, :d, :e, :f);

julia> @inferred named_tuple(MyType(1, 1.0, 1, 1.0, 1, 1.0))
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
```
"""
@inline function named_tuple(data)
    to_Names(propertynames(data))(data)
end
export named_tuple

@inline function haskey(data::Some{Named}, name::Name)
    in_unrolled(name, map_unrolled(key, data)...)
end
@inline function haskey(data, ::Name{name}) where {name}
    hasproperty(data, name)
end

"""
    remove(data, old_names...)

Remove `old_names` from `data`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = @name (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))

julia> @name @inferred remove(data, :c, :f)
((`a`, 1), (`b`, 1.0), (`d`, 1.0), (`e`, 1))
```
"""
@inline function remove(data, old_names...)
    diff_unrolled(old_names, map_unrolled(key, data)...)(data)
end

export remove

"""
    transform(data, assignments...)

Merge `assignments` into `data`, overwriting old values.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = @name (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))

julia> @name @inferred transform(data, c = 2.0, f = 2, g = 1, h = 1.0)
((`a`, 1), (`b`, 1.0), (`d`, 1.0), (`e`, 1), (`c`, 2.0), (`f`, 2), (`g`, 1), (`h`, 1.0))
```
"""
@inline function transform(data, assignments...)
    remove(data, map_unrolled(key, assignments)...)..., assignments...
end
export transform

@inline function merge_2(data1, data2)
    transform(data1, data2...)
end
@inline function merge(datas::Some{Named}...)
    reduce_unrolled(merge_2, datas...)
end

@inline function rename_one(data, (new_name, old_name))
    new_name, old_name(data)
end
"""
    rename(data, new_name_old_names...)

Rename `data`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = @name (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))

julia> @name @inferred rename(data, c2 = :c, f2 = :f)
((`a`, 1), (`b`, 1.0), (`d`, 1.0), (`e`, 1), (`c2`, 1), (`f2`, 1.0))
```
"""
@inline function rename(data, new_name_old_names...)
    remove(data, map_unrolled(value, new_name_old_names)...)...,
    partial_map(rename_one, data, new_name_old_names)...
end
export rename

"""
    gather(data, new_name_old_names...)

For each `new_name, old_names` pair in `new_name_old_names`, gather the `old_names` into a single `new_name`. Inverse of [`spread`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> data = @name (a = 1, b = 1.0, c = 1, d = 1.0, e = 1, f = 1.0)
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))

julia> @name @inferred gather(data, g = (:b, :e), h = (:c, :f))
((`a`, 1), (`d`, 1.0), (`g`, ((`b`, 1.0), (`e`, 1))), (`h`, ((`c`, 1), (`f`, 1.0))))
```
"""
@inline function gather(data, new_name_old_names...)
    remove(
        data,
        flatten_unrolled(map_unrolled(value, new_name_old_names)...)...
    )...,
    partial_map(rename_one, data, new_name_old_names)...
end
export gather

"""
    spread(data, some_names...)

Unnest nested named tuples. Inverse of [`gather`](@ref).

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> gathered = @name (a = 1, d = 1.0, g = (b = 1.0, e = 1), h = (c = 1, f = 1.0))
((`a`, 1), (`d`, 1.0), (`g`, ((`b`, 1.0), (`e`, 1))), (`h`, ((`c`, 1), (`f`, 1.0))))

julia> @name @inferred spread(gathered, :g, :h)
((`a`, 1), (`d`, 1.0), (`b`, 1.0), (`e`, 1), (`c`, 1), (`f`, 1.0))
```
"""
@inline function spread(data, some_names...)
    remove(data, some_names...)...,
    flatten_unrolled(map_unrolled(value, some_names(data))...)...
end
export spread

"""
    struct Apply{Names}

Apply [`Name`](@ref)s to unnamed values.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @name @inferred Apply((:a, :b, :c, :d, :e, :f))((1, 1.0, 1, 1.0, 1, 1.0))
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
```
"""
struct Apply{Names <: Some{Name}}
    names::Names
end
export Apply

@inline function (apply::Apply)(them)
    map_unrolled(tuple, apply.names, them)
end

struct InRow{name, type, position}
end

@pure function InRow(name, type, position)
    InRow{name, type, position}()
end
@inline function (::InRow{name, type, position})(row::Row) where {name, type, position}
    getcell(getfile(row), type, position, getrow(row))::type
end
@inline function get_pair(row::Row, column::InRow{name}) where {name}
    name, column(row)
end
@inline function (columns::Some{InRow})(data)
    partial_map(get_pair, data, columns)
end

@inline function InRow_at(::Schema{Names, Types}, index) where {Names, Types}
    InRow(Name{Names[index]}(), fieldtype(Types, index), index)
end

"""
    row_info(::Tables.Schema)

Get row info for the schema. Can be used as a type stable selector function.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> using CSV: File

julia> using Tables: schema

julia> test = File("test.csv")
CSV.File("test.csv"):
Size: 1 x 6
Tables.Schema:
 :a  Int64
 :b  Float64
 :c  Int64
 :d  Float64
 :e  Int64
 :f  Float64

julia> template = @inferred row_info(schema(test))
(LightQuery.InRow{`a`,Int64,1}(), LightQuery.InRow{`b`,Float64,2}(), LightQuery.InRow{`c`,Int64,3}(), LightQuery.InRow{`d`,Float64,4}(), LightQuery.InRow{`e`,Int64,5}(), LightQuery.InRow{`f`,Float64,6}())

julia> @inferred template(first(test))
((`a`, 1), (`b`, 1.0), (`c`, 1), (`d`, 1.0), (`e`, 1), (`f`, 1.0))
```
"""
@inline function row_info(a_schema::Schema{Names}) where {Names}
    ntuple(let a_schema = a_schema
        @inline InRow_at_capture(index) = InRow_at(a_schema, index)
    end, Val{length(Names)}())
end

export row_info


text:
using Timeout
using Test

@testset "Timeout.jl" begin
    @test timeout(()->sleep(1), 2)() == nothing
    @test_throws TimeoutException timeout(()->sleep(2), 1)()
end


text:
using Jacobi
using Polynomials

# Script for writing the legerdre_RefLine file with Legendre polynomials up to
# a maximum degree of max_degree.
max_degree = 50

setprecision(BigFloat, 512)

function weights(points)
    N = length(points)
    return [ 2/((1-x^2)*(dlegendre(x, N)^2)) for x in points ]
end

open("legendre_line.jl", "w") do io
    for i = 1:max_degree
        # Legendre = Jacobi with α = β = 0
        x = jacobi_zeros(i, 0, 0, BigFloat)
        w = weights(x)/2
        println("Error in sum of the weights for degree $i: $(1-sum(w))")
        x = (x .+ 1)/2
        println(io, "@generated function gauss_quadrature(form::Val{:legendre},")
        println(io, "                                     shape::RefLine,")
        println(io, "                                     degree::Val{$i},")
        println(io, "                                     type::Type{T}) where {T}")
        println(io, """    weights = SVector(:(\$(T(big"$(w[1])"))),""")
        for j in 2:i
        println(io, """                      :(\$(T(big"$(w[j])"))),""")
        end
        println(io, "                       )")
        println(io, """    points = SVector(:(\$(NTuple{1,T}(big"$(x[1])"))),""")
        for j in 2:i
        println(io, """                     :(\$(NTuple{1,T}(big"$(x[j])"))),""")
        end
        println(io, "                      )")
        println(io, "    return weights, points")
        println(io, "end")
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#========================================================================================#
"""
	Casinos

Module Casinos: A repository for fast access to arrays random numbers.

Author: Niall Palfreyman, 9/02/22
"""
module Casinos

test = 5

end		# ... of module Casinos


text:
"""
The main WaylandScanner module provides a way to scan any Wayland protocol files and generate Julia code for using the protocols.
This can either be done by generating expressions into the module scope at compile time or outputting a Julia script file with the definitions.
"""
module WaylandScanner

export tostructid, wlparse, genprotocol, generate_script, output_script

using WaylandCore

import Unicode, Base.read, Base.write
using Base: finalizer

# Parsing
include("Parsing.jl")
using .Parsing

# Utils
tostructid(identifier) = Symbol(replace(Unicode.titlecase(string(identifier)), "_" => "")) # wl_display to WlDisplay

# Generation
include("Generation.jl")
using .Generation

end  # module WaylandScanner


text:
module SimulacoesTemplate

  f(x)=2*x + 1
  export f
end


text:
function read_ace_header(filename)
    open(filename) do f
        idx = 1
        line = readline(f)
        if length(split(line)) == 3
            readline(f)
            line = readline(f)
        end
        HZ, AW, TZ, HD = split(line, " ", keepempty = false)
        suffix = HZ[end-2 : end]
        if suffix == "00c"
            T = 293.6
        elseif suffix == "01c"
            T = 600.0
        elseif suffix == "02c"
            T = 900.0
        elseif suffix == "03c"
            T = 1200.0
        elseif suffix == "04c"
            T = 2500.0
        elseif suffix == "05c"
            T = 0.1
        elseif suffix == "06c"
            T = 250.0
        end
        return HZ, T, parse(Float64, AW)
    end
end

xs_dir = "C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\Lib80x\\Lib80x"

xs_paths = []
for (root, dirs, files) in walkdir(xs_dir)
    for filename in joinpath.(root, files)
        HZ, T, AW = read_ace_header(filename)
        push!(xs_paths, [HZ, AW, T, filename])
        #println(HZ)
    end
end

xs_dir_filename = "C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\Lib80x\\xs_dir.txt"
using DelimitedFiles
writedlm(xs_dir_filename, xs_paths, ',')

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
module TestWithprogress
include("preamble.jl")

@testset "traits" begin
    @test Base.IteratorSize(withprogress(1:10)) == Base.HasShape{1}()
    @test length(withprogress(1:10)) == 10
    @test Base.IteratorEltype(withprogress(1:10)) == Base.HasEltype()
    @test eltype(withprogress(1:10)) == Int
end

end  # module


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_PrimitiveType <: ProtoEnum
    PRIMITIVE_TYPE_INVALID::Int32
    PRED::Int32
    S8::Int32
    S16::Int32
    S32::Int32
    S64::Int32
    U8::Int32
    U16::Int32
    U32::Int32
    U64::Int32
    F16::Int32
    F32::Int32
    BF16::Int32
    F64::Int32
    C64::Int32
    TUPLE::Int32
    OPAQUE::Int32
    TOKEN::Int32
    __enum_PrimitiveType() = new(0,1,2,3,4,5,6,7,8,9,10,11,16,12,15,13,14,17)
end #struct __enum_PrimitiveType
const PrimitiveType = __enum_PrimitiveType()

struct __enum_Format <: ProtoEnum
    INVALID_FORMAT::Int32
    DENSE::Int32
    SPARSE::Int32
    __enum_Format() = new(0,1,2)
end #struct __enum_Format
const Format = __enum_Format()

struct __enum_FftType <: ProtoEnum
    FFT::Int32
    IFFT::Int32
    RFFT::Int32
    IRFFT::Int32
    __enum_FftType() = new(0,1,2,3)
end #struct __enum_FftType
const FftType = __enum_FftType()

struct __enum_RandomDistribution <: ProtoEnum
    RNG_INVALID::Int32
    RNG_UNIFORM::Int32
    RNG_NORMAL::Int32
    __enum_RandomDistribution() = new(0,1,2)
end #struct __enum_RandomDistribution
const RandomDistribution = __enum_RandomDistribution()

mutable struct PaddingConfig_PaddingConfigDimension <: ProtoType
    edge_padding_low::Int64
    edge_padding_high::Int64
    interior_padding::Int64
    PaddingConfig_PaddingConfigDimension(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct PaddingConfig_PaddingConfigDimension

mutable struct PaddingConfig <: ProtoType
    dimensions::Base.Vector{PaddingConfig_PaddingConfigDimension}
    PaddingConfig(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct PaddingConfig

mutable struct Tile <: ProtoType
    dimensions::Base.Vector{Int64}
    Tile(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Tile
const __pack_Tile = Symbol[:dimensions]
meta(t::Type{Tile}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_Tile, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct Layout <: ProtoType
    format::Int32
    minor_to_major::Base.Vector{Int64}
    max_sparse_elements::Int64
    tiles::Base.Vector{Tile}
    element_size_in_bits::Int64
    Layout(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Layout
const __fnum_Layout = Int[4,1,5,6,7]
const __pack_Layout = Symbol[:minor_to_major]
meta(t::Type{Layout}) = meta(t, ProtoBuf.DEF_REQ, __fnum_Layout, ProtoBuf.DEF_VAL, true, __pack_Layout, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ShapeProto <: ProtoType
    element_type::Int32
    dimensions::Base.Vector{Int64}
    tuple_shapes::Base.Vector{ShapeProto}
    layout::Layout
    ShapeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ShapeProto
const __fnum_ShapeProto = Int[2,3,4,5]
const __pack_ShapeProto = Symbol[:dimensions]
meta(t::Type{ShapeProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ShapeProto, ProtoBuf.DEF_VAL, true, __pack_ShapeProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ProgramShapeProto <: ProtoType
    parameters::Base.Vector{ShapeProto}
    result::ShapeProto
    parameter_names::Base.Vector{AbstractString}
    ProgramShapeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ProgramShapeProto

mutable struct ComputationStats <: ProtoType
    flop_count::Float64
    transcendental_count::Float64
    ComputationStats(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ComputationStats

mutable struct OpMetadata <: ProtoType
    op_type::AbstractString
    op_name::AbstractString
    source_file::AbstractString
    source_line::Int32
    OpMetadata(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpMetadata

mutable struct ExecutionProfile <: ProtoType
    compilation_cache_hit::Bool
    compile_time_ms::Int64
    compute_cycle_count::Int64
    compute_time_ns::Int64
    compute_and_transfer_time_ns::Int64
    executable_size_in_bytes::Int64
    ExecutionProfile(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecutionProfile

mutable struct ExecutionHandle <: ProtoType
    handle::Int64
    ExecutionHandle(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ExecutionHandle

mutable struct GlobalDataHandle <: ProtoType
    handle::Int64
    GlobalDataHandle(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GlobalDataHandle

mutable struct DeviceHandle <: ProtoType
    handle::Int64
    device_count::Int64
    DeviceHandle(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DeviceHandle

struct __enum_ChannelHandle_ChannelType <: ProtoEnum
    CHANNEL_TYPE_INVALID::Int32
    DEVICE_TO_DEVICE::Int32
    DEVICE_TO_HOST::Int32
    HOST_TO_DEVICE::Int32
    __enum_ChannelHandle_ChannelType() = new(0,1,2,3)
end #struct __enum_ChannelHandle_ChannelType
const ChannelHandle_ChannelType = __enum_ChannelHandle_ChannelType()

mutable struct ChannelHandle <: ProtoType
    handle::Int64
    _type::Int32
    ChannelHandle(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ChannelHandle

mutable struct DeviceAssignmentProto_ComputationDevice <: ProtoType
    replica_device_ids::Base.Vector{Int32}
    DeviceAssignmentProto_ComputationDevice(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DeviceAssignmentProto_ComputationDevice
const __pack_DeviceAssignmentProto_ComputationDevice = Symbol[:replica_device_ids]
meta(t::Type{DeviceAssignmentProto_ComputationDevice}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_DeviceAssignmentProto_ComputationDevice, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct DeviceAssignmentProto <: ProtoType
    replica_count::Int32
    computation_count::Int32
    computation_devices::Base.Vector{DeviceAssignmentProto_ComputationDevice}
    DeviceAssignmentProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DeviceAssignmentProto

mutable struct LiteralProto <: ProtoType
    shape::ShapeProto
    preds::Base.Vector{Bool}
    s8s::Array{UInt8,1}
    u8s::Array{UInt8,1}
    s32s::Base.Vector{Int32}
    s64s::Base.Vector{Int64}
    u32s::Base.Vector{UInt32}
    u64s::Base.Vector{UInt64}
    f32s::Base.Vector{Float32}
    f64s::Base.Vector{Float64}
    c64s::Base.Vector{Float32}
    tuple_literals::Base.Vector{LiteralProto}
    f16s::Array{UInt8,1}
    bf16s::Array{UInt8,1}
    u16s::Array{UInt8,1}
    s16s::Array{UInt8,1}
    sparse_indices::Base.Vector{Int64}
    LiteralProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct LiteralProto
const __fnum_LiteralProto = Int[1,2,15,3,4,5,6,7,8,9,12,10,11,13,16,17,14]
const __pack_LiteralProto = Symbol[:preds,:s32s,:s64s,:u32s,:u64s,:f32s,:f64s,:c64s,:sparse_indices]
meta(t::Type{LiteralProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_LiteralProto, ProtoBuf.DEF_VAL, true, __pack_LiteralProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct WindowDimension <: ProtoType
    size::Int64
    stride::Int64
    padding_low::Int64
    padding_high::Int64
    window_dilation::Int64
    base_dilation::Int64
    window_reversal::Bool
    WindowDimension(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct WindowDimension

mutable struct Window <: ProtoType
    dimensions::Base.Vector{WindowDimension}
    Window(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Window

mutable struct GatherDimensionNumbers <: ProtoType
    offset_dims::Base.Vector{Int64}
    collapsed_slice_dims::Base.Vector{Int64}
    start_index_map::Base.Vector{Int64}
    index_vector_dim::Int64
    GatherDimensionNumbers(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GatherDimensionNumbers
const __pack_GatherDimensionNumbers = Symbol[:offset_dims,:collapsed_slice_dims,:start_index_map]
meta(t::Type{GatherDimensionNumbers}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_GatherDimensionNumbers, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ScatterDimensionNumbers <: ProtoType
    update_window_dims::Base.Vector{Int64}
    inserted_window_dims::Base.Vector{Int64}
    scatter_dims_to_operand_dims::Base.Vector{Int64}
    index_vector_dim::Int64
    ScatterDimensionNumbers(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ScatterDimensionNumbers
const __pack_ScatterDimensionNumbers = Symbol[:update_window_dims,:inserted_window_dims,:scatter_dims_to_operand_dims]
meta(t::Type{ScatterDimensionNumbers}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_ScatterDimensionNumbers, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ConvolutionDimensionNumbers <: ProtoType
    input_batch_dimension::Int64
    input_feature_dimension::Int64
    input_spatial_dimensions::Base.Vector{Int64}
    kernel_input_feature_dimension::Int64
    kernel_output_feature_dimension::Int64
    kernel_spatial_dimensions::Base.Vector{Int64}
    output_batch_dimension::Int64
    output_feature_dimension::Int64
    output_spatial_dimensions::Base.Vector{Int64}
    ConvolutionDimensionNumbers(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ConvolutionDimensionNumbers
const __fnum_ConvolutionDimensionNumbers = Int[7,8,11,3,4,6,9,10,12]
const __pack_ConvolutionDimensionNumbers = Symbol[:input_spatial_dimensions,:kernel_spatial_dimensions,:output_spatial_dimensions]
meta(t::Type{ConvolutionDimensionNumbers}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ConvolutionDimensionNumbers, ProtoBuf.DEF_VAL, true, __pack_ConvolutionDimensionNumbers, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct DotDimensionNumbers <: ProtoType
    lhs_contracting_dimensions::Base.Vector{Int64}
    rhs_contracting_dimensions::Base.Vector{Int64}
    lhs_batch_dimensions::Base.Vector{Int64}
    rhs_batch_dimensions::Base.Vector{Int64}
    DotDimensionNumbers(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DotDimensionNumbers
const __pack_DotDimensionNumbers = Symbol[:lhs_contracting_dimensions,:rhs_contracting_dimensions,:lhs_batch_dimensions,:rhs_batch_dimensions]
meta(t::Type{DotDimensionNumbers}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_DotDimensionNumbers, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

struct __enum_OpSharding_Type <: ProtoEnum
    REPLICATED::Int32
    MAXIMAL::Int32
    TUPLE::Int32
    OTHER::Int32
    __enum_OpSharding_Type() = new(0,1,2,3)
end #struct __enum_OpSharding_Type
const OpSharding_Type = __enum_OpSharding_Type()

mutable struct OpSharding <: ProtoType
    _type::Int32
    tile_shape::ShapeProto
    tile_assignment_dimensions::Base.Vector{Int64}
    tile_assignment_devices::Base.Vector{Int64}
    tuple_shardings::Base.Vector{OpSharding}
    OpSharding(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OpSharding
const __pack_OpSharding = Symbol[:tile_assignment_dimensions,:tile_assignment_devices]
meta(t::Type{OpSharding}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_OpSharding, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct ReplicaGroup <: ProtoType
    replica_ids::Base.Vector{Int64}
    ReplicaGroup(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ReplicaGroup
const __pack_ReplicaGroup = Symbol[:replica_ids]
meta(t::Type{ReplicaGroup}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_ReplicaGroup, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct SourceTarget <: ProtoType
    source::Int64
    target::Int64
    SourceTarget(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct SourceTarget

struct __enum_PrecisionConfig_Precision <: ProtoEnum
    DEFAULT::Int32
    HIGH::Int32
    HIGHEST::Int32
    __enum_PrecisionConfig_Precision() = new(0,1,2)
end #struct __enum_PrecisionConfig_Precision
const PrecisionConfig_Precision = __enum_PrecisionConfig_Precision()

mutable struct PrecisionConfig <: ProtoType
    operand_precision::Base.Vector{Int32}
    PrecisionConfig(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct PrecisionConfig
const __pack_PrecisionConfig = Symbol[:operand_precision]
meta(t::Type{PrecisionConfig}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, __pack_PrecisionConfig, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export PrimitiveType, Format, FftType, RandomDistribution, PaddingConfig_PaddingConfigDimension, PaddingConfig, Tile, Layout, ShapeProto, ProgramShapeProto, ComputationStats, OpMetadata, ExecutionProfile, ExecutionHandle, GlobalDataHandle, DeviceHandle, ChannelHandle_ChannelType, ChannelHandle, DeviceAssignmentProto_ComputationDevice, DeviceAssignmentProto, LiteralProto, WindowDimension, Window, GatherDimensionNumbers, ScatterDimensionNumbers, ConvolutionDimensionNumbers, DotDimensionNumbers, OpSharding_Type, OpSharding, ReplicaGroup, SourceTarget, PrecisionConfig_Precision, PrecisionConfig


text:
############################
###### Boundary Types ######
############################

abstract type AbstractProblem  end
abstract type BoundaryID       end
abstract type BoundaryFunction end
abstract type BoundaryArray end

struct DirichletID <: BoundaryID id::Vector{Symbol} end
struct NeumannID   <: BoundaryID id::Vector{Symbol} end
struct RobinID     <: BoundaryID id::Vector{Symbol} end

struct DirichletFunction <: BoundaryFunction f::Function end
struct NeumannFunction   <: BoundaryFunction f::Function end
struct RobinFunction     <: BoundaryFunction f::Function end
struct ForcingFunction   <: BoundaryFunction f::Function end

struct DirichletArray <: BoundaryArray v::Vector end
struct NeumannArray   <: BoundaryArray v::Vector end
struct RobinArray     <: BoundaryArray v::Vector end
struct ForcingArray   <: BoundaryArray v::Vector end

#### constructors ####

function Dirichlet(s::Symbol ...)
  return DirichletID([i for i in s])
end

function Neumann(s::Symbol ...)
  return NeumannID([i for i in s])
end

function Robin(s::Symbol ...)
  return RobinID([i for i in s])
end

function Dirichlet(f::Function)
  return DirichletFunction(f)
end
function Dirichlet(v::Vector)
  return DirichletArray(v)
end

function Neumann(f::Function)
  return NeumannFunction(f)
end
function Neumann(v::Vector)
  return NeumannArray(v)
end

function Robin(f::Function)
  return RobinFunction(f)
end
function Robin(v::Vector)
  return RobinArray(v)
end

function Forcing(f::Function)
  return ForcingFunction(f)
end
function Forcing(v::Vector)
  return ForcingArray(v)
end

#### Methods #####

isDirichletID(x)       = (typeof(x)==DirichletID ? true : false)
isNeumannID(x)         = (typeof(x)==NeumannID   ? true : false)
isRobinID(x)           = (typeof(x)==RobinID     ? true : false)

isDirichletFunction(x) = (typeof(x)==DirichletFunction ? true : false)
isNeumannFunction(x)   = (typeof(x)==NeumannFunction   ? true : false)
isRobinFunction(x)     = (typeof(x)==RobinFunction     ? true : false)
isForcingFunction(x)   = (typeof(x)==ForcingFunction   ? true : false)

isDirichletArray(x) = (typeof(x)==DirichletArray ? true : false)
isNeumannArray(x)   = (typeof(x)==NeumannArray   ? true : false)
isRobinArray(x)     = (typeof(x)==RobinArray     ? true : false)
isForcingArray(x)   = (typeof(x)==ForcingArray   ? true : false)

############################
###### Problem Types #######
############################

mutable struct Problem <: AbstractProblem
  OperatorType::Symbol
  bcid::DictList{Int}
  bcval::DictList{Float64}
end

######  Constructors for Problem  ##########

# Assumes that only 1 dNodes and nNodes and rNodes (maximally ...)
function Problem(mesh::ScalarMesh,Nodes,bcfun,OpType)
  bcIDNodes::Vector{Vector{Int}} = []
  bcIDNamesSym::Vector{Symbol} = []
  bcValNodes::Vector{Vector{Float64}} = []
  bcValNamesSym::Vector{Symbol} = []

  # derive boundary node ID information
  dNodesSym::Vector{DirichletID} = filter(isDirichletID,Nodes)
  if length(dNodesSym)>0
    dNodes::Vector{Int} = vcat([mesh.bdry[s] for s in dNodesSym[1].id]...)
    push!(bcIDNodes,dNodes)
    push!(bcIDNamesSym,:dNodes)
  end

  nNodesSym::Vector{NeumannID} = filter(isNeumannID,Nodes)
  if length(nNodesSym)>0
    nNodes::Vector{Int} = vcat([mesh.bdry[s] for s in nNodesSym[1].id] ...)
    push!(bcIDNodes,nNodes)
    push!(bcIDNamesSym,:nNodes)
  end

  rNodesSym = filter(isRobinID,Nodes)
  if length(rNodesSym)>0
    rNodes::Vector{Int} = vcat([mesh.bdry[s] for s in rNodesSym[1].id] ...)
    push!(bcIDNodes,rNodes)
    push!(bcIDNamesSym,:rNodes)
  end

  # derive boundary value information
  # dirichlet
  dNodesFun::Vector{DirichletFunction} = filter(isDirichletFunction,bcfun)
  if length(dNodesFun)>0
    dBCarr = [dNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i in dNodes]
    push!(bcValNodes,dBCarr)
    push!(bcValNamesSym,:dBC)
  end
  dNodesArr::Vector{DirichletArray} = filter(isDirichletArray,bcfun)
  if length(dNodesArr)>0
    dBCarr = [dNodesArr[1].v[i] for i in dNodes]
    push!(bcValNodes,dBCarr)
    push!(bcValNamesSym,:dBC)
  end

  # neumann
  nNodesFun::Vector{NeumannFunction} = filter(isNeumannFunction,bcfun)
  nNodesArr::Vector{NeumannArray} = filter(isNeumannArray,bcfun)
  if (length(nNodesFun)>0) || (length(nNodesArr)>0)
    if length(nNodesFun)>0
      nBCarr = [nNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                                for i=1:length(nNodes)]
      push!(bcValNodes,nBCarr)
      push!(bcValNamesSym,:nBC)
    elseif length(nNodesArr)>0
      nBCarr = [nNodesArr[1].v[i] for i=1:length(nNodes)]
      push!(bcValNodes,nBCarr)
      push!(bcValNamesSym,:nBC)
    end
  else
    if :nNodes in bcIDNamesSym
      nBCarr = zeros(Float64,length(nNodes))
      push!(bcValNodes,nBCarr)
      push!(bcValNamesSym,:nBC)
    end
  end

  # robin
  rNodesFun::Vector{RobinFunction} = filter(isRobinFunction,bcfun)
  if length(rNodesFun)>0
    rBCarr = [rNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i in rNodes]
    push!(bcValNodes,rBCarr)
    push!(bcValNamesSym,:rBC)
  end
  rNodesArr::Vector{DirichletArray} = filter(isRobinArray,bcfun)
  if length(rNodesArr)>0
    rBCarr = [rNodesArr[1].v[i] for i in rNodes]
    push!(bcValNodes,rBCarr)
    push!(bcValNamesSym,:rBC)
  end

  # forcing
  fNodesFun::Vector{ForcingFunction} = filter(isForcingFunction,bcfun)
  fNodesArr::Vector{ForcingArray} = filter(isForcingArray,bcfun)
  if (length(fNodesFun)>0) || (length(fNodesArr)>0)
    if length(fNodesFun)>0
      farr = [fNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                                for i=1:length(mesh.xy)]
      push!(bcValNodes,farr)
      push!(bcValNamesSym,:forcing)
    elseif length(fNodesArr)>0
      farr = [fNodesArr[1].v[i] for i=1:length(mesh.xy)]
      push!(bcValNodes,farr)
      push!(bcValNamesSym,:forcing)
    end
  else
    farr = zeros(Float64,length(mesh.xy))
    push!(bcValNodes,farr)
    push!(bcValNamesSym,:forcing)
  end

  # convert all these into appropriate types for problem object
  bcIDNames = Index(bcIDNamesSym)
  bcid  = DictList{Int}(bcIDNodes,bcIDNames)

  bcValNames = Index(bcValNamesSym)
  bcval  = DictList{Float64}(bcValNodes,bcValNames)

  # return FluidProblem object
  return Problem(OpType,bcid,bcval)
end

# assumes that dUNodes = dVNodes and anything that isn't dirichlet is neumann
function Problem(mesh::FluidMesh,Nodes,bcfun,OperatorType)
  bcIDNodes::Vector{Vector{Int}} = []
  bcIDNamesSym::Vector{Symbol} = []
  bcValNodes::Vector{Vector{Float64}} = []
  bcValNamesSym::Vector{Symbol} = []

  # derive boundary node ID information
  dNodesSym::Vector{DirichletID} = filter(isDirichletID,Nodes)
  if length(dNodesSym)>0
    dUNodes::Vector{Int} = unique(vcat([mesh.bdry[s]
                                  for s in dNodesSym[1].id]...))
    push!(bcIDNodes,dUNodes)
    push!(bcIDNamesSym,:dUNodes)

    if length(dNodesSym)==2
      dVNodes::Vector{Int} = unique(vcat([mesh.bdry[s]
                                    for s in dNodesSym[2].id]...))
    else
      dVNodes = copy(dUNodes)
    end

    push!(bcIDNodes,dVNodes)
    push!(bcIDNamesSym,:dVNodes)
  end

  nNodesSym::Vector{NeumannID} = filter(isNeumannID,Nodes)
  if length(nNodesSym)==1
    nNodes::Vector{Int} = unique(vcat([mesh.bdry[s]
                                 for s in nNodesSym[1].id] ...))
    push!(bcIDNodes,nNodes)
    push!(bcIDNamesSym,:nNodes)
  end

  # derive boundary value information
  dNodesFun::Vector{DirichletFunction} = filter(isDirichletFunction,bcfun)
  if length(dNodesFun)==2
    dUBCarr = [dNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i in dUNodes]
    push!(bcValNodes,dUBCarr)
    push!(bcValNamesSym,:dUBC)

    dVBCarr = [dNodesFun[2].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i in dVNodes]
    push!(bcValNodes,dVBCarr)
    push!(bcValNamesSym,:dVBC)
  end
  dNodesArr::Vector{DirichletFunction} = filter(isDirichletArray,bcfun)
  if length(dNodesArr)==2
    dUBCarr = [dNodesArr[1].v[i] for i in dUNodes]
    push!(bcValNodes,dUBCarr)
    push!(bcValNamesSym,:dUBC)

    dVBCarr = [dNodesArr[2].v[i] for i in dVNodes]
    push!(bcValNodes,dVBCarr)
    push!(bcValNamesSym,:dVBC)
  end

  nNodesFun::Vector{NeumannFunction} = filter(isNeumannFunction,bcfun)
  if length(nNodesFun)==1
    nBCarr::Vector{Float64} = [nNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i in nNodes]
    push!(bcValNodes,nBCarr)
    push!(bcValNamesSym,:nBC)
  else
    if :nNodes in bcIDNamesSym
      nBCarr = zeros(Float64,length(nNodes))
      push!(bcValNodes,nBCarr)
      push!(bcValNamesSym,:nBC)
    end
  end

  fNodesFun::Vector{ForcingFunction} = filter(isForcingFunction,bcfun)
  fNodesArr::Vector{ForcingArray} = filter(isForcingArray,bcfun)
  if length(fNodesFun)==2
    farrX = [fNodesFun[1].f(mesh.xy[i].x,mesh.xy[i].y)
                               for i=1:length(mesh.xy)]
    farrY = [fNodesFun[2].f(node.x,node.y)
                               for node in mesh.xy]
    push!(bcValNodes,farrX)
    push!(bcValNamesSym,:forcingX)
    push!(bcValNodes,farrY)
    push!(bcValNamesSym,:forcingY)
  elseif length(fNodesArr)==2
    farrX = [fNodesArr[1].v[i] for i=1:length(mesh.xy)]
    farrY = [fNodesArr[2].v[i] for i=1:length(mesh.xy)]
    push!(bcValNodes,farrX)
    push!(bcValNamesSym,:forcingX)
    push!(bcValNodes,farrY)
    push!(bcValNamesSym,:forcingY)
  else
    farrX = zeros(Float64,length(mesh.xy))
    farrY = zeros(Float64,length(mesh.xy))
    push!(bcValNodes,farrX)
    push!(bcValNamesSym,:forcingX)
    push!(bcValNodes,farrY)
    push!(bcValNamesSym,:forcingY)
  end

  # convert all these into appropriate types for problem object
  bcIDNames = Index(bcIDNamesSym)
  bcid  = DictList{Int}(bcIDNodes,bcIDNames)

  bcValNames = Index(bcValNamesSym)
  bcval  = DictList{Float64}(bcValNodes,bcValNames)

  # return Problem object
  return Problem(OperatorType,bcid,bcval)
end

####### Problem Methods #######

bc(prob::Problem) = indexlist(prob.bcid)


text:
using .JSort
import PyPlot; const plt = PyPlot
import PyCall
#using Interactive

function GatePicker(matrix, x, y)
    pushfirst!(PyVector(pyimport("sys")."path"), "/home/erdos/gits/JSort/src/")
    interactivegating = pyimport("interactivegating")
    GP = interactivegating.GatePicker
    return GP(matrix, x, y)
end

function GatePicker(det::ParticleDetector; nbins=1000)
    (matrix, (x, y)) = histogram(det, nbins=nbins)
    gp = GatePicker(matrix, x, y)
    return () -> Gate2D(gp.extent...)
end

@enum Action Continue Quit Align
struct ActionEvent
    action::Action
    data::Vector{Any}
end

function ActionEvent(input::AbstractString)
    if input == "q"
        return ActionEvent(Quit, [])
    end
    tokens = split(input, r"\s+")
    if length(tokens) ≥ 2
        try
            data = parse.(Float64, tokens[1:2])
            for token in tokens[3:end]
                key, val = split(token, r"\s*=\s*")
                push!(data, Symbol(key) => parse(Float64, val))
            end
            return ActionEvent(Align, data)
        catch
            return ActionEvent(Continue, [])
        end
    end
    return ActionEvent(Continue, [])
end

quit(event::ActionEvent) = event.action == Quit
align(event::ActionEvent) = event.action == Align

function pickpgate(det::ParticleDetector, parameters::Parameters; nbins=1000)
    quick = QuickReader(parameters)
    read(quick)
    ungated = ChannelDetector(deepcopy(quick))
    char = ""
    default = GammaCalibrator([-913.5, 5.6])
    default(ungated)
    while true
        g = GatePicker(det, nbins=nbins)
        # Pause to let Python write to g
        event = ActionEvent(readline())
        aligner = (x...) -> ()
        if quit(event)
            break
        elseif align(event)
            aligner = (gated) -> begin
                calib = ialign(ungated, event.data[1]:event.data[2], event.data[3:end]...)
                calib(gated)
            end
        end
        @show event

        read(quick, g())
        det2 = ChannelDetector(quick)
        default(det2)
        aligner(det2)
        det2
        plot(det2)
        plot(det2 |> sum)
        plot(det2[1])
    end
end

function ialign(detector::ChannelDetector, region; nbins=5000, kwargs...)
    full, edges = histogram(detector, nbins=nbins)
    chns = goodchannels(detector)

    signal = pyimport("scipy.signal")
    smoother(x) = signal.savgol_filter(x, 51, 4)

    coeff, _ = alignspectra(full, edges, lowregion=region, lowsmoother=smoother,
                            lownumregions=get(kwargs, :nregions, 2),
                            lowregionwidth=get(kwargs, :swidth, 200),
                            lowsearchwidth=get(kwargs, :lwidth, 60),
                            plot=get(kwargs, :plot, false), order=1)
    coefficients = repeat([0.0 1.0], outer=(33, 1))
    coefficients[chns, :] .= hcat(coeff...)'
    calib = GammaCalibrator(coefficients)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using TupleTools
export i_einsum!

"""
	i_einsum!(ixs, xs, iy, y::AbstractArray{T})

A naive reversible implementation of `i_einsum` function for tropical numbers.
    * `ixs`: input tensor indices,
    * `xs`: input tensors,
    * `iy`: output tensor indices,
    * `y`: accumulated tensor, notice it is initialized to 0 as output!

# NOTE: this function is general purposed and slow!
"""
@i function i_einsum!(ixs, xs, iy, y::AbstractArray{T}) where {T<:Tropical}
	@routine @invcheckoff begin
	    # outer legs and inner legs
	    outer_indices ← unique(iy)
	    inner_indices ← setdiff(TupleTools.vcat(ixs...), outer_indices)

	    # find size for each leg
	    all_indices ← TupleTools.vcat(ixs..., iy)
	    all_sizes ← TupleTools.vcat(size.(xs)..., size(y))
	    outer_sizes ← [map(i->all_sizes[i], indexin(outer_indices, [all_indices...]))...]
	    inner_sizes ← [map(i->all_sizes[i], indexin(inner_indices, [all_indices...]))...]

	    # cartesian indices for outer and inner legs
	    outer_ci ← CartesianIndices((outer_sizes...,))
	    inner_ci ← CartesianIndices((inner_sizes...,))

	    # for indexing tensors (leg binding)
	    indices ← (outer_indices..., inner_indices...)
	    locs_xs ← map(ix->map(i->findfirst(isequal(i), indices), ix), ixs)
    	locs_y ← map(i->findfirst(isequal(i), outer_indices), iy)
	end
    i_loop!(locs_xs, xs, locs_y, y, outer_ci, inner_ci)
	~@routine
end

"""take an index subset from `ind`"""
index_map(ind::CartesianIndex, locs::Tuple) = CartesianIndex(TupleTools.getindices(Tuple(ind), locs))

"""
loop and accumulate products to y, the GPU version, the CPU version.
"""
@i function i_loop!(locs_xs::NTuple{N}, xs::NTuple{N, AbstractArray}, locs_y, y::AbstractArray{T}, outer_ci::CartesianIndices, inner_ci::CartesianIndices) where {N, T<:Tropical}
    @invcheckoff @inbounds for i in outer_ci
		@routine begin
            ind_y ← outer_ci[i]
            iy ← index_map(ind_y, locs_y)
            branch_keeper ← zeros(Bool, length(inner_ci))
            pl ← ones(T, length(inner_ci))
            el ← zero(T)
            for ind_x in inner_ci
                @routine ind_xy ← CartesianIndex(TupleTools.vcat(ind_y.I, ind_x.I))
                #y[iy] += map_prod(T, xs, ind_xy, locs_xs)
                for I=1:N
                    pl[ind_x] *= xs[I][index_map(ind_xy, locs_xs[I])]
                end
                if (el.n < pl[ind_x].n, branch_keeper[ind_x])
                    FLIP(branch_keeper[ind_x])
                    NiLang.SWAP(el, pl[ind_x])
                end
                ~@routine
            end
        end
        y[iy] *= el
        ~@routine
    end
end


text:
#!/usr/bin/env julia
# File              : .julia/dev/VSFlow/src/VortexElement.jl
# Author            : Denis N Dumoulin <denis.dumoulin@uclouvain.be>
# Date              : 13.03.2019
# Last Modified Date: 17.09.2021
include("Utils.jl")

const pi2_inv = .5 / pi
const sq3_inv = 1 / sqrt(3.)
const w14 = (18 - sqrt(30)) / 36
const w23 = (18 + sqrt(30)) / 36
const X14 = sqrt((15 + 2sqrt(30)) / 35)
const X23 = sqrt((15 - 2sqrt(30)) / 35)

"""
    VortexPoint

Structure representing a point vortex.

# Fields
 - `X`: X-position of the point vortex in the inertial frame.
 - `Y`: Y-position of the point vortex in the inertial frame.
 - `Γ`: Circulation of the point vortex.
 - `dv`: Velocity correction in case of lumping.
"""
mutable struct VortexPoint
	X::Float64
	Y::Float64
	Γ::Float64
	dv::Vector{Float64}
end

VortexPoint(X, Y, Γ) = VortexPoint(X, Y, Γ, [0., 0.])

"""
    blobkernel(X, Y, vpX, vpY, δ, isgaussian)

Return the velocity at point `X`, `Y` induced by a vortex point located at `vpX`,`vpY`.
"""
function blobkernel(X, Y, vpX, vpY, δ, isgaussian)
	## Low order algebraic i.e. Blob
    !isgaussian && (coeff = pi2_inv / norm2(X - vpX, Y - vpY, δ))

    ## Gaussian regularization kernel
    isgaussian && begin
        x2 = norm2(X - vpX, Y - vpY)
        abs(x2) < 1e-12 && return zeros(2)
        coeff = pi2_inv * (1 - exp(-x2 / δ^2)) / x2
    end

	return [(vpY - Y), (X - vpX)] .* coeff

end

"""
    getinducedvelocity(vp::VortexPoint, X, Y, δ)

Return the velocity intuced by vortex `vp` at location `X`, `Y`.
"""
function getinducedvelocity(vp::VortexPoint, X, Y, δ, isgaussian)
	return vp.Γ * blobkernel(X, Y, vp.X, vp.Y, δ, isgaussian)
end

"""
    getcloudvelocities(vortex_points, δ, isgaussian)

Return the velocity induced by all other vortices in the cloud.
"""
function getcloudvelocities(vortex_points, δ, isgaussian)
	Nv = length(vortex_points)
	U = zeros(Nv, Nv)
	V = zeros(Nv, Nv)
	for i = 1:Nv
		@inbounds vi = vortex_points[i]
		@inbounds for j = i + 1:Nv
			@inbounds vj = vortex_points[j]
	 		@inbounds U[i, j], V[i, j] = blobkernel(vi.X, vi.Y, vj.X, vj.Y, δ, isgaussian)
			@inbounds U[j, i], V[j, i] = -U[i, j], -V[i, j]
		end
	end
	Γs = map(vp -> vp.Γ, vortex_points)
	return map((u, v) -> [u, v], U * Γs, V * Γs)
end

"""
    getvelocitygradient(vp::VortexPoint, X, Y, θ)

Computes the gradient of induced velocity at position `X`, `Y`
due to a change of position of `vp`.

# Arguments
 - `vp` : sensitive point vortex
 - `X`  : X position where velocity is computed
 - `Y`  : Y position where velocity is computed
 - `θ`  : rotation angle
 - `δ`  : Gaussian kernel cut-off width
"""
function getvelocitygradient(vp::VortexPoint, X, Y, θ, δ)
	coeff = -pi2_inv / norm2(X - vp.X, Y - vp.Y, δ)^2
	bx = (((X - vp.X)^2 - norm2(Y - vp.Y, δ))cos(θ) - 2(vp.X - X) * (vp.Y - Y)sin(θ))coeff
	by = ((norm2(X - vp.X, δ) - (Y - vp.Y)^2)sin(θ) + 2(vp.X - X) * (vp.Y - Y)cos(θ))coeff
	return [bx, by]
end

abstract type Panel
end

"""
    BuildPanel(T, X1, Y1, X2, Y2, profID)

Build a panel of type `T` starting at point `X1`, `Y1` and ending at point `X2`, `Y2`.
"""
function BuildPanel(T, X1, Y1, X2, Y2, profID)
	Xc = .5 * (X1 + X2)
	Yc = .5 * (Y1 + Y2)
	b = .5norm([X2 - X1, Y2 - Y1])
	θ = atan(Y2 - Y1, X1 - X2)
	return T(X1, Y1, X2, Y2, Xc, Yc, b, θ, profID)
end

"""
    inpanelframe(p::Panel, X, Y)

Return the point `X`, `Y` in the frame of panel `p`.
"""
function inpanelframe(p::Panel, X, Y)
	X, Y = rotate(X - p.Xc, Y - p.Yc, π + p.θ)
end

"""
    inreferenceframe(p::Panel, X, Y)

Return the point `x`, `y` in the reference frame.
"""
function inreferenceframe(p::Panel, x, y)
	X, Y = rotate(x, y, -(π + p.θ))
	return X + p.Xc, Y + p.Yc
end

"""
    cropegde!(p::Panel, ϵ, isend)

Shortens the panel by a length `ϵ` at its end if `isend == true`
or at its start otherwise.
"""
function cropedge!(p::Panel, ϵ, isend)
	p.b -= .5ϵ
	if isend
		p.X2 = p.X1 + 2(p.b * cos(pi - p.θ))
		p.Y2 = p.Y1 + 2(p.b * sin(pi - p.θ))
	else
		p.X1 = p.X2 + 2(p.b * cos(p.θ))
		p.Y1 = p.Y2 - 2(p.b * sin(p.θ))
	end
	p.Xc = .5(p.X1 + p.X2)
	p.Yc = .5(p.Y1 + p.Y2)
end

    """
    setangle!(p::Panel, α, pivotlocation)

Rotate panel `p` by an angle `α` with rotation center located at `pivotlocation`.
"""
function setangle!(p::Panel, α, pivotlocation)
	x1 = p.X1 - pivotlocation[1]
	y1 = p.Y1 - pivotlocation[2]
	x2 = p.X2 - pivotlocation[1]
	y2 = p.Y2 - pivotlocation[2]
	p.X1, p.Y1 = pivotlocation + rotate(x1, y1, -α)
	p.X2, p.Y2 = pivotlocation + rotate(x2, y2, -α)
	p.Xc = .5(p.X1 + p.X2)
	p.Yc = .5(p.Y1 + p.Y2)
	p.θ = atan(p.Y2 - p.Y1, p.X1 - p.X2)
end

"""
    translate!(p::Panel, Xt, Yt)

Translate panel `p`.
"""
function translate!(p::Panel, Xt, Yt)
	p.X1 += Xt
	p.Y1 += Yt
	p.X2 += Xt
	p.Y2 += Yt
	p.Xc += Xt
	p.Yc += Yt
end

"""
    ConstantPanel

Structure representing a vortex panel of uniform strength.

# Fields
 - `X1`: X-position of the start of the panel in the inertial frame.
 - `Y1`: Y-position of the start of the panel in the inertial frame.
 - `X2`: X-position of the end of the panel in the inertial frame.
 - `Y2`: Y-position of the end of the panel in the inertial frame.
 - `Xc`: X-position of the center of the panel in the inertial frame.
 - `Yc`: Y-position of the center of the panel in the inertial frame.
 - `b`: half-length of the panel.
 - `θ`: panel orientation with respect to the X axis of the inertial frame.
 - `profID`: ID of the profile the panel belongs to.
"""
mutable struct ConstantPanel <: Panel
	X1::Float64
	Y1::Float64
	X2::Float64
	Y2::Float64
	Xc::Float64
	Yc::Float64
	b::Float64
	θ::Float64
	profID::String
end

"""
    getparams(p::ConstantPanel, X, Y)

Return the velocity induced by panel `p` if it had unit strength at location `X`, `Y`.
"""
function getparams(p::ConstantPanel, X, Y)
	x, y = inpanelframe(p, X, Y)
	u = (atan((x - p.b) / y) - atan((x + p.b) / y))pi2_inv
	v = .5pi2_inv * log(norm2(x + p.b, y) / norm2(x - p.b, y))
	return u, v
end

"""
    LinearPanel

Structure representing a vortex panel of linearly varying strength.

# Fields
 - `X1`: X-position of the start of the panel in the inertial frame.
 - `Y1`: Y-position of the start of the panel in the inertial frame.
 - `X2`: X-position of the end of the panel in the inertial frame.
 - `Y2`: Y-position of the end of the panel in the inertial frame.
 - `Xc`: X-position of the center of the panel in the inertial frame.
 - `Yc`: Y-position of the center of the panel in the inertial frame.
 - `b`: half-length of the panel.
 - `θ`: panel orientation with respect to the X axis of the inertial frame.
 - `profID`: ID of the profile the panel belongs to.
"""
mutable struct LinearPanel <: Panel
	X1::Float64
	Y1::Float64
	X2::Float64
	Y2::Float64
	Xc::Float64
	Yc::Float64
	b::Float64
	θ::Float64
	profID::String
end

"""
    getparams(p::LinearPanel, X, Y)

Obtain the parameters to multiply with γL, γR in order to get the
velocity at point `X`, `Y` induced by linear vortex panel `p`.
"""
function getparams(p::LinearPanel, X, Y)
	(p.b == 0) && return [zeros(2), zeros(2)]
	x, y = inpanelframe(p, X, Y)
	b = p.b
	xb = x / b
	if abs(y) != 0
		y4b = y / 4b
		logstuff = log(norm2(x - b, y) / norm2(x + b, y))
		atanstuff =  atan((x - b) / y) - atan((x + b) / y)
		u1 = (y4b * logstuff - .5atanstuff * (xb - 1.))pi2_inv
		u2 = -(y4b * logstuff - .5atanstuff * (xb + 1.))pi2_inv
		v1 = (1. + 2y4b * atanstuff + .25(xb - 1.) * logstuff)pi2_inv
		v2 = -(1. + 2y4b * atanstuff + .25(xb + 1.) * logstuff)pi2_inv
	else
		logstuff = log((b - x)^2 / (b + x)^2)
		v1 = (1. + .25(xb - 1.)logstuff)pi2_inv
		v2 = -(1. + .25(xb + 1.)logstuff)pi2_inv
		u1 = 0
		u2 = 0
	end
	return [[u1, u2], [v1, v2]]
end

"""
    getmotionkelvin(p::LinearPanel, ω, Ut, Vt, X0, Y0)

Get the circulation induced by the motion of panel `p`.

# Arguments
 - `p`  : the linear panel to integrate on
 - `ω`  : angular velocity in the z direction
 - `Ut` : translational horizontal body velocity
 - `Vt` : translational vertical body velocity
 - `X0` : position X of the center of rotation of the body
 - `Y0` : position Y of the center of rotation of the body
"""
function getmotionkelvin(p::LinearPanel, ω, Ut, Vt, X0, Y0)
	x0, y0 = inpanelframe(p, X0, Y0)
	ut, vt = rotate(Ut, Vt, π + p.θ)
	return 2p.b * (ut + ω * y0)
end

"""
    getmotionparams(p::LinearPanel, X, Y)

Get the velocity parameters at point `X`, `Y` due to the motion of panel `p`.
Integrals are performed in the panel frame of reference then re-cast
in the inertial frame through a rotation.

# Arguments
 - `p`  : the linear panel to integrate on
 - `X`  : position X where we want the velocity induced by the body motion
 - `Y`  : position Y where we want the velocity induced by the body motion
 - `ω`  : angular velocity in the z direction
 - `Ut` : translational horizontal body velocity
 - `Vt` : translational vertical body velocity
 - `X0` : position X of the center of rotation of the body
 - `Y0` : position Y of the center of rotation of the body
 - `vs` : stretching velocity
"""
function getmotionparams(p::LinearPanel, X, Y, ω, Ut, Vt, X0, Y0, vs=0)
	(p.b == 0) && return zeros(2)
	x, y = inpanelframe(p, X, Y)
	x0, y0 = inpanelframe(p, X0, Y0)
	ut, vt = rotate(Ut, Vt, π + p.θ)
	k = ut + ω * (y + y0) + .5vs * x / p.b
	l = vt + ω * (x - x0) - .5vs * y / p.b
	atanstuff = 0
    logstuff = .5log(norm2(x - p.b, y) / norm2(x + p.b, y))
    y != 0 && (atanstuff = atan((x - p.b) / y) - atan((p.b + x) / y))
	u = (atanstuff * k + l * logstuff + 2ω * p.b)pi2_inv
	v = (l * atanstuff - k * logstuff)pi2_inv
	return rotate(u, v, -π - p.θ)
end

"""
    getinducedvelocity(p::LinearPanel, XY, k, γL, γR, ω, Ut, Vt, Xp, Yp, vs)

Computes the velocity due to panel `p` at location `XY`
`k` specifies which integral is required:
 - `0` full integral; vortex
 - `1` full integral; source
 - `2` half-left integral; vortex
 - `3` half-left integral; source
 - `4` half-right integral; vortex
 - `5` half-right integral; source
"""
function getinducedvelocity(p::LinearPanel, XY, k, γL, γR, ω, Ut, Vt, Xp, Yp, vs=0)
    @assert 0 ≤ k < 6
	x, y = inpanelframe(p, XY...)
	xp, yp = inpanelframe(p, Xp, Yp)
	ut, vt = rotate(Ut, Vt, π + p.θ)
    u = ut + ω * yp
    v = vt - ω * xp
    function p1(α, β, ξ)
        atanstuff = 0
        y != 0 && (atanstuff = atan((x - ξ) / y))
        logstuff = log(norm2(x - ξ, y))
        return [(α * x + β) * atanstuff - .5α * y * logstuff;
                -.5(α * x + β) * logstuff + α * ((x - ξ) - y * atanstuff)] .* pi2_inv
    end
    iseven(k) && begin
        α = .5(γR - γL + vs) / p.b
        β = .5(γR + γL + 2u)
        k == 0 && return rotate((p1(α, β, p.b) - p1(α, β, -p.b))..., -π - p.θ)
        x0 = -β / α
        k == 2 && return rotate((p1(α, β, x0) - p1(α, β, -p.b))..., -π - p.θ)
        k == 4 && return rotate((p1(α, β, p.b) - p1(α, β, x0))..., -π - p.θ)
    end
    p2(α, β, ξ) = [0 -1;1 0] * p1(α, β, ξ)
    k == 1 && return rotate((p2(ω, v, p.b) - p2(ω, v, -p.b))..., -π - p.θ)
    x0 = -v / ω
    k == 3 && return rotate((p2(ω, v, x0) - p2(ω, v, -p.b))..., -π - p.θ)
    k == 5 && return rotate((p2(ω, v, p.b) - p2(ω, v, x0))..., -π - p.θ)
end

"""
    getimpulse(p::LinearPanel, γL, γR, ω, Ut, Vt, X0, Y0, Xp, Yp, nc, vs=0)

Compute the impulse by integrating X^γ and X^X^γ analytically on the linear panel.
You can specify if you need to add the integral of X^(n^Vb) where n
points towards the fluid.
The integral is performed in the panel frame and the result is then
rotated in the inertial frame of reference.
Inertial frame: `X`,
Panel frame:   `-x + x'`

# Arguments
 - `p`  : the linear panel to integrate on
 - `γL` : vortex strength on the left side of the panel (in its frame)
 - `γR` : vortex strength on the right side of the panel
 - `ω`  : angular velocity in the z direction
 - `Ut` : translational horizontal body velocity
 - `Vt` : translational vertical body velocity
 - `X0` : position X of origin of the inertial frame
 - `Y0` : position Y of origin of the inertial frame
 - `Xp` : position X of the center of rotation of the body
 - `Yp` : position Y of the center of rotation of the body
 - `nc` : need body velocity contribution?
 - `vs` : stretching velocity
"""
function getimpulse(p::LinearPanel, γL, γR, ω, Ut, Vt, X0, Y0, Xp, Yp, nc, vs=0)
	Θ = π + p.θ
	x, y = inpanelframe(p, X0, Y0)
	_, yp = inpanelframe(p, Xp, Yp)
	ut, _ = rotate(Ut, Vt, Θ)

	b3 = p.b / 3
	xy2 = norm2(x, y)
	Σ = p.b * (γR + γL)
	Δ = p.b * (γR - γL)
	m = -.5b3 * (p.b * Σ - 2x * Δ) - .5Σ * xy2
	Δ *= b3
	!nc && return [rotate(-y * Σ, x * Σ + Δ, -Θ)..., m]

	vb = (ut + ω * yp) * p.b
	b23 = p.b * b3
	vsb3 = vs * b23
	Σ += 2vb
	return [rotate(-y * Σ, x * Σ - vsb3 + Δ, -Θ)..., m - vb * (xy2 + b23) - x * vsb3]
end

"""
    getmonopolevortex(p::LinearPanel, id, γL, γR, ω, Ut, Vt, Xp, Yp, vs=0)

Compute the monopole of vorticity of a linear panel and its location in the reference frame.
If vortex strength crosses the x axis, then the panel is split into to smaller
panels of positive and negative strengths.
In the case it is split, then the id is modified to be able to retrace it later.

# Arguments
 - `p`  : the linear panel to integrate on
 - `id` : the id of the panel
 - `γL` : vortex strength on the left side of the panel (in its frame)
 - `γR` : vortex strength on the right side of the panel
"""
function getmonopolevortex(p::LinearPanel, id, γL, γR, ω, Ut, Vt, Xp, Yp, vs=0)
	Θ = π + p.θ
	yp, yp = inpanelframe(p, Xp, Yp)
	ut, _ = rotate(Ut, Vt, Θ)
    u = ut + ω * yp
    γL += u - .5vs
    γR += u + .5vs

    γL * γR ≥ 0 && return [[getcentroid(p, γL, γR)..., 1im * p.b * (γR + γL), id]]
    #= println("vortex panel $id is split") =#

    x0 = p.b * (γL + γR) / (γL - γR)
    xc1 = (x0^3 - 3x0 * p.b^2 - 2p.b^3) / 3(p.b + x0)^2
    xc2 = (x0^3 - 3x0 * p.b^2 + 2p.b^3) / 3(p.b - x0)^2
    m1 = [inreferenceframe(p, xc1, 0)..., .5im * γL * (p.b + x0), -id]
    m2 = [inreferenceframe(p, xc2, 0)..., .5im * γR * (p.b - x0), -id * 1im]
    return [m1, m2]
end

"""
    getmonopolesource(p::LinearPanel, id, ω, Ut, Vt, Xp, Yp)

Compute the monopole of source of a linear panel and its location in the reference frame.
If vortex strength crosses the x axis, then the panel is split into to smaller
panels of positive and negative strengths.
In the case it is split, then the id is modified to be able to retrace it later.

# Arguments
 - `p`  : the linear panel to integrate on
 - `id` : the id of the panel
"""
function getmonopolesource(p::LinearPanel, id, ω, Ut, Vt, Xp, Yp)
    @assert id < 1e9
	Θ = π + p.θ
	xp, _ = inpanelframe(p, Xp, Yp)
	_, vt = rotate(Ut, Vt, Θ)
    v = vt - ω * xp
    x0 = -v / ω

    abs(x0) ≥ p.b &&
        return [[inreferenceframe(p, ω * p.b^2 / 3v, 0)..., 2p.b * v, id]]

    #= println("source panel $id is split") =#
    xc1 = (x0^3 - 3x0 * p.b^2 - 2p.b^3) / 3(p.b + x0)^2
    xc2 = (x0^3 - 3x0 * p.b^2 + 2p.b^3) / 3(p.b - x0)^2
    m1 = [inreferenceframe(p, xc1, 0)..., .5 * (p.b + x0) * (v - ω * p.b), -id]
    m2 = [inreferenceframe(p, xc2, 0)..., .5 * (p.b - x0) * (v + ω * p.b), -id * 1im]
    return [m1, m2]
end

"""
    getcentroid(p::LinearPanel, γL, γR)

Return the absolute position (in reference frame) of the vorticity centroid of panel `p`.
"""
function getcentroid(p::LinearPanel, γL, γR)
	xc = p.b * (γR - γL) / 3(γR + γL)
	return inreferenceframe(p, xc, 0)
end

"""
    getpressureforce(p::LinearPanel, p1, p2, Xp, Yp)

Return the force on panel `p` due to pressure distribution.
"""
function getpressureforce(p::LinearPanel, p1, p2, Xp, Yp)
	x, _ = inpanelframe(p, Xp, Yp)
	[rotate(0, (p1 + p2) * p.b, -π - p.θ)..., p.b * ((p2 - p1) * p.b / 3 - (p1 + p2) * x)]
end

"""
    getpotential(p::LinearPanel, γL, γR, ω, Ut, Xp, Yp)

Compute the potential through the integral of tangent velocity
along the panel `p`.
"""
function getpotential(p::LinearPanel, γL, γR, ω, Ut, Vt, Xp, Yp)
	_, yp = inpanelframe(p, Xp, Yp)
	ut, _ = rotate(Ut, Vt, π + p.θ)

	vv = ut + ω * yp
	return p.b * ((γL + γR) + 2vv)
end

"""
    getenergyflux(p::LinearPanel, γL, γR, ω, Ut, Vt, X0, Y0, Xp, Yp, vs=0)

Return the flux of energy at panel `p` to obtain the force/moment with a compact
formulation depending on velocity and its derivatives.
"""
function getenergyflux(p::LinearPanel, γL, γR, ω, Ut, Vt, X0, Y0, Xp, Yp, vs=0)
	xp, yp = inpanelframe(p, Xp, Yp)
	x, y = inpanelframe(p, X0, Y0)
	ut, vt = rotate(Ut, Vt, π + p.θ)

	αb = .5(γR - γL + vs)
	uu = .5(γL + γR) + ut + ω * yp
	vv = vt - ω * xp
	ωb = ω * p.b

	αωb2 = 2αb * ωb
	Δb = (ωb^2 - αb^2)
	uv2 = 2vv * uu

	I1 = p.b * (αωb2 / 3 + uv2)
	I2 = p.b * (Δb / 3 + (vv + uu) * (vv - uu))
	I3 = p.b * ((αωb2 * y - Δb * x + 2p.b * (ωb * vv - αb * uu)) / 3 + y * uv2 - x * (vv + uu) * (vv - uu))
	return [rotate(I1, I2, -π - p.θ)..., I3]
end

"""
    getangularimpulsecorrection(p::LinearPanel, ω, Up, Vp, Ut, Vt, X0, Y0, Xp, Yp)

Return the angular impulse correction at panel `p` to obtain the moment with a compact
formulation depending on velocity and its derivatives.
This correction is necessary if integrals are computed in the body reference frame.
"""
function getangularimpulsecorrection(p::LinearPanel, ω, Up, Vp, Ut, Vt, X0, Y0, Xp, Yp)
	xp, _ = inpanelframe(p, Xp, Yp)
	x, y = inpanelframe(p, X0, Y0)
	_, vp = rotate(Up, Vp, π + p.θ)
	ut, vt = rotate(Ut, Vt, π + p.θ)

	vv = vp - ω * xp
	return -2p.b * (p.b * p.b * vt * ω / 3 + (vt * x - y * ut) * vv)
end

"""
    getbodymomentofinertia(p::LinearPanel, X0, Y0)

Return the contribution of panel `p` on the moment of inertia of
the body.
"""
function getbodymomentofinertia(p::LinearPanel, X0, Y0)
	x, y = inpanelframe(p, X0, Y0)
	return .5 * y * p.b * ((x^2 + y^2) + p.b^2 / 3)
end

"""
    getbodyvolume(p::LinearPanel, X0, Y0)

Return the contribution of panel `p` on the body volume.
"""
function getbodyvolume(p::LinearPanel, X0, Y0)
	_, y = inpanelframe(p, X0, Y0)
	return y * p.b
end

"""
    getbodycentroid(p::LinearPanel, X0, Y0)

Returns the contribution of panel `p` on the position of the
centroid of the body with respect to position `X0`, `Y0`,
multiplied by the body volume.
"""
function getbodycentroid(p::LinearPanel, X0, Y0)
	x, y = inpanelframe(p, X0, Y0)
	return rotate(0, -p.b * ((x^2 + y^2) + p.b^2 / 3), -π - p.θ)
end

"""
    getlinearspeedcoeffs(p::LinearPanel, o::LinearPanel, vertical)

Return the velocity parameters induced by linear panel `o` at center of panel `p`.
If `vertical == true`, then only the normal components are returned.
"""
function getlinearspeedcoeffs(p::LinearPanel, o::LinearPanel, vertical)
	u, v = getparams(o, p.Xc, p.Yc)
	uv1 = rotate(u[1], v[1], p.θ - o.θ)
	uv2 = rotate(u[2], v[2], p.θ - o.θ)
	vertical && return uv1[2], uv2[2]
	return uv1[1], uv2[1]
end

"""
    getconstantspeedcoeffs(p::LinearPanel, o::ConstantPanel, vertical)

Return the velocity parameters induced by uniform panel `o` at center of panel `p`.
If `vertical == true`, then only the normal components are returned.
"""
function getconstantspeedcoeffs(p::LinearPanel, o::ConstantPanel, vertical)
	u, v = getparams(o, p.Xc, p.Yc)
	uv = rotate(u, v, p.θ - o.θ)
	vertical && return uv[2]
	return uv[1]
end

"""
    getnormal(p::LinearPanel)

Get the normal of the panel in its frame of reference (inside of the body).
"""
function getnormal(p::LinearPanel)
	return [-sin(p.θ), -cos(p.θ)]
end

"""
    getγ(p::LinearPanel, γ1, γ2, X, Y)

Return the vortex strength of panel `p` at point `X`, `Y`.
This point must be part of the panel.
"""
function getγ(p::LinearPanel, γ1, γ2, X, Y)
	x, y = inpanelframe(p, X, Y)
	@assert	abs(y) < 1e-12
	return (γ2 - γ1) * x / 2p.b + .5(γ2 + γ1)
end

"""
    get∫fdV(pa::LinearPanel, XY, f)

Computes the integral of `f` on the triangle created by the
panel and point `XY` using a Gauss-Legendre quadrature.
`f` can return a vector.
"""
function get∫fdV(pa::LinearPanel, XY, f)
	J =	(XY[1] - pa.X2) * (pa.Y1 - pa.Y2) - (pa.X1 - pa.X2) * (XY[2] - pa.Y2)
	x(u, v)	= pa.X2 + u * (XY[1] - pa.X2) + v * (pa.X1 - pa.X2)
	y(u, v)	= pa.Y2 + u * (XY[2] - pa.Y2) + v * (pa.Y1 - pa.Y2)
	Is = gausslegendretriangle(f, x, y)
	return abs(J) .* Is
end


text:
"""
    UniversalFallback

The `UniversalFallback` can be applied on a `ModelLike` `model` to create the
model `UniversalFallback(model)` supporting *any* constaints and attributes.
This allows to have a specialized implementation in `model` for performance
critical constraints and attributes while still supporting other attributes
with a small performance penalty. Note that `model` is unaware of constraints
and attributes stored by `UniversalFallback` so this is not appropriate if
`model` is an optimizer (for this reason, `optimize!` has not been
implemented). In that case, optimizer bridges should be used instead.
"""
mutable struct UniversalFallback{MT} <: MOI.ModelLike
    model::MT
    constraints::Dict{Tuple{DataType, DataType}, Dict} # See https://github.com/JuliaOpt/JuMP.jl/issues/1152
    nextconstraintid::Int64
    con_to_name::Dict{CI, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
    optattr::Dict{MOI.AbstractOptimizerAttribute, Any}
    modattr::Dict{MOI.AbstractModelAttribute, Any}
    varattr::Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}
    conattr::Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}
    function UniversalFallback{MT}(model::MOI.ModelLike) where {MT}
        new{typeof(model)}(model,
                           Dict{Tuple{DataType, DataType}, Dict}(),
                           0,
                           Dict{CI, String}(),
                           nothing,
                           Dict{MOI.AbstractOptimizerAttribute, Any}(),
                           Dict{MOI.AbstractModelAttribute, Any}(),
                           Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}(),
                           Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}())
    end
end
UniversalFallback(model::MOI.ModelLike) = UniversalFallback{typeof(model)}(model)

MOI.is_empty(uf::UniversalFallback) = MOI.is_empty(uf.model) && isempty(uf.constraints) && isempty(uf.optattr) && isempty(uf.modattr) && isempty(uf.varattr) && isempty(uf.conattr)
function MOI.empty!(uf::UniversalFallback)
    MOI.empty!(uf.model)
    empty!(uf.constraints)
    uf.nextconstraintid = 0
    empty!(uf.con_to_name)
    uf.name_to_con = nothing
    empty!(uf.optattr)
    empty!(uf.modattr)
    empty!(uf.varattr)
    empty!(uf.conattr)
end
function MOI.copy_to(uf::UniversalFallback, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(uf, src; kws...)
end
function supports_default_copy_to(uf::UniversalFallback, copy_names::Bool)
    return supports_default_copy_to(uf.model, copy_names)
end

# References
MOI.is_valid(uf::UniversalFallback, idx::VI) = MOI.is_valid(uf.model, idx)
function MOI.is_valid(uf::UniversalFallback, idx::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.is_valid(uf.model, idx)
    else
        haskey(uf.constraints, (F, S)) && haskey(uf.constraints[(F, S)], idx)
    end
end
function MOI.delete(uf::UniversalFallback, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.delete(uf.model, ci)
    else
        if !MOI.is_valid(uf, ci)
            throw(MOI.InvalidIndex(ci))
        end
        delete!(uf.constraints[(F, S)], ci)
        if haskey(uf.con_to_name, ci)
            delete!(uf.con_to_name, ci)
        end
        uf.name_to_con = nothing
    end
    for d in values(uf.conattr)
        delete!(d, ci)
    end
end
function MOI.delete(uf::UniversalFallback, vi::VI)
    MOI.delete(uf.model, vi)
    for d in values(uf.varattr)
        delete!(d, vi)
    end
    for (FS, constraints) in uf.constraints
        for (ci, constraint) in constraints
            f, s = constraint
            if f isa MOI.SingleVariable
                if f.variable == vi
                    delete!(constraints, ci)
                end
            else
                constraints[ci] = removevariable(f, s, vi)
            end
        end
    end
end

# Attributes
_get(uf, attr::MOI.AbstractOptimizerAttribute)          = uf.optattr[attr]
_get(uf, attr::MOI.AbstractModelAttribute)              = uf.modattr[attr]
function _get(uf, attr::MOI.AbstractVariableAttribute, vi::VI)
    return get(uf.varattr[attr], vi, nothing)
end
function _get(uf, attr::MOI.AbstractConstraintAttribute, ci::CI)
    return get(uf.conattr[attr], ci, nothing)
end
function MOI.get(uf::UniversalFallback,
                 attr::Union{MOI.AbstractOptimizerAttribute,
                             MOI.AbstractModelAttribute})
    if !MOI.is_copyable(attr) || MOI.supports(uf.model, attr)
        MOI.get(uf.model, attr)
    else
        _get(uf, attr)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::Union{MOI.AbstractVariableAttribute,
                             MOI.AbstractConstraintAttribute}, idx::MOI.Index)
    if MOI.supports(uf.model, attr, typeof(idx))
        MOI.get(uf.model, attr, idx)
    else
        _get(uf, attr, idx)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.NumberOfConstraints{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, attr)
    else
        return length(get(uf.constraints, (F, S), Dict{CI{F, S}, Tuple{F, S}}()))
    end
end
function MOI.get(uf::UniversalFallback,
                 listattr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, listattr)
    else
        collect(keys(get(uf.constraints, (F, S), Dict{CI{F, S}, Tuple{F, S}}())))
    end
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfConstraints)
    list = MOI.get(uf.model, listattr)
    for (FS, constraints) in uf.constraints
        if !isempty(constraints)
            push!(list, FS)
        end
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfOptimizerAttributesSet)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.optattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfModelAttributesSet)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.modattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfVariableAttributesSet)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.varattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfConstraintAttributesSet{F, S}) where {F, S}
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.conattr)
        push!(list, attr)
    end
    return list
end

# Name
# The names of constraints not supported by `uf.model` need to be handled
function MOI.set(uf::UniversalFallback, attr::MOI.ConstraintName, ci::CI{F, S}, name::String) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, attr, ci, name)
    else
        uf.con_to_name[ci] = name
        uf.name_to_con = nothing # Invalidate the name map.
    end
    return
end
function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintName, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, attr, ci)
    else
        return get(uf.con_to_name, ci, EMPTYSTRING)
    end
end

MOI.get(uf::UniversalFallback, ::Type{VI}, name::String) = MOI.get(uf.model, VI, name)
function MOI.get(uf::UniversalFallback, ::Type{CI{F, S}}, name::String) where {F, S}
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end

    if MOI.supports_constraint(uf.model, F, S)
        ci = MOI.get(uf.model, CI{F, S}, name)
        if ci !== nothing && haskey(uf.name_to_con, name)
            error("Multiple constraints have the name $name.")
        end
        return ci
    else
        ci = get(uf.name_to_con, name, nothing)
        if ci == CI{Nothing, Nothing}(-1)
            error("Multiple constraints have the name $name.")
        elseif ci isa CI{F, S}
            return ci
        else
            return nothing
        end
    end
end
function MOI.get(uf::UniversalFallback, ::Type{CI}, name::String)
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end

    ci = MOI.get(uf.model, CI, name)
    if ci === nothing
        uf_ci = get(uf.name_to_con, name, nothing)
        if uf_ci == CI{Nothing, Nothing}(-1)
            error("Multiple constraints have the name $name.")
        else
            return uf_ci
        end
    else
        if haskey(uf.name_to_con, name)
            error("Multiple constraints have the name $name.")
        end
        return ci
    end
end

_set(uf, attr::MOI.AbstractOptimizerAttribute, value) = uf.optattr[attr] = value
_set(uf, attr::MOI.AbstractModelAttribute, value)     = uf.modattr[attr] = value
function _set(uf, attr::MOI.AbstractVariableAttribute, vi::VI, value)
    if !haskey(uf.varattr, attr)
        uf.varattr[attr] = Dict{VI, Any}()
    end
    uf.varattr[attr][vi] = value
end
function _set(uf, attr::MOI.AbstractConstraintAttribute, ci::CI, value)
    if !haskey(uf.conattr, attr)
        uf.conattr[attr] = Dict{CI, Any}()
    end
    uf.conattr[attr][ci] = value
end
MOI.supports(::UniversalFallback, ::Union{MOI.AbstractModelAttribute, MOI.AbstractOptimizerAttribute}) = true
function MOI.set(uf::UniversalFallback, attr::Union{MOI.AbstractOptimizerAttribute, MOI.AbstractModelAttribute}, value)
    if MOI.supports(uf.model, attr)
        return MOI.set(uf.model, attr, value)
    else
        return _set(uf, attr, value)
    end
end
MOI.supports(::UniversalFallback, ::Union{MOI.AbstractVariableAttribute, MOI.AbstractConstraintAttribute}, ::Type{<:MOI.Index}) = true
function MOI.set(uf::UniversalFallback, attr::MOI.AbstractVariableAttribute, idx::VI, value)
    if MOI.supports(uf.model, attr, typeof(idx))
        return MOI.set(uf.model, attr, idx, value)
    else
        return _set(uf, attr, idx, value)
    end
end
function MOI.set(uf::UniversalFallback, attr::MOI.AbstractConstraintAttribute, idx::CI{F, S}, value) where {F, S}
    if MOI.supports_constraint(uf.model, F, S) && MOI.supports(uf.model, attr, CI{F, S})
        return MOI.set(uf.model, attr, idx, value)
    else
        return _set(uf, attr, idx, value)
    end
end

# Constraints
MOI.supports_constraint(uf::UniversalFallback, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet} = true
function MOI.add_constraint(uf::UniversalFallback, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    F = typeof(f)
    S = typeof(s)
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.add_constraint(uf.model, f, s)
    else
        constraints = get!(uf.constraints, (F, S)) do
            Dict{CI{F, S}, Tuple{F, S}}()
        end::Dict{CI{F, S}, Tuple{F, S}}
        uf.nextconstraintid += 1
        ci = CI{F, S}(uf.nextconstraintid)
        constraints[ci] = (f, s)
        return ci
    end
end
function MOI.modify(uf::UniversalFallback, ci::CI{F, S}, change::MOI.AbstractFunctionModification) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.modify(uf.model, ci, change)
    else
        (f, s) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (modifyfunction(f, change), s)
    end
end

function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, attr, ci)
    else
        uf.constraints[(F, S)][ci][1]
    end
end
function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, attr, ci)
    else
        uf.constraints[(F, S)][ci][2]
    end
end
function MOI.set(uf::UniversalFallback, ::MOI.ConstraintFunction, ci::CI{F,S}, func::F) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, MOI.ConstraintFunction(), ci, func)
    else
        (_, s) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (func, s)
    end
end
function MOI.set(uf::UniversalFallback, ::MOI.ConstraintSet, ci::CI{F,S}, set::S) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, MOI.ConstraintSet(), ci, set)
    else
        (f, _) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (f, set)
    end
end

# Objective
MOI.modify(uf::UniversalFallback, obj::MOI.ObjectiveFunction, change::MOI.AbstractFunctionModification) = MOI.modify(uf.model, obj, change)

# Variables
MOI.add_variable(uf::UniversalFallback) = MOI.add_variable(uf.model)
MOI.add_variables(uf::UniversalFallback, n) = MOI.add_variables(uf.model, n)


text:
struct ZygoteAD <: ADBackend end
ADBackend(::Val{:zygote}) = ZygoteAD
function setadbackend(::Val{:zygote})
    ADBACKEND[] = :zygote
end

function gradient_logp(
    backend::ZygoteAD,
    θ::AbstractVector{<:Real},
    vi::VarInfo,
    model::Model,
    sampler::AbstractSampler = SampleFromPrior(),
)
    T = typeof(getlogp(vi))

    # Specify objective function.
    function f(θ)
        new_vi = VarInfo(vi, sampler, θ)
        return getlogp(runmodel!(model, new_vi, sampler))
    end

    # Compute forward and reverse passes.
    l::T, ȳ = Zygote.pullback(f, θ)
    ∂l∂θ::typeof(θ) = ȳ(1)[1]

    return l, ∂l∂θ
end

Zygote.@nograd DynamicPPL.updategid!