text:
# PROBABILISTIC CIRCUITS LIBRARY ROOT

module ProbabilisticCircuits

# USE EXTERNAL MODULES

using Reexport

include("Utils/Utils.jl")

@reexport using .Utils

# INCLUDE CHILD MODULES
include("Probabilistic/Probabilistic.jl")
include("Logistic/Logistic.jl")
include("IO/IO.jl")
include("StructureLearner/StructureLearner.jl")
include("Reasoning/Reasoning.jl")

# USE CHILD MODULES (in order to re-export some functions)
@reexport using .Probabilistic
@reexport using .Logistic
@reexport using .IO
@reexport using .StructureLearner
@reexport using .Reasoning

end


text:
using MusicManipulations, Test, Statistics

@testset "humanize" begin

#test of humanize functionality
notes = randomnotes(100)
for n in notes; n.velocity = 100; end
hnotes = humanize(notes, :velocity, 10)
@test any(i -> notes[i].velocity ≠ hnotes[i].velocity, 1:length(notes))
σ = std(velocities(hnotes))
@test 9.5 ≤ σ ≤ 10.5

end

ext: jl


text:
@testset "Test Type Constructor" begin
    I = [1, 4, 3, 5]; J = [4, 7, 18, 9];
    lin = LinearSparse(18, 5, I, J)
    @test lin.ninputs == 18
    @test lin.noutputs == 5
    @test lin.nparams == 4
    @test lin.i == I
    @test lin.j == J
end

@testset "Test dx Buffer" begin
    I = [1, 4, 3, 5]; J = [4, 7, 18, 9];
    lin = LinearSparse(18, 5, I, J)
    dx = dx_buffer(lin)
    dx_expected = sparse([1, 4, 3, 5], [1, 4, 3, 5], zeros(4), 5, 18)
    @test typeof(dx) == typeof(dx_expected)
    @test dx.m  == dx_expected.m
    @test dx.n  == dx_expected.n
    @test dx.colptr == dx.colptr
    @test dx.rowval == dx.rowval
end

@testset "Test dΘ Buffer" begin
    I = [1, 4, 3, 5]; J = [4, 7, 18, 9];
    lin = LinearSparse(18, 5, I, J)
    dΘ = dΘ_buffer(lin)
    dΘ_expected = sparse([1, 4, 3, 5], [1, 2, 3, 4], zeros(4),
                         5, 4)
    @test typeof(dΘ) == typeof(dΘ_expected)
    @test dΘ.m  == dΘ_expected.m
    @test dΘ.n  == dΘ_expected.n
    @test dΘ.colptr == dΘ.colptr
    @test dΘ.rowval == dΘ.rowval
end

@testset "Test LinearSparse Evaluation" begin
    I = [1, 2, 3, 4]; J = [4, 3, 2, 1]; K = [1, 1, 1, 1]
    lin = LinearSparse(4, 4, I, J)
    z = z_buffer(lin)
    evaluate!(lin, [2, 4, 6, 5], K, z)
    @test z == [5, 6, 4, 2]
    evaluate!(lin, [1, 1, 1, 5], K, z)
    @test z == [5, 1, 1, 1]
    evaluate!(lin, [3, 3, 3, 5], K, z)
    @test z == [5, 3, 3, 3]
end

@testset "Test Jacobian Matrix Computation" begin
    @testset "Test Sparse" begin
        I = [1, 2, 3, 4]; J = [4, 3, 2, 1];
        lin = LinearSparse(4, 4, I, J)
        x = [10, 20, 30, 40]
        Θ = [1, 2, 3, 4]

        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)
        evaluate!(lin, x, Θ, z, dx=dx, dΘ=dΘ)
        @test dx  == [0 0 0 1;
                      0 0 2 0;
                      0 3 0 0;
                      4 0 0 0]
        @test dΘ == [40 0 0 0;
                     0 30 0 0;
                     0 0 20 0;
                     0 0 0 10]
    end

    @testset "Test Finite Diferences" begin
        I = [1, 2, 3, 4, 1, 2, 4]; J = [4, 3, 2, 1, 3, 2, 2];
        lin = LinearSparse(4, 4, I, J)
        x = [10.0, 20.0, 30.0, 40]
        Θ = [1, 2, 3, 4, 4.5, 5.5, 6.1]

        function cal_x(y)
            z = z_buffer(lin)
            evaluate!(lin, y, Θ, z)
            return z
        end
        function cal_Θ(theta)
            z = z_buffer(lin)
            evaluate!(lin, x, theta, z)
            return z
        end

        dx_numeric = Calculus.finite_difference_jacobian(cal_x, x)
        dΘ_numeric = Calculus.finite_difference_jacobian(cal_Θ, Θ)

        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)
        evaluate!(lin, x, Θ, z, dx=dx, dΘ=dΘ)
        @test dx ≈ dx_numeric
        @test dΘ ≈ dΘ_numeric
    end
    @testset "Test Argument Check" begin
        I = [1, 2, 3, 4, 1, 2, 4]; J = [4, 3, 2, 1, 3, 2, 2];
        lin = LinearSparse(4, 4, I, J)
        x = [10.0, 20.0, 30.0, 40]
        Θ = [1, 2, 3, 4, 4.5, 5.5, 6.1]

        z = z_buffer(lin)
        dx = dx_buffer(lin)
        dΘ = dΘ_buffer(lin)

        @test_throws ArgumentError evaluate!(lin, [1, 2, 3],
                                             Θ,
                                             z,
                                             dx=dx,
                                             dΘ=dΘ)
        @test_throws ArgumentError evaluate!(lin, x,
                                             [1, 2, 3, 4, 5],
                                             z,
                                             dx=dx,
                                             dΘ=dΘ)
        @test_throws ArgumentError evaluate!(lin, x, Θ,
                                             Vector{Float64}(10),
                                             dx=dx,
                                             dΘ=dΘ)
        @test_throws ArgumentError evaluate!(lin, x, Θ, z,
                                             dx=Matrix{Float64}(3, 4),
                                             dΘ=dΘ)
        @test_throws ArgumentError evaluate!(lin, x, Θ, z,
                                             dx=Matrix{Float64}(4, 3),
                                             dΘ=dΘ)
        @test_throws ArgumentError evaluate!(lin, x, Θ, z,
                                             dx=dx,
                                             dΘ=Matrix{Float64}(5, 7))
        @test_throws ArgumentError evaluate!(lin, x, Θ, z,
                                             dx=dx,
                                             dΘ=Matrix{Float64}(4, 2))
    end
end


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput05.txt")
function findRowColumn(string; first = 0, last = 127, upStr = 'B', lowStr = 'F')
    planeRows = first:last
    cStr = string[1]
    for s in string
        if s == upStr
            first = planeRows.start + Int(length(planeRows)/2)
            last = planeRows.stop
            planeRows = first:last
        elseif s == lowStr
            first = planeRows.start
            last = first + Int(length(planeRows)/2-1)
            planeRows = first:last
        end
        cStr = s
    end
    cStr == lowStr ? planeRows.start : planeRows.stop
end
function maxIdBoardingPass(data)
    idBoardingPass = 0
    for s in data
        posRow = findRowColumn(s)
        posCol = findRowColumn(s; first=0, last = 7, upStr = 'R', lowStr = 'L')
        idBoadPass = posRow *8 + posCol
        idBoardingPass =  idBoadPass > idBoardingPass ? idBoadPass : idBoardingPass
    end
    idBoardingPass
end
function IdsBoardingPass(data)
    idsBPass = []
    for s in data
        posRow = findRowColumn(s)
        posCol = findRowColumn(s; first=0, last = 7, upStr = 'R', lowStr = 'L')
        idBoadPass = posRow *8 + posCol
        push!(idsBPass, idBoadPass)
    end
    idsBPass
end
function findSeat(IdsBP)
    for s_id in 1:maximum(IdsBP)
        ismyseat = s_id ∉ IdsBP
        isplus1  = s_id + 1 ∈ IdsBP
        isminus1 = s_id - 1 ∈ IdsBP
        if ismyseat && isplus1 && isminus1
            return s_id
        end
    end
end
maximum(IdsBoardingPass(data))
findSeat(IdsBoardingPass(data))


text:
# Extract the matrix diagonal
Base.diag(dm::DataMatrix{T}) where {T} = DataArray(diag(dm.data), diag(dm.na))


text:
function genlib(cpucompiler::String, gpucompiler::String)
# Examples: genlib("g++","");
#           genlib("g++","nvcc");

if length(cpucompiler)==0
    error("cpucompiler is empty");
end

rm("cpuCore.a");
rm("libcpuCore.so");
str = `$cpucompiler -fPIC -O3 -c cpuCore.cpp`;
run(str);
str = `$cpucompiler --shared cpuCore.o -o libcpuCore.so`;
run(str);
str = `ar rvs cpuCore.a cpuCore.o`;
run(str);
rm("cpuCore.o");

if length(gpucompiler)>0
    rm("gpuCore.a");
    rm("libgpuCore.so");
    str = `$gpucompiler -D_FORCE_INLINES -O3 -c --compiler-options '-fPIC' -w gpuCore.cu`;
    run(str);
    str = `$cpucompiler --shared gpuCore.o -o libgpuCore.so`;
    run(str);
    str = `ar rvs gpuCore.a gpuCore.o`;
    run(str);
    rm("gpuCore.o");
end

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#######################################################################
# Green function
#######################################################################
abstract type AbstractGreensSolver end

struct GreensFunction{S<:AbstractGreensSolver,L,B<:NTuple{L,Union{Int,Missing}},H<:Hamiltonian}
    solver::S
    h::H
    boundaries::B
end

"""
    greens(h::Hamiltonian, solveobject; boundaries::NTuple{L,Integer} = missing)

Construct the Green's function `g::GreensFunction` of `L`-dimensional Hamiltonian `h` using
the provided `solveobject`. Currently valid `solveobject`s are

- `Schur1D()` (single-shot generalized eigenvalue approach for 1D Hamiltonians)

If a `boundaries = (n₁, n₂, ...)` is provided, a reflecting boundary is assumed for each
non-missing `nᵢ` perpendicular to Bravais vector `i` at a cell distance `nᵢ` from the
origin.

    h |> greens(h -> solveobject(h), args...)

Curried form equivalent to the above, giving `greens(h, solveobject(h), args...)`.

    g(ω, cells::Pair)
    g(ω)[cells::Pair]

From a constructed `g::GreensFunction`, obtain the retarded Green's function matrix at
frequency `ω` between unit cells `src` and `dst`, where `src, dst` are `::NTuple{L,Int}` or
`SVector{L,Int}`. If allowed by the used `solveobject`, `g0=g(ω)` builds an solution object
that can efficiently produce the Greens function between different cells at fixed `ω` with
`g0[cells]` without repeating cell-independent parts of the computation.

# Examples

```jldoctest
julia> g = LatticePresets.square() |> hamiltonian(hopping(-1)) |> unitcell((1,0), region = r->0<r[2]<3) |> greens(Schur1D())
GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method
  Flat matrix size      : 2 × 2
  Flat deflated size    : 2 × 2
  Original element type : scalar (ComplexF64)
  Boundaries            : (missing,)

julia> g(0.2, 3=>2) ≈ g(0.2)[3=>2]
true
```

# See also
    `Schur1D`
"""
greens(h::Hamiltonian{<:Any,L}, solverobject; boundaries = filltuple(missing, Val(L))) where {L} =
    GreensFunction(greensolver(solverobject, h), h, boundaries)
greens(s; kw...) = h -> greens(h, greensolver(s, h); kw...)

# fallback
greensolver(s::AbstractGreensSolver, h) = s
greensolver(s::Function, h) = s(h)

sanitize_cells((cell0, cell1)::Pair{<:Number,<:Number}) =
    SA[Int(cell0)] => SA[Int(cell1)]
sanitize_cells((cell0, cell1)::Pair{<:NTuple{L,Number},<:NTuple{L,Number}}) where {L} =
    SVector(Int.(cell0)) => SVector(Int.(cell1))
sanitize_cells(cells) =
    throw(ArgumentError("Cells should be of the form `cᵢ => cⱼ`, with each `c` an `Integer` or an `NTuple{L,Integer}`, got $cells"))

Base.size(g::GreensFunction, args...) = size(g.h, args...)
Base.eltype(g::GreensFunction) = eltype(g.h)

flatsize(g::GreensFunction, args...) = flatsize(g.h, args...)
blockeltype(g::GreensFunction) = blockeltype(g.h)
blocktype(g::GreensFunction) = blocktype(g.h)
orbitaltype(g::GreensFunction) = orbitaltype(g.h)
orbitalstructure(g::GreensFunction) = orbitalstructure(g.h)

