text:
using MPI
using p4est_wrapper
using Test

el_num = Cint(8)
el_size = Cint(32)

# Initialize MPI if not initialized yet
if !MPI.Initialized()
    MPI.Init()
end

mpicomm = p4est_wrapper.P4EST_ENABLE_MPI ? MPI.COMM_WORLD : Cint(0)

# SC library initialization
sc_init(mpicomm, Cint(true), Cint(true), C_NULL, p4est_wrapper.SC_LP_DEFAULT)
@test typeof(sc_is_root()) == Cint

# SC package related API
id = sc_package_register(C_NULL, p4est_wrapper.SC_LP_INFO, "p4est", "A forest of octrees")
@test sc_package_is_registered(id) == 1
sc_package_lock(id)
sc_package_unlock(id)
sc_package_set_verbosity(id, p4est_wrapper.SC_LP_DEBUG)
sc_package_print_summary(id)

# SC memory allocation/deallocation related API
@test sc_memory_status(id) == 0 && sc_memory_check(id) == nothing
mem_ptr = sc_malloc(id, el_size)
@test sc_memory_status(id) == 1
mem_ptr = sc_realloc(id, mem_ptr, 2*el_size)
@test sc_memory_status(id) == 1
@test sc_free(id, mem_ptr) == nothing
@test sc_memory_status(id) == 0 && sc_memory_check(id) == nothing

@test sc_memory_status(id) == 0 && sc_memory_check(id) == nothing
mem_ptr = sc_calloc(id, el_num, el_size)
@test sc_memory_status(id) == 1
str_ptr = sc_strdup(id, "p4est")
@test sc_memory_status(id) == 2
@test sc_free(id, mem_ptr) == nothing
@test sc_free(id, Ptr{Cvoid}(str_ptr)) == nothing
@test sc_memory_status(id) == 0 && sc_memory_check(id) == nothing

# SC number comparison related API
@test sc_int_compare([Int(1)],        [Int(1)])     == 0
@test sc_int8_compare([Int8(2)],      [Int8(2)])    == 0
@test sc_int16_compare([Int16(3)],    [Int16(3)])   == 0
@test sc_int32_compare([Int32(4)],    [Int32(4)])   == 0
@test sc_int64_compare([Int64(5)],    [Int64(5)])   == 0
@test sc_double_compare([Float64(6)], [Float64(6)]) == 0

# SC_ARRAY related API
sc_array_ptr = sc_array_new(el_size)
@test sc_array_init(sc_array_ptr, el_size) == nothing
@test sc_array_memset(sc_array_ptr,0) == nothing
@test sc_array_reset(sc_array_ptr) == nothing
@test sc_array_truncate(sc_array_ptr) == nothing
@test sc_array_rewind(sc_array_ptr,2*el_num) == nothing
@test sc_array_resize(sc_array_ptr,2*el_num) == nothing
@test sc_array_destroy(sc_array_ptr) == nothing

sc_array_ptr = sc_array_new_count(el_size, el_num)
#@test sc_array_init_size(sc_array_ptr, el_size, el_num) == nothing
#@test sc_array_init_count(sc_array_ptr, el_size, el_num) == nothing
#sc_array_data_ptr = sc_array_new_data(sc_array_ptr, el_size, el_num)
#sc_array_view_ptr = sc_array_new_view(sc_array_ptr, 0, 0)
#@test sc_array_destroy(sc_array_data_ptr) == nothing
#@test sc_array_destroy(sc_array_view_ptr) == nothing
@test sc_array_destroy(sc_array_ptr) == nothing

# SC_MEMPOOL related API
sc_mempool_ptr = sc_mempool_new(el_size)
#@test sc_mempool_init(sc_mempool_ptr, el_size) == nothing
#@test sc_mempool_truncate(sc_mempool_ptr) == nothing
#@test sc_mempool_reset(sc_mempool_ptr) == nothing
#@test sc_mempool_destroy(sc_mempool_ptr) == nothing
#sc_mempool_ptr = sc_mempool_new_zero_and_persist(el_size)
@test sc_mempool_destroy(sc_mempool_ptr) == nothing

# SC library finalization
@test sc_package_unregister(id) == nothing
@test sc_finalize() == nothing

# Finalize MPI if initialized and session is not interactive
if (MPI.Initialized() && !isinteractive())
    MPI.Finalize()
end


text:
using Distances

export Normalised

struct Normalised{T<:Metric} <: Metric
    d::T
end

function (d_n::Normalised{T})(x,y) where {T<:Metric}
    d = d_n.d
    d_tmp = d(x,y)
    return 2 * d_tmp / (d(x,nothing) + d(y,nothing) + d_tmp)
end

(d_n::Normalised{T} where {T<:Metric})(x::Nothing, y::Nothing) = 0.0 # To avoid NaN

ext: jl


text:
#
using SpectralElements.Spectral

nr = 8
ns = 12

u = rand(nr,ns) |> Field
v = rand(nr,ns) |> Field

## DiagonalOp
d = rand(nr,ns) |> Field
D = DiagonalOp(d)

@test mul!(v,D,u) ≈ d .* u

## TensorProductOp2D

A = rand(nr,nr)
B = rand(ns,ns)
T = TensorProductOp2D(A, B)