#######################################################################
# Schur1DGreensSolver
#######################################################################
"""
    Schur1D()

Return a Greens function solver using the generalized eigenvalue approach, whereby given the
energy `ω`, the eigenmodes of the infinite 1D Hamiltonian, and the corresponding infinite
and semi-infinite Greens function can be computed by solving the generalized eigenvalue
equation

    A⋅φχ = λ B⋅φχ
    A = [0 I; -h₊ ω-h₀]
    B = [I 0; 0 h₋]

This is the matrix form of the problem `λ(ω-h₀)φ - h₊φ - λ²h₋φ = 0`, where `φχ = [φ; λφ]`,
and `φ` are `ω`-energy eigenmodes, with (possibly complex) momentum `q`, and eigenvalues are
`λ = exp(-iqa₀)`. The algorithm assumes the Hamiltonian has only `dn = (0,)` and `dn = (±1,
)` Bloch harmonics (`h₀`, `h₊` and `h₋`), and will error otherwise instructing the user to
grow the unit cell. Bound states in the spectrum will yield delta functions in the density
of states that can be resolved by adding a broadening in the form of a small positive
imaginary part to `ω`. If `ω::Real`, a small imaginary part will be added automatically.

For performace, the eigenvalue equation may be `deflated' and `stabilized', i.e. singular
solutions `λ=0,∞` will be removed, and an inverse-free algorithm is used to preserve
precision even in the presence of singularities.

# Examples
```jldoctest
julia> using LinearAlgebra

julia> h = LP.honeycomb() |> hamiltonian(hopping(1)) |> unitcell((1,-1), (10,10)) |> wrap(2);

julia> g = greens(h, Schur1D(), boundaries = (0,))
GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method
  Flat matrix size      : 40 × 40
  Flat deflated size    : 20 × 20
  Original element type : scalar (ComplexF64)
  Boundaries            : (0,)

julia> tr(g(0.3, 1=>1))
-32.193416071797216 - 3.4400038418349084im
```

# See also
    `greens`
"""
struct Schur1D end

greensolver(s::Schur1D, h) = Schur1DGreensSolver(s, h)

#### Schur1DGreensSolver ###################################################################

struct Schur1DWorkspace{T}
    nn1::Matrix{T}
    nn2::Matrix{T}
    nr1::Matrix{T}
    nr2::Matrix{T}
    rr1::Matrix{T}
    rr2::Matrix{T}
    rr3::Matrix{T}
    rr4::Matrix{T}
    r2r21::Matrix{T}
    r2r22::Matrix{T}
    n2r1::Matrix{T}
    n2r2::Matrix{T}
end

Schur1DWorkspace(R::AbstractMatrix{T}) where {T} = Schur1DWorkspace{T}(size(R)...)

Schur1DWorkspace{T}(n, r) where {T} = Schur1DWorkspace(Matrix{T}.(undef,
    ((n, n), (n, n), (n, r), (n, r), (r, r), (r, r), (r, r), (r, r), (2r, 2r), (2r, 2r), (n+2r, n+2r), (n+2r, n+2r)))...)

struct Deflator{T,S}
    L::Matrix{T}        # h₊ = L*R'
    R::Matrix{T}        # h₋ = R*L'
    iG::Matrix{T}       # Matrix(-h₀ + (LL' + RR') * im)
    ωshifter::S         # metadata to aid in ω-shifting iG
end

struct Schur1DGreensSolver{D<:Deflator,T,M} <: AbstractGreensSolver
    h0::M
    hp::M
    hm::M
    effmat::EffectiveMatrix{T}
    deflator::D
    tmp::Schur1DWorkspace{T}
end

function Schur1DGreensSolver(s::Schur1D, h)
    latdim(h) == 1 || throw(ArgumentError("Cannot use a Schur1D Green function solver with an $(latdim(h))-dimensional Hamiltonian"))
    H = maybeflatten(h)
    maxdn = maximum(har -> abs(first(har.dn)), h.harmonics)
    maxdn > 1 && throw(ArgumentError("The Hamiltonian has next-nearest unitcell hoppings. Please enlarge the unit cell with `unitcell(h, $maxdn)` to reduce to nearest-cell couplings."))
    h₋, h₀, h₊ = H[(-1,)], H[(0,)], H[(1,)]
    deflator = Deflator(h₊, h₀, h₋)
    L, R = deflator.L, deflator.R
    effmat = EffectiveMatrix(h₀, L, R)
    tmp = Schur1DWorkspace(R)
    return Schur1DGreensSolver(h₀, h₊, h₋, effmat, deflator, tmp)
end

function Deflator(h₊, h₀, h₋)
    h₊ ≈ h₋' || throw(ArgumentError("Deflation requires mutually adjoint intercell h₊ = h₋'. If you intended to build a non-Hermitian Hamiltonian, please use the undeflated method `Schur(deflation = nothing)`."))
    Ls, Rs   = svd_sparse(h₊)
    L, R     = Matrix(Ls), Matrix(Rs)
    G⁻¹      = Matrix(-h₀ + (Ls * Ls' + Rs * Rs') * im)
    ωshifter = diag(G⁻¹)
    return Deflator(L, R, G⁻¹, ωshifter)
end

function Schur1DWorkspace(d::Deflator, h₀::AbstractMatrix{T}) where {T}
    n = size(h₀, 1)
    r = size(d.R, 2)
    return Schur1DWorkspace{T}(n, r)
end

function Base.show(io::IO, g::GreensFunction{<:Schur1DGreensSolver})
    print(io, summary(g), "\n",
"  Flat matrix size      : $(size(g.solver.h0, 1)) × $(size(g.solver.h0, 2))
  Flat deflated size    : $(deflated_size_text(g))
  Original element type : $(displayelements(g.h))
  Boundaries            : $(g.boundaries)")
end

undeflated_size_text(g::GreensFunction) = undeflated_size_text(g.solver.deflator.R)
undeflated_size_text(R::AbstractMatrix) = "$(size(R, 1)) × $(size(R, 1))"

deflated_size_text(g::GreensFunction) = deflated_size_text(g.solver.deflator.R)
deflated_size_text(R::AbstractMatrix) = "$(size(R, 2)) × $(size(R, 2))"

Base.summary(g::GreensFunction{<:Schur1DGreensSolver}) =
    "GreensFunction{Schur1DGreensSolver}: Green's function using the Schur1D method"

Base.size(s::Schur1DGreensSolver, args...) = size(s.deflator, args...)

#### Deflation and modes ###################################################################

# Return deflated Adef and Bdef
function deflated_pencil!(tmp::Schur1DWorkspace, d::Deflator{T}) where {T}
    luiG = lu!(copy!(tmp.nn1, d.iG))
    GR = ldiv!(luiG, copy!(tmp.nr1, d.R))
    GL = ldiv!(luiG, copy!(tmp.nr2, d.L))
    R = d.R
    L = d.L
    r = size(R, 2)
    i1, i2 = 1:r, r+1:2r
    Adef, Bdef = tmp.r2r21, tmp.r2r22
    fill!(Adef, zero(T))
    fill!(Bdef, zero(T))
    one!(view(Adef, i1, i1))
    mul!(view(Adef, i1, i1), L', GL, -im, 1)
    mul!(view(Adef, i1, i2), L', GL, -1, 1)
    mul!(view(Adef, i2, i1), R', GL, im, 1)
    mul!(view(Adef, i2, i2), R', GL, 1, 1)
    one!(view(Bdef, i2, i2))
    mul!(view(Bdef, i1, i1), L', GR, 1, 1)
    mul!(view(Bdef, i1, i2), L', GR, im, 1)
    mul!(view(Bdef, i2, i1), R', GR, -1, 1)
    mul!(view(Bdef, i2, i2), R', GR, -im, 1)
    return Adef, Bdef
end

function shiftω!(d::Deflator, ω)
    diagG⁻¹ = d.ωshifter
    for (n, v) in enumerate(diagG⁻¹)
        d.iG[n, n] = ω + v
    end
    return d
end

function mode_subspaces(s::Schur1DGreensSolver{<:Deflator}, ω)
    d = s.deflator
    shiftω!(d, ω)
    A, B = deflated_pencil!(s.tmp, d)
    return mode_subspaces(s, A, B, imag(ω))
end

# returns invariant subspaces of retarded and advanced eigenmodes
function mode_subspaces(s::Schur1DGreensSolver, A::AbstractArray{T}, B::AbstractArray{T}, imω) where {T}
    ZrL, ZrR, ZaL, ZaR = s.tmp.rr1, s.tmp.rr2, s.tmp.rr3, s.tmp.rr4
    r = size(A, 1) ÷ 2
    if !iszero(r)
        sch = schur!(A, B)
        # Retarded modes
        whichmodes = Vector{Bool}(undef, length(sch.α))
        retarded_modes!(whichmodes, sch, imω)
        ordschur!(sch, whichmodes)
        copy!(ZrL, view(sch.Z, 1:r, 1:sum(whichmodes)))
        copy!(ZrR, view(sch.Z, r+1:2r, 1:sum(whichmodes)))
        # Advanced modes
        advanced_modes!(whichmodes, sch, imω)
        ordschur!(sch, whichmodes)
        copy!(ZaL, view(sch.Z, 1:r, 1:sum(whichmodes)))
        copy!(ZaR, view(sch.Z, r+1:2r, 1:sum(whichmodes)))
    end
    return ZrL, ZrR, ZaL, ZaR
end

# need this barrier for type-stability (sch.α and sch.β are finicky)
function retarded_modes!(whichmodes, sch, imω)
    whichmodes .= abs.(sch.α) .< abs.(sch.β)
    sum(whichmodes) == length(whichmodes) ÷ 2 || throw_imω(imω)
    return whichmodes
end

function advanced_modes!(whichmodes, sch, imω)
    whichmodes .= abs.(sch.β) .< abs.(sch.α)
    sum(whichmodes) == length(whichmodes) ÷ 2 || throw_imω(imω)
    return whichmodes
end

throw_imω(imω) =
    throw(ArgumentError("Couldn't separate advanced from retarded modes. Consider adding a larger positive imaginary part to ω, currently $imω"))

#### g(ω) ##################################################################################

(g::GreensFunction{<:Schur1DGreensSolver})(ω; kw...) =
    Schur1DGreensSolution(g, ensurecomplex(ω); kw...)

ensurecomplex(ω::T) where {T<:Real} = ω + im * default_tol(T)
ensurecomplex(ω::Complex) = ω

#### Schur1DGreensSolution #################################################################

struct Schur1DGreensSolutionWorkspace{T}
    rs1::Matrix{T}
    rs2::Matrix{T}
    ns1::Matrix{T}
    ns2::Matrix{T}
    dr::Matrix{T}
    rr::Matrix{T}
end

Schur1DGreensSolutionWorkspace(srcmat, dstmat, R::AbstractMatrix{T}) where {T} =
    Schur1DGreensSolutionWorkspace{T}(size(srcmat, 2), size(dstmat, 2), size(R)...)

Schur1DGreensSolutionWorkspace{T}(s, d, n, r) where {T} =
    Schur1DGreensSolutionWorkspace(Matrix{T}.(undef, ((r, s), (r, s), (n, s), (n, s), (d, r), (r, r)))...)

struct Schur1DGreensSolution{B<:Union{Int,Missing},T,S<:SubArray{T},O<:OrbitalStructure}
    ω::T
    L::Matrix{T}
    R::Matrix{T}
    dstmat::Matrix{T}
    G∞S::S      # G^∞ * source
    GLR::S      # G^L * R
    GRL::S      # G^R * L
    boundary::B
    orbstruct::O
    tmp::Schur1DGreensSolutionWorkspace{T}
end

function Base.show(io::IO, s::Schur1DGreensSolution{B}) where {B}
    print(io, summary(s), "\n",
"  ω             : $(s.ω)
  Matrix size   : $(undeflated_size_text(s.R))
  Deflated size : $(deflated_size_text(s.R))
  Element type  : $(blocktype(s.orbstruct))
  Boundary      : $(s.boundary)")
end

Base.summary(s::Schur1DGreensSolution) = "Schur1DGreensSolution : Schur1D solution of a GreensFunction at fixed ω"

function Schur1DGreensSolution(g, ω; source = I, dest = I)
    s = g.solver
    L, R, em, G⁻¹, G⁻¹backup = s.deflator.L, s.deflator.R, s.effmat, s.tmp.n2r1, s.tmp.n2r2
    modes = mode_subspaces(s, ω)
    effective_matrix!(G⁻¹, em, ω, modes)
    copy!(G⁻¹backup, G⁻¹)
    srcmat = flat_source_dest_matrix(source, g)
    G∞S = padded_ldiv(G⁻¹, srcmat, em, Val(:LR))
    copy!(G⁻¹, G⁻¹backup)
    GRL = padded_ldiv(G⁻¹, L, em, Val(:R))
    copy!(G⁻¹, G⁻¹backup)
    GLR = padded_ldiv(G⁻¹, R, em, Val(:L))
    dstmat = flat_source_dest_matrix(dest, g)
    boundary = only(g.boundaries)
    tmp = Schur1DGreensSolutionWorkspace(srcmat, dstmat, R)
    orbstruct = orbitalstructure(g.h)

    # Σ = R * L' * GRL * R'
    # @show sum(abs.(Σ - s.hm * ((ω*I - s.h0 - Σ) \ Matrix(s.hp))))
    # Σ = L * R' * GLR * L'
    # @show sum(abs.(Σ - s.hp * ((ω*I - s.h0 - Σ) \ Matrix(s.hm))))

    return Schur1DGreensSolution(ω, L, R, dstmat, G∞S, GLR, GRL, boundary, orbstruct, tmp)
end

function flat_source_dest_matrix(id::UniformScaling, g)
    n = flatsize(g, 1)
    return Matrix(one(blockeltype(g)) * id, n, n)
end

flat_source_dest_matrix(sourcemat::AbstractMatrix, g) =
    flatten(sourcemat, orbitalstructure(g))

function flat_source_dest_matrix(sources, g)
    n = length(sources)
    m = size(g, 1)
    T = blocktype(g)
    sourcemat = zeros(T, m, n)
    for (i, s) in enumerate(sources)
        sourcemat[s, i] = one(T)
    end
    return flat_source_dest_matrix(sourcemat, g)
end

flat_source_dest_dim(s::AbstractArray, g) = size(s, 1) * blockdim(g)
flat_source_dest_dim(::UniformScaling, g) = flatsize(g, 1)

## Generic path (indexing g(ω; kw...)[cells])

# Infinite:
# G∞ₙₙ = G∞₀₀ = (ω*I - h0 - ΣR - ΣL)⁻¹
# G∞ₙₘ = (G₁₁h₊)ⁿ⁻ᵐ G∞₀₀ = G₁₁L (R'G₁₁L)ⁿ⁻ᵐ⁻¹ R'G∞₀₀     for n-m > 0
# G∞ₙₘ = (G₋₁₋₁h₋)ᵐ⁻ⁿ G∞₀₀ = G₋₁₋₁R(L'G₋₁₋₁R)ᵐ⁻ⁿ⁻¹L'G∞₀₀ for n-m < 0
function Base.getindex(s::Schur1DGreensSolution{Missing,T}, cells) where {T}
    src, dst = sanitize_cells(cells)
    dist = only(dst) - only(src)
    G = zeros(T, size(s.dstmat, 2), size(s.G∞S, 2))
    add_G∞!(G, s, dist)
    G´ = maybe_unflatten_blocks(G, s.orbstruct)
    return G´
end

# Semiinifinite:
# Gₙₘ = (Ghⁿ⁻ᵐ - GhⁿGh⁻ᵐ)G∞₀₀ = G∞ₙₘ - GhⁿG∞₀ₘ
# Gₙₘ = G∞ₙₘ - G₁₁L(R'G₁₁L)ⁿ⁻¹ R'G∞₀ₘ       for n > 1
# Gₙₘ = G∞ₙₘ - G₋₁₋₁R(L'G₋₁₋₁R)⁻ⁿ⁻¹L'G∞₀ₘ   for n < -1
function Base.getindex(s::Schur1DGreensSolution{Int,T}, cells) where {T}
    cells´ = sanitize_cells(cells)
    G = zeros(T, size(s.dstmat, 2), size(s.G∞S, 2))
    if !is_across_boundary(cells´, s.boundary)
        m, n = dist_to_boundary.(cells´, s.boundary)  # m, n = dist_src, dist_dest
        add_G∞!(G, s, n-m, 0, 1)
        add_G∞!(G, s, n, -m, -1)
    end
    G´ = maybe_unflatten_blocks(G, s.orbstruct)
    return G´
end

# G += α dstmat' * Gh^dist * Gh^dist´ * G∞S
function add_G∞!(G, s::Schur1DGreensSolution, dist, dist´ = 0, α = 1)
    TG∞ = Ghⁿmul!(s.tmp.ns1, s.G∞S, s, dist´)
    TG∞ = Ghⁿmul!(s.tmp.ns2, TG∞, s, dist)
    mul!(G, s.dstmat', TG∞, α, 1)
    return G
end

# Gs´ = Ghⁿ * Gs = GX * (X´GX)ⁿ⁻¹ * X´ * Gs
function Ghⁿmul!(Gs´, Gs, s::Schur1DGreensSolution, n)
    if n == 0
        Gs´ = Gs  # no need to copy!(Gs´, Gs)
    else
        isplus = n > 0
        GX  = ifelse(isplus, s.GRL,  s.GLR)
        X´ = ifelse(isplus, s.R', s.L')
        X´GX = mul!(s.tmp.rr, X´, GX)
        X´Gs = mul!(s.tmp.rs1, X´, Gs)
        mul!(Gs´, GX, mul!(s.tmp.rs2, (X´GX)^(abs(n)-1), X´Gs))
    end
    return Gs´
end

## Fast-paths

function (g::GreensFunction)(ω, cells; kw...)
    cells´ = sanitize_cells(cells)
    ω´ = ensurecomplex(ω)
    b = only(g.boundaries)
    if is_infinite_local(cells´, b)
        G = infinite_local_fastpath(g, ω´; kw...)
    elseif is_across_boundary(cells´, b)
        G = zero_fastpath(g; kw...)
    elseif is_from_surface(cells´, b)
        G = from_surface_fastpath(g, ω´, dist_to_boundary(cells´, b); kw...)
    else # general form
        G = g(ω´; kw...)[cells]
    end
    G´ = maybe_unflatten_blocks(G, orbitalstructure(g.h))  # only if non-scalar eltype
    return G´
end

is_infinite_local((src, dst), b::Missing) = only(src) == only(dst)
is_infinite_local(cells, b) = false

is_from_surface((src, dst), b::Int) = abs(dist_to_boundary(src, b)) == 1
is_from_surface(cells, b) = false

is_across_boundary((src, dst), b::Int) =
    sign(dist_to_boundary(src, b)) != sign(dist_to_boundary(dst, b)) ||
    dist_to_boundary(src, b) == 0 || dist_to_boundary(dst, b) == 0
is_across_boundary(cells, b) = false

dist_to_boundary(cell, b) = only(cell) - b
dist_to_boundary((src, dst)::Pair, b) = dist_to_boundary(src, b), dist_to_boundary(dst, b)

# Infinite:
# G∞ₙₙ = G∞₀₀ = (ω*I - h0 - ΣR - ΣL)⁻¹
function infinite_local_fastpath(g, ω; source = I, dest = I)
    em, G⁻¹ = g.solver.effmat, g.solver.tmp.n2r1
    modes = mode_subspaces(g.solver, ω)
    effective_matrix!(G⁻¹, em, ω, modes)
    srcmat = flat_source_dest_matrix(source, g)
    dstmat = flat_source_dest_matrix(dest, g)
    G = dstmat' * padded_ldiv(G⁻¹, srcmat, em, Val(:LR))
    return G
end

# Surface-bulk semi-infinite:
# G₁₁ = (ω*I - h0 - ΣR)⁻¹, G₋₁₋₁ = (ω*I - h0 - ΣL)⁻¹
# Gₙ₀ = G₀ₙ = 0
# h₊ = LR', h₋ = RL'
# Gₙ₁ = (G₁₁h₊)ⁿ⁻¹G₁₁       = G₁₁L(R'G₁₁L)ⁿ⁻²R'G₁₁         for n ≥ 2
# Gₙ₁ = (G₋₁₋₁h₋)⁻ⁿ⁻¹G₋₁₋₁  = G₋₁₋₁R(L'G₋₁₋₁R)⁻ⁿ⁻²L'G₋₁₋₁  for n ≤ -2
function from_surface_fastpath(g, ω, (dsrc, ddst); source = I, dest = I)
    dist = ddst - dsrc
    em, G⁻¹ = g.solver.effmat, g.solver.tmp.n2r1
    modes = mode_subspaces(g.solver, ω)
    effective_matrix!(G⁻¹, em, ω, modes)
    isplus = dsrc > 0
    side = ifelse(isplus, Val(:R), Val(:L))
    luG⁻¹ = padded_lu!(G⁻¹, em, side)
    srcmat = flat_source_dest_matrix(source, g)
    dstmat = flat_source_dest_matrix(dest, g)

    if dist == 0
        G = dstmat' * padded_ldiv(luG⁻¹, srcmat, em, side)
    else
        X  = ifelse(isplus, em.L, em.R)
        X´ = ifelse(isplus, em.R', em.L')
        X´Gs = X´ * padded_ldiv(luG⁻¹, srcmat, em, side)
        GX = padded_ldiv(luG⁻¹, X, em, side)
        dGX = dstmat' * GX
        if abs(dist) > 1
            X´GX = mul!(g.solver.tmp.rr1, X´, GX)
            G = dGX * (X´GX)^(abs(dist)-1) * X´Gs
        else
            G = dGX*X´Gs
        end
    end
    return G
end

zero_fastpath(g; source = I, dest = I) =
    zeros(blockeltype(g), flat_source_dest_dim.((source, dest), Ref(g))...)

#######################################################################
# BandGreensSolver
#######################################################################
# struct SimplexData{D,E,T,C<:SMatrix,DD,SA<:SubArray}
#     ε0::T
#     εmin::T
#     εmax::T
#     k0::SVector{D,T}
#     Δks::SMatrix{D,D,T,DD}     # k - k0 = Δks * z
#     volume::T
#     zvelocity::SVector{D,T}
#     edgecoeffs::NTuple{E,Tuple{T,C}} # s*det(Λ)/w.w and Λc for each edge
#     dωzs::NTuple{E,NTuple{2,SVector{D,T}}}
#     defaultdη::SVector{D,T}
#     φ0::SA
#     φs::NTuple{D,SA}
# end

# struct BandGreensSolver{P<:SimplexData,E,H<:Hamiltonian} <: AbstractGreensSolver
#     simplexdata::Vector{P}
#     indsedges::NTuple{E,Tuple{Int,Int}} # all distinct pairs of 1:V, where V=D+1=num verts
#     h::H
# end

# function Base.show(io::IO, g::GreensFunction{<:BandGreensSolver})
#     print(io, summary(g), "\n",
# "  Matrix size    : $(size(g.h, 1)) × $(size(g.h, 2))
#   Element type   : $(displayelements(g.h))
#   Band simplices : $(length(g.solver.simplexdata))")
# end

# Base.summary(g::GreensFunction{<:BandGreensSolver}) =
#     "GreensFunction{Bandstructure}: Green's function from a $(latdim(g.h))D bandstructure"

# function greensolver(b::Bandstructure{D}, h) where {D}
#     indsedges = tuplepairs(Val(D))
#     v = [SimplexData(simplex, band, indsedges) for band in bands(b) for simplex in band.simplices]
#     return BandGreensSolver(v,  indsedges, h)
# end

# edges_per_simplex(L) = binomial(L,2)

# function SimplexData(simplex::NTuple{V}, band, indsedges) where {V}
#     D = V - 1
#     vs = ntuple(i -> vertices(band)[simplex[i]], Val(V))
#     ks = ntuple(i -> SVector(Base.front(Tuple(vs[i]))), Val(V))
#     εs = ntuple(i -> last(vs[i]), Val(V))
#     εmin, εmax = extrema(εs)
#     ε0 = first(εs)
#     k0 = first(ks)
#     Δks = hcat(tuple_minus_first(ks)...)
#     zvelocity = SVector(tuple_minus_first(εs))
#     volume = abs(det(Δks))
#     edgecoeffs = edgecoeff.(indsedges, Ref(zvelocity))
#     dωzs = sectionpoint.(indsedges, Ref(zvelocity))
#     defaultdη = dummydη(zvelocity)
#     φ0 = vertexstate(first(simplex), band)
#     φs = vertexstate.(Base.tail(simplex), Ref(band))
#     return SimplexData(ε0, εmin, εmax, k0, Δks, volume, zvelocity, edgecoeffs, dωzs, defaultdη, φ0, φs)
# end

# function edgecoeff(indsedge, zvelocity::SVector{D}) where {D}
#     basis = edgebasis(indsedge, Val(D))
#     othervecs = Base.tail(basis)
#     edgevec = first(basis)
#     cutvecs = (v -> dot(zvelocity, edgevec) * v - dot(zvelocity, v) * edgevec).(othervecs)
#     Λc = hcat(cutvecs...)
#     Λ = hcat(zvelocity, Λc)
#     s = sign(det(hcat(basis...)))
#     coeff = s * (det(Λ)/dot(zvelocity, zvelocity))
#     return coeff, Λc
# end

# function edgebasis(indsedge, ::Val{D}) where {D}
#     inds = ntuple(identity, Val(D+1))
#     swappedinds = tupleswapfront(inds, indsedge) # places the two edge vertindices first
#     zverts = (i->unitvector(SVector{D,Int}, i-1)).(swappedinds)
#     basis = (z -> z - first(zverts)).(Base.tail(zverts)) # first of basis is edge vector
#     return basis
# end

# function sectionpoint((i, j), zvelocity::SVector{D,T}) where {D,T}
#     z0, z1 = unitvector(SVector{D,Int}, i-1), unitvector(SVector{D,Int}, j-1)
#     z10 = z1 - z0
#     # avoid numerical cancellation errors due to zvelocity perpendicular to edge
#     d = chop(dot(zvelocity, z10), maximum(abs.(zvelocity)))
#     dzdω = z10 / d
#     dz0 = z0 - z10 * dot(zvelocity, z0) / d
#     return dzdω, dz0   # The section z is dω * dzdω + dz0
# end

# # A vector, not parallel to zvelocity, and with all nonzero components and none equal
# function dummydη(zvelocity::SVector{D,T}) where {D,T}
#     (D == 1 || iszero(zvelocity)) && return SVector(ntuple(i->T(i), Val(D)))
#     rng = MersenneTwister(0)
#     while true
#         dη = rand(rng, SVector{D,T})
#         isparallel = dot(dη, zvelocity)^2 ≈ dot(zvelocity, zvelocity) * dot(dη, dη)
#         isvalid = allunique(dη) && !isparallel
#         isvalid && return dη
#     end
#     throw(error("Unexpected error finding dummy dη"))
# end

# function vertexstate(ind, band)
#     ϕind = 1 + band.dimstates*(ind - 1)
#     state = view(band.states, ϕind:(ϕind+band.dimstates-1))
#     return state
# end

# ## Call API

# function greens!(matrix, g::GreensFunction{<:BandGreensSolver,L}, ω::Number, (src, dst)::SVectorPair{L}) where {L}
#     fill!(matrix, zero(eltype(matrix)))
#     dn = dst - src
#     for simplexdata in g.solver.simplexdata
#         g0, gjs = green_simplex(ω, dn, simplexdata, g.solver.indsedges)
#         addsimplex!(matrix, g0, gjs, simplexdata)
#     end
#     return matrix
# end

# function green_simplex(ω, dn, data::SimplexData{L}, indsedges) where {L}
#     dη = data.Δks' * dn
#     phase = cis(dot(dn, data.k0))
#     dω = ω - data.ε0
#     gz = simplexterm.(dω, Ref(dη), Ref(data), data.edgecoeffs, data.dωzs, indsedges)
#     g0z, gjz = first.(gz), last.(gz)
#     g0 = im^(L-1) * phase * sum(g0z)
#     gj = -im^L * phase * sum(gjz)
#     return g0, gj
# end

# function simplexterm(dω, dη::SVector{D,T}, data, coeffs, (dzdω, dz0), (i, j)) where {D,T}
#     bailout = Complex(zero(T)), Complex.(zero(dη))
#     z = dω * dzdω + dz0
#     # Edges with divergent sections do not contribute
#     all(isfinite, z) || return bailout
#     z0 = unitvector(SVector{D,T},i-1)
#     z1 = unitvector(SVector{D,T},j-1)
#     coeff, Λc = coeffs
#     # If dη is zero (DOS) use a precomputed (non-problematic) simplex-constant vector
#     dη´ = iszero(dη) ? data.defaultdη : dη
#     d = dot(dη´, z)
#     d0 = dot(dη´, z0)
#     d1 = dot(dη´, z1)
#     # Skip if singularity in formula
#     (d ≈ d0 || d ≈ d1) && return bailout
#     s = sign(dot(dη´, dzdω))
#     coeff0 = coeff / prod(Λc' * dη´)
#     coeffj = isempty(Λc) ? zero(dη) : (Λc ./ ((dη´)' * Λc)) * sumvec(Λc)
#     params = s, d, d0, d1
#     zs = z, z0, z1
#     g0z = iszero(dη) ? g0z_asymptotic(D, coeff0, params) : g0z_general(coeff0, params)
#     gjz = iszero(dη) ? gjz_asymptotic(D, g0z, coeffj, coeff0, zs, params) :
#                        gjz_general(g0z, coeffj, coeff0, zs, params)
#     return g0z, gjz
# end

# sumvec(::SMatrix{N,M,T}) where {N,M,T} = SVector(ntuple(_->one(T),Val(M)))

# g0z_general(coeff0, (s, d, d0, d1)) =
#     coeff0 * cis(d) * ((cosint_c(-s*(d0-d)) + im*sinint(d0-d)) - (cosint_c(-s*(d1-d)) + im*sinint(d1-d)))

# gjz_general(g0z, coeffj, coeff0, (z, z0, z1), (s, d, d0, d1)) =
#     g0z * (im * z - coeffj) + coeff0 * ((z0-z) * cis(d0) / (d0-d) - (z1-z) * cis(d1) / (d1-d))

# g0z_asymptotic(D, coeff0, (s, d, d0, d1)) =
#     coeff0 * (cosint_a(-s*(d0-d)) - cosint_a(-s*(d1-d))) * (im*d)^(D-1)/factorial(D-1)

# function gjz_asymptotic(D, g0z, coeffj, coeff0, (z, z0, z1), (s, d, d0, d1))
#     g0z´ = g0z
#     for n in 1:(D-1)
#         g0z´ += coeff0 * im^n * (im*d)^(D-1-n)/factorial(D-1-n) *
#                 ((d0-d)^n - (d1-d)^n)/(n*factorial(n))
#     end
#     gjz = g0z´ * (im * z - im * coeffj * d / D) +
#         coeff0 * ((z0-z) * (im*d0)^D / (d0-d) - (z1-z) * (im*d1)^D / (d1-d)) / factorial(D)
#     return gjz
# end

# cosint_c(x::Real) = ifelse(iszero(abs(x)), zero(x), cosint(abs(x))) + im*pi*(x<=0)

# cosint_a(x::Real) = ifelse(iszero(abs(x)), zero(x), log(abs(x))) + im*pi*(x<=0)

# function addsimplex!(matrix, g0, gjs, simplexdata)
#     φ0 = simplexdata.φ0
#     φs = simplexdata.φs
#     vol = simplexdata.volume
#     for c in CartesianIndices(matrix)
#         (row, col) = Tuple(c)
#         x = g0 * (φ0[row] * φ0[col]')
#         for (φ, gj) in zip(φs, gjs)
#             x += (φ[row]*φ[col]' - φ0[row]*φ0[col]') * gj
#         end
#         matrix[row, col] += vol * x
#     end
#     return matrix
# end


text:
using GraphQLParser
using Documenter

DocMeta.setdocmeta!(GraphQLParser, :DocTestSetup, :(using GraphQLParser); recursive=true)

makedocs(;
    modules=[GraphQLParser],
    authors="Mal Miller and contributors",
    repo="https://github.com/mmiller-max/GraphQLParser.jl/blob/{commit}{path}#{line}",
    sitename="GraphQLParser.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://mmiller-max.github.io/GraphQLParser.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Library" => [
            "Public" => "public.md",
            "Private" => "private.md",
        ]
    ],
)

deploydocs(;
    repo="github.com/mmiller-max/GraphQLParser.jl",
    devbranch="main",
)


text:
using Test
using JET
using LinearAlgebra, IntervalSets

import ClimaCore:
    Geometry, Fields, Domains, Topologies, Meshes, Spaces, Operators

# We need to pull these broadcasted expressions out as
# toplevel functions due to how broadcast expressions are
# lowered so JETTest can have a single callsite to analyze.

function opt_Interpolate(I, field)
    return I.(field)
end

function opt_Restrict(R, Ifield)
    return R.(Ifield)
end

function opt_Gradient(field)
    grad = Operators.Gradient()
    return grad.(field)
end

function opt_WeakGradient(field)
    wgrad = Operators.WeakGradient()
    return wgrad.(field)
end

function opt_Curl(field)
    curl = Operators.Curl()
    return curl.(field)
end

function opt_WeakCurl(field)
    wcurl = Operators.WeakCurl()
    return wcurl.(field)
end

function opt_CurlCurl(field)
    curl = Operators.Curl()
    return curl.(curl.(field))
end

function opt_Divergence(field)
    div = Operators.Divergence()
    return div.(field)
end

function opt_WeakDivergence(field)
    wdiv = Operators.WeakDivergence()
    return wdiv.(field)
end

function opt_ScalarDSS(field)
    grad = opt_Gradient(field)
    Spaces.weighted_dss!(grad)
    return grad
end

function opt_VectorDss_Curl(field)
    return Spaces.weighted_dss!(opt_Curl(field))
end

function opt_VectorDss_DivGrad(field)
    sdiv = Operators.Divergence()
    wgrad = Operators.WeakGradient()
    return Spaces.weighted_dss!(@. wgrad(sdiv(field)))
end

function opt_ScalarHyperdiffusion(field)
    grad = Operators.Gradient()
    wdiv = Operators.WeakDivergence()
    χ = Spaces.weighted_dss!(@. wdiv(grad(field)))
    ∇⁴field = Spaces.weighted_dss!(@. wdiv(grad(χ)))
    return ∇⁴field
end

function opt_VectorHyperdiffusion(field)
    curl = Operators.Curl()
    wcurl = Operators.WeakCurl()

    sdiv = Operators.Divergence()
    wgrad = Operators.WeakGradient()

    χ = Spaces.weighted_dss!(
        @. Geometry.UVVector(wgrad(sdiv(field))) -
           Geometry.UVVector(wcurl(Geometry.Covariant3Vector(curl(field))))
    )
    ∇⁴field = Spaces.weighted_dss!(
        @. Geometry.UVVector(wgrad(sdiv(χ))) -
           Geometry.UVVector(wcurl(Geometry.Covariant3Vector(curl(χ))))
    )
    return ∇⁴field
end

@static if @isdefined(var"@test_opt")

    function test_operators(field, vfield)
        @test_opt opt_Gradient(field)
        opt_WeakGradient(field)

        @test_opt opt_Curl(vfield)
        @test_opt opt_WeakCurl(vfield)
        @test_opt opt_CurlCurl(vfield)

        @test_opt opt_Divergence(vfield)
        @test_opt opt_WeakDivergence(vfield)

        @test_opt opt_ScalarDSS(field)
        @test_opt opt_VectorDss_Curl(vfield)
        @test_opt opt_VectorDss_DivGrad(vfield)

        @test_opt opt_ScalarHyperdiffusion(field)

        # TODO: Work on getting vector hyperdiffusion to optimize
        # after curl operator changes
        #@test_opt opt_VectorHyperdiffusion(vfield)
    end

end

# Test that Julia ia able to optimize spectral element operations v1.7+
@static if @isdefined(var"@test_opt")

    @testset "Spectral Element 2D Field optimizations" begin
        for FT in (Float64, Float32)
            domain = Domains.RectangleDomain(
                Geometry.XPoint{FT}(-pi) .. Geometry.XPoint{FT}(pi),
                Geometry.YPoint{FT}(-pi) .. Geometry.YPoint{FT}(pi);
                x1periodic = true,
                x2periodic = true,
            )

            Nq = 3
            quad = Spaces.Quadratures.GLL{Nq}()
            mesh = Meshes.RectilinearMesh(domain, 3, 3)

            topology = Topologies.Topology2D(mesh)
            space = Spaces.SpectralElementSpace2D(topology, quad)

            coords = Fields.coordinate_field(space)

            field = ones(FT, space)
            vfield =
                Geometry.UVVector.(
                    sin.(coords.x .+ 2 .* coords.y),
                    cos.(3 .* coords.x .+ 4 .* coords.y),
                )

            INq = 4
            Iquad = Spaces.Quadratures.GLL{INq}()
            Ispace =
                Spaces.SpectralElementSpace2D(Spaces.topology(space), Iquad)

            I = Operators.Interpolate(Ispace)
            @test_opt opt_Interpolate(I, field)

            Ifield = opt_Interpolate(I, field)
            R = Operators.Restrict(space)
            @test_opt opt_Restrict(R, Ifield)

            test_operators(field, vfield)
        end
    end

    @testset "Spectral Element 3D Hybrid Field optimizations" begin
        for FT in (Float64, Float32)
            xelem = 3
            yelem = 3
            velem = 5
            npoly = 3

            vertdomain = Domains.IntervalDomain(
                Geometry.ZPoint{FT}(0),
                Geometry.ZPoint{FT}(1000);
                boundary_tags = (:bottom, :top),
            )
            vertmesh = Meshes.IntervalMesh(vertdomain, nelems = velem)
            vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

            xdomain = Domains.IntervalDomain(
                Geometry.XPoint{FT}(-500) .. Geometry.XPoint{FT}(500),
                periodic = true,
            )
            ydomain = Domains.IntervalDomain(
                Geometry.YPoint{FT}(-100) .. Geometry.YPoint{FT}(100),
                periodic = true,
            )

            horzdomain = Domains.RectangleDomain(xdomain, ydomain)
            horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem)
            horztopology = Topologies.Topology2D(horzmesh)

            quad = Spaces.Quadratures.GLL{npoly + 1}()
            horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

            hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
                horzspace,
                vert_center_space,
            )
            hv_face_space =
                Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)

            for space in (hv_center_space, hv_face_space)
                coords = Fields.coordinate_field(space)
                field = ones(FT, space)
                vfield =
                    Geometry.UVWVector.(
                        sin.(coords.x .+ 2 .* coords.y),
                        cos.(3 .* coords.x .+ 4 .* coords.y),
                        zero(FT),
                    )
                test_operators(field, vfield)
            end
        end
    end
end


text:
"""
---
Quadrature Amplitude Modulation (QAM) function
   Apply symbol mapping to a input binary sequence (of size 1xL) with constellation size M.
	Output is a vector (1xN) with N = L / log2(M)
	Conventional gray mapping is used. Output constellation is casted in float, with unitary average power
 Supported constellation
* QPSK
* 16-QAM
* 64-QAM
* 256-QAM
# --- Syntax
  bitMappingQAM!(qamMat,M,bitSeq)
# --- Input parameters
- qamMat	: Complex Vector to populate of size length(bitSeq) / log2(M) [Array{Complex{Float64}}]
- M			: Modulation size (i.e from 4 to 256) such as bit per symbol is log2(M) [Int]
- bitSeq	: Binary sequence to be transformed into QPSK symbols [Array{UInt8}]
# --- Output parameters
- []
"""
function bitMappingQAM!(qamMat,M, bitSeq)
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	# --- Defining scaling factor
	# Constellation output should have an unitary average power
	# σ^2 = ∑ p(a_i)  ⃒ a_i ⃒ ^2 = 1
	scalingFactor	= sqrt(2/3*(M-1));
	nbSymb			= length(qamMat);
	# ----------------------------------------------------
	# --- Switch on modulation order
	# ----------------------------------------------------
	if M == 2
		# --- BPSK modulator
		# No need to rescale data, binary is direcly inherited from binary stream
		qamMat	.= 2 .*bitSeq[:] .- 1;
	elseif M == 4
		## --- QPSK modulator
		# QPSK mapping
		# Odd bit are in I; Even bits are in Q = [I1 Q1 ...]
		# -------------------------#
		#   |      |
		#   0      1
		#  -1     -1
		# 0 -> -1
		# 1 -> 1
		#nbSymb	  = Int64(length(bitSeq)/2);
		#qamMat	  = zeros(Complex{Float64},nbSymb);
		alphabet  = [-1 1];
		for i = 1 : 1 : nbSymb
			## Binary to Symbol conversion
			# Conversion from 2 bit to one symbol (I path):
			classI	  = bitSeq[2(i-1)+1];
			# Conversion from 2 bit to one symbol (Q path)
			classQ	  = bitSeq[2(i-1)+2];
			## Gray encoding
			outI  = alphabet[1+classI];
			outQ  = alphabet[1+classQ];
			# Setting output
			qamMat[i] = (outI + 1im*outQ)/scalingFactor;
		end
	elseif M == 16
		# 16-QAM mapping
		# Odd bit are in I; Even bits are in Q = [I1 Q1 I2 Q2 ...]
		# I1 :> MSB || I2:> LSB
		# -------------------------#
		#   |      |     |     |
		#  10     11    01    00
		#  -3     -1     1     3
		# 0 -> 3
		# 1 -> 1
		# 2 -> -3
		# 3 -> -1
		nbSymb	  = Int64(length(bitSeq)/4);
		alphabet  = [3 1 -3 -1];
		for i = 1 : 1 : nbSymb
			## Binary to Symbol conversion
			# Conversion from 2 bit to one symbol (I path): {0,1,2,3}
			classI	  = 2bitSeq[4(i-1)+1] + bitSeq[4(i-1)+3];
			# Conversion from 2 bit to one symbol (Q path)
			classQ	  = 2bitSeq[4(i-1)+2] + bitSeq[4(i-1)+4];
			## Gray encoding
			outI  = alphabet[1+classI];
			outQ  = alphabet[1+classQ];
			# Setting output
			qamMat[i] = (outI + 1im*outQ)/scalingFactor;
		end
	elseif M == 64
		# 64-QAM mapping
		# Odd bit are in I; Even bits are in Q = [I1 Q1 I2 Q2  I3 Q3 ...]
		# I1 :> MSB || I4:> LSB
		# ------------------------------------------------#
		#   |      |     |     |     |     |      |      |
		#  100    101   111   110   010   011    001    000
		#   -7     -5    -3    -1    +1    +3     +5     +7
		# 0 -> 000 -> +7
		# 1 -> 001 -> +5
		# 2 -> 010 -> +1
		# 3 -> 011 -> +3
		# 4 -> 100 -> -7
		# 5 -> 101 -> -5
		# 6 -> 110 -> -1
		# 7 -> 111 -> -3
		#nbSymb	  = Int64(length(bitSeq)/6);
		#qamMat	  = zeros(Complex{Float64},nbSymb);
		alphabet  = [7 5 1 3 -7 -5 -1 -3];
		for i = 1 : 1 : nbSymb
			## Binary to Symbol conversion
			# Conversion from 2 bit to one symbol (I path): {0,1,2,3}
			classI	  = 4bitSeq[6(i-1)+1] + 2bitSeq[6(i-1)+3] + bitSeq[6(i-1)+5];
			# Conversion from 2 bit to one symbol (Q path)
			classQ	  = 4bitSeq[6(i-1)+2] + 2bitSeq[6(i-1)+4] + bitSeq[6(i-1)+6];
			## Gray encoding
			outI  = alphabet[1+classI];
			outQ  = alphabet[1+classQ];
			# Setting output
			qamMat[i] = (outI + 1im*outQ)/scalingFactor;
		end
	elseif M == 256
		# 256-QAM mapping
		# Inherited from 64-QAM with classic tree genereation
		# -> [0 Sequence(1:end)] [1 Sequence(end:-1:1)]
		# Odd bit are in I; Even bits are in Q = [I1 Q1 I2 Q2  I3 Q3 I4 Q4...]
		# I1 :> MSB || I4:> LSB
		# ------------------------------------------------#
		#	   |       |         |        |        |        |        |        |
		#	0100    0101 .    0111 .   0110 .   0010 .   0011 .   0001 .   0000
		#	 -15     -13 .     -11 .     -9 .     -7 .     -5 .     -3 .     -1
		#	   |       |         |        |        |        |        |        |
		#   1000 .   1001 .   1011 .   1010 .   1110 .   1111 .   1101 .   1100
		#     +1 .     +3 .     +5 .     +7 .     +9 .    +11 .    +13      +15.
		# 0 ->  0000 -> -1
		# 1 ->  0001 -> -3
		# 2 ->  0010 -> -7
		# 3 ->  0011 -> -5
		# 4 ->  0100 -> -15
		# 5 ->  0101 -> -13
		# 6 ->  0110 -> -9
		# 7 ->  0111 -> -11
		# 8 ->  1000 -> +1
		# 9 ->  1001 -> +3
		# 10 -> 1010 -> +7
		# 11 -> 1011 -> +5
		# 12 -> 1100 -> +15
		# 13 -> 1101 -> +13
		# 14 -> 1110 -> +9
		# 15 -> 1111 -> +11
		#nbSymb	  = Int64(length(bitSeq)/8);	#
		#qamMat	  = zeros(Complex{Float64},nbSymb);
		alphabet  = [-1 -3 -7 -5 -15 -13 -9 -11 +1 +3 +7 +5 +15 +13 +9 +11];
		for i = 1 : 1 : nbSymb
			## Binary to Symbol conversion
			# Conversion from 2 bit to one symbol (I path): {0,1,2,3}
			classI	  = 8bitSeq[8(i-1)+1] + 4bitSeq[8(i-1)+3] + 2bitSeq[8(i-1)+5] + bitSeq[8(i-1)+7];
			# Conversion from 2 bit to one symbol (Q path)
			classQ	  = 8bitSeq[8(i-1)+2] + 4bitSeq[8(i-1)+4] + 2bitSeq[8(i-1)+6] + bitSeq[8(i-1)+8];
			## Gray encoding
			outI  = alphabet[1+classI];
			outQ  = alphabet[1+classQ];
			# Setting output
			qamMat[i] = (outI + 1im*outQ)/scalingFactor;
		end
	end
end

# ----------------------------------------------------
# --- Allocating function
# ----------------------------------------------------
"""
---
Quadrature Amplitude Modulation (QAM) function
   Apply symbol mapping to a input binary sequence (of size 1xL) with constellation size M.
	Output is a vector (1xN) with N = L / log2(M)
	Conventional gray mapping is used. Output constellation is casted in float, with unitary average power
 Supported constellation
* QPSK
* 16-QAM
* 64-QAM
* 256-QAM
# --- Syntax
qamMat =  bitMappingQAM(M,bitSeq)
# --- Input parameters
- M			: Modulation size (i.e from 4 to 256) such as bit per symbol is log2(M) [Int]
- bitSeq	: Binary sequence to be transformed into QPSK symbols [Array{UInt8}]
# --- Output parameters
- qamMat	: Complex Vector to populate of size length(bitSeq) / log2(M) [Array{Complex{Float64}}]
"""
function bitMappingQAM(M,bitSeq)
	nbBits	= length(bitSeq);
	nbSymb	= Int( nbBits ÷ log2(M));
	buffer	= zeros(Complex{Float64},nbSymb);
	bitMappingQAM!(buffer,M,bitSeq);
	return buffer;
end

# ----------------------------------------------------
# --- Multiple dispatch handling
# ----------------------------------------------------
# --- MD: String case for modulation order
function bitMappingQAM!(qamMat,M::String, bitSeq)
	# --- Casting modulation order to int8
	if M == "QPSK" || M == "4-QAM" || M == "QAM-4" || M=="4QAM" || M == "QAM4"
		bitMappingQAM!(qamMat,4,bitSeq);
	elseif M == "16-QAM" || M == "QAM-16" || M=="16QAM" || M == "QAM16"
		bitMappingQAM!(qamMat,16,bitSeq);
	elseif M == "64-QAM" || M == "QAM-64" || M=="64QAM" || M == "QAM64"
		bitMappingQAM!(qamMat,64,bitSeq);
	elseif M == "256-QAM" || M == "QAM-256" || M=="256QAM" || M == "QAM256"
		bitMappingQAM!(qamMat,256,bitSeq);
	end
end
function bitMappingQAM(M::String, bitSeq)
	# --- Casting modulation order to int8
	if M == "QPSK" || M == "4-QAM" || M == "QAM-4" || M=="4QAM" || M == "QAM4"
		bitMappingQAM(4,bitSeq);
	elseif M == "16-QAM" || M == "QAM-16" || M=="16QAM" || M == "QAM16"
		bitMappingQAM(16,bitSeq);
	elseif M == "64-QAM" || M == "QAM-64" || M=="64QAM" || M == "QAM64"
		bitMappingQAM(64,bitSeq);
	elseif M == "256-QAM" || M == "QAM-256" || M=="256QAM" || M == "QAM256"
		bitMappingQAM(256,bitSeq);
	end
end


text:
using SparseArrays;

"""
  function naive_petsc_read(filename, int_type=Int32, float_type=ComplexF64)

Loading petsc data into julia arrays and vector
The function does not support all Petsc class types.
It only supports sparse arrays and vectors.

File format documentation:
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatLoad.html
https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecLoad.html
"""

function naive_petsc_read(filename, int_type=Int32, float_type=ComplexF64)
    open(filename) do io
        class_id = ntoh(read(io, int_type))

        MAT_FILE_CLASSID=1211216;
        VEC_FILE_CLASSID=1211214;

        if (class_id != MAT_FILE_CLASSID  && class_id != VEC_FILE_CLASSID)
            @show class_id
            error("Unsupported class_id. This function can only load sparse arrays and vectors: ")
        end

        if (class_id == MAT_FILE_CLASSID)
            # It's a sparse array

            # Get size
            rows = ntoh(read(io, int_type)) # From documentation: number of rows
            cols = ntoh(read(io, int_type)) # From documentation: number of columns

            # Number of nnz elements
            nnz = ntoh(read(io, int_type))  # From documentation: total number of nonzeros

            # Prep for reading
            row_ptr = Vector{int_type}(undef, rows+1)
            row_ptr[1] = 1

            ## Read data

            # row lengths
            # From documentation: number nonzeros in each row
            row_ptr[2:end] = ntoh.(read!(io, Vector{int_type}(undef, rows)))
            cumsum!(row_ptr, row_ptr)

            # columns
            # From documentation: column indices of all nonzeros
            colvals = ntoh.(read!(io, Vector{int_type}(undef, nnz)))
            colvals .+= int_type(1) # Start index is 0 in PETsc and 1 in Julia

            # From documentation: values of all nonzeros
            vals = ntoh.(read!(io, Vector{float_type}(undef, nnz)))

            mat = SparseMatrixCSC(cols, rows, row_ptr, colvals, vals)
            return SparseMatrixCSC(transpose(mat))
        else
            # It's a vec

            # From documentation: number of rows
            rows = ntoh(read(io, int_type))

            # From documentation: values of all entries
            vec=ntoh.(read!(io, Vector{float_type}(undef, rows)))
        end

    end

end


text:
# using GLMakie
# using CairoMakie
using Plots

function vfplot(v,g;kwargs...)
x=g.x
y=g.y
u =getindex.(v,1)
v =getindex.(v,2)
x, y, u, v = vec.([x, y, u, v])
quiver(x, y, quiver = (u, v), c = :blue)
end
# function vis(x::EquivConv)
#     r=LinRange(0.,x.rmax,64)
#     plots = [Plots.plot(r,p.op.radfunc.(r)) for p in x.paths]
#     Plots.plot(plots..., layout = length(plots))
# end


text:
# ---
# title: 1155. Number of Dice Rolls With Target Sum
# id: problem1155
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-dice-rolls-with-target-sum/description/>
# hidden: true
# ---
#
# You have `d` dice, and each die has `f` faces numbered `1, 2, ..., f`.
#
# Return the number of possible ways (out of `fd` total ways) **modulo`10^9 +
# 7`** to roll the dice so the sum of the face up numbers equals `target`.
#
#
#
# **Example 1:**
#
#
#
#     Input: d = 1, f = 6, target = 3
#     Output: 1
#     Explanation:
#     You throw one die with 6 faces.  There is only one way to get a sum of 3.
#
#
# **Example 2:**
#
#
#
#     Input: d = 2, f = 6, target = 7
#     Output: 6
#     Explanation:
#     You throw two dice, each with 6 faces.  There are 6 ways to get a sum of 7:
#     1+6, 2+5, 3+4, 4+3, 5+2, 6+1.
#
#
# **Example 3:**
#
#
#
#     Input: d = 2, f = 5, target = 10
#     Output: 1
#     Explanation:
#     You throw two dice, each with 5 faces.  There is only one way to get a sum of 10: 5+5.
#
#
# **Example 4:**
#
#
#
#     Input: d = 1, f = 2, target = 3
#     Output: 0
#     Explanation:
#     You throw one die with 2 faces.  There is no way to get a sum of 3.
#
#
# **Example 5:**
#
#
#
#     Input: d = 30, f = 30, target = 500
#     Output: 222616187
#     Explanation:
#     The answer must be returned modulo 10^9 + 7.
#
#
#
#
# **Constraints:**
#
#   * `1 <= d, f <= 30`
#   * `1 <= target <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Sanity checking for the CoNLL-X format parsing module.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-07

include("../hdr.jl")

require("conllx.jl")

using CoNLLX

respath = string(dirname(source_path()), "/../../res/")
debugpath = string(respath, "debug.conllx")

# Check that reading in and then writing the data causes no corruption.
lines = String[]
open(debugpath, "r") do debug_f
    for line in eachline(debug_f)
        if line[1] == '#'
            continue
        end

        push!(lines, line)
    end
end

before = join(lines, "")

parsed = collect(Sentence, conllxparse(IOBuffer(before)))
buff = IOBuffer()
print(buff, parsed)
seekstart(buff)
after = strip(readall(buff))

cmp = strip(before)
if after != cmp
    before_f = tempname()
    after_f = tempname()
    try
        run(`echo $cmp` |> before_f)
        run(`echo $after` |> after_f)
        run(ignorestatus(`diff -u $before_f $after_f`))
    finally
        rm(after_f)
        rm(before_f)
    end
    @assert false
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Internal function to keep track of string state when iterating through the
#   characters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _process_string_state(c::Char, state::Symbol = :text)

Return the current state of the string given the new character `c` and the
previous state `state`.

The following states are possible:

- `:text`: The character is part of the printable text.
- `:escape_state_begin`: Beginning of an ANSI escape sequence.
- `:escape_state_opening`: Matches the `[`.
- `:escape_state_1`: First state of an ANSI escape sequence.
- `:escape_state_2`: Second state of an ANSI escape sequence.
- `:escape_state_3`: Third state of an ANSI escape sequence.
- `:escape_state_end`: End of an ANSI escape sequence.
"""
function _process_string_state(c::Char, state::Symbol = :text)
    if state == :text
        # Here, we need to check if an escape sequence is found.
        if c == '\x1b'
            state = :escape_state_begin
        end

    elseif state == :escape_state_begin
        if (c == '[')
            state = :escape_state_opening

        elseif ('@' ≤ c ≤ 'Z') || ('\\' ≤ c ≤ '_')
            state = :escape_state_1
        else
            state = :text
        end

    elseif state == :escape_state_opening
        state = :escape_state_1
        return _process_string_state(c, state)

    elseif state == :escape_state_1
        if !('0' ≤ c ≤ '?')
            state = :escape_state_2
            return _process_string_state(c, state)
        end

    elseif state == :escape_state_2
        if !(' ' ≤ c ≤ '/')
            state = :escape_state_3
            return _process_string_state(c, state)
        end

    elseif state == :escape_state_3
        if ('@' ≤ c ≤ '~')
            state = :escape_state_end
        else
            state = :text
        end
    elseif state == :escape_state_end
        state = :text

        # We need to recall this function because the next character can be the
        # beginning of a new ANSI escape sequence.
        return _process_string_state(c, state)
    end

    return state
end


text:
const MODULE_GLOBALS = Dict{Module,Py}()

function _pyeval_args(globals, locals)
    if globals isa Module
        globals_ = get!(pydict, MODULE_GLOBALS, globals)
    elseif ispy(globals)
        globals_ = globals
    else
        ArgumentError("globals must be a module or a Python dict")
    end
    if locals === nothing
        locals_ = Py(globals_)
    elseif ispy(locals)
        locals_ = Py(locals)
    else
        locals_ = pydict(locals)
    end
    return (globals_, locals_)
end

"""
    pyeval([T=Py], code, globals, locals=nothing)

Evaluate the given Python `code`, returning the result as a `T`.

If `globals` is a `Module`, then a persistent `dict` unique to that module is used.

By default the code runs in global scope (i.e. `locals===globals`). To use a temporary
local scope, set `locals` to `()`, or to a `NamedTuple` of variables to include in the
scope.

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

# Examples

The following computes `1.1+2.2` in the `Main` module as a `Float64`:
```
pyeval(Float64, "x+y", Main, (x=1.1, y=2.2))  # returns 3.3
```
"""
function pyeval(::Type{T}, code, globals, locals=nothing) where {T}
    globals_, locals_ = _pyeval_args(globals, locals)
    ans = pybuiltins.eval(code, globals_, locals_)
    pydel!(locals_)
    return T == Py ? ans : pyconvert_and_del(T, ans)
end
pyeval(code, globals, locals=nothing) = pyeval(Py, code, globals, locals)
export pyeval

_pyexec_ans(::Type{Nothing}, globals, locals) = nothing
@generated function _pyexec_ans(::Type{NamedTuple{names, types}}, globals, locals) where {names, types}
    # TODO: use precomputed interned strings
    # TODO: try to load from globals too
    n = length(names)
    code = []
    vars = Symbol[]
    for i in 1:n
        v = Symbol(:ans, i)
        push!(vars, v)
        push!(code, :($v = pyconvert_and_del($(types.parameters[i]), pygetitem(locals, $(string(names[i]))))))
    end
    push!(code, :(return $(NamedTuple{names, types})(($(vars...),))))
    return Expr(:block, code...)
end

"""
    pyexec([T=Nothing], code, globals, locals=nothing)

Execute the given Python `code`.

If `globals` is a `Module`, then a persistent `dict` unique to that module is used.

By default the code runs in global scope (i.e. `locals===globals`). To use a temporary
local scope, set `locals` to `()`, or to a `NamedTuple` of variables to include in the
scope.

If `T==Nothing` then returns `nothing`. Otherwise `T` must be a concrete `NamedTuple` type
and the corresponding items from `locals` are extracted and returned.

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

# Examples

The following computes `1.1+2.2` in the `Main` module as a `Float64`:
```
pyexec(@NamedTuple{ans::Float64}, "ans=x+y", Main, (x=1.1, y=2.2))  # returns (ans = 3.3,)
```

Marking variables as `global` saves them into the module scope, so that they are available
in subsequent invocations:
```
pyexec("global x; x=12", Main)
pyeval(Int, "x", Main)  # returns 12
```
"""
function pyexec(::Type{T}, code, globals, locals=nothing) where {T}
    globals_, locals_ = _pyeval_args(globals, locals)
    pydel!(pybuiltins.exec(code, globals_, locals_))
    ans = _pyexec_ans(T, globals_, locals_)
    pydel!(locals_)
    return ans
end
pyexec(code, globals, locals=nothing) = pyexec(Nothing, code, globals, locals)
export pyexec

function _pyeval_macro_code(arg)
    if arg isa String
        return arg
    elseif arg isa Expr && arg.head === :macrocall && arg.args[1] == :(`foo`).args[1]
        return arg.args[3]
    else
        return nothing
    end
end

function _pyeval_macro_args(arg, filename, mode)
    # separate out inputs => code => outputs (with only code being required)
    if @capture(arg, inputs_ => code_ => outputs_)
        code = _pyeval_macro_code(code)
        code === nothing && error("invalid code")
    elseif @capture(arg, lhs_ => rhs_)
        code = _pyeval_macro_code(lhs)
        if code === nothing
            code = _pyeval_macro_code(rhs)
            code === nothing && error("invalid code")
            inputs = lhs
            outputs = nothing
        else
            inputs = nothing
            outputs = rhs
        end
    else
        code = _pyeval_macro_code(arg)
        code === nothing && error("invalid code")
        inputs = outputs = nothing
    end
    # precompile the code
    codestr = code
    codeobj = pynew()
    codeready = Ref(false)
    code = quote
        if !$codeready[]
            $pycopy!($codeobj, $pybuiltins.compile($codestr, $filename, $mode))
            $codeready[] = true
        end
        $codeobj
    end
    # convert inputs to locals
    if inputs === nothing
        locals = ()
    else
        if inputs isa Expr && inputs.head === :tuple
            inputs = inputs.args
        else
            inputs = [inputs]
        end
        locals = []
        for input in inputs
            if @capture(input, var_Symbol)
                push!(locals, var => var)
            elseif @capture(input, var_Symbol = ex_)
                push!(locals, var => ex)
            else
                error("invalid input: $input")
            end
        end
        locals = :(($([:($var = $ex) for (var,ex) in locals]...),))
    end
    # done
    return locals, code, outputs
end

"""
    @pyeval [inputs =>] code [=> T]

Evaluate the given `code` in a new local scope and return the answer as a `T`.

The global scope is persistent and unique to the current module.

The `code` must be a literal string or command.

The `inputs` is a tuple of inputs of the form `v=expr` to be included in the local scope.
Only `v` is required, `expr` defaults to `v`.

# Examples

The following computes `1.1+2.2` and returns a `Float64`:
```
@pyeval (x=1.1, y=2.2) => `x+y` => Float64  # returns 3.3
```
"""
macro pyeval(arg)
    locals, code, outputs = _pyeval_macro_args(arg, "$(__source__.file):$(__source__.line)", "eval")
    if outputs === nothing
        outputs = Py
    end
    esc(:($pyeval($outputs, $code, $__module__, $locals)))
end
export @pyeval

"""
    @pyexec [inputs =>] code [=> outputs]

Execute the given `code` in a new local scope.

The global scope is persistent and unique to the current module.

The `code` must be a literal string or command.

The `inputs` is a tuple of inputs of the form `v=expr` to be included in the local scope.
Only `v` is required, `expr` defaults to `v`.

The `outputs` is a tuple of outputs of the form `x::T=v`, meaning that `v` is extracted from
locals, converted to `T` and assigned to `x`. Only `x` is required: `T` defaults to `Py`
and `v` defaults to `x`.

# Examples

The following computes `1.1+2.2` and assigns its value to `ans` as a `Float64`:
```
@pyexec (x=1.1, y=2.2) => `ans=x+y` => ans::Float64  # returns 3.3
```

Marking variables as `global` saves them into the module scope, so that they are available
in subsequent invocations:
```
@pyexec `global x; x=12`
@pyeval `x` => Int  # returns 12
```
"""
macro pyexec(arg)
    locals, code, outputs = _pyeval_macro_args(arg, "$(__source__.file):$(__source__.line)", "exec")
    if outputs === nothing
        outputs = Nothing
        esc(:($pyexec(Nothing, $code, $__module__, $locals)))
    else
        if outputs isa Expr && outputs.head === :tuple
            oneoutput = false
            outputs = outputs.args
        else
            oneoutput = true
            outputs = [outputs]
        end
        pyvars = Symbol[]
        jlvars = Symbol[]
        types = []
        for output in outputs
            if @capture(output, lhs_ = rhs_)
                rhs isa Symbol || error("invalid output: $output")
                output = lhs
                pyvar = rhs
            else
                pyvar = missing
            end
            if @capture(output, lhs_ :: rhs_)
                outtype = rhs
                output = lhs
            else
                outtype = Py
            end
            output isa Symbol || error("invalid output: $output")
            if pyvar === missing
                pyvar = output
            end
            push!(pyvars, pyvar)
            push!(jlvars, output)
            push!(types, outtype)
        end
        outtype = :($NamedTuple{($(map(QuoteNode, pyvars)...),), Tuple{$(types...),}})
        ans = :($pyexec($outtype, $code, $__module__, $locals))
        if oneoutput
            ans = :($(jlvars[1]) = $ans[1])
        else
            if pyvars != jlvars
                outtype2 = :($NamedTuple{($(map(QuoteNode, jlvars)...),), Tuple{$(types...),}})
                ans = :($outtype2($ans))
            end
            ans = :(($(jlvars...),) = $ans)
        end
        esc(ans)
    end
end
export @pyexec


text:
import Clp

include("LPBenchmark.jl")

function run_clp(fname::String)
    clp = Clp.Optimizer()
    cache = MOIU.UniversalFallback(MOIU.Model{Float64}())
    cached = MOIU.CachingOptimizer(cache, clp)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)

    load_problem!(bridged, fname)

    # Set some parameters
    MOI.set(bridged, MOI.TimeLimitSec(), 10_000)
    MOI.set(bridged, MOI.RawParameter("SolveType"), 4)  # barrier, no crossover

    # Solve the problem
    run(bridged)

    # Display timing info
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_clp(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
using DelimitedFiles;
using BioSequences;
using BioAlignments;
using CSV;
using DataFrames;

data = readdlm("/Users/stefftaelman/Downloads/Documents_local/PROJECTS/PAPER2019/PhaLP_unique_seqs.csv", ',', '\n')
seqs = data[:, 1]

### started around 10:05pm on 29/04
sim_matrix = zeros(length(seqs), length(seqs))
for i in 1:length(seqs)
    for j in i:length(seqs)
        tmp = pairalign(LocalAlignment(), seqs[i], seqs[j], AffineGapScoreModel(BLOSUM62, gap_open=-10, gap_extend=-1); score_only=true)
        sim_matrix[i, j] = score(tmp)
        sim_matrix[j, i] = score(tmp)
    end
end  #started 30/12 around 14:45

CSV.write("/Users/stefftaelman/Downloads/Documents_local/PROJECTS/PAPER2019/similarity_matrix_010120.csv",  DataFrame(sim_matrix), writeheader=false)

max_stars_repo_name: bjcriel/PhaLP

max_issues_repo_name: bjcriel/PhaLP

max_forks_repo_name: bjcriel/PhaLP


text:
function get_next_gridpoints(trade_grid, price)
    k = sort(string.(keys(trade_grid)))
    index = findall(x->x == price, k)[1]
    result = Dict()
    if index > 1
        result["BUY"] = k[index-1]
    end
    if index < length(k)
        result["SELL"] = k[index+1]
    end
    return result
end
export get_next_gridpoints


text:
module NetworkDismantle

using LightGraphs
using DataStructures

export
HDA,
CoreHD,
WeakNei,
high_degree_adaptive,
decore,
weak_neighbor_decore,
tree_break!,
reverse_greedy!,
giant_component_size,
recover_add_nodes

include("dismantle.jl")
include("high_degree_adaptive.jl")
include("decore.jl")
include("weak_neighbor_decore.jl")
include("tree_break.jl")
include("reverse_greedy.jl")
include("utils.jl")

end # module


text:
@testset "Herm" begin
  include("Herm/Spaces.jl")
  include("Herm/Genus.jl")
  include("Herm/GenusRep.jl")
end


text:
function search_subgraph_ss(g::Graph, ks::Vector{Int})
  for k=ks
    tmp_result = search_subgraph_ss(g, k)
    if !isempty(tmp_result)
      return tmp_result
    end
  end
  return k
end

function search_subgraph_ss(g::Graph, k::Int)
  @assert 2 <= k <= nv(g)-1 "subgraphs order must be at least 2 and at most order of graph minus 1"
  for vertices = subsets(1:nv(g), k)
    ind_g = induced_subgraph(g, vertices)[1]
    if is_connected(ind_g) && have_ss_connected(g, vertices)
      return vertices
    end
  end
  return Int[]
end

function search_subgraph_ss(g::Graph)
  for vertices = subsets(1:nv(g))
    if 2 <= length(vertices) <= nv(g)-1
      ind_g = induced_subgraph(g, vertices)[1]
      if is_connected(ind_g) && have_ss_connected(g, vertices)
        return vertices
      end
    end
  end
  return Int[]
end

function search_subgraph_noss(g::Graph, ks::Vector{Int})
  for k=ks
    tmp_result = search_subgraph_noss(g, k)
    if !isempty(tmp_result)
      return tmp_result
    end
  end
  return k
end

function search_subgraph_noss(g::Graph, k::Int)
  @assert 2 <= k <= nv(g)-1 "subgraphs order must be at least 2 and at most order of graph minus 1"
  for vertices = subsets(1:nv(g), k)
    ind_g = induced_subgraph(g, vertices)[1]
    if is_connected(ind_g) && !have_ss_connected(g, vertices)
      return vertices
    end
  end
  return Int[]
end

function search_subgraph_noss(g::Graph)
  for vertices = subsets(1:nv(g))
    if 2 <= length(vertices) <= nv(g)-1
      ind_g = induced_subgraph(g, vertices)[1]
      if is_connected(ind_g) && !have_ss_connected(g, vertices)
        return vertices
      end
    end
  end
  return Int[]
end

max_stars_repo_name: QuantumWalks/ExceptionalConfiguration

max_issues_repo_name: QuantumWalks/ExceptionalConfiguration

max_forks_repo_name: QuantumWalks/ExceptionalConfiguration


text:
"""
    function addIntegratorMediumLax_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions related with Medium Update.
"""
function addIntegratorMediumLax_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    if "UpdateMedium" in keys(abm.declaredUpdates)

        #Check updated
        up = symbols_(abm,abm.declaredUpdates["UpdateMedium"])
        up = up[Bool.((up[:,"placeDeclaration"].==:Model) .* Bool.((up[:,"assigned"].==true) .+ (up[:,"updated"].==true))),:]

        #Construct functions

            #Make function to compute everything inside the volume

        f = abm.declaredUpdates["UpdateMedium"]
        for (i,j) in enumerate(abm.declaredSymbols["Medium"]) # Change symbol for update

            if abm.dims == 1
                add = :((mediumV[ic1_+1,$i] + mediumV[ic1_+1,$i])/2)
            elseif abm.dims == 2
                add = :((mediumV[ic1_-1,ic2_,$i] + mediumV[ic1_+1,ic2_,$i]
                        +mediumV[ic1_,ic2_-1,$i] + mediumV[ic1_,ic2_+1,$i])/4)
            elseif abm.dims == 3
                add = :((mediumV[ic1_-1,ic2_,ic3_,$i] + mediumV[ic1_+1,ic2_,ic3_,$i]
                        +mediumV[ic1_,ic2_-1,ic3_,$i] + mediumV[ic1_,ic2_+1,ic3_,$i]
                        +mediumV[ic1_,ic2_,ic3_-1,$i] + mediumV[ic1_,ic2_,ic3_+1,$i])/6)
            end

            s = Meta.parse(string(MEDIUMSYMBOL,j))
            f = postwalk(x -> @capture(x,ss_=v_) && ss == s ? :($j = $add + $v*dt) : x, f)
        end
        f = postwalk(x->@capture(x,s_(v_)) ? :($(adaptOperatorsMediumLax_(v,s,abm,space,p))) : x, f) # Adapt
        f = postwalk(x->@capture(x,s_) ? adaptSymbolsMediumLax_(s,abm,space,p) : x, f) # Adapt

        f = vectorizeMedium_(abm,f,p)

        f = simpleGridLoopWrapInFunction_(platform,:mediumInnerStep_!, f, abm.dims)

        #Remove count over boundaries
        if platform == "cpu"
            f = postwalk(x->@capture(x,1:Nx_) && Nx == :Nx_ ? :(2:Nx_-1) : x, f)
            f = postwalk(x->@capture(x,1:Ny_) && Ny == :Ny_ ? :(2:Ny_-1) : x, f)
            f = postwalk(x->@capture(x,1:Nz_) && Nz == :Nz_ ? :(2:Nz_-1) : x, f)
        elseif platform == "gpu"
            f = postwalk(x->@capture(x,indexX_:strideX_:Nx_) && Nx == :Nx_ ? :(indexX_+1:strideX_:Nx_-1) : x, f)
            f = postwalk(x->@capture(x,indexY_:strideY_:Ny_) && Ny == :Ny_ ? :(indexY_+1:strideY_:Ny_-1) : x, f)
            f = postwalk(x->@capture(x,indexZ_:strideZ_:Nz_) && Nz == :Nz_ ? :(indexZ_+1:strideZ_:Nz_-1) : x, f)
        end

            #Make function to compute the boundaries
        code = quote end

        lUpdate = length(abm.declaredSymbols["Medium"])
        for i in 1:abm.dims #Adapt boundary updates depending on the boudary type
            ic = Meta.parse(string("ic",i,"_"))
            nm = Meta.parse(string("N",["x","y","z"][i],"_"))
            subcode = quote end

            if abm.dims == 1
                v = :(mediumVCopy[ic1_,ic4_])
            elseif abm.dims == 2
                v = :(mediumVCopy[ic1_,ic2_,ic4_])
            elseif abm.dims == 3
                v = :(mediumVCopy[ic1_,ic2_,ic3_,ic4_])
            end

            if space.medium[i].minBoundaryType == "Periodic"

                #Lower boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)
                vAss = postwalk(x->@capture(x, t_) && t == ic ? :($nm-1) : x, v)

                push!(subcode.args, :($vUp=$vAss))

                #Upper boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? :($nm) : x, v)
                vAss = postwalk(x->@capture(x, t_) && t == ic ? 2 : x, v)

                push!(subcode.args, :($vUp=$vAss))
            end

            if space.medium[i].minBoundaryType == "Newmann"

                #Lower boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)
                vAss = postwalk(x->@capture(x, t_) && t == ic ? :(2) : x, v)

                push!(subcode.args, :($vUp=$vAss))
            end

            if space.medium[i].maxBoundaryType == "Newmann"

                #Lower boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? nm : x, v)
                vAss = postwalk(x->@capture(x, t_) && t == ic ? :($nm-1) : x, v)

                push!(subcode.args, :($vUp=$vAss))
            end

            if space.medium[i].minBoundaryType == "Dirichlet"

                #Lower boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? 1 : x, v)

                push!(subcode.args, :($vUp=0.))
            end

            if space.medium[i].maxBoundaryType == "Dirichlet"

                #Lower boundary
                vUp = postwalk(x->@capture(x, t_) && t == ic ? nm : x, v)

                push!(subcode.args, :($vUp=0.))
            end

            subcode = :(for ic4_ in 1:$lUpdate
                            $subcode
                        end
                    )

            subcode = simpleGridLoop_(platform,subcode,abm.dims-1, indexes = [1,2,3][findall([1,2,3].!=i)])

            push!(code.args,subcode)
        end

        f2 = wrapInFunction_(:mediumBoundaryStep_!,code)

        if "UpdateMediumInteraction" in keys(p.agent.declaredUpdates)
            fWrap = wrapInFunction_(:mediumStep_!,
                                    :(begin
                                        @platformAdapt mediumInnerStep_!(ARGS_)
                                        @platformAdapt mediumBoundaryStep_!(ARGS_)
                                        @platformAdapt mediumInteractionStep_!(ARGS_)
                                    end)
                                    )
        else
            fWrap = wrapInFunction_(:mediumStep_!,
                                    :(begin
                                        @platformAdapt mediumInnerStep_!(ARGS_)
                                        @platformAdapt mediumBoundaryStep_!(ARGS_)
                                    end)
                                    )
        end

        # println(prettify(code))

        push!(p.declareF.args,f,f2,fWrap)

        push!(p.execInloop.args,
                :(mediumStep_!(ARGS_))
            )
    end

    return nothing
end

function adaptSymbolsMediumLax_(s,abm,space,p)

    if s == :xₘ

        return :(ic1_*dxₘ_+$(space.box[1].min))

    elseif s == :yₘ

        return :(ic2_*dyₘ_+$(space.box[2].min))

    elseif s == :zₘ

        return :(ic1_*dzₘ_+$(space.box[3].min))

    end

    return s

end

function adaptOperatorsMediumLax_(f,op,abm,space,p)

    f = vectorizeMedium_(abm,f,p)

    if op == :δx

        if space.medium[1].minBoundaryType == "Periodic"
            f = :(δMedium_(round(Int,($f-$(space.box[1].min))/dxₘ_)+1,ic1_))
        else
            f = :(δMedium_(round(Int,($f-$(space.box[1].min))/dxₘ_),ic1_))
        end

    elseif op == :δy

        if abm.dims < 2
            error("Operator δy cannot be used in Agent with dimension 1.")
        end

        if space.medium[2].minBoundaryType == "Periodic"
            f = :(δMedium_(round(Int,($f-$(space.box[2].min))/dyₘ_)+1,ic2_))
        else
            f = :(δMedium_(round(Int,($f-$(space.box[2].min))/dyₘ_),ic2_))
        end

    elseif op == :δz

        if abm.dims < 3
            error("Operator δz cannot be used in Agent with dimension 1 or 2.")
        end

        if space.medium[3].minBoundaryType == "Periodic"
            f = :(δMedium_(round(Int,($f-$(space.box[3].min))/dzₘ_)+1,ic3_))
        else
            f = :(δMedium_(round(Int,($f-$(space.box[3].min))/dzₘ_),ic3_))
        end

    elseif op == :Δ
        if abm.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif abm.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2+($addy+$subsy-2*$f)/dyₘ_^2)
        elseif abm.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2+($addy+$subsy-2*$f)/dyₘ_^2+($addz+$subsz-2*$f)/dzₘ_^2)

        end
    elseif op == :Δx
        if abm.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif abm.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        elseif abm.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)

            f = :(($add+$subs-2*$f)/dxₘ_^2)
        end
    elseif op == :Δy
        if abm.dims == 1
            error("Operator Δy cannot exist in Agent with dimension 1.")
        elseif abm.dims == 2
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($addy+$subsy-2*$f)/dyₘ_^2)
        elseif abm.dims == 3
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)

            f = :(($addy+$subsy-2*$f)/dyₘ_^2)
        end
    elseif op == :Δz
        if abm.dims == 1
            error("Operator Δz cannot exist in Agent with dimension 1.")
        elseif abm.dims == 2
            error("Operator Δy cannot exist in Agent with dimension 2.")
        elseif abm.dims == 3
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($addz+$subsz-2*$f)/dzₘ_^2)
        end
    elseif op == :∇
        if abm.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif abm.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_)+($addy-$subsy)/(2*dyₘ_))
        elseif abm.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_)+($addy-$subsy)/(2*dyₘ_)+($addz-$subsz)/(2*dzₘ_))
        end
    elseif op == :∇x
        if abm.dims == 1
            add = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif abm.dims == 2
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        elseif abm.dims == 3
            add = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_+1,ic2_,ic3_,$j]) : x,f)
            subs = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_-1,ic2_,ic3_,$j]) : x,f)

            f = :(($add-$subs)/(2*dxₘ_))
        end
    elseif op == :∇y
        if abm.dims == 1
            error("Operator Δy cannot exist in Agent with dimension 1.")
        elseif abm.dims == 2
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,$j]) : x,f)

            f = :(($addy-$subsy)/(2*dyₘ_))
        elseif abm.dims == 3
            addy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_+1,ic3_,$j]) : x,f)
            subsy = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_-1,ic3_,$j]) : x,f)

            f = :(($addy-$subsy)/(2*dyₘ_))
        end
    elseif op == :∇z
        if abm.dims == 1
            error("Operator Δz cannot exist in Agent with dimension 1.")
        elseif abm.dims == 2
            error("Operator Δy cannot exist in Agent with dimension 2.")
        elseif abm.dims == 3
            addz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_+1,$j]) : x,f)
            subsz = postwalk(x->@capture(x,s_[ic1_,ic2_,ic3_,j_]) && s == :mediumV ? :(mediumV[ic1_,ic2_,ic3_-1,$j]) : x,f)

            f = :(($addz-$subsz)/(2*dzₘ_))
        end
    end

    return f

end

δMedium_(i1,i2) = if i1==i2 1. else 0. end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
# test sonar functions directly

using RoME
using DistributedFactorGraphs
using Distributions
using Manifolds: TranslationGroup
using Test
##

@testset "bearing range elevation functions" begin

##

meas = LinearRangeBearingElevation((3.0,3e-4),(0.2,3e-4))

N=100
Xc = [0.01*randn(6) for _=1:N]
pts = getPoint.(Pose3, Xc)

# pre-emptively populate the measurements, kept separate since nlsolve calls fp(x, res) multiple times
# measurement = getSample(meas, N)

fg = initfg()
X0 = addVariable!(fg, :x0, Pose3)
initManual!(fg, :x0, pts)
X1 = addVariable!(fg, :x1, Point3)
addFactor!(fg, [:x0;:x1], meas, graphinit=false)

pts = approxConv(fg, :x0x1f1, :x1)

##

@warn "still need to insert kld(..) test to ensure this is working"

p1 = manikde!(TranslationGroup(3), pts);

println("Test back projection from ")

meas = LinearRangeBearingElevation((3.0,3e-4),(0.2,3e-4))

N = 100
L = [[3.0; 0.65; 0] for _=1:N]

fg = initfg()
X0 = addVariable!(fg, :x0, Pose3)
X1 = addVariable!(fg, :x1, Point3)
initManual!(fg, :x1, L)
addFactor!(fg, [:x0;:x1], meas, graphinit=false)

pts = approxConv(fg, :x0x1f1, :x0)

p2 = manikde!(SpecialEuclidean(3), pts);

end

# using Gadfly
#
# axis = [[1.5;3.5]';[-1.25;1.25]';[-1.0;1.0]']
# draw( PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
#       plotKDE( p1, dimLbls=["x";"y";"z"], axis=axis)  )
# #
#
# axis = [[-2.5;2.5]';[-2.5;2.5]';[-2.5;2.5]';[-2pi;2pi]';[-2pi;2pi]';[-2pi;2pi]']
# draw( PDF("/home/dehann/Desktop/test.pdf",30cm,20cm),
#       plotKDE( p2, dimLbls=["x";"y";"z"; "roll"; "pitch"; "yaw"], axis=axis)  )

##

@testset "LinearRangeBearingElevation in solve" begin

##

# do test directly in factor graph
fg = initfg()
# @show fg.registeredModuleFunctions
N = 100
@warn "Breaks if not set to 100"

initCov = Matrix{Float64}(LinearAlgebra.I, 6,6)
[initCov[i,i] = 0.01 for i in 4:6];
odoCov = deepcopy(initCov)

v1 = addVariable!(fg,:x1, Pose3,  N=N)
initPosePrior = PriorPose3( MvNormal(zeros(6), initCov) )
f1  = addFactor!(fg,[v1], initPosePrior)

# implemented in SensorModels.jl
meas = LinearRangeBearingElevation((3.0,3e-4),(0.2,3e-4))

@time X = getVal(v1)

v2 = addVariable!(fg, :l1, Point3, N=N)
f2 = addFactor!(fg, [:x1;:l1], meas) #, threadmodel=MultiThreaded)

# L1pts = approxConv(fg, :x1l1f1, :l1)
# X1pts = approxConv(fg, :x1l1f1, :x1)

solveTree!(fg);

##

end


text:
using LazySequences

# Cons
c = cons(1, [42])
# Test first(s::Cons)
@assert first(c) == 1
# Test rest(s::Cons)
@assert first(rest(c)) == 42

ct = cat([1], [42])
# Test first(s::Cat)
@assert first(ct) == 1
# Test rest(s::Cat)

# Test getindex implementation
fibs = cat([0, 1], @lazyseq map(+, rest(fibs), fibs))
@assert fibs[1] == 0
@assert fibs[2] == 1
@assert fibs[3] == 1
@assert fibs[4] == 2


text:
type Hemorrhage
    hID::Array{Int64,1}
    tID::Array{Int64,1}
    Amax::Array{Float64,1}
    Amin::Array{Float64,1}
    bmax::Array{Float64,1}
    bmin::Array{Float64,1}
    Vloss::Float64
    Vlossinit::Float64
    totalloss::Float64
    Ph::Float64
    ttn::Float64
    injured::Vector{Bool}

    function Hemorrhage()
        this = new()
        this.hID = [40]; # hepatic artery ID is 14, femoral artery is 40
        this.tID = [39,40]; # apply tourniquet to femoral & profundis arteries
        this.Amax = Array{Float64,1}[];
        this.Amin = Array{Float64,1}[];
        this.bmax = Array{Float64,1}[];
        this.bmin = Array{Float64,1}[];
        this.Vloss = 0;
        this.Vlossinit = 0;
        this.totalloss = 550*cm3Tom3;
        this.Ph = 0*mmHgToPa;
        this.ttn = 10;
        this.injured = Vector{String}[];
        return this
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
@testset "Perspective transformation" begin
    @test PerspectiveMap()([2.0, -1.0, 0.5]) ≈ [4.0, -2.0]

    if VERSION >= v"0.6.0-dev.948"
        @test cameramap(2.0)                  ≈ LinearMap(UniformScaling(2.0)) ∘ PerspectiveMap()
        @test cameramap((1.1,2.2))            ≈ LinearMap([1.1 0; 0 2.2]) ∘ PerspectiveMap()
        @test cameramap(1.1,       (3.3,4.4)) ≈ Translation([-3.3,-4.4]) ∘ LinearMap(UniformScaling(1.1)) ∘ PerspectiveMap()
        @test cameramap((1.1,2.2), (3.3,4.4)) ≈ Translation([-3.3,-4.4]) ∘ LinearMap([1.1 0; 0 2.2]) ∘ PerspectiveMap()
    else
        # `isapprox` for `UniformScaling` not defined. Just put a few points through and
        # test that we get the same results
        points = [[1.0,1.0,1.0], [-1.0,1.0,1.0], [-1.0,-1.0,1.0], [1.0,-1.0,1.0]]

        trans1 = cameramap(2.0)
        trans2 = LinearMap(UniformScaling(2.0)) ∘ PerspectiveMap()
        @test all(isapprox.(trans1.(points), trans2.(points)))

        trans1 = cameramap((1.1,2.2))
        trans2 = LinearMap([1.1 0; 0 2.2]) ∘ PerspectiveMap()
        @test all(isapprox.(trans1.(points), trans2.(points)))

        trans1 = cameramap(1.1,       (3.3,4.4))
        trans2 = Translation([-3.3,-4.4]) ∘ LinearMap(UniformScaling(1.1)) ∘ PerspectiveMap()
        @test all(isapprox.(trans1.(points), trans2.(points)))

        trans1 = cameramap((1.1,2.2), (3.3,4.4))
        trans2 = Translation([-3.3,-4.4]) ∘ LinearMap([1.1 0; 0 2.2]) ∘ PerspectiveMap()
        @test all(isapprox.(trans1.(points), trans2.(points)))

    end
end