@test mul!(v,T,u) ≈ Field(A * u.array * B')
@test mul!(v,T,u) ≈ Field(A * u.array * B')

## ZeroOp

## IdentityOp

## AffineOp

#### +,-(op,op)
#### +,-(λ, op)
#### *,/(λ, op)

## ComposeOp

# caching mechanism is broken. fix later.

## InverseOp


text:
"""
    inverse(f)

Return the inverse of function `f`.

This is an internal function to avoid type piracies such as `inv(::typeof(exp))`. At some
point `inv` might support standard Julia functions which would allow us to remove `inverse`.
"""
inverse(f)

inverse(::typeof(exp)) = log
inverse(::typeof(log)) = exp

inverse(::typeof(log1pexp)) = logexpm1
inverse(::typeof(logexpm1)) = log1pexp

inverse(::typeof(inv)) = inv

square(x) = x^2
inverse(::typeof(sqrt)) = square
inverse(::typeof(square)) = sqrt

inverse(::typeof(logit)) = logistic
inverse(::typeof(logistic)) = logit

inverse(::typeof(normcdf)) = norminvcdf
inverse(::typeof(norminvcdf)) = normcdf


text:
#= ------------------------------------------------------------------------

    Region

------------------------------------------------------------------------ =#

import Base: copy, -, union, ==, show
export Region
export isempty, ==, copy, invert, -, translate, translate!, contains, ∈
export left, top, right, bottom, bounds
export complement, complement!
export union, intersection, difference
export region_from_box
export region_to_image

"""
    Region

A region is a discrete set of coordinates in two-dimensional euclidean space.

A region consists of zero or more runs, which are sorted in ascending order.

```jldoctest reg
julia> using Regions

julia> Region() # create an empty region
Region(Run[], false)

julia> Region([Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)]) # create a region with 3 runs
Region(Run[Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)], false)

julia> Region([Run(y, x1:x2) for y=-10:10, x1=-10, x2=10]); # create a region with many runs using comprehension
```
"""
struct Region
    runs::Vector{Run}
    complement::Bool
end

Region(runs::Vector{Run}) = Region(runs, false)
Region() = Region(Run[], false)

"""
    isempty(x::Region)

Discover whether the region is empty.

```jldoctest
julia> using Regions

julia> isempty(Region(Run[]))
true

julia> isempty(Region([Run(2, 1:1)]))
false
```
"""
isempty(x::Region) = isempty(x.runs)

"""
    ==(a::Region, b::Region)

Equality operator for two regions. Two regions are equal, if both their runs and their
complement flags are equal.

```jldoctest
julia> using Regions

julia> a = Region([Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)]);

julia> b = Region([Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)]);

julia> a == b
true
```
"""
==(a::Region, b::Region) = a.runs == b.runs && a.complement == b.complement

"""
    copy(x::Region)

Create a copy of a region.

```jldoctest
julia> using Regions

julia> a = Region([Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)]);

julia> b = copy(a)
Region(Run[Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)], false)
```
"""
copy(x::Region) = Region(copy(x.runs), x.complement)

"""
    invert(x::Region)
    -(x::Region)

Invert a region. Inversion mirrors a region at the origin. A region is inverted
by inverting each of its runs. Since the runs of a region are sorted by their row and
column coordinates, the order of the runs is inversed as well.

```jldoctest
julia> using Regions

julia> a = Region([Run(0, 0:2), Run(1, 0:2), Run(2, 0:2)]);

julia> b = invert(a)
Region(Run[Run(-2, -2:0), Run(-1, -2:0), Run(0, -2:0)], false)
```
"""
function invert(x::Region)
    result = Region(Run[], x.complement)
    # iterating backwards maintains the correct sort order of the runs
    for i in length(x.runs):-1:1
        push!(result.runs, -x.runs[i])
    end
    return result
end
-(x::Region) = invert(x)

"""
    translate(r::Region, x::Integer, y::Integer)
    translate(r::Region, a::Vector{Int64})

Translate a region. Translation moves a region. A region is translated by translating each
of its runs.

```jldoctest
julia> using Regions

julia> a = Region([Run(0, 0:2), Run(1, 0:2), Run(2, 0:2)]);

julia> b = translate(a, -1, -1)
Region(Run[Run(-1, -1:1), Run(0, -1:1), Run(1, -1:1)], false)
```
"""
function translate(r::Region, x::Integer, y::Integer)
    return translate!(copy(r), x, y)
end
translate(r::Region, d::Vector{Int64}) = translate(r, d[1], d[2])
+(x::Region, y::Vector{Int64}) = translate(x, y[1], y[2])
+(x::Vector{Int64}, y::Region) = translate(y, x[1], x[2])
-(x::Region, y::Vector{Int64}) = translate(x, -y[1], -y[2])
function translate!(r::Region, x::Integer, y::Integer)
    for i in 1:length(r.runs)
        r.runs[i] = translate(r.runs[i], x, y)
    end
    return r
end
translate!(r::Region, d::Vector{Int64}) = translate!(r, d[1], d[2])

"""
    contains(r::Region, x::Integer, y::Integer)
    contains(r::Region, a::Array{Int64, 1})

Test if region r contains position (x, y).
"""
contains(r::Region, x::Integer, y::Integer) = r.complement ? !any(run -> contains(run, x, y), r.runs) : any(run -> contains(run, x, y), r.runs)
contains(r::Region, a::Array{Int64, 1}) = contains(r, a[1], a[2])
∈(a::Array{Int64, 1}, r::Region) = contains(r, a)

"""
    left(x::Region)

Calculates the leftmost region coordinate.

This function works for non-complement and non-empty regions only.
"""
function left(r::Region)
    @assert !r.complement "cannot calculate for infinite (complement) regions"
    @assert length(r.runs)>0 "cannot calculate for empty regions"
    v = r.runs[1].columns.start
    for i in r.runs
        v = min(v, i.columns.start)
    end
    return v
end

"""
    top(x::Region)

Calculates the topmost region coordinate.

This function works for non-complement and non-empty regions only.
"""
function top(r::Region)
    @assert !r.complement "cannot calculate for infinite (complement) regions"
    @assert length(r.runs)>0 "cannot calculate for empty regions"
    v = r.runs[1].row
    for i in r.runs
        v = max(v, i.row)
    end
    return v
end

"""
    right(x::Region)

Calculates the rightmost region coordinate.

This function works for non-complement and non-empty regions only.
"""
function right(r::Region)
    @assert !r.complement "cannot calculate for infinite (complement) regions"
    @assert length(r.runs)>0 "cannot calculate for empty regions"
    v = r.runs[1].columns.stop
    for i in r.runs
        v = max(v, i.columns.stop)
    end
    return v
end

"""
    bottom(x::Region)

Calculates the bottommost region coordinate.

This function works for non-complement and non-empty regions only.
"""
function bottom(r::Region)
    @assert !r.complement "cannot calculate for infinite (complement) regions"
    @assert length(r.runs)>0 "cannot calculate for empty regions"
    v = r.runs[1].row
    for i in r.runs
        v = min(v, i.row)
    end
    return v
end

"""
    bounds(x::Region)

Calculates the left, top, right and bottom region coordinate and returns them
as a tuple.

This function works for non-complement and non-empty regions only.
"""
function bounds(region::Region)
    @assert !region.complement "cannot calculate for infinite (complement) regions"
    @assert length(region.runs)>0 "cannot calculate for empty regions"
    l = region.runs[1].columns.start
    t = region.runs[1].row
    r = region.runs[1].columns.stop
    b = region.runs[1].row
    for run in region.runs
        l = min(l, run.columns.start)
        t = max(t, run.row)
        r = max(r, run.columns.stop)
        b = min(b, run.row)
    end
    return l, t, r, b
end

"""
    complement(x::Region)

Calculates the set-theoretic complement of a region.

With a non-complemented region, the runs specify the contained pixels, i.e. they specify
what is included within the region. With a complemented region, the runs specify the
non-contained pixels, i.e. they specify what is not included within the region.
"""
complement(x::Region) = Region(copy(x.runs), !x.complement)
complement!(x::Region) = x.complement =! x.complement

"""
    merge(a::Vector{Run}, b::Vector{Run})

Merge sorted vectors `a` and `b`. Assumes that `a` and `b` are sorted
and does not check whether `a` or `b` are sorted.

merge is not exported, since its basic usage is within this file and it conflicts with
a definition in Base.
"""
function merge(a::Vector{Run}, b::Vector{Run})
    res = Run[]
    i = 1
    j = 1
    while i <= length(a)
        if j > length(b)
            while i <= length(a)
                push!(res, a[i])
                i = i+1
            end
            return res
        end
        if (b[j] < a[i])
            push!(res, b[j])
            j = j+1
        else
            push!(res, a[i])
            i = i+1
        end
    end
    while j <= length(b)
        push!(res, b[j])
        j = j+1
    end
    return res
end

"""
    sort(a::Vector{Run})

Variant of sort that returns a sorted copy of a leaving a itself unmodified. This ensures that
runs are sorted after an operation that might have destroyed the sort order, such as
downsampling.
"""
sort(a::Vector{Run}) = sort(a)

"""
    sort!(a::Vector{Run})

Sort the vector a in place. This ensures that runs are sorted after an operation that might
have destroyed the sort order, such as downsampling.
"""
sort!(a::Vector{Run}) = sort!(a)

"""
    pack!(a::Vector{Run})

Packs runs together. pack! is not exported, since its basic usage is within this file as a
building block for union.
"""
function pack!(a::Vector{Run})
    read = 1
    write = 1

    while read <= length(a)
        a[write] = a[read]
        read += 1

        while read <= length(a) && a[write].row == a[read].row && a[write].columns.stop + 1 >= a[read].columns.start
            if a[read].columns.stop > a[write].columns.stop
                a[write] = Run(a[read].row, a[write].columns.start:a[read].columns.stop)
            end
            read += 1
        end
        write += 1
    end
    deleteat!(a, write:length(a))
end

"""
    union(a::Vector{Run}, b::Vector{Run})

Calculates the union of two sorted arrays of runs. The function assumes that the runs are sorted
but does not check this.
"""
function union(a::Vector{Run}, b::Vector{Run})
    res = merge(a, b)
    pack!(res)
    return res
end

"""
    union(a::Region, b::Region)

Calculates the union of two regions. This function supports complement regions and uses
DeMorgan's rules to eliminate the complement.
"""
function union(a::Region, b::Region)
    if a.complement && b.complement
        return Region(intersection(a.runs, b.runs), true)
    elseif a.complement
        return Region(difference(a.runs, b.runs), true)
    elseif b.complement
        return Region(difference(b.runs, a.runs), true)
    else
        return Region(union(a.runs, b.runs), false)
    end
end

"""
    intersect!(a::Vector{Run})

Intersects runs. intersect! is not exported, since its basic usage is within this file as a
building block for intersection.
"""
function intersect!(a::Vector{Run})
    read = 1
    if read > length(a)
        return
    end
    next = read + 1
    write = read
    while next <= length(a)
        if a[read].row != a[next].row
            read = next
            next += 1
        else
            if a[next].columns.start > a[read].columns.stop
                read = next
                next += 1
            else
                a[write] = Run(a[read].row, a[next].columns.start:min(a[read].columns.stop, a[next].columns.stop))
                if a[next].columns.stop < a[read].columns.stop
                    a[next] = a[read]
                end
                read = next
                next += 1
                write += 1
            end
        end
    end
    deleteat!(a, write:length(a))
end

"""
    intersection(a::Vector{Run}, b::Vector{Run})

Calculates the intersection of two sorted arrays of runs. The function assumes that the runs are sorted
but does not check this.
"""
function intersection(a::Vector{Run}, b::Vector{Run})
    res = merge(a, b)
    intersect!(res)
    return res
end

"""
    intersection(a::Region, b::Region)

Calculates the intersection of two regions. This function supports complement regions and uses
DeMorgan's rules to eliminate the complement.
"""
function intersection(a::Region, b::Region)
    if a.complement && b.complement
        return Region(union(a.runs, b.runs), true)
    elseif a.complement
        return Region(difference(b.runs, a.runs), false)
    elseif b.complement
        return Region(difference(a.runs, b.runs), false)
    else
        return Region(intersection(a.runs, b.runs), false)
    end
end

"""
    difference(a::Vector{Run}, b::Vector{Run})

Calculates the difference of two sorted vectors of runs. The function assumes that the runs are sorted
but does not check this.
"""
function difference(a::Vector{Run}, b::Vector{Run})
    if length(a) == 0
        return Run[]
    end

    if length(b) == 0
        return a
    end

    res = Run[]

    # first_b and last_b form a range of chords in b that are in the same row
    first_b = findfirst(x -> x.row >= a[1].row, b)
    last_b = first_b
    if first_b != nothing
        last_b = findlast(x -> x.row == b[first_b].row, view(b, first_b:length(b)))
    end

    for a_index in 1:length(a)
        if first_b != nothing && a[a_index].row > b[first_b].row
            # update the range
            first_b = findfirst(x -> x.row >= a[a_index].row, b)
            last_b = first_b
            if first_b != nothing
                last_b = findlast(x -> x.row == b[first_b].row, view(b, first_b:length(b)))
            end
        end
        if first_b == nothing || a[a_index].row != b[first_b].row
            push!(res, a[a_index])
        else
            for i in first_b:last_b
                if isoverlapping(a[a_index], b[i])
                    # total overlap, erase all of a_chord and break
                    if b[i].columns.start <= a[a_index].columns.start && b[i].columns.stop >= a[a_index].columns.stop
                        a[a_index] = Run(a[a_index].row, a[a_index].columns.start:a[a_index].columns.start-1)
                        break;
                    # overlap at left only, shorten a_chord at left and continue
                    elseif b[i].columns.start <= a[a_index].columns.start
                        a[a_index] = Run(a[a_index].row, b[i].columns.stop+1:a[a_index].columns.stop)
                    # overlap at right only, shorten a_chord at right and continue
                    elseif b[i].columns.stop >= a[a_index].columns.stop
                        a[a_index] = Run(a[a_index].row, a[a_index].columns.start:b[i].columns.start-1)
                    # overlap in the middle, split a_chord into two and continue
                    else
                        push!(res, Run(a[a_index].row, a[a_index].columns.start:b[i].columns.start-1))
                        a[a_index] = Run(a[a_index].row, b[i].columns.stop+1:a[a_index].columns.stop)
                    end
                end
            end
            if !isempty(a[a_index])
                push!(res, a[a_index])
            end
        end
    end

    return res
end

"""
    difference(a::Region, b::Region)

Calculates the union of two regions. This function supports complement regions and uses
DeMorgan's rules to eliminate the complement.
"""
function difference(a::Region, b::Region)
    if a.complement && b.complement
        return Region(difference(b.runs, a.runs), false)
    elseif a.complement
        return Region(union(a.runs, b.runs), true)
    elseif b.complement
        return Region(intersection(a.runs, b.runs), false)
    else
        return Region(difference(a.runs, b.runs), false)
    end
end

"""
    region_from_box(left::Integer, top::Integer, right::Integer, bottom::Integer)

Create a region given box coordinates. The region consists of all coordinates within the box.
"""
function region_from_box(left::Integer, top::Integer, right::Integer, bottom::Integer)
    @assert bottom < top "bottom must be smaller than top"
    @assert left < right "left must be smaller than right"

    region = Region(Run[])
    for row in bottom:top
        push!(region.runs, Run(row, left:right))
    end
    return region
end

"""
    region_to_image(r::Region, color=Gray(true))

Converts a region to an image. The function determines the bounds of the
region and then renders the region into the image.

The background of the image is filled with zeroes, the region pixels are
colored with the passed in color.

Some examples of colors that you can pass:
Gray(0.5) : mid gray value
RGB(1, 0, 0) : brightest red value
RGBA(0, 0.5, 0, 0.5) : half transparent mid green value
"""
function region_to_image(region::Region, color=Gray(true))
    (l, t, r, b) = bounds(region)
    img = zeros(typeof(color), t-b+1, r-l+1)
    for run in region.runs
        for column in run.columns
            img[run.row-b+1, column-l+1] = color
        end
    end
    return img
end

"""
    Base.show(io, mime::MIME"image/png", r::Region)

Shows a rich graphical display of a region. The region is displayed in a half
transparent blue color.
"""
function Base.show(io::IO, mime::MIME"image/png", region::Region)
    if !isempty(region)
        Base.show(io, mime, region_to_image(region, RGBA(0,0,1,0.5)))
    end
end


text:
"note that this uses the non-dimensional X = k*depth"
function discrete_wave_system(ω::T, X::AbstractVector{T}, source::PlaneSource{T,2,1,Acoustic{T,2}}, material::Material{2,Halfspace{T,2}};
        tol::T = 1e-6,
        scheme::Symbol = :trapezoidal,
        basis_order::Int = 2,
        kws...
    ) where T<:AbstractFloat

    specie = material.species[1]
    t_vec = t_matrix(specie, source.medium, ω, basis_order)

    k = real(ω / source.medium.c)
    a12k = specie.exclusion_distance * T(2)*real(k * outer_radius(specie));
    M = basis_order;

    θin = transmission_angle(source,material)

    J = length(X) - 1
    h = X[2] - X[1]

    PQ_quad = intergrand_kernel(X, a12k; M = M, θin = θin, scheme=scheme);

    MM_quad = [
        (number_density(specie) / (k^2))*t_vec[m+M+1,m+M+1]*PQ_quad[l,m+M+1,j,n+M+1] - ( (m==n && j==l) ? 1.0+0.0im : 0.0+0.0im)
    for  l=1:(J+1), m=-M:M, j=1:(J+1), n=-M:M];

    b_mat = [
        -t_vec[m+M+1,m+M+1]*exp(im*X[l]*cos(θin))*exp(im*m*(pi/2.0 - θin))
    for l = 1:(J+1), m = -M:M]

    return (MM_quad,b_mat)
end

"Returns x the mesh used to discretise the integral equations."
function x_mesh(wave_eff_long::EffectivePlaneWaveMode{T}, wave_eff_short::EffectivePlaneWaveMode{T} = wave_eff_long;
        tol::T = T(1e-5),  a12::T = T(Inf),
        max_size::Int = 1000,
        min_size::Int = 5,
        max_x::T = -log(tol) / abs(imag(wave_eff_long.wavenumber))
        # max_x::T = -log(tol) / abs(imag(wave_eff_short.wavenumber))
        # max_x::T = (-log(tol))/abs(cos(wave_eff_long.θ_eff)*imag(wave_eff_long.wavenumber))
) where T<:AbstractFloat

    #= The default max_x results in:
        abs(exp(im*max_x*cos(θ_effs[end])*k_effs[end]/k)) < tol
    =#
    # estimate a reasonable derivative based on more rapidly varying wave_eff_short.
    df = abs(wave_eff_short.wavenumber)
    # NOTE was: df = abs(wave_eff_short.wavenumber * cos(wave_eff_short.θ_eff))

    # Based on Simpson's rule
        # dX  = (tol*90 / (df^4))^(1/5)
    # Based on trapezoidal integration
        dx  = (tol * T(30) / (df^2))^(1/3)

    # prioritise a mesh fine enough to give accurate discrete integrals
    if dx > a12 dx = a12 end

    if max_x/dx + 1 > max_size
        if dx == a12
            max_x = (max_size - 1)*dx
            @warn("The mesh max_size = $max_size which was too small for tol = $tol. Will shrink meshed region.")
        else # otherwise priortise shrink max_x and increase dx proportionally
            a = sqrt(max_x/((max_size - 1)*dx))
            dx = min(dx*a,a12)
            max_x = (max_size - 1)*dx
            @warn("The mesh max_size = $max_size which was too small for tol = $tol. Will make a smaller and coarser mesh.")
        end
    elseif max_x/dx + 1 < min_size
        max_x = (min_size - 1)*dx
    end

    # if whole correction length a12k was given, then make dX/a12k = integer
    if a12 != T(Inf) && dx < a12
        n = ceil(a12 / dx)
        dx = a12/n
    end

    return zero(T):dx:max_x
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
    Shift(N, [k=1])

Shift operator. Optionally specify the distance to shift (default: 1).

# Fields
- `N::Int`: total number of spins
- `amount::Int`: shift by how many places
"""
struct Shift <: AbstractSymmetry
    N::Int
    amount::Int
    Shift(N, k=1) = new(N, mod(k, N))
end

_order(s::Shift) = div(s.N, gcd(s.N, s.amount))

_mask(k) = (2^k)-1 # mask out the k lowest bits
_roll_bits(N, i, amount) = ((i & _mask(amount)) << (N-amount)) | (i >> amount)
(s::Shift)(index) = _roll_bits(s.N, index, s.amount)

"""
    Flip(positions)
    Flip(system_size)

Spin Flip operator on N spins. Either provide the specific locations to flip or the system's
size to flip them all.

# Fields
- `mask::Int`: encodes the positions to flip in binary
"""
struct Flip <: AbstractSymmetry
    mask::Int
    """
        Flip(positions)
        Flip(system_size)

    Flip the spins at the positions or all if the system's size is given.
    """
    Flip(positions) = new(_compute_mask(positions))
end
Flip(system_size::Int) = Flip(1:system_size)
_compute_mask(inds) = sum(1 .<< (inds .- 1))

_order(::Flip) = 2

(f::Flip)(index)  = xor(index, f.mask)

"""
    Swap(pos1, pos2)

Swaps the 2 specified spins.

# Fields
- `ind1::Int`: Position of the first spin to flip
- `ind2::Int`: Position of the second spin to flip
Positions should be given in the range 1:N.
"""
struct Swap <: AbstractSymmetry
    ind1::Int
    ind2::Int
    function Swap(pos1, pos2)
        pos1 == pos2 && @warn "Swap with identical positions is nonsensical and won't work properly!"
        new(pos1, pos2)
    end
end

_order(::Swap) = 2

function _swap_bits(x, ind1, ind2)
    # compare bits, swap if unequal
    sw = ((x >> ind1) & 1) ⊻ ((x >> ind2) & 1) # xor == 1 if unequal
    return xor(x, (sw << ind1) | (sw << ind2)) # xor(bit, 1) -> flips bit
end

# convert spin positions to bitstring position
(s::Swap)(index) = _swap_bits(index, s.ind1-1, s.ind2-1)

"""
    SpatialReflection(N)

Reflects the whole chain in space.

# Fields
- `N::Int`: Number of spins
"""
struct SpatialReflection <: AbstractSymmetry
    N::Int
    function SpatialReflection(N)
        N == 1 && @warn "SpatialReflection with N=1 is nonsensical and won't work properly!"
        new(N)
    end
end

_order(::SpatialReflection) = 2

# approx 10times faster than the simpler
# parse(Int, string(x; base=2, pad=N)[N:-1:1]; base=2)
# Also note that the former does not account for cases where there are actually more than N
# spins
function _reflect_bits(N, bits)
    for i in 0:div(N,2)-1
        bits = _swap_bits(bits, i, N-1-i)
    end
    return bits
end

(sr::SpatialReflection)(index) = _reflect_bits(sr.N, index)

"""
    GenericSymmetry(f, L)

Allows a user-defined symmetry. Just provide a function to transform a binary basis index
and the order of the symmetry.

# Fields
- `f::Function`: transformation of the indices given in binary representation (0:2^N-1)
- `L::Int`: order of the symmetry, i.e. the smallest positive number such that f^L = identity
"""
struct GenericSymmetry <: AbstractSymmetry
    f::Function
    L::Int
end

_order(gs::GenericSymmetry) = gs.L

(gs::GenericSymmetry)(index) = gs.f(index)


text:
abstract type Continuity{T} end

function raviartthomas(mesh, ::Type{Continuity{:none}})

    @assert dimension(mesh) == 2

    P = vertextype(mesh)
    S = Shape{coordtype(mesh)}
    F = Vector{S}

    nf = 3*numcells(mesh)
    functions = Vector{F}(undef, nf)
    positions = Vector{P}(undef, nf)

    for i in 1:numcells(mesh)
        functions[3*(i-1)+1] = [S(i,1,+1.0)]
        functions[3*(i-1)+2] = [S(i,2,+1.0)]
        functions[3*(i-1)+3] = [S(i,3,+1.0)]

        ctr = sum(mesh.vertices[mesh.faces[i]])/3
        positions[3*(i-1)+1] = ctr
        positions[3*(i-1)+2] = ctr
        positions[3*(i-1)+3] = ctr
    end

    return RTBasis(mesh, functions, positions)
end


text:
c = ingleton(4,1,2,3,4)
cut = nonnegative(4, 1234)
n = 4
h = polymatroidcone(Float64, 4)
newcut = :AddImmediately

using CutPruners
using StructDualDynProg

@testset "Ingleton score" begin
    @testset "n = $n, m = $m, nits=$nits, objval=$objval" for (n, m, nnodes, nits, objval, ncuts, ncuts_nonred) in ((4, 1, 1, 1, -1/4,  0, 0), # Shannon cone without any adhesivity
                                                                                (5, 1, 4, 1, -1/4,  0, 0), # Adhesivity adding up to n=5 but paths of length m=1 so they are not encountered
                                                                                (5, 2, 4, 2, -1/6, 24, 4)) # n=5, m=2 so they are encountered and we have 6 new cuts that improves the Ingleton bound to -1/6
        (sp, allnodes) = StructDualDynProg.SOI.stochasticprogram(m, c, h, lp_solver, n, cut, newcut,
                                                                 AvgCutPruningAlgo.([-1,-1,-1,-1,-1,-1,-1]))
        algo = StructDualDynProg.SDDP.Algorithm(K = -1)
        info = StructDualDynProg.SOI.optimize!(sp, algo, StructDualDynProg.SOI.IterLimit(nits))
        @test length(allnodes) == nnodes
        res = StructDualDynProg.SOI.last_result(info)
        @test res.status == :Optimal
        @test res.lowerbound ≈ objval
        @test res.upperbound ≈ objval
        @test CutPruners.ncuts(sp.data[1].nlds.FCpruner) == ncuts
        CutPruners.exactpruning!(sp.data[1].nlds.FCpruner, lp_solver)
        @test CutPruners.ncuts(sp.data[1].nlds.FCpruner) == ncuts_nonred
    end
end


text:
export fixed_field, subfields

add_verbose_scope(:Subfields)

# Compute basis for the subfield of K that is generated by the elements of as.
function _subfield_basis(K::S, as::Vector{T}) where {
    S <: Union{AnticNumberField, Hecke.NfRel},
    T <: Union{nf_elem, Hecke.NfRelElem}
   }
  if isempty(as)
    return [gen(K)]
  end

  # Notation: k base field, K the ambient field, F the field generated by as

  k = base_field(K)

  d = degree(K)
  Kvs = VectorSpace(k,d)
  # We transition the coefficients of a in reverse order, so that the
  # first vector in the row reduced echelon form yields the highest
  # degree among all elements of Fas.
  (Fvs,phivs) = sub(Kvs, [Kvs([coeff(a,n) for n in d-1:-1:0])
                          for a in as])
  dF = length(Fvs.gens) # dim(Fvs)
  bs = as
  while !isempty(bs)
    nbs = elem_type(K)[]
    for b in bs
      abs = [a*b for a in as]
      abvs,_ = sub(Kvs, [Kvs([coeff(ab,n) for n in d-1:-1:0])
                         for ab in abs])
      (Fvs,phivs) = sub(Kvs, typeof(Fvs)[Fvs, abvs])
      if dF != length(Fvs.gens) # dim(Fvs)
        dF = length(Fvs.gens) # dim(Fvs)
        append!(nbs, abs)
      end
    end
    bs = nbs
  end

  kx = parent(K.pol)
  return [let Kv = phivs(v)
            K(kx([Kv[n] for n in d:-1:1]))
          end
          for v in gens(Fvs)]
end

function _improve_subfield_basis(K, bas)
  # First compute the maximal order of <bas> by intersecting and saturating
  # Then B_Ok = N * B_LLL_OK
  # Then B' defined as lllN * B_LLL_OK will hopefully be small
  OK = maximal_order(K)
  OKbmatinv = basis_mat_inv(OK, copy = false)
  basinOK = bas * fmpq_mat(OKbmatinv.num) * fmpq(1, OKbmatinv.den)
  deno = fmpz(1)
  for i in 1:nrows(basinOK)
    for j in 1:ncols(basinOK)
      deno = lcm(deno, denominator(basinOK[i, j]))
    end
  end
  S = saturate(map_entries(FlintZZ, basinOK * deno))
  SS = S * basis_matrix(OK, copy = false)
  lllOK = lll(OK)
  N = (SS * basis_mat_inv(lllOK)).num
  lllN = lll(N)
  maybesmaller = lllN * basis_matrix(lllOK)
  return maybesmaller
end

function _improve_subfield_basis_no_lll(K, bas)
  OK = maximal_order(K)
  OKbmatinv = basis_mat_inv(OK, copy = false)
  basinOK = bas * fmpq_mat(OKbmatinv.num) * fmpq(1, OKbmatinv.den)
  deno = fmpz(1)
  for i in 1:nrows(basinOK)
    for j in 1:ncols(basinOK)
      deno = lcm(deno, denominator(basinOK[i, j]))
    end
  end
  S = saturate(map_entries(FlintZZ, basinOK * deno))
  SS = S * basis_matrix(OK, copy = false)
  return SS
end

# Compute a primitive element given a basis of a subfield
function _subfield_primitive_element_from_basis(K::S, as::Vector{T}) where {
    S <: Union{AnticNumberField, Hecke.NfRel},
    T <: Union{nf_elem, Hecke.NfRelElem}
   }
  if isempty(as)
    return gen(K)
  end

  d = length(as)

  # First check basis elements
  i = findfirst(a -> degree(minpoly(a)) == d, as)
  if i <= d
    return as[i]
  end

  k = base_field(K)

  # Notation: cs the coefficients in a linear combination of the as, ca the dot
  # product of these vectors.
  cs = fmpz[zero(ZZ) for n in 1:d]
  cs[1] = one(ZZ)
  while true
    ca = sum(c*a for (c,a) in zip(cs,as))
    if degree(minpoly(ca)) == d
      return ca
    end

    # increment the components of cs
    cs[1] += 1
    let i = 2
      while i <= d && cs[i-1] > cs[i]+1
        cs[i-1] = zero(ZZ)
        cs[i] += 1
        i += 1
      end
    end
  end
end

#As above, but for AnticNumberField type
#In this case, we can use block system to find if an element is primitive.
function _subfield_primitive_element_from_basis(K::AnticNumberField, as::Vector{nf_elem})
  if isempty(as)
    return gen(K)
  end

  dsubfield = length(as)

  @vprint :Subfields 1 "Sieving for primitive elements\n"
  # First check basis elements
  Zx = PolynomialRing(FlintZZ, "x", cached = false)[1]
  f = Zx(K.pol*denominator(K.pol))
  p, d = _find_prime(f)
  #First, we search for elements that are primitive using block systems
  F = FlintFiniteField(p, d, "w", cached = false)[1]
  Ft = PolynomialRing(F, "t", cached = false)[1]
  ap = zero(Ft)
  fit!(ap, degree(K)+1)
  rt = roots(f, F)
  indices = Int[]
  for i = 1:length(as)
    b = _block(as[i], rt, ap)
    if length(b) == dsubfield
      push!(indices, i)
    end
  end

  @vprint :Subfields 1 "Found $(length(indices)) primitive elements in the basis\n"
  #Now, we select the one of smallest T2 norm
  if !isempty(indices)
    a = as[indices[1]]
    I = t2(a)
    for i = 2:length(indices)
      t2n = t2(as[indices[i]])
      if t2n < I
        a = as[indices[i]]
        I = t2n
      end
    end
    @vprint :Subfields 1 "Primitive element found\n"
    return a
  end

  @vprint :Subfields 1 "Trying combinations of elements in the basis\n"
  # Notation: cs the coefficients in a linear combination of the as, ca the dot
  # product of these vectors.
  cs = fmpz[rand(FlintZZ, -2:2) for n in 1:dsubfield]
  k = 0
  s = 1
  first = true
  a = one(K)
  I = t2(a)
  while true
    s += 1
    ca = sum(c*a for (c,a) in zip(cs,as))
    b = _block(ca, rt, ap)
    if length(b) == dsubfield
      t2n = t2(ca)
      if first
        a = ca
        I = t2n
        first = false
      elseif t2n < I
        a = ca
        I = t2n
      end
      k += 1
      if k == 5
        @vprint :Subfields 1 "Primitive element found\n"
        return a
      end
    end

    # increment the components of cs
    bb = div(s, 10)+1
    for n = 1:dsubfield
      cs[n] = rand(FlintZZ, -bb:bb)
    end
    #=
    cs[1] += 1
    i = 2
    while i <= dsubfield && cs[i-1] > cs[i]+1
      cs[i-1] = zero(ZZ)
      cs[i] += 1
      i += 1
    end
    =#
  end
end

################################################################################
#
#  Subfield
#
################################################################################

@doc Markdown.doc"""
    subfield(L::NumField, elt::Vector{<: NumFieldelem};
                          isbasis::Bool = false) -> NumField, Map

The simple number field $k$ generated by the elements of `elt` over the base
field $K$ of $L$ together with the embedding $k \to L$.

If isbasis is true, it is assumed that `elt` holds a $K$-basis of $k$.
"""
function subfield(K::NumField, elt::Vector{<:NumFieldElem}; isbasis::Bool = false)
  if length(elt) == 1
    return _subfield_from_primitive_element(K, elt[1])
  end

  if isbasis
    s = _subfield_primitive_element_from_basis(K, elt)
  else
    bas = _subfield_basis(K, elt)
    s = _subfield_primitive_element_from_basis(K, bas)
  end

  return _subfield_from_primitive_element(K, s)
end

function _subfield_from_primitive_element(K::AnticNumberField, s::nf_elem)
  Qx = QQ["x"][1]
  if ismaximal_order_known(K)
    OK = maximal_order(K)
    @vtime :Subfields 1 f = Qx(minpoly(representation_matrix(OK(s, false))))
  else
    @vtime :Subfields 1 f = minpoly(Qx, s)
  end
  L, _ = NumberField(f, cached = false)
  return L, hom(L, K, s, check = false)
end

function _subfield_from_primitive_element(K, s)
  @vtime :Subfields 1 f = minpoly(s)
  L, _ = NumberField(f, cached = false)
  return L, hom(L, K, s, check = false)
end

################################################################################
#
#  Fixed field
#
################################################################################

@doc Markdown.doc"""
    fixed_field(K::SimpleNumField,
                sigma::Map;
                simplify::Bool = true) -> NumberField, NfToNfMor

Given a number field $K$ and an automorphisms $\sigma$ of $K$, this function
returns the fixed field of $\sigma$ as a pair $(L, i)$ consisting of a number
field $L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::SimpleNumField, sigma::T; simplify::Bool = true) where {T <: Union{NfToNfMor, NfRelToNfRelMor}}
  return fixed_field(K, T[sigma], simplify = simplify)
end

@doc Markdown.doc"""
    fixed_field(K::SimpleNumField, A::Vector{NfToNfMor}) -> NumberField, NfToNfMor

Given a number field $K$ and a set $A$ of automorphisms of $K$, this function
returns the fixed field of $A$ as a pair $(L, i)$ consisting of a number field
$L$ and an embedding of $L$ into $K$.

By default, the function tries to find a small defining polynomial of $L$. This
can be disabled by setting `simplify = false`.
"""
function fixed_field(K::AnticNumberField, A::Vector{NfToNfMor}; simplify::Bool = true)

  autos = small_generating_set(A)
  if length(autos) == 0
    return K, id_hom(K)
  end

  if ismaximal_order_known(K)
    OK = maximal_order(K)
    if isdefined(OK, :lllO)
      k, mk = fixed_field1(K, A)
      return k, mk
    end
  end

  a = gen(K)
  n = degree(K)
  ar_mat = Vector{fmpq_mat}()
  v = Vector{nf_elem}(undef, n)
  for i in 1:length(autos)
    domain(autos[i]) !== codomain(autos[i]) && throw(error("Maps must be automorphisms"))
    domain(autos[i]) !== K && throw(error("Maps must be automorphisms of K"))
    o = one(K)
    # Compute the image of the basis 1,a,...,a^(n - 1) under autos[i] and write
    # the coordinates in a matrix. This is the matrix of autos[i] with respect
    # to 1,a,...a^(n - 1).
    as = autos[i](a)
    if a == as
      continue
    end
    v[1] = o
    for j in 2:n
      o = o * as
      v[j] = o
    end
    bm = basis_matrix(v, FakeFmpqMat)
    # We have to be a bit careful (clever) since in the absolute case the
    # basis matrix is a FakeFmpqMat

    m = fmpq_mat(bm.num)
    for j in 1:n
      m[j, j] = m[j, j] - bm.den # This is autos[i] - identity
    end

    push!(ar_mat, m)
  end

  if length(ar_mat) == 0
    return K, id_hom(K)
  else
    bigmatrix = hcat(ar_mat)
    k, Ker = kernel(bigmatrix, side = :left)
    bas = Vector{elem_type(K)}(undef, k)
    if simplify
      KasFMat = _improve_subfield_basis(K, Ker)
      for i in 1:k
        bas[i] = elem_from_mat_row(K, KasFMat.num, i, KasFMat.den)
      end
    else
      #KasFMat = _improve_subfield_basis_no_lll(K, Ker)
      KasFMat = FakeFmpqMat(Ker)
      Ksat = saturate(KasFMat.num)
      Ksat = lll(Ksat)
      onee = one(fmpz)
      for i in 1:k
        #bas[i] = elem_from_mat_row(K, KasFMat.num, i, KasFMat.den)
        bas[i] = elem_from_mat_row(K, Ksat, i, onee)
      end
    end
  end
  return subfield(K, bas, isbasis = true)
end

function fixed_field(K::NfRel, A::Vector{T}; simplify::Bool = true) where {T <: NfRelToNfRelMor}
  autos = A

  # Everything is fixed by nothing :)
  if length(autos) == 0
    return K, id_hom(K)
  end

  F = base_field(K)
  a = gen(K)
  n = degree(K)
  ar_mat = Vector{dense_matrix_type(elem_type(F))}()
  v = Vector{elem_type(K)}(undef, n)
  for i in 1:length(autos)
    domain(autos[i]) !== codomain(autos[i]) && throw(error("Maps must be automorphisms"))
    domain(autos[i]) !== K && throw(error("Maps must be automorphisms of K"))
    o = one(K)
    # Compute the image of the basis 1,a,...,a^(n - 1) under autos[i] and write
    # the coordinates in a matrix. This is the matrix of autos[i] with respect
    # to 1,a,...a^(n - 1).
    as = autos[i](a)
    if a == as
      continue
    end
    v[1] = o
    for j in 2:n
      o = o * as
      v[j] = o
    end

    bm = basis_matrix(v)
    # In the generic case just subtract the identity
    m = bm - identity_matrix(F, degree(K))
    push!(ar_mat, m)
  end

  if length(ar_mat) == 0
    return K, id_hom(K)
  else
    bigmatrix = hcat(ar_mat)
    k, Ker = kernel(bigmatrix, side = :left)
    bas = Vector{elem_type(K)}(undef, k)
    for i in 1:k
      bas[i] = elem_from_mat_row(K, Ker, i)
    end
  end
  return subfield(K, bas, isbasis = true)
end

function fixed_field1(K::AnticNumberField, auts::Vector{NfToNfMor})

	auts_new = small_generating_set(auts)
  orderG = _order(auts)
  degree_subfield = divexact(degree(K), orderG)
  #TODO: Experiments to see if this is helpful
  #=
  if length(auts_new) == 1 && isprime_power(degree_subfield)
    #In this case, one of the coefficients of the minpoly of gen(K)
    #over the subfield is a generator for the subfield.
    #if the given generator was not too large, also this element will be ok
		gens = auts
		if orderG != length(auts)
		  gens = closure(auts, orderG)
    end
		conjs = nf_elem[x.prim_img for x in gens]
		prim_el = sum(conjs)
		def_pol = minpoly(prim_el)
    if degree(def_pol) != degree_subfield
			conjs1 = copy(conjs)
      while degree(def_pol) != degree_subfield
				for i = 1:length(conjs)
          conjs1[i] *= conjs[i]
				end
				prim_el = sum(conjs1)
       	def_pol = minpoly(prim_el)
			end
		end
    subK = number_field(def_pol, cached = false)[1]
    mp = hom(subK, K, prim_el, check = false)
    return subK, mp
	end
  =#
  OK = maximal_order(K)
  if isdefined(OK, :lllO)
    OK = lll(OK)
  end
  M = zero_matrix(FlintZZ, degree(K), degree(K)*length(auts_new))
  v = Vector{nf_elem}(undef, degree(K))
  MOK = basis_matrix(OK, copy = false)
  MOKinv = basis_mat_inv(OK, copy = false)
  for i = 1:length(auts_new)
		v[1] = one(K)
    v[2] = auts_new[i].prim_img
    for j = 3:degree(K)
      v[j] = v[j-1]*v[2]
		end
    B = basis_matrix(v, FakeFmpqMat)
    mul!(B, B, MOKinv)
    mul!(B, MOK, B)
    @assert isone(B.den)
    for i = 1:degree(K)
      B.num[i, i] -= 1
    end
		_copy_matrix_into_matrix(M, 1, (i-1)*degree(K)+1, B.num)
	end
	@vtime :Subfields 1 rk, Ker = kernel(M, side = :left)
  @assert rk == degree_subfield
  Ker = view(Ker, 1:rk, 1:degree(K))
  @vtime :Subfields 1 Ker = lll(Ker)
	#The kernel is the maximal order of the subfield.
  bas = Vector{nf_elem}(undef, degree_subfield)
	for i = 1:degree_subfield
    bas[i] = elem_from_mat_row(OK, Ker, i).elem_in_nf
	end
  return subfield(K, bas, isbasis = true)
end

################################################################################
#
#  Fixed field as relative extension
#
################################################################################

function fixed_field(K::AnticNumberField, auts::Vector{NfToNfMor}, ::Type{NfRel{nf_elem}}; simplify_subfield::Bool = true)
  F, mF = fixed_field(K, auts)
  if simplify_subfield
    F, mF1 = simplify(F, cached = false)
    mF = mF1*mF
  end
  all_auts = closure(auts, div(degree(K), degree(F)))
  Kx, x = PolynomialRing(K, "x", cached = false)
  p = prod(x-y.prim_img for y in all_auts)
  def_eq = map_coeffs(x -> haspreimage(mF, x)[2], p)
  L, gL = number_field(def_eq, cached = false, check = false)
  iso = hom(K, L, gL, mF.prim_img, gen(K))
  #I also set the automorphisms...
  autsL = Vector{NfRelToNfRelMor{nf_elem, nf_elem}}(undef, length(all_auts))
  for i = 1:length(autsL)
    autsL[i] = hom(L, L, iso(all_auts[i].prim_img))
  end
  Hecke._set_automorphisms_nf_rel(L, autsL)
  return L, iso
end


text:
"""

`SPGBoxResult`: Data structure that contains the results of the `spgbox`
optimization method.

  `x`: Array containing final point (solution or best point found).

  `f`: Final objective function value.

  `gnorm`: Norm of the projected gradient on the bounds.

  `nit`: Number of iterations performed.

  `nfeval`: Number of function evaluations.

  `ierr`: Status of execution:

         0. Converged successfully to solution.

         1. Maximum number of iterations achieved.

         2. Maximum number of function evaluations achieved.

"""
struct SPGBoxResult{T,F,G}
    x::T
    f::F
    gnorm::G
    nit::Int64
    nfeval::Int64
    ierr::Int64
end

function Base.show(io::IO, R::SPGBoxResult)
    println("")
    println(" SPGBOX RESULT: ")
    println("")
    R.ierr == 0 && println(" Convergence achieved. ")
    R.ierr == 1 && println(" Maximum number of iterations (nitmax) reached.")
    R.ierr == 2 && println(" Maximum number of function evaluations (nfevalmax) reached.")
    println("")
    println(" Final objective function value = ", R.f)
    print(" Sample of best point = $(typeof(R.x))[ ")
    for i = 1:min(length(R.x) - 1, 3)
        print("$(R.x[i]), ")
    end
    if length(R.x) > 4
        print("..., ")
    end
    println(R.x[length(R.x)], "]")
    println(" Projected gradient norm = ", R.gnorm)
    println("")
    println(" Number of iterations = ", R.nit)
    println(" Number of function evaluations = ", R.nfeval)
end


text:
"""
    group(by, iter)

Group the elements `x` of the iterable `iter` into groups labeled by `by(x)`.

# Example

```jldoctest
julia> group(iseven, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [1, 3, 5, 7, 9]
  true  => [2, 4, 6, 8, 10]
```
"""
group(by, iter) = group(by, identity, iter)

group(iter) = group(identity, identity, iter) # Not particularly useful, yet consistent with `groupinds`

"""
    group(by, f, iter...)

Sorts the elements `x` of the iterable `iter` into groups labeled by `by(x)`,
transforming each element to `f(x)`. If multiple collections (of the same length)
are provided, the transformations occur elementwise.

# Example

```jldoctest
julia> group(iseven, x -> x ÷ 2, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [0, 1, 2, 3, 4]
  true  => [1, 2, 3, 4, 5]
```
"""
function group(by, f, iter)
    # TODO Do this inference-free, like comprehensions...
    T = eltype(iter)
    K = promote_op(by, T)
    V = promote_op(f, T)

    out = Dict{K, Vector{V}}()
    for x ∈ iter
        key = by(x)
        push!(get!(Vector{V}, out, key), f(x))
    end
    return out
end

group(by, f, iter1, iter2, iters...) = group((x -> by(x...)), (x -> f(x...)), zip(iter1, iter2, iters...))

# For arrays we follow a different algorithm
group(by, f, a::AbstractArray) = group2(mapview(by, a), mapview(f, a))

function group2(groups, values)
    # TODO: assert that keys(groups) match up to keys(values)
    out = Dict{eltype(groups), typeof(empty(values))}()
    @inbounds for i ∈ keys(groups)
        group = groups[i]
        value = values[i]
        push!(get!(() -> empty(values), out, group), value)
    end
    return out
end

"""
    groupinds(by, iter...)

Sorts the indices `i` of `iter` into groups labeled by `by(iter[i])`. If multiple
collections (with matching indices) are provided, the groups are formed elementwise.

# Example

```jldoctest
julia> groupinds(iseven, [3,4,2,6,5,8])
Dict{Bool,Array{Int64,1}} with 2 entries:
  false => [1, 5]
  true  => [2, 3, 4, 6]
```
"""
function groupinds(by, iter)
    T = eltype(iter)
    K = promote_op(by, T)
    inds = keys(iter)
    V = eltype(inds)

    out = Dict{K, Vector{V}}()
    for i ∈ inds
        @inbounds x = iter[i]
        key = by(x)
        push!(get!(Vector{V}, out, key), i)
    end
    return out
end

groupinds(iter) = groupinds(identity, iter)

function groupinds(by, a::AbstractArray)
    _groupinds(mapview(by, a))
end

function _groupinds(a::AbstractArray)
    K = eltype(a)
    inds = keys(a)
    V = eltype(inds)

    out = Dict{K, Vector{V}}()
    @inbounds for i ∈ inds
        x = a[i]
        push!(get!(Vector{V}, out, x), i)
    end
    return out
end

# Semi-lazy grouping container
struct Groups{K, V, T, Inds} <: AbstractDict{K, V}
    parent::T
    inds::Inds
end

Base.keys(g::Groups) = keys(g.inds)
Base.@propagate_inbounds function Base.getindex(g::Groups, k)
    inds = g.inds[k]
    return view(g.parent, inds)
end
Base.length(g::Groups) = length(keys(g))

function Base.iterate(g::Groups, state...)
    i = iterate(g.inds, state...)
    if i === nothing
        return nothing
    else
        ((key, inds), newstate) = i
        return (key => @inbounds(g.parent[inds]), newstate)
    end
end

"""
    groupview(by, iter)

Like `group`, but each grouping is a view of the input collection `iter`.

# Example

```jldoctest
julia> v = [3,4,2,6,5,8]
6-element Array{Int64,1}:
 3
 4
 2
 6
 5
 8

julia> groups = groupview(iseven, v)
SAC.Groups{Bool,Any,Array{Int64,1},Dict{Bool,Array{Int64,1}}} with 2 entries:
  false => [3, 5]
  true  => [4, 2, 6, 8]

julia> groups[false][:] = 99
99

julia> v
6-element Array{Int64,1}:
 99
  4
  2
  6
 99
  8
```
"""
function groupview(by, iter)
    inds = groupinds(by, iter)
    V = promote_op(view, iter, eltype(inds))
    return Groups{keytype(inds), V, typeof(iter), typeof(inds)}(iter, inds)
end

groupview(iter) = groupview(identity, iter)

"""
    groupreduce(by, f, op, iter...; [init])

Applies a `mapreduce`-like operation on the groupings labeled by passing the elements of
`iter` through `by`. Mostly equivalent to `map(g -> reduce(op, g), group(by, f, iter))`,
but designed to be more efficient. If multiple collections (of the same length) are
provided, the transformations `by` and `f` occur elementwise.
"""
function groupreduce(by, f, op, iter; kw...)
    # TODO Do this inference-free, like comprehensions...
    nt = kw.data
    T = eltype(iter)
    K = promote_op(by, T)
    T2 = promote_op(f, T)
    V = promote_op(op, T2, T2)

    out = Dict{K, V}()
    for x ∈ iter
        key = by(x)
        dict_index = ht_keyindex2!(out, key)
        if dict_index > 0
            @inbounds out.vals[dict_index] = op(out.vals[dict_index], f(x))
        else
            if nt isa NamedTuple{()}
                Base._setindex!(out, convert(V, f(x)), key, -dict_index)
            elseif nt isa NamedTuple{(:init,)}
                Base._setindex!(out, convert(V, op(nt.init, f(x))), key, -dict_index)
            else
                throw(ArgumentError("groupreduce doesn't support the keyword arguments $(setdiff(keys(nt), (:init,)))"))
            end
        end
    end
    return out
end

groupreduce(by, f, op, iter1, iter2, iters...; kw...) = groupreduce((x -> by(x...)), (x -> f(x...)), op, zip(iter1, iter2, iters...); kw...)

"""
    groupreduce(by, op, iter; [init])

Like `groupreduce(by, identity, op, iter; init=init)`.

# Example

```jldoctest
julia> groupreduce(iseven, +, 1:10)
Dict{Bool,Int64} with 2 entries:
false => 25
true  => 30
```
"""
groupreduce(by, op, iter; kw...) = groupreduce(by, identity, op, iter; kw...)
groupreduce(op, iter; kw...) = groupreduce(identity, identity, op, iter; kw...)

# Special group operators

@deprecate grouplength(by, iter) groupcount(by, iter)
export grouplength

"""
    groupcount([by], iter)

Determine the number of elements of `iter` belonging to each group.
"""
groupcount(by, iter) = groupreduce(by, x->1, +, iter)
groupcount(iter) = groupcount(identity, iter) # A handy extension of `unique`

"""
    groupsum(by, [f], iter)

Sum the elements of `iter` belonging to different groups, optionally mapping by `f`.
"""
groupsum(by, iter) = groupreduce(by, identity, +, iter)
groupsum(by, f, iter) = groupreduce(by, f, +, iter)
groupsum(iter) = groupreduce(identity, identity, +, iter) # For consistency

"""
    groupprod(by, [f], iter)

Multiply the elements of `iter` belonging to different groups, optionally mapping by `f`.
"""
groupprod(by, iter) = groupreduce(by, identity, *, iter)
groupprod(by, f, iter) = groupreduce(by, f, *, iter)
groupprod(iter) = groupreduce(identity, identity, *, iter) # For consistency


text:
import Dates

const monitor_default = Vector{Vector{Symbol}}([])
const thin_default = Int[]

"""
Partitions a list `xs` into `(good, bad)` according
to a condition.
"""
function partition(condition::Function, xs::Vector{T}) where T
  good = T[]
  bad = T[]

  for x in xs
    if condition(x)
      append!(good, [x])
    else
      append!(bad, [x])
    end
  end

  return good, bad
end

"""
Checks if a field is a subtype.
If a field has the form "a__b", the name
of the field is `a`, and it has a field `b`.
"""
issubtype(x::Symbol)::Bool = occursin("__", String(x))

function deepcopyFields(state::T, fields::Vector{Symbol}) where T
  substate = Dict{Symbol, Any}()

  # Partition fields into subtypes and regular fields
  subtypefields, regfields = partition(issubtype, fields)

  # Get level1 fields
  for field in regfields
    substate[field] = deepcopy(getfield(state, field))
  end

  # Get level2 fields
  for field in subtypefields
    statename, _field = split(String(field), "__")
    _state = getfield(state, Symbol(statename))
    substate[field] = deepcopy(getfield(_state, Symbol(_field)))
  end

  return substate
end

showtime() = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")

"""
TODO...
"""
function gibbs(init::T,
               update::Function;
               monitors::Vector{Vector{Symbol}}=deepcopy(monitor_default),
               thins::Vector{Int}=deepcopy(thin_default),
               nmcmc::Int64=1000, nburn::Int=0,
               printFreq::Int=0, loglike=missing,
               printlnAfterMsg::Bool=true) where T

  @assert printFreq >= -1
  if printFreq == 0
    numPrints = 10
    printFreq = Int(ceil((nburn + nmcmc) / numPrints))
  end

  state = deepcopy(init)

  # Checking number of monitors.
  numMonitors = length(monitors)
  println("Number of monitors: $(numMonitors)"); flush(stdout)
  @assert numMonitors == length(thins)

  # Check monitor
  if numMonitors == 0
    println("Using default monitor."); flush(stdout)
    fnames = [ fname for fname in fieldnames(typeof(init)) ]
    append!(monitors, [fnames])
    append!(thins, 1)
    numMonitors = 1
  end

  # Number of Samples for each Monitor
  numSamps = [ div(nmcmc, thins[i]) for i in 1:numMonitors ]

  if printFreq > 0
    println("Preallocating memory..."); flush(stdout)
  end

  # Create object to return
  @time out = [fill(deepcopyFields(state, monitors[i]), numSamps[i])
               for i in 1:numMonitors]

  function printMsg(i::Int)
    if (printFreq > 0) && (i % printFreq == 0) && (i > 1)
      loglikeMsg = ismissing(loglike) ? "" : "-- loglike: $(last(loglike))"
      print("$(showtime()) -- $(i)/$(nburn+nmcmc) $loglikeMsg"); flush(stdout)

      if printlnAfterMsg
        println(); flush(stdout)
      end

      flush(stdout)
    end
  end

  # burn in
  for i in 1:nburn
    update(state, i, out)
    printMsg(i)
  end

  counters = zeros(Int, numMonitors)
  # Gibbs loop
  for i in 1:nmcmc
    update(state, i + nburn, out)
    printMsg(i + nburn)

    for j in 1:numMonitors
      if i % thins[j] == 0
        substate = deepcopyFields(state, monitors[j])
        counters[j] += 1
        out[j][counters[j]] = substate
      end
    end
  end

  return (out, state)
end


text:
# ---
# title: 28. Implement strStr()
# id: problem28
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Two Pointers, String
# link: <https://leetcode.com/problems/implement-strstr/description/>
# hidden: true
# ---
#
# Implement [strStr()](http://www.cplusplus.com/reference/cstring/strstr/).
#
# Return the index of the first occurrence of needle in haystack, or `-1` if
# `needle` is not part of `haystack`.
#
# **Clarification:**
#
# What should we return when `needle` is an empty string? This is a great
# question to ask during an interview.
#
# For the purpose of this problem, we will return 0 when `needle` is an empty
# string. This is consistent to C's
# [strstr()](http://www.cplusplus.com/reference/cstring/strstr/) and Java's
# [indexOf()](https://docs.oracle.com/javase/7/docs/api/java/lang/String.html#indexOf\(java.lang.String\)).
#
#
#
# **Example 1:**
#
#
#
#     Input: haystack = "hello", needle = "ll"
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: haystack = "aaaaa", needle = "bba"
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: haystack = "", needle = ""
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `0 <= haystack.length, needle.length <= 5 * 104`
#   * `haystack` and `needle` consist of only lower-case English characters.
#
#
## @lc code=start
using LeetCode

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


text:
# Integer conic problems

function intsoc1test(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    #@test MOI.supportsproblem(model, MOI.ScalarAffineFunction{Float64},
    #    [(MOI.VectorAffineFunction{Float64},MOI.Zeros),
    #     (MOI.SingleVariable,MOI.ZeroOne),
    #     (MOI.VectorOfVariables,MOI.SecondOrderCone)])
    @testset "INTSOC1" begin

        # Problem SINTSOC1
        # min 0x - 2y - 1z
        #  st  x            == 1
        #      x >= ||(y,z)||
        #      (y,z) binary

        MOI.empty!(model)
        @test MOI.isempty(model)

        MOI.canaddvariable(model)
        x,y,z = MOI.addvariables!(model, 3)

        @test MOI.canset(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
        MOI.set!(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([y,z],[-2.0,-1.0],0.0))
        @test MOI.canset(model, MOI.ObjectiveSense())
        MOI.set!(model, MOI.ObjectiveSense(), MOI.MinSense)

        @test MOI.canaddconstraint(model, MOI.VectorAffineFunction{Float64}, MOI.Zeros)
        ceq = MOI.addconstraint!(model, MOI.VectorAffineFunction([1],[x],[1.0],[-1.0]), MOI.Zeros(1))
        @test MOI.canaddconstraint(model, MOI.VectorOfVariables, MOI.SecondOrderCone)
        csoc = MOI.addconstraint!(model, MOI.VectorOfVariables([x,y,z]), MOI.SecondOrderCone(3))

        @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.Zeros}()) == 1
        @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorOfVariables,MOI.SecondOrderCone}()) == 1
        loc = MOI.get(model, MOI.ListOfConstraints())
        @test length(loc) == 2
        @test (MOI.VectorAffineFunction{Float64},MOI.Zeros) in loc
        @test (MOI.VectorOfVariables,MOI.SecondOrderCone) in loc

        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(y)), typeof(MOI.ZeroOne))
        bin1 = MOI.addconstraint!(model, MOI.SingleVariable(y), MOI.ZeroOne)
        @test MOI.canaddconstraint(model, typeof(MOI.SingleVariable(z)), typeof(MOI.ZeroOne))
        bin2 = MOI.addconstraint!(model, MOI.SingleVariable(z), MOI.ZeroOne)

        if config.solve
            MOI.optimize!(model)

            @test MOI.canget(model, MOI.TerminationStatus())
            @test MOI.get(model, MOI.TerminationStatus()) == MOI.Success

            @test MOI.canget(model, MOI.PrimalStatus())
            @test MOI.get(model, MOI.PrimalStatus()) == MOI.FeasiblePoint

            @test MOI.canget(model, MOI.ObjectiveValue())
            @test MOI.get(model, MOI.ObjectiveValue()) ≈ -2 atol=atol rtol=rtol

            @test MOI.canget(model, MOI.VariablePrimal(), MOI.VariableIndex)
            @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 atol=atol rtol=rtol
            @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1 atol=atol rtol=rtol
            @test MOI.get(model, MOI.VariablePrimal(), z) ≈ 0 atol=atol rtol=rtol
        end

    end
end

function intsoctests(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    intsoc1test(model, atol=atol, rtol=rtol)
end

function intconictests(model::MOI.ModelLike; atol=Base.rtoldefault(Float64), rtol=Base.rtoldefault(Float64))
    intsoctests(model, atol=atol, rtol=rtol)
end


text:
################################################################################
# MPI-specific communication methods

# Execute a command on all MPI ranks
# This uses MPI as communication method even if @everywhere uses TCP
function mpi_do(mgr::Union{MPIManager,MPIWorkerManager}, expr)
  !mgr.initialized && wait(mgr.cond_initialized)
  jpids = keys(mgr.j2mpi)
  refs = Array{Any}(undef, length(jpids))
  for (i,p) in enumerate(Iterators.filter(x -> x != myid(), jpids))
      refs[i] = remotecall(expr, p)
  end
  # Execution on local process should be last, since it can block the main
  # event loop
  if myid() in jpids
      refs[end] = remotecall(expr, myid())
  end

  # Retrieve remote exceptions if any
  @sync begin
      for r in refs
          @async begin
              resp = remotecall_fetch(r.where, r) do rr
                  wrkr_result = rr[]
                  # Only return result if it is an exception, i.e. don't
                  # return a valid result of a worker computation. This is
                  # a mpi_do and not mpi_callfetch.
                  isa(wrkr_result, Exception) ? wrkr_result : nothing
              end
              isa(resp, Exception) && throw(resp)
          end
      end
  end
  nothing
end

macro mpi_do(mgr, expr)
  quote
      # Evaluate expression in Main module
      thunk = () -> (Core.eval(Main, $(Expr(:quote, expr))); nothing)
      mpi_do($(esc(mgr)), thunk)
  end
end


text:
using JuliaBlackBoard

# used to create .txt file with same basename and same directory as script
OPEN = JuliaBlackBoard.OPEN(@__FILE__)

OPEN() do io

    ## Questions can be simple strings
    ## The text is in markdown
    q = "What is *1 + 1*?"
    question(io, NUM, q, 1+1)

    ## Typically though questions use Mustache templates
    ## substitution can be done by named arguments
    q = mt"""
# Addition

What is *{{:x}} + {{:y}}*?
"""
    x, y = 2, 3
    question(io, NUM, q(x=x, y=y), x + y)

    ## Math Markup
    ##
    ## BlackBoard's question editor has a formula editor for math
    ## dispaly.  The math editor is written by WIRIS and based on
    ## standards like MathML for internal representation and the PNG
    ## image format for displaying formulas. The math editor is based
    ## on Javascript and runs on any browser and operating system,
    ## including smartphones and tablets.
    ##
    ## This is not directly accessible from a file upload, thoug the
    ## last example shows how it can be. Rather, from within this framework,
    ## we can use LaTeX to mark up questions in different ways

    ## Using TtH for markup
    ## LaTeX-formatted equations are rendered in HTML through TtH
    ## This has limitations for some formula (e.g. \bar{x} or \frac{x}{y}), but
    ## is okay.
    q = mt"""
# Division

Compute $z$ where

$$
z = \frac{ {{:x}} }{ {{:y}} }.
$$
"""
    x, y = 4, 2
    question(io, NUM, q(x=x, y=y), x/y)

    ## using LaTeX. The entire question can be a png image from
    ## LaTeX output. The `LaTeX` command formats the question.
    ## Alternatively, using `ltx"""` in place of `mt"""` and LaTeX
    ## will be called without specifying it, as below:
    q = mt"""
\noindent\textbf{Fractions}

Compute $z$ where
$$
z = \frac{1}{ {{:x}} } + \frac{1}{ {{:y}} }.
$$
"""
    x, y = 2, 3
    str = q(x=x, y=y)
    question(io, NUM, LaTeX(str), 1/x + 1/y)

    ## Using markdown and using CommonMark to convert to latex, then running LaTeX
    ## This is done by using `LaTeX′` (`\prime[tab]`) or, as below, the `mdltx"""`
    ## string macro
    q = mdltx"""
# exponents

Compute $z$ where
$$
z = {{:n}}^{{:m}}
$$

"""
    n,m = 2, 3
    str = q(n=n, m=m)
    question(io, NUM, str, n^m)

    ## Mix and match. It seems silly, but ala y2k's LaTeX2HTML
    ## The `LaTeX(str)` command produces an image that can be included within
    ## a markdown-formatted question, as follows:
    q = mt"""
# subtraction

Compute {{:latex_z}}  where

{{:latex_x_y}}
"""
    x, y = 5, 3
    latex_z = LaTeX(mt"$z$")
    latex_x_y = LaTeX(mt"$$z = {{:x}} - {{:y}}$$"(x=x, y=y))
    question(io, NUM, q(latex_z=latex_z, latex_x_y=latex_x_y), x - y)

    ## Finally, in the question editor on BlackBoard, LaTeX markup can be
    ## used. (Inline markup is indicated with `$$...$$`.) *However*,
    ## this only renders after the question is submitted for
    ## saving. (At which point, an image is generated that is
    ## displayed).  If the environment variable `USE_MATHJAX` is set
    ## to `true`, then the markup will include the dollar signs. The
    ## question must be displayed in the BlackBoard editor and
    ## submitted for the image to be generated (though not
    ## modified). (https://blackboard.secure.force.com/publickbarticleview?id=kA339000000L6QH)

    ENV["USE_MATHJAX"] = true

    q = mt"""
# Division

Compute $z$ where

$$
z = \frac{ {{:x}} }{ {{:y}} }.
$$
"""
    x, y = 7, 2
    question(io, NUM, q(x=x, y=y), x/y)

    ENV["USE_MATHJAX"] = false # restore default

end


text:
using MonteCarloX
using Test
using Random
using HypothesisTests
using Distributions
import Distributions.pdf
import Distributions.cdf

@testset "MonteCarloX.jl" begin
    # basics and helper
    include("test_utils.jl")
    @test test_histogram_set_get()

    include("test_event_handler.jl")
    @test test_event_handler_rate("ListEventRateSimple")
    @test test_event_handler_rate("ListEventRateActiveMask")

    # equilibrium
    include("test_equilibrium.jl")
    @test test_sweep_random_element()
    @test test_unimodal_metropolis()
    @test test_2D_unimodal_metropolis()
    @test test_unimodal_sweep()

    include("test_ising.jl")
    @test test_ising_reweighting()
    @test test_ising_metropolis()
    @test test_ising_cluster()

    # non-equilibrium
    include("test_inhomogenous_poisson.jl")
    @test test_poisson_single()
    @test test_poisson_constant()
    @test test_poisson_sin_wave()

    # test_kinetic_monte_carlo.jl
    include("test_kinetic_monte_carlo.jl")
    @test test_kmc_next()
    include("test_gillespie.jl")
    # @test test_gillespie()

end


text:
module PolaronMakie

export viewfunction

# using AbstractPlotting.MakieLayout
using AbstractPlotting
using GLMakie

include("../polaronmakie/function.jl")

end # PolaronMakie


text:
module ColorBrewer

using Colors
using JSON

# The colorbrewer.json file contains the schemes,
# and was obtained from colorbrewer2.org.
# The contents of that file are covered under the Apache License 2.0, included
# in this distribution under data/.
# Here, dirname is used to accommodate non-standard package directory
# configurations (as well as standard ones, obviously).
colorSchemes = JSON.parsefile(
    dirname(@__FILE__) * "/../data/colorbrewer.json");

# Creates a color palette based on the colorbrewer2 color schemes.
#
# Args:
#   schemeName: The name of the scheme. The schemes are as follows:
#       - sequential: OrRd, PuBu, BuPu, Oranges, BuGn, YlOrBr, YlGn, Reds,
#           RdPu, Greens, YlGnBu, Purples, GnBu, Greys, YlOrRd, PuRd, Blues,
#           PuBuGn
#       - diverging: Spectral, RdYlGn, RdBu, PiYG, PRGn, RdYlBu, BrBG, RdGy,
#           PuOr
#       - qualitative: Set1, Set2, Set3, Accent, Dark2, Paired, Pastel1, Pastel2
#   n: The number of colors in the palette. Must be within the proper range.
#       The range for each type is as follows:
#       - sequential: between 3 and 9, inclusive
#       - diverging: between 3 and 11, inclusive
#       - qualitative: between 3 and 12 for Paired, Set3
#                      between 3 and 9 for Pastel1, Set1
#                      between 3 and 8 for Set2, Accent, Dark2, Pastel2
#
# Throws:
#   ArgumentError: If the scheme name is invalid.
#   ArgumentError: If the number of colors in the palette is invalid.
#
# Returns: An array of RGBs (from the Colors package).
function palette(schemeName::AbstractString, n::Integer)
    if !haskey(colorSchemes, schemeName)
        throw(ArgumentError("Scheme $schemeName is undefined."));
    end #Close if statement validating scheme name.

    if !haskey(colorSchemes[schemeName], string(n))
        throw(ArgumentError(
            "Scheme $schemeName is undefined for $n colors."));
    end #Close if statement validating n for the scheme name.

    return [parse(Colorant, s) for s in colorSchemes[schemeName][string(n)]];
end

export palette

end # module


text:
using Images, Docile

makedocs()

cp(joinpath(dirname(@__FILE__), "src", "img"), joinpath(dirname(@__FILE__), "build", "img"), remove_destination=true)
cp(joinpath(dirname(@__FILE__), "..", "LICENSE.md"), joinpath(dirname(@__FILE__), "build", "LICENSE.md"), remove_destination=true)

# When ready, deploy from top level Images dir with:
#   mkdocs gh-deploy --clean


text:
# intension
# for performance test, I want to have a big formula
# genform(dep, wid, syms::Array{Symbol}) makes such thing.
# wid is width size, dep is depth size.
# It may be a rectangle.
# the symbols are choosed from syms
# rand(syms) may be ok.
#

function baseff(var, fun, n)
  exp = :(f())
  exp.args = fill(var, n+1)
  exp.args[1] = fun
  exp
end

function replace1(var, base, arg)
  rexp = deepcopy(base)
  for ix in 1:length(base.args)
    if var == rexp.args[ix]
      rexp.args[ix] = arg
    end
  end
  return rexp
end

function makedd(var, exp, sexp, ndep, nwid)
  rexp = exp
  if ndep <= 0
    return replace1(var, exp, sexp)
  else
    sexp = makedd(var, exp, sexp, ndep-1, nwid)
    rep = replace1(var, exp, sexp)
    return rep
  end
end

function makeff(var, fsym, ssym, ndep, nwid)
  bexp =  baseff(var, fsym, nwid)
  if ndep <= 0; return bexp end
  base = baseff(var, ssym, nwid)
  narg = deepcopy(base)
  narg = makedd(var, base, narg, ndep-2, nwid)
  rexp = replace1(var, bexp, narg)
  return rexp

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
#=
Very similar to examples/plasticityFFtest.jl
Here I consider triplet rule, two high frequency neurons and two low
frequency ones (all Poisson).  Desired outcome: weight increases between the
high firing neurons, decreases in all other cases.
=#
push!(LOAD_PATH, abspath(@__DIR__,".."))

using Test,ProgressMeter
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using BenchmarkTools

using Random ; Random.seed!(1)

##
const dt = 0.1E-3
const Ttot = 20.0
const highrate = 80.0
const lowrate = 5.0
const nneus = 4

nt,ps = let spkgen = S.SGPoissonMulti([highrate,highrate,lowrate,lowrate])
  # spkgen = S.SGPoissonMulti([highrate,highrate])
  sker = S.SKExp(-1E6) # this will never be used
  vrev = -3000.0 # same
  nt = S.NTInputConductance(spkgen,sker,vrev)
  ps = S.PSInputConductance(nt,nneus)
  (nt,ps)
end

# plasticity and initial weights

wstart = let w = 100.0
  wm = fill(w,(nneus,nneus))
  wm[diagind(wm)] .= 0.0
  sparse(wm)
end

triplets_plasticity = let (n_post,n_pre) = size(wstart),
  τplus = 17E-3 # 16.8
  τminus = 34E-3 # 33.7
  τx = 100E-3 # 101
  τy = 120E-3 # 125
  plast_eps = 1E-3
  A2plus = 7.5E-7*plast_eps
  A3plus = 9.3*plast_eps
  A2minus = 7.0*plast_eps
  A3minus = 2.3E-1*plast_eps
  S.PlasticityTriplets(τplus,τminus,τx,τy,A2plus,A3plus,
    A2minus,A3minus,n_post,n_pre)
end

conn = S.ConnectionPlasticityTest(wstart,triplets_plasticity)

pop = S.Population(ps,(conn,ps))

# that's it, let's make the network
ntw = S.RecurrentNetwork(dt,pop)

## recorders : spikes
krec = 1
rec_spikes = S.RecSpikes(ps,highrate,Ttot)
rec_weights = S.RecWeights(conn,krec,dt,Ttot)

## Run

times = (0.0:ntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([ps,rec_spikes,rec_weights,conn.plasticities[1]])
# initial conditions

@time @showprogress for (k,t) in enumerate(times)
  rec_spikes(t,k,ntw)
  rec_weights(t,k,ntw)
  S.dynamics_step!(t,ntw)
end

spk_out = S.get_spiketimes_dictionary(rec_spikes)
tw,weights_all = S.rebuild_weights(rec_weights)
##
# check spikes
_ = let plt=plot(leg=false)
  for i in 1:nneus
   scatter!(plt, _->i, spk_out[i];marker=:vline)
  end
  plt
 end

##

_ = let plt=plot(legend=:topleft,xlabel="time",ylabel="weight")
  x = rec_weights.times
  n = size(weights_all[1],1)
  pre = 1
  for post in 1:n
    if post != pre
      y = getindex.(weights_all,post,pre)
      plot!(plt,x,y;linewidth=2,label="$(post) <- $(pre)")
    end
  end
  plot(plt)
end

##

##############################
# Let's make a function to test it

function test_plasticity(lowrate::R,highrate::R,
    myplasticity::S.PlasticityRule;(Ttot::R)=50.0,(dt::R)=0.1E-3) where R
  nneus = 2
  nt,ps = let spkgen = S.SGPoissonMulti([highrate,lowrate])
    sker = S.SKExp(-1E6) # this will never be used
    vrev = -3000.0 # same
    nt = S.NTInputConductance(spkgen,sker,vrev)
    ps = S.PSInputConductance(nt,nneus)
    (nt,ps)
  end
  wstart = let w = 100.0
    wm = fill(w,(nneus,nneus))
    wm[diagind(wm)] .= 0.0
    sparse(wm)
  end
  conn = S.ConnectionPlasticityTest(wstart,myplasticity)
  pop = S.Population(ps,(conn,ps))
  ntw = S.RecurrentNetwork(dt,pop)

  ## recorders : spikes
  krec = 1
  rec_spikes = S.RecSpikes(ps,highrate,Ttot)
  rec_weights = S.RecWeights(conn,krec,dt,Ttot)

  ## Run
  times = (0.0:ntw.dt:Ttot)
  nt = length(times)
  # clean up
  S.reset!.([ps,rec_spikes,rec_weights,conn.plasticities[1]])
  for (k,t) in enumerate(times)
    rec_spikes(t,k,ntw)
    rec_weights(t,k,ntw)
    S.dynamics_step!(t,ntw)
  end
  # FROM HIGH TO LOW
  idx_low_high = findfirst(idx->Tuple(idx)==(2,1),rec_weights.idx_save)
  w_low_high = rec_weights.weights_now[idx_low_high,:]
  # FROM LOW TO HIGH
  idx_high_low = findfirst(idx->Tuple(idx)==(1,2),rec_weights.idx_save)
  w_high_low = rec_weights.weights_now[idx_high_low,:]
  _dostuff = function(v::Vector{R})
    vd=diff(v)
    filter!(!=(0.0),vd)
    sum(vd)/Ttot
  end
  return _dostuff.((w_low_high,w_high_low))
end
##
triplets_plasticity = let (n_post,n_pre) = size(wstart),
  τplus = 17E-3 # 16.8
  τminus = 34E-3 # 33.7
  τx = 100E-3 # 101
  τy = 120E-3 # 125
  plast_eps = 1E-3
  A2plus = 7.5E-7*plast_eps
  A3plus = 9.3*plast_eps
  A2minus = 7.0*plast_eps
  A3minus = 2.3E-1*plast_eps
  S.PlasticityTriplets(τplus,τminus,τx,τy,A2plus,A3plus,
    A2minus,A3minus,n_post,n_pre)
end

booh = test_plasticity(10.0,80.0,triplets_plasticity)

##
nrats = 25
ratestest = range(0.5,60;length=nrats)
testret = fill(NaN,nrats,nrats)

for i in 1:nrats, j in i:nrats
  # j is high, i is low or equal
  @info "now processing $i,$j (out of $nrats)"
  (low_high,high_low) = test_plasticity(ratestest[i],
    ratestest[j], triplets_plasticity ; Ttot=25.0)
  testret[i,j] = low_high
  testret[j,i] = high_low
end
##

heatmap(ratestest,ratestest,testret;color=:vik)

function plus_minus_rescale(mat::Array{R}) where R
  down,up = extrema(mat)
  @assert (down < 0.0) && (up > 0.0)
  pluspart = mat ./ up
  minuspart = mat ./ abs(down)
  ret = similar(mat)
  for i in eachindex(mat)
    if mat[i] > 0
      ret[i] = pluspart[i]
    else
      ret[i] = minuspart[i]
    end
  end
  return ret
end

heatmap(ratestest,ratestest,testret;
  color=:seismic,ratio=1)

heatmap(ratestest,ratestest,plus_minus_rescale(testret);
  color=:seismic,ratio=1)


text:
export dnorm, dgumbel, dgamma, dlognorm

"""
    dnorm(x::Real, μ::Real, σ::Real)

Normal distribution density function
"""
dnorm(x::Real, μ::Real, σ::Real) = exp(-0.5 * ((x - μ)/σ)^2)/(σ*sqrt2π)

"""
    dgumbel(x::Real, μ::Real, σ::Real)

Gumbel density, parameterized by mean (μ) and SD (σ) of x.
"""
function dgumbel(x::Real, μ::Real, σ::Real)
    β = (√6 * σ) / π
    α = μ - (β*γ)
    z = (x - α) / β
    exp(-(z + exp(-z)))
end

"""
    dgamma(x::Real, μ::Real, σ::Real)

Gamma density, parameterized by mean (μ) and sigma(σ) of x.
"""
function dgamma(x::Real, μ::Real, σ::Real)
    α = μ^2 / σ^2
    β = μ / σ^2
    (x^(α - 1) * exp(-β * x))  * (β^α / Γ(α))
end

"""
    dlognorm(x::Real, μ::Real, σ::Real)

Lognormal density, parameterized by mean (μ) and SD (σ) of x, NOT of log(x).
"""
function dlognorm(x::Real, μ::Real, σ::Real)
    lμ = log(μ^2 / sqrt(1 + σ^2 / μ^2))
    lσ = √(log(1 + σ^2 / μ^2))
    dnorm(log(x), lμ, lσ) / x
end


text:
# url to remote organization repo
const REMOTE_URL = "https://github.com/JuliaPOMDP/"

# supported solvers
const SUPPORTED_SOLVERS = Set{AbstractString}(
                          ["DiscreteValueIteration",
                           "MCTS",
                           "QMDP",
                           "SARSOP",
                           "POMCP",
                           "POMDPSolve"])

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# Optimize AD loads via one-step MILP after the shiftable load's decisions

function optimize_AD_loads(h::AbstractDict, ρ::AbstractFloat, PPV::AbstractFloat,
    P_shiftable::AbstractVector; K::Integer = 20, verbose=false)
    ls = h["AD"]
    g = h["UtilityGrid"]
    @assert length(h["SU"]) + length(h["SI"]) == length(P_shiftable)
    model = JuMP.Model(Cbc.Optimizer)
    set_optimizer_attribute(model, "logLevel", 0)
    # only consider the current time t
    for load in ls
        add_AD_load_t!(model, load, ρ; K)
    end
    setup_power_exchange_t!(model, g, ρ)
    # power balance
    P_other = sum(P_shiftable) + sum(l["P"] for l in h["CR"])   # critical and shiftable load consumption
    @constraint(model, cst_power_balance, P_other + sum(model[Symbol("P_", l["id"])]
                                                        for l in ls) == model[:Pb] - model[:Ps] + PPV)
    ρs = g["β"] .* ρ
    @objective(model, Min, sum(model[Symbol("C_", l["id"])] for l in ls) + (ρ - ρs) * model[:Ps])

    optimize!(model)
    if termination_status(model) == OPTIMAL
        # compute total cost
        # @show objective_value(model)
    else
        @warn "Not OPTIMAL"
    end
    if verbose
        @show P_other value(sum(model[Symbol("P_", l["id"])] for l in ls)) value(model[:Pb]) value(model[:Ps]) PPV
    end
    return model
end

function setup_power_exchange_t!(model, g, ρ)
    Pmax = g["Pmax"] # any large value is OK
    β = g["β"]
    @assert 0 <= β <= 1
    @variable(model, zb, Bin)
    @variable(model, zs, Bin)
    @variable(model, 0 <= Pb <= Pmax)
    @variable(model, 0 <= Ps <= Pmax)
    @constraint(model, Pb <= zb * Pmax)
    @constraint(model, Ps <= zs * Pmax)
    @constraint(model, zb + zs <= 1)
    ρs = β .* ρ
    @expression(model, Cele, ρ .* Pb - ρs .* Ps)
    return model
end

function add_AD_load_t!(model::JuMP.Model, load::AbstractDict, ρ::Float64; K::Integer = 20)
    @assert K >= 3
    Pmin, Pmax = load["Pmin"], load["Pmax"]
    id = load["id"]
    model[Symbol("P_", id)] = P = @variable(model, lower_bound = Pmin, upper_bound = Pmax)
    model[Symbol("C_", id)] = C = @variable(model)
    f(P) = ρ * P + load["α"] * (P - Pmax)^2
    xs, ks, bs = linearize(f, Pmin, Pmax, K)
    @constraint(model, C .>= ks .* P .+ bs)
    model[Symbol("C_", id, "_dc")] = @expression(model, load["α"] .* (P .- Pmax) .^ 2)
    return model
end

"""
    manage_loads(s::DataFrame, h::AbstractDict, agents::AbstractDict; K::Integer = 20) -> DataFrame

Given the `agents` managing shiftable loads, perform load management for home `h` in scenario `s`.
`K` is the number of segments in linear approximation of the quadratic cost.

The management strategy is returned in a new data frame along with the original scenario info.
"""
function manage_loads(s::DataFrame, h::AbstractDict, agents::AbstractDict; K::Integer = 20)
    s = s[:, [:ρ, :PPV]]    # a copy due to :
    ρ, PPV = s.ρ, s.PPV
    T = size(s)[1]
    # insert relevant columns
    for l in h["AD"]
        s[!, Symbol("P_", l["id"])] .= fill(NaN, T)
        s[!, Symbol("C_", l["id"])] .= fill(NaN, T)
        s[!, Symbol("C_", l["id"], "_dc")] .= fill(NaN, T)
    end
    for l in Iterators.Flatten((h["SU"], h["SI"]))
        s[!, Symbol("z_", l["id"])] .= fill(NaN, T)
        s[!, Symbol("P_", l["id"])] .= fill(NaN, T)
        s[!, Symbol("C_", l["id"])] .= fill(NaN, T)
    end
    for col in (:zb, :zs, :Pb, :Ps)
        s[!, col] = fill(NaN, T)
    end
    # schedule shiftable loads
    for l in h["SU"]
        schedule_SU_load!(s, l, agents[l["id"]])
    end
    for l in h["SI"]
        schedule_SI_load!(s, l, agents[l["id"]])
    end
    # AD loads
    info_AD = Dict(l["id"] => (l["α"], l["Pmax"]) for l in h["AD"])
    for t = 1:T
        P_shiftable = [s[t, Symbol("P_", l["id"])] for l in Iterators.Flatten((h["SU"], h["SI"]))]
        model = optimize_AD_loads(h, ρ[t], PPV[t], P_shiftable; K, verbose = false)
        Ps = Dict(l["id"] => value(model[Symbol("P_", l["id"])]) for l in h["AD"])
        for (id, P) in Ps
            s[t, Symbol("P_", id)] = P
            α, Pmax = info_AD[id]
            s[t, Symbol("C_", id)] = ρ[t] * P + α * (P - Pmax)^2
            s[t, Symbol("C_", id, "_dc")] = α * (P - Pmax)^2
        end
        # write also the grid power exchange
        for col in (:zb, :zs, :Pb, :Ps)
            s[t, col] = value(model[col])
        end
    end
    β = h["UtilityGrid"]["β"]
    s.Cele = ρ .* s.Pb .- β * ρ .* s.Ps
    #step-wise cost
    s.C = push!([s[!, Symbol("C_", l["id"], "_dc")] for l in h["AD"]], s.Cele) |> sum
    return s
end

function schedule_SU_load!(s::DataFrame, l, agent)
    id, ts, tf = l["id"], l["ts"], l["tf"]
    P, L = l["P"], l["L"]
    s[!, Symbol("z_", id)] .= 0
    s[!, Symbol("P_", id)] .= 0.0
    s[!, Symbol("C_", id)] .= 0.0

    function switch_on(t)
        s[t, Symbol("z_", id)] = 1
        for dt = 0:L-1
            s[t+dt, Symbol("P_", id)] = P
            s[t+dt, Symbol("C_", id)] = s.ρ[t+dt] * P
        end
    end
    # action can be 1 only inside this time range
    for t = ts:tf-L+1
        state = [t, s.ρ[t], s.PPV[t]]
        action = only(agent(state)) >= 0.5
        if action == 1 || t == tf - L + 1  # tf - L + 1 is the last chance
            switch_on(t)
            break  # only a single 1 is allowed
        end
    end
    @assert sum(s[!, Symbol("z_", id)]) == 1
end

function schedule_SI_load!(s::DataFrame, l, agent)
    id, ts, tf = l["id"], l["ts"], l["tf"]
    P, L = l["P"], l["L"]
    s[!, Symbol("z_", id)] .= 0
    s[!, Symbol("P_", id)] .= 0.0
    s[!, Symbol("C_", id)] .= 0.0

    function keep_on(t)
        s[t, Symbol("z_", id)] = 1
        s[t, Symbol("P_", id)] = P
        s[t, Symbol("C_", id)] = s.ρ[t] * P
    end

    n_already_on = 0
    for t = ts:tf
        if tf - t + 1 == L - n_already_on   # all remaining steps must be ON
            for tr = t:tf
                keep_on(tr)
            end
            break
        end
        state = [t, s.ρ[t], s.PPV[t], n_already_on]
        action = only(agent(state)) >= 0.5
        if action == 1
            n_already_on += 1
            keep_on(t)
            if n_already_on == L
                break
            end
        end
    end
    @assert sum(s[!, Symbol("z_", id)]) == L
end

max_stars_repo_name: ShuhuaGao/HybridHEMS

max_issues_repo_name: ShuhuaGao/HybridHEMS

max_forks_repo_name: ShuhuaGao/HybridHEMS


text:
using DiskArrayTools: InterpolatedDiskArray, DiskArrayTools
using DiskArrayTools.Interpolations: Linear, Flat, Constant, NoInterp
using DiskArrays: eachchunk, GridChunks, approx_chunksize, grid_offset
"""
  spatialinterp(c,newlons::AbstractRange,newlats::AbstractRange;order=Linear(),bc = Flat())
"""
function spatialinterp(c,newlons::AbstractRange,newlats::AbstractRange;order=Linear(),bc = Flat())
  interpolatecube(c,Dict("Lon"=>newlons, "Lat"=>newlats), order = Dict("Lon"=>order,"Lat"=>order))
end
spatialinterp(c,newlons::CubeAxis,newlats::CubeAxis;kwargs...)=
  spatialinterp(c,newlons.values,newlats.values;kwargs...)

function getsteprat(newax::AbstractRange, oldax::AbstractRange)
    step(oldax)/step(newax)
end

sorted(x,y,rev=false) = rev ? sorted(y,x,false) : x<=y ? (x,y) : (y,x)

function getsteprat(newax, oldax)
  val1,val2 = sorted(first(newax),last(newax))
  rev = last(oldax) < first(oldax)
  i1,i2 = searchsortedfirst(oldax,val1,rev=rev), searchsortedlast(oldax,val2,rev=rev)
  return length(newax)/(i2-i1+1)
end

  function expandchunksizes(c,newaxdict)
    k = collect(keys(newaxdict))
    chunkold = eachchunk(c.data)
    axold = caxes(c)
    axinds  = map(i->findAxis(i,c),k)
    steprats = map((inew,iold)->getsteprat(newaxdict[inew], axold[iold].values),k,axinds)
    newcs = ntuple(ndims(c)) do i
        ii = findfirst(isequal(i),axinds)
        round(Int,approx_chunksize(chunkold.chunks[i]) * (ii==nothing ? 1 : steprats[ii]))
    end
    newco = ntuple(ndims(c)) do i
        ii = findfirst(isequal(i),axinds)
        round(Int,grid_offset(chunkold.chunks[i]) * (ii==nothing ? 1 : steprats[ii]))
    end
    return GridChunks(size(c),newcs, offset = newco)
  end

function interpolatecube(c,
  newaxes::Dict,
  newchunks = expandchunksizes(c,newaxes);
  order=Dict()
  )
  ii = map(axn->(axn,findAxis(axn,c)),collect(keys(newaxes)))
  oo = ntuple(ndims(c)) do i
    ai = findfirst(j->j[2]==i,ii)
    if ai === nothing
      nothing,NoInterp(),nothing
    else
      oldvals = caxes(c)[i].values
      newvals = newaxes[ii[ai][1]]
      getinterpinds(oldvals, newvals),get(order,ii[ai][1],Constant()),newvals
    end
  end
  newinds = getindex.(oo,1)
  intorder = getindex.(oo,2)
  @show intorder
  ar = InterpolatedDiskArray(c.data,newchunks,newinds..., order = intorder)
  newvals = getindex.(oo,3)
  newax = map(caxes(c),newvals) do ax,val
    if val === nothing
      ax
    else
      RangeAxis(axname(ax),val)
    end
  end
  YAXArray(newax,ar,c.properties, cleaner = c.cleaner)
end

function getinterpinds(oldvals::AbstractRange, newvals::AbstractRange)
  (newvals.-first(oldvals))./step(oldvals).+1
end
function getinterpinds(r1,r2)
  rev = issorted(r1) ? false : issorted(r1,rev=true) ? true : error("Axis values are not sorted")
  map(r2) do ir
    ii = searchsortedfirst(r1,ir,rev=rev)
    ii1 = max(min(ii-1,length(r1)),1)
    ii2 = max(min(ii,length(r1)),1)
    ind = if ii1 == ii2
      Float64(ii1)
    else
      ii1+(ir-r1[ii1])/(r1[ii2]-r1[ii1])
    end
    ind
  end
end


text:
using Morbit
using Test

#%%
#test_files = ["optimization.jl", "batch_objectives.jl",]

test_files = ["taylor_models.jl", "lagrange_models.jl", "model_derivatives.jl"]
for fn in test_files
    include( joinpath(@__DIR__, fn) )   # prefixe with path for PackageCompiler
end


text:
module FrostHelperEntityMover

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/EntityMover" EntityMoverFrostHelper(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, types::String="", blacklist::Bool=false, easing::String="CubeInOut", moveDuration::Number=1.0, mustCollide::Bool=true, relativeMovementMode::Bool=false, pauseTimeLength::Number=0.0, onEndSFX::String="", startPauseTimeLength::Number=0.0, nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "Entity Mover (Frost Helper)" => Ahorn.EntityPlacement(
        EntityMoverFrostHelper,
        "rectangle"
    ),
)

Ahorn.nodeLimits(entity::EntityMoverFrostHelper) = 1, 1
Ahorn.minimumSize(entity::EntityMoverFrostHelper) = 8, 8
Ahorn.resizable(entity::EntityMoverFrostHelper) = true, true

function Ahorn.selection(entity::EntityMoverFrostHelper)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    nodes = get(entity.data, "nodes", ())
    if isempty(nodes)
        return Ahorn.Rectangle(x, y, width, height)

    else
        nx, ny = Int.(nodes[1])
        return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
    end
end

edgeColor = (38, 0, 0, 128) ./ 255
centerColor = (64, 0, 0, 10) ./ 255

function renderSpaceJam(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number)
    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    Ahorn.drawRectangle(ctx, x, y, width, height, edgeColor, centerColor)

    Ahorn.restore(ctx)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::EntityMoverFrostHelper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        renderSpaceJam(ctx, nx, ny, width, height)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::EntityMoverFrostHelper, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderSpaceJam(ctx, 0, 0, width, height)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper