text:
import Base: convert, ==

# import JuLIP: energy, virial, site_energies, forces
# import JuLIP.Potentials: evaluate

using JuLIP: JVec, JMat, chemical_symbols

export OneBody, MOneBody

"""
`mutable struct OneBody{T}  <: NBodyFunction{1}`

this should not normally be constructed by a user, but instead E0 should be
passed to the relevant lsq functions, which will construct it.
"""
mutable struct OneBody{T} <: AbstractCalculator
   E0::T
end

@pot OneBody

evaluate(V::OneBody) = V.E0

site_energies(V::OneBody, at::AbstractAtoms) = fill(V(), length(at))

energy(V::OneBody, at::AbstractAtoms; domain = 1:length(at)) = length(domain) * V()
forces(V::OneBody, at::AbstractAtoms{T}; kwargs...) where {T} = zeros(JVec{T}, length(at))
virial(V::OneBody, at::AbstractAtoms{T}; kwargs...) where {T} = zero(JMat{T})
site_energy(V::OneBody, at::AbstractAtoms, i0::Integer) = V()
site_energy_d(V::OneBody, at::AbstractAtoms{T}, i0::Integer) where {T} =
      zeros(JVec{T}, length(at))

Dict(V::OneBody) = Dict("__id__" => "OneBody", "E0" => V.E0)
OneBody(D::Dict) = OneBody(D["E0"])
convert(::Val{:OneBody}, D::Dict) = OneBody(D)

==(V1::OneBody, V2::OneBody) = (V1.E0 == V2.E0)

import Base: *
*(c::Real, V::OneBody) = OneBody(V.E0 * c)
*(V::OneBody, c::Real) = c * V

"""
`mutable struct MOneBody{T}  <: NBodyFunction{1}`

this should not normally be constructed by a user, but instead E0 should be
passed to the relevant lsq functions, which will construct it.
This structure deals with multi-species configurations.
"""
mutable struct MOneBody{T <: AbstractFloat} <: AbstractCalculator
   E0::Dict{Symbol, T}
end

@pot MOneBody

MOneBody(args...) = MOneBody(Dict(args...))

evaluate(V::MOneBody, sp) = V.E0[sp]

function site_energies(V::MOneBody, at::AbstractAtoms{T}) where {T}
   E = zeros(T, length(at))
   for i in 1:length(at)
      E[i] = V(chemical_symbols(at)[i])
   end
   return E
end

energy(V::MOneBody, at::AbstractAtoms) = sum(site_energies(V, at))

forces(V::MOneBody, at::AbstractAtoms{T}) where {T} = zeros(JVec{T}, length(at))

virial(V::MOneBody, at::AbstractAtoms{T}) where {T} = zero(JMat{T})

Dict(V::MOneBody) = Dict("__id__" => "MOneBody", "E0" => V.E0)

convert_str_2_symb(D::Dict{Symbol,T}) where {T} = D #already in the correct form

# convert the E0 Dict{String} read from JSON into a Dict{Symbol}
MOneBody(D::Dict{String}, T=Float64) =
   Dict( Symbol(key) => T(val) for (key, val) in D ) |> MOneBody

convert(::Val{:MOneBody}, D::Dict) = MOneBody(D["E0"])

==(V1::MOneBody, V2::MOneBody) = (V1.E0 == V2.E0)


text:
struct WedderburnDecomposition{B,iV,DS<:DirectSummand,Hom}
    basis::B
    invariants::iV
    Uπs::Vector{DS}
    hom::Hom
end

"""
    WedderburnDecomposition([T::Type, ]G::Group, action::Action, basis_full, basis_half[, S; semisimple])
Compute `WedderburnDecomposition` related to `G` acting on `basis_half`, representing a form of Wedderburn-Artin decomposition.

This object is intended to be used to simplify problems of positive-semidefinite optimization.
# Arguments
 * `basis_full` corresponds to the basis (or index set) for the objective functional (the set of constraints);
 * `basis_half` corresponds to the basis (or index set) of the PSD constraint.

For description of the remaining arguments see [symmetry_adapted_basis](@ref).

Return a `wd::WedderburnDecomposition` object which defines:
 * `basis(wd)`: the original `basis_full`.
 * `invariant_vectors(wd)`: a basis for the subspace `⟨basis_full⟩^G` invariant under the action of `G`.
 * `direct_summands(wd)`: a vector of `DirectSummands` defining a map from `⟨basis_half⟩` to a direct (orthogonal) sum of subspaces
 * `diagonalize(M::AbstractMatrix, wd::WedderburnDecomposition)`: a map implementing the Wedderburn-Artin decomposition for matrices `M ∈ End(⟨basis_half⟩)` (i.e. with rows and columns indexed by the elements of `basis_half`).

See also: [symmetry_adapted_basis](@ref).
"""
function WedderburnDecomposition(
    T::Type,
    G::Group,
    action::Action,
    basis_full,
    basis_half,
    S = Rational{Int};
    semisimple = false,
)
    tbl = CharacterTable(S, G)
    ehom = CachedExtensionHomomorphism(G, action, basis_half, precompute = true)

    Uπs = symmetry_adapted_basis(T, tbl, ehom; semisimple = semisimple)

    basis = StarAlgebras.Basis{UInt32}(basis_full)
    invariants = invariant_vectors(tbl, action, basis)

    return WedderburnDecomposition(basis, invariants, Uπs, ehom)
end

function Base.show(io::IO, wbdec::SymbolicWedderburn.WedderburnDecomposition)
    ds = direct_summands(wbdec)
    simple = all(issimple.(ds))
    dims = size.(ds, 1)
    norbs = length(invariant_vectors(wbdec))

    print(io, "Wedderburn Decomposition into $norbs orbits and $(length(ds))")
    all(simple) && print(io, " simple")
    println(io, " summands of sizes")
    return print(io, dims)
end

invariant_vectors(wbdec::WedderburnDecomposition) = wbdec.invariants
StarAlgebras.basis(wbdec::WedderburnDecomposition) = wbdec.basis
direct_summands(wbdec::WedderburnDecomposition) = wbdec.Uπs
Base.eltype(wbdec::WedderburnDecomposition) =
    eltype(eltype(direct_summands(wbdec)))

_tmps(wbdec::WedderburnDecomposition) =
    zeros.(eltype(wbdec), size.(direct_summands(wbdec)))

_fillfactor(M::AbstractMatrix) = count(!iszero, M) / length(M)
_fillfactor(M::AbstractSparseMatrix) = nnz(M) / length(M)

function diagonalize(
    M::AbstractMatrix,
    wbdec::WedderburnDecomposition,
    tmps = _tmps(wbdec),
)
    # return [degree(Uπ).*(Uπ*(M*Uπ')) for Uπ in summands(wbdec)]

    T = eltype(eltype(direct_summands(wbdec)))
    Mπs = [(d = size(Uπ, 1); zeros(T, d, d)) for Uπ in direct_summands(wbdec)]
    return diagonalize!(Mπs, M, direct_summands(wbdec), tmps)
end

function diagonalize!(
    Mπs,
    M::AbstractMatrix,
    wbdec::WedderburnDecomposition,
    tmps = _tmps(wbdec),
)
    return diagonalize!(Mπs, M, direct_summands(wbdec), tmps)
end

function diagonalize!(
    Mπs::AbstractVector{<:AbstractMatrix},
    M::AbstractMatrix,
    Uπs::AbstractVector{<:DirectSummand},
    tmps,
)
    @assert length(Mπs) == length(Uπs)

    for (π, Uπ) in enumerate(Uπs)
        imUπ = convert(Matrix, image_basis(Uπ)) # Base.Matrix to allow BLAS paths below
        LinearAlgebra.mul!(tmps[π], imUπ, M)
        LinearAlgebra.mul!(Mπs[π], tmps[π], imUπ')
        zerotol!(Mπs[π], atol = eps(eltype(imUπ)) * max(size(imUπ)...))
        Mπs[π] .*= degree(Uπ)
    end
    return Mπs
end

function invariant_vectors(
    tbl::Characters.CharacterTable,
    act::Action,
    basis::StarAlgebras.Basis,
)
    triv_χ = Characters.Character{Rational{Int}}(Characters.trivial_character(tbl))
    ehom =
        CachedExtensionHomomorphism(parent(tbl), act, basis, precompute = true)
    # ehom = ExtensionHomomorphism(act, basis)

    mpr = matrix_projection_irr(ehom, triv_χ)
    mpr, pivots = row_echelon_form!(mpr)
    img = mpr[1:length(pivots), :]

    # change the format of invariant_vectors to image_basis(ehom, trχ)
    return sparsevec.(eachrow(img))
end

function invariant_vectors(
    tbl::Characters.CharacterTable,
    act::ByPermutations,
    basis::StarAlgebras.Basis{T,I},
) where {T,I}
    G = parent(tbl)
    tovisit = trues(length(basis))
    invariant_vs = Vector{SparseVector{Rational{Int}}}()

    ordG = order(Int, G)
    elts = collect(G)
    orbit = zeros(I, ordG)

    for i in eachindex(basis)
        if tovisit[i]
            bi = basis[i]
            Threads.@threads for j in eachindex(elts)
                orbit[j] = basis[action(act, elts[j], bi)]
            end
            tovisit[orbit] .= false
            push!(
                invariant_vs,
                sparsevec(orbit, fill(1 // ordG, ordG), length(basis)),
            )
        end
    end
    return invariant_vs
end

ext: jl


text:
"""
    @with_kw mutable struct RemoteImage
        imageDigest::Union{Missing, String} = missing
        imageTag::Union{Missing, String} = missing
        ecr_repo::Union{Missing, ECRRepo} = missing
        exists::Bool = true
    end

Represents a docker image stored in an AWS ECR Repository. The `exists` attribute indicates whether
the RemoteImage still exists.
"""
@with_kw mutable struct RemoteImage <: LambdaComponent
  imageDigest::Union{Missing, String} = missing
  imageTag::Union{Missing, String} = missing
  ecr_repo::Union{Missing, ECRRepo} = missing
  exists::Bool = true
end
StructTypes.StructType(::Type{RemoteImage}) = StructTypes.Mutable()
Base.:(==)(a::RemoteImage, b::RemoteImage) = a.imageDigest == b.imageDigest

"""
    function get_all_remote_images(jot_generated_only::Bool = true)::Vector{RemoteImage}

Returns all remote images stored on AWS ECR. By default, filters for jot-generated images only.
"""
function get_all_remote_images(jot_generated_only::Bool = true)::Vector{RemoteImage}
  repos = get_all_ecr_repos()
  remote_images = [img for repo in repos for img in get_remote_images(repo)]
  jot_generated_only ? filter(is_jot_generated, remote_images) : remote_images
end

function get_remote_images(repo::ECRRepo)::Vector{RemoteImage}
  @debug repo
  get_images_script = get_images_in_ecr_repo_script(repo)
  @debug get_images_script
  images_json = readchomp(`bash -c $get_images_script`)
  images = JSON3.read(images_json, Dict{String, Vector{RemoteImage}})["imageIds"]
  map(images) do img
    @set img.ecr_repo = repo
  end
end

function delete!(r::RemoteImage)
  r.exists || error("Remote image does not exist")
  delete_script = get_delete_remote_image_script(r)
  output = readchomp(`bash -c $delete_script`)
  r.exists = false
  @debug get_remote_images(r.ecr_repo)
  if length(get_remote_images(r.ecr_repo)) == 0
    delete!(r.ecr_repo)
  end
  nothing
end

"""
    get_remote_image(local_image::LocalImage)::Union{Nothing, RemoteImage}

Queries AWS and returns a `RemoteImage` instance corresponding to the given `local_image`.

If multiple valid images exist, this will return the first only. If none exists, returns `nothing`.
"""
function get_remote_image(local_image::LocalImage)::Union{Nothing, RemoteImage}
  all_remote_images = get_all_remote_images()
  @debug local_image.Digest
  @debug all_remote_images
  index = findfirst(remote_image -> matches(local_image, remote_image), all_remote_images)
  @debug index
  isnothing(index) ? nothing : all_remote_images[index]
end

"""
    get_remote_image(identity::AbstractString)::Union{Nothing, RemoteImage}

Queries AWS and returns a `RemoteImage` instance corresponding to the given `identity` string.

The identity string will attempt to match on the name of the remote image, or the image's Digest.

If multiple valid images exist, this will return the first only. If none exists, returns `nothing`.
"""
function get_remote_image(identity::AbstractString)::Union{Nothing, RemoteImage}
  all_remote_images = get_all_remote_images()
  index = findfirst(all_remote_images) do ri
    name_matches = get_lambda_name(ri) == identity
    hash_matches = begin
      search_hash = split(identity, ":") |> last
      ri_hash = split(ri.imageDigest, ":") |> last
      comparison_length = minimum([length(search_hash), length(ri_hash)])
      ri_hash[begin:comparison_length] == search_hash[begin:comparison_length]
    end
    name_matches || hash_matches
  end
  isnothing(index) ? nothing : all_remote_images[index]
end


text:
@testset "Automatic Convert" begin
    x = rand()
    v = rand(3)
    @test convert(PriorMean, x) isa ConstantMean
    @test convert(PriorMean, v) isa EmpiricalMean
end


text:
const Rembedded = Ref{Bool}(false)

@static if Sys.iswindows()
    import WinReg

    function ask_yes_no_cancel(prompt::Ptr{Cchar})
        println(String(prompt))
        query = readline(STDIN)
        c = uppercase(query[1])
        local r::Cint
        r = (c=='Y' ? 1 : c=='N' ? -1 : 0)
        return r
    end

    """
        RStart

    This type mirrors `structRstart` in `R_ext/RStartup.h`. It is used to initialize the R engine.
    """
    mutable struct RStart # mirror structRstart in R_ext/RStartup.h
        R_Quiet::Cint
        R_Slave::Cint
        R_Interactive::Cint
        R_Verbose::Cint
        LoadSiteFile::Cint
        LoadInitFile::Cint
        DebugInitFile::Cint
        RestoreAction::Cint
        SaveAction::Cint
        vsize::Csize_t
        nsize::Csize_t
        max_vsize::Csize_t
        max_nsize::Csize_t
        ppsize::Csize_t
        NoRenviron::Cint
        rhome::Ptr{Cchar}
        home::Ptr{Cchar}
        ReadConsole::Ptr{Cvoid}
        WriteConsole::Ptr{Cvoid}
        CallBack::Ptr{Cvoid}
        ShowMessage::Ptr{Cvoid}
        YesNoCancel::Ptr{Cvoid}
        Busy::Ptr{Cvoid}
        CharacterMode::Cint
        WriteConsoleEx::Ptr{Cvoid}
    end
    RStart() = RStart(0,0,0,0,0,
                      0,0,0,0,0,
                      0,0,0,0,0,
                      C_NULL,C_NULL,
                      C_NULL,C_NULL,C_NULL,C_NULL,
                      C_NULL,C_NULL,2,C_NULL)

end

"""
    initEmbeddedR()

This initializes an embedded R session. It should only be called when R is not already running (e.g. if Julia is running inside an R session)
"""
function initEmbeddedR()

    # disable R signal handling
    unsafe_store!(cglobal((:R_SignalHandlers,RCall.libR),Cint),0)

    @static if Sys.iswindows()
        # TODO: Use direct Windows interface, see §8.2.2 "Calling R.dll directly"
        # of "Writing R Extensions" (aka R-exts)

        Ruser_ptr = ccall((:getRUser,libR),Ptr{Cchar},())
        Ruser = unsafe_string(Ruser_ptr)

        ccall(:_wputenv,Cint,(Cwstring,),"PATH="*ENV["PATH"]*";"*dirname(libR))
        ccall(:_wputenv,Cint,(Cwstring,),"R_USER="*Ruser)

        # otherwise R will set it itself, which can be wrong on Windows
        if !("HOME" in keys(ENV))
            ccall(:_wputenv,Cint,(Cwstring,),"HOME="*homedir())
        end

        argv = ["REmbeddedJulia","--silent","--no-save"]
        i = ccall((:Rf_initEmbeddedR,libR),Cint,(Cint,Ptr{Ptr{Cchar}}),length(argv),argv)
        if i == 0
            error("Could not start embedded R session.")
        end

        rs = RStart()
        ccall((:R_DefParams,libR),Nothing,(Ptr{RStart},), Ref(rs))

        rs.rhome          = ccall((:get_R_HOME,libR), Ptr{Cchar}, ())
        rs.home           = Ruser_ptr
        rs.ReadConsole    = cglobal((:R_ReadConsole, libR), Nothing)
        rs.CallBack       = @cfunction($event_callback, Nothing, ()).ptr
        rs.ShowMessage    = cglobal((:R_ShowMessage, libR), Nothing)
        rs.YesNoCancel    = @cfunction($ask_yes_no_cancel, Cint, (Ptr{Cchar},)).ptr
        rs.Busy           = cglobal((:R_Busy, libR), Nothing)
        rs.WriteConsole   = C_NULL
        rs.WriteConsoleEx = @cfunction($write_console_ex, Nothing, (Ptr{UInt8}, Cint, Cint)).ptr

        ccall((:R_SetParams,libR),Nothing,(Ptr{RStart},), Ref(rs))
    end

    @static if Sys.isunix()
        # set necessary environmental variables
        ENV["R_HOME"] = Rhome
        ENV["R_DOC_DIR"] = joinpath(Rhome,"doc")
        ENV["R_INCLUDE_DIR"] = joinpath(Rhome,"include")
        ENV["R_SHARE_DIR"] = joinpath(Rhome,"share")

        # initialize library
        argv = ["REmbeddedJulia","--silent","--no-save"]
        i = ccall((:Rf_initEmbeddedR,libR),Cint,(Cint,Ptr{Ptr{Cchar}}),length(argv),argv)
        if i == 0
            error("Could not start embedded R session.")
        end

        ptr_write_console_ex = @cfunction($write_console_ex,Nothing,(Ptr{UInt8},Cint,Cint)).ptr
        unsafe_store!(cglobal((:ptr_R_WriteConsole,libR),Ptr{Cvoid}), C_NULL)
        unsafe_store!(cglobal((:ptr_R_WriteConsoleEx,libR),Ptr{Cvoid}), ptr_write_console_ex)
        unsafe_store!(cglobal((:R_Consolefile,libR),Ptr{Cvoid}), C_NULL)
        unsafe_store!(cglobal((:R_Outputfile,libR),Ptr{Cvoid}), C_NULL)
        ptr_polled_events = @cfunction($polled_events,Nothing,()).ptr
        unsafe_store!(cglobal((:R_PolledEvents,libR),Ptr{Cvoid}), ptr_polled_events)
    end

    Rembedded[] = true
    atexit(endEmbeddedR)

end

"""
    endEmbeddedR()

Close embedded R session.
"""
function endEmbeddedR()
    if Rembedded[]
        ccall((:Rf_endEmbeddedR, libR),Nothing,(Cint,),0)
        Rembedded[] = false
    end
end

# for validate_libR
include(joinpath(dirname(@__FILE__),"..","deps","setup.jl"))

function __init__()
    validate_libR(libR)

    # Check if R already running
    # for some reaons, cglobal((:R_NilValue, libR)) doesn't work on rstudio/linux
    # https://github.com/Non-Contradiction/JuliaCall/issues/34
    Rinited, from_libR = try
        unsafe_load(cglobal(:R_NilValue, Ptr{Cvoid})) != C_NULL, false
    catch
        unsafe_load(cglobal((:R_NilValue, libR), Ptr{Cvoid})) != C_NULL, true
    end

    if !Rinited
        initEmbeddedR()
    end

    ip = ccall((:Rf_ScalarInteger, libR),Ptr{Cvoid},(Cint,),0)

    Const.load(from_libR)

    # set up function callbacks
    setup_callbacks()

    if !Rinited
        # print warnings as they arise
        # we can't use Rf_PrintWarnings as not exported on all platforms.
        rcall_p(:options,warn=1)

        # R gui eventloop
        isinteractive() && rgui_init()
    end

    @require AxisArrays="39de3d68-74b9-583c-8d2d-e117c070f3a9" begin
        include("convert/axisarray.jl")
    end

    # R REPL mode
    isdefined(Base, :active_repl) &&
        isinteractive() && typeof(Base.active_repl) != REPL.BasicREPL &&
            !RPrompt.repl_inited(Base.active_repl) && RPrompt.repl_init(Base.active_repl)

    # # IJulia hooks
    isdefined(Main, :IJulia) && Main.IJulia.inited && ijulia_init()
end


text:
using ImageChipper
using Test

@testset "Image Chipping" begin
    faux_image = rand( 256, 256 )
    chip_dims = ( 64, 64 )
    chipped = chip( faux_image, chip_dims, ( 0.0, 0.0 ) )
    #checkoutput dimensionality
    @test all( size(chipped)[1:2] .== chip_dims )
    @test size(chipped)[3] == 16
    @test all( chipped[:,:,1] .== faux_image[ 1:chip_dims[1], 1:chip_dims[2] ] )
    @test all( chipped[:,:,end] .== faux_image[   (end - chip_dims[1] + 1):end,
                                        (end - chip_dims[2] + 1):end ] )
    #checkoutput dimensionality
    chipped = chip( faux_image, ( 60, 60 ), ( 0.0, 0.0 ) )
    @test size(chipped)[3] .==  25
    #checkoutput dimensionality
    chipped = chip( faux_image, ( 32, 32 ), ( 0.0, 0.0 ) )
    @test size(chipped)[3] .==  64
end

@testset "Image Chipping with Boxes" begin
    using ObjectDetectionStats
    faux_image = rand( 256, 256 )
    chip_dims = ( 128, 128 )

    boxes = [   ObjectDetectionStats.Box(   1,   1,  22,  22 ),
                ObjectDetectionStats.Box(  55,  55, 222, 222 ),
                ObjectDetectionStats.Box( 222, 222, 264, 264 )   ]

    chips, boxassignments = chip_with_boxes( faux_image, boxes, chip_dims, ( 0.0, 0.0 ), 0.0 )

    @test( length(boxassignments[ 1 ]) == 2 )
    @test( length(boxassignments[ 2 ]) == 1 )
    @test( length(boxassignments[ 3 ]) == 1 )
    @test( length(boxassignments[ 4 ]) == 2 )
end


text:
# Example driver for the geometric multigrid program, GMG.jl, to solve a Poisson equation on a 2D grid using
#  second order finite difference method.
#  Dirichlet Boundary conditions

include("GMG.jl")
using .GMG: V_cycle!
using Formatting: printfmt

# Define the problem being solved. It is on a 1x1 square domain with
# homogeneous dirichlet boundary conditions
function U_exact(x,y)
  return (x^3-x)*(y^3-y)
end

function source(x,y)
  return 6*x*y*(x^2 + y^2 -2)
end

max_cycles = 10              # max number of V cycles to execute
nlevels    = 10              # number of grid levels. 1 means no multigrid, 2 means one coarse grid. etc
nx         = 1*2^(nlevels-1) # Nx and Ny are given as function of grid levels
ny         = 1*2^(nlevels-1) #

const u_ex = zeros(Float64,(nx+2,ny+2))
const u = zeros(Float64,(nx+2,ny+2))
const f = zeros(Float64,(nx+2,ny+2))
const res = zeros(Float64,(nx+2,ny+2))
const error = zeros(Float64,(nx,ny))

#calcualte the RHS and exact solution
DX=1.0/nx
DY=1.0/ny
xc=range(0.5*DX,1-0.5*DX,length=nx)
yc=range(0.5*DY,1-0.5*DY,length=ny)

#set the exact solution and the source term corresponding to it
for i= 1:nx, j=1:ny;
  u_ex[i+1,j+1] = U_exact(xc[i],yc[j])
  f[i+1,j+1]    = source(xc[i],yc[j])
end

println("-------------------------New Run------------------------------")
println("           Poisson equation on $(nx)X$(ny) grid")

for it=1:max_cycles

  u[:,:],res[:,:]=V_cycle!(nx,ny,nlevels,u,f,1)

  error[:,:]=u_ex[2:nx+1,2:ny+1]-u[2:nx+1,2:ny+1]

  printfmt( "|V cycle: {:2d}| res. norm: {:16.9e}| error norm: {:16.9e}|\n", it,findmax(abs.(res))[1],findmax(abs.(error))[1] )
end

max_stars_repo_name: AbhilashReddyM/Multigrid-in-Julia

max_issues_repo_name: AbhilashReddyM/Multigrid-in-Julia

max_forks_repo_name: AbhilashReddyM/Multigrid-in-Julia


text:
"""
See [`simulate!`](@ref)
"""
function simulate end

function simulate(rng::AbstractRNG, m::MixedModel{T}, newdata; kwargs...) where {T}
    dat = Tables.columntable(newdata)
    y = zeros(T, length(first(dat)))
    return simulate!(rng, y, m, newdata; kwargs...)
end

function simulate(rng::AbstractRNG, m::MixedModel; kwargs...)
    return simulate!(rng, similar(response(m)), m; kwargs...)
end

function simulate(m::MixedModel, args...; kwargs...)
    return simulate(Random.GLOBAL_RNG, m, args...; kwargs...)
end

"""
    simulate!(rng::AbstractRNG, m::MixedModel{T}; β=m.β, σ=m.σ, θ=T[])
    simulate!(m::MixedModel; β=m.β, σ=m.σ, θ=m.θ)

Overwrite the response (i.e. `m.trms[end]`) with a simulated response vector from model `m`.

This simulation includes sampling new values for the random effects.

!!! note
    Note that `simulate!` methods with a `y::AbstractVector` as the first argument
    (besides the RNG) and `simulate` methods return the simulated response. This is
    in contrast to `simulate!` methods with a `m::MixedModel` as the first argument,
    which modify the model's response and return the entire modified model.
"""
function simulate!(
    rng::AbstractRNG, m::LinearMixedModel{T}; β=coef(m), σ=m.σ, θ=T[]
) where {T}
    # XXX should we add support for doing something with weights?
    simulate!(rng, m.y, m; β, σ, θ)
    return unfit!(m)
end

function simulate!(
    rng::AbstractRNG, m::GeneralizedLinearMixedModel{T}; β=coef(m), σ=m.σ, θ=T[]
) where {T}
    # note that these m.resp.y and m.LMM.y will later be sychronized in (re)fit!()
    # but for now we use them as distinct scratch buffers to avoid allocations

    # the noise term is actually in the GLM and not the LMM part so no noise
    # at the LMM level
    η = fill!(copy(m.LMM.y), zero(T))  # ensure that η is a vector - needed for GLM.updateμ! below
    # A better approach is to change the signature for updateμ!
    y = m.resp.y

    _simulate!(rng, y, η, m, β, σ, θ, m.resp)

    return unfit!(m)
end

"""
    _rand(rng::AbstractRNG, d::Distribution, location, scale=missing, n=1)

A convenience function taking a draw from a distrbution.

Note that `d` is specified as an existing distribution, such as
from the `GlmResp.d` field. This isn't vectorized nicely because
for distributions where the scale/dispersion is dependent on the
location (e.g. Bernoulli, Binomial, Poisson), it's not really
possible to avoid creating multiple `Distribution` objects.

Note that `n` is the `n` parameter for the Binomial distribution,
*not* the number of draws from the RNG. It is then used to change the
random draw (an integer in [0, n]) into a probability (a float in [0,1]).
"""
function _rand(rng::AbstractRNG, d::Distribution, location, scale=NaN, n=1)
    if !ismissing(scale)
        throw(ArgumentError("Families with a dispersion parameter not yet supported"))
    end

    if d isa Binomial
        dist = Binomial(Int(n), location)
    else
        dist = typeof(d)(location)
    end

    return rand(rng, dist) / n
end

function simulate!(m::MixedModel{T}; β=coef(m), σ=m.σ, θ=T[]) where {T}
    return simulate!(Random.GLOBAL_RNG, m; β, σ, θ)
end

"""
    simulate!([rng::AbstractRNG,] y::AbstractVector, m::MixedModel{T}[, newdata];
                    β = coef(m), σ = m.σ, θ = T[], wts=m.wts)
    simulate([rng::AbstractRNG,] m::MixedModel{T}[, newdata];
                    β = coef(m), σ = m.σ, θ = T[], wts=m.wts)

Simulate a new response vector, optionally overwriting a pre-allocated vector.

New data can be optionally provided in tabular format.

This simulation includes sampling new values for the random effects. Thus in
contrast to `predict`, there is no distinction in between "new" and
"old" / previously observed random-effects levels.

Unlike `predict`, there is no `type` parameter for `GeneralizedLinearMixedModel`
because the noise term in the model and simulation is always on the response
scale.

The `wts` argument is currently ignored except for `GeneralizedLinearMixedModel`
models with a `Binomial` distribution.

!!! warning
    Models are assumed to be full rank.

!!! note
    Note that `simulate!` methods with a `y::AbstractVector` as the first argument
    (besides the RNG) and `simulate` methods return the simulated response. This is
    in contrast to `simulate!` methods with a `m::MixedModel` as the first argument,
    which modify the model's response and return the entire modified model.
"""
function simulate!(
    rng::AbstractRNG,
    y::AbstractVector,
    m::LinearMixedModel,
    newdata::Tables.ColumnTable;
    β=m.β,
    σ=m.σ,
    θ=m.θ,
)
    # the easiest thing here is to just assemble a new model and
    # pass that to the other simulate methods....
    # this can probably be made much more efficient
    # (for one thing, this still allocates for the model's response)
    # note that the contrasts get copied over with the formula
    # (as part of the applied schema)
    # contr here are the fast Grouping contrasts
    f, contr = _abstractify_grouping(m.formula)
    mnew = LinearMixedModel(f, newdata; contrasts=contr)
    # XXX why not do simulate!(rng, y, mnew; β=β, σ=σ, θ=θ)
    # instead of simulating the model and then copying?
    # Well, it turns out that the call to randn!(rng, y)
    # gives different results at the tail end of the array
    # for y <: view(::Matrix{Float64}, :, 3) than y <: Vector{Float64}
    # I don't know why, but this doesn't actually incur an
    # extra computation and gives consistent results at the price
    # of an allocationless copy
    simulate!(rng, mnew; β, σ, θ)
    return copy!(y, mnew.y)
end

function simulate!(
    rng::AbstractRNG, y::AbstractVector, m::LinearMixedModel{T}; β=m.β, σ=m.σ, θ=m.θ
) where {T}
    length(β) == length(fixef(m)) ||
        length(β) == length(coef(m)) ||
        throw(ArgumentError("You must specify all (non-singular) βs"))

    β = convert(Vector{T}, β)
    σ = T(σ)
    θ = convert(Vector{T}, θ)
    isempty(θ) || setθ!(m, θ)

    if length(β) ≠ length(coef(m))
        padding = length(coef(m)) - length(β)
        for ii in 1:padding
            push!(β, -0.0)
        end
    end

    # initialize y to standard normal
    randn!(rng, y)

    # add the unscaled random effects
    for trm in m.reterms
        unscaledre!(rng, y, trm)
    end

    # scale by σ and add fixed-effects contribution
    return mul!(y, m.X, β, one(T), σ)
end

function simulate!(
    rng::AbstractRNG,
    y::AbstractVector,
    m::GeneralizedLinearMixedModel,
    newdata::Tables.ColumnTable;
    β=m.β,
    σ=m.σ,
    θ=m.θ,
)
    # the easiest thing here is to just assemble a new model and
    # pass that to the other simulate methods....
    # this can probably be made much more efficient
    # (for one thing, this still allocates for the model's response)
    # note that the contrasts get copied over with the formula
    # (as part of the applied schema)
    # contr here are the fast Grouping contrasts
    f, contr = _abstractify_grouping(m.formula)
    mnew = GeneralizedLinearMixedModel(f, newdata, m.resp.d, Link(m.resp); contrasts=contr)
    # XXX why not do simulate!(rng, y, mnew; β, σ, θ)
    # instead of simulating the model and then copying?
    # Well, it turns out that the call to randn!(rng, y)
    # gives different results at the tail end of the array
    # for y <: view(::Matrix{Float64}, :, 3) than y <: Vector{Float64}
    # I don't know why, but this doesn't actually incur an
    # extra computation and gives consistent results at the price
    # of an allocationless copy
    simulate!(rng, mnew; β, σ, θ)
    return copy!(y, mnew.y)
end

function simulate!(
    rng::AbstractRNG,
    y::AbstractVector,
    m::GeneralizedLinearMixedModel{T};
    β=m.β,
    σ=m.σ,
    θ=m.θ,
) where {T}
    # make sure both scratch arrays are init'd to zero
    η = zeros(T, size(y))
    copyto!(y, η)
    return _simulate!(rng, y, η, m, β, σ, θ)
end

function _simulate!(
    rng::AbstractRNG,
    y::AbstractVector,
    η::AbstractVector,
    m::GeneralizedLinearMixedModel{T},
    β,
    σ,
    θ,
    resp=nothing,
) where {T}
    length(β) == length(fixef(m)) ||
        length(β) == length(coef(m)) ||
        throw(ArgumentError("You must specify all (non-singular) βs"))

    dispersion_parameter(m) ||
        ismissing(σ) ||
        throw(
            ArgumentError(
                "You must not specify a dispersion parameter for model families without a dispersion parameter",
            ),
        )

    β = convert(Vector{T}, β)
    if σ !== missing
        σ = T(σ)
    end
    θ = convert(Vector{T}, θ)

    d = m.resp.d

    if length(β) ≠ length(coef(m))
        padding = length(coef(m)) - length(β)
        for ii in 1:padding
            push!(β, -0.0)
        end
    end

    fast = (length(m.θ) == length(m.optsum.final))
    setpar! = fast ? setθ! : setβθ!
    params = fast ? θ : vcat(β, θ)
    setpar!(m, params)

    lm = m.LMM

    # assemble the linear predictor

    # add the unscaled random effects
    # note that unit scaling may not be correct for
    # families with a dispersion parameter
    @inbounds for trm in m.reterms
        unscaledre!(rng, η, trm)
    end

    # add fixed-effects contribution
    # note that unit scaling may not be correct for
    # families with a dispersion parameter
    mul!(η, lm.X, β, one(T), one(T))

    μ = resp === nothing ? linkinv.(Link(m), η) : GLM.updateμ!(resp, η).mu

    # convert to the distribution / add in noise
    @inbounds for (idx, val) in enumerate(μ)
        n = isempty(m.wt) ? 1 : m.wt[idx]
        y[idx] = _rand(rng, d, val, σ, n)
    end

    return y
end

function simulate!(rng::AbstractRNG, y::AbstractVector, m::MixedModel, newdata; kwargs...)
    return simulate!(rng, y, m, Tables.columntable(newdata); kwargs...)
end
function simulate!(y::AbstractVector, m::MixedModel, newdata; kwargs...)
    return simulate!(Random.GLOBAL_RNG, y, m, Tables.columntable(newdata); kwargs...)
end

"""
    unscaledre!(y::AbstractVector{T}, M::ReMat{T}) where {T}
    unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, M::ReMat{T}) where {T}

Add unscaled random effects simulated from `M` to `y`.

These are unscaled random effects (i.e. they incorporate λ but not σ) because
the scaling is done after the per-observation noise is added as a standard normal.
"""
function unscaledre! end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,S}) where {T,S}
    return mul!(y, A, vec(lmul!(A.λ, randn(rng, S, nlevs(A)))), one(T), one(T))
end

function unscaledre!(rng::AbstractRNG, y::AbstractVector{T}, A::ReMat{T,1}) where {T}
    return mul!(y, A, lmul!(first(A.λ), randn(rng, nlevs(A))), one(T), one(T))
end

unscaledre!(y::AbstractVector, A::ReMat) = unscaledre!(Random.GLOBAL_RNG, y, A)


text:
@testset "Common" begin
    @test get_code_center_frequency_ratio(GPSL1()) ≈ 1/1540
end

@testset "Code spectra $(get_system_string(system))" for system = [GPSL1(), GPSL5(), BOCcos(GPSL1(),0,15)]
    @test get_code_spectrum(system, 0) ≈ 1.0Hz/get_code_frequency(system)
    @testset "Test $(m). zero" for m = 1:10
        @test get_code_spectrum(system, m * get_code_frequency(system)) == 0
        @test get_code_spectrum(system, -m * get_code_frequency(system)) == 0
    end
    @test sum(get_code_spectrum.(system, -1e12:1e4:1e12)) * 1e4 ≈ 1 rtol = 1e-5
end

@testset "Code spectra $(get_system_string(system))" for system = [GalileoE1B(), BOCcos(GPSL1(),15,2.5)]
    @test get_code_spectrum(system, 0) == 0.0
    @test sum(get_code_spectrum.(system, -1e12:1e4:1e12)) * 1e4 ≈ 1 rtol = 1e-5
end

@testset "Code generation $(get_system_string(system))" for system = [GalileoE1B(), GPSL1(), GPSL5()]
    sampling_rate = 25e6Hz
    samples = 1000
    code = zeros(Int16, samples)
    code = gen_code!(code, system, 1, sampling_rate, get_code_frequency(system), 0)
    phase = (0:length(code) - 1) * get_code_frequency(system) / sampling_rate
    @test code == get_code.(system, phase, 1)
    @test code == gen_code(samples, system, 1, sampling_rate, get_code_frequency(system), 0)
end

@testset "Code generation with start_phase bigger than code_length" begin
    system = GPSL1()
    sampling_rate = 2.5e6Hz
    samples = 2500
    code = zeros(Int16, samples)
    code = gen_code!(code, system, 1, sampling_rate, get_code_frequency(system), 2065)
    phase = (0:samples - 1) * get_code_frequency(system) / sampling_rate .+ 2065
    @test code == get_code.(system, phase, 1)
    @test code == gen_code(samples, system, 1, sampling_rate, get_code_frequency(system), 2065)
end

@testset "Code generation $(get_system_string(system)) with different index" for system = [GalileoE1B(), GPSL1(), GPSL5()]
    sampling_rate = 25e6Hz
    samples = 4002
    code = zeros(Int16, samples)
    code = gen_code!(code, system, 1, sampling_rate, get_code_frequency(system), 0.0, -1)
    phase = (-1:4000) * get_code_frequency(system) / sampling_rate
    @test code == get_code.(system, phase, 1)
    @test code == gen_code(samples, system, 1, sampling_rate, get_code_frequency(system), 0.0, -1)
end


text:
########
# Set up parameters
########
parameters = (
    ρₒ = 1.0,  # reference density
    cₛ = 1e-2, # sound speed
    ℓᵐ = 10,   # jet thickness, (larger is thinner)
    ℓ = 20,   # perturbation thickness, (larger is thinner)
    m = 2,    # number of sign changes on equator for perturbation
    ϕᵖ = π / 2 * 0.05, # of centerdness of perturbation
    ϵ = 0.3,  # perturbation amplitude
    vˢ = 5e-4, # velocity scale
    α = 2e-4,
    Ω = 1e-3,
    λ_coupler = (500 / 60 / 86400), #(L_airsea / τ_airsea)
)

########
# Set up inital condition
########
uᵐ(p, λ, ϕ, r) = p.ℓᵐ * sech(p.ℓᵐ * ϕ)^2
vᵐ(p, λ, ϕ, r) = 0.0
hᵐ(p, λ, ϕ, r) = 0.0

u1(p, λ, ϕ, r) = p.ℓ * 2 * (ϕ - p.ϕᵖ) * exp(-p.ℓ * (ϕ - p.ϕᵖ)^2) * cos(ϕ) * cos(2 * (ϕ - p.ϕᵖ)) * sin(p.m * λ)
u2(p, λ, ϕ, r) = exp(-p.ℓ * (ϕ - p.ϕᵖ)^2) * sin(ϕ) * cos(2 * (ϕ - p.ϕᵖ)) * sin(p.m * λ)
u3(p, λ, ϕ, r) = 2 * exp(-p.ℓ * (ϕ - p.ϕᵖ)^2) * cos(ϕ) * sin(2 * (ϕ - p.ϕᵖ)) * sin(p.m * λ)
uᵖ(p, λ, ϕ, r) = u1(p, λ, ϕ, r) + u2(p, λ, ϕ, r) + u3(p, λ, ϕ, r)
vᵖ(p, λ, ϕ, r) = p.m * exp(-p.ℓ * (ϕ - p.ϕᵖ)^2) * cos(2 * (ϕ - p.ϕᵖ)) * cos(p.m * λ)
hᵖ(p, λ, ϕ, r) = 0.0

ρ₀(p, λ, ϕ, r) = p.ρₒ
ρuʳᵃᵈ(p, λ, ϕ, r) = 0
ρuˡᵃᵗ(p, λ, ϕ, r) = p.vˢ * ρ₀(p, λ, ϕ, r) * (p.ϵ * vᵖ(p, λ, ϕ, r))
ρuˡᵒⁿ(p, λ, ϕ, r) = p.vˢ * ρ₀(p, λ, ϕ, r) * (uᵐ(p, λ, ϕ, r) + p.ϵ * uᵖ(p, λ, ϕ, r))
ρθ₀B(p, λ, ϕ, r) = ρ₀(p, λ, ϕ, r) * tanh(p.ℓᵐ * ϕ)
ρθ₀A(p, λ, ϕ, r) = 0

# Cartesian Representation for initialization (boiler plate really)
ρ₀ᶜᵃʳᵗ(p, x...) = ρ₀(p, lon(x...), lat(x...), rad(x...))
ρu⃗₀ᶜᵃʳᵗ(p, x...) = (
    ρuʳᵃᵈ(p, lon(x...), lat(x...), rad(x...)) * r̂(x...) +
    ρuˡᵃᵗ(p, lon(x...), lat(x...), rad(x...)) * ϕ̂(x...) +
    ρuˡᵒⁿ(p, lon(x...), lat(x...), rad(x...)) * λ̂(x...)
)
ρθ₀ᶜᵃʳᵗA(p, x...) = ρθ₀A(p, lon(x...), lat(x...), rad(x...))
ρθ₀ᶜᵃʳᵗB(p, x...) = ρθ₀B(p, lon(x...), lat(x...), rad(x...))

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module LinearMaps

export LinearMap
export ⊗, kronsum, ⊕

using LinearAlgebra
using SparseArrays

if VERSION < v"1.2-"
    import Base: has_offset_axes
    require_one_based_indexing(A...) = !has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    import Base: require_one_based_indexing
end

abstract type LinearMap{T} end

const MapOrMatrix{T} = Union{LinearMap{T},AbstractMatrix{T}}

Base.eltype(::LinearMap{T}) where {T} = T

abstract type MulStyle end

struct FiveArg <: MulStyle end
struct ThreeArg <: MulStyle end

MulStyle(::FiveArg, ::FiveArg) = FiveArg()
MulStyle(::ThreeArg, ::FiveArg) = ThreeArg()
MulStyle(::FiveArg, ::ThreeArg) = ThreeArg()
MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg()
MulStyle(::LinearMap) = ThreeArg() # default
@static if VERSION ≥ v"1.3.0-alpha.115"
    MulStyle(::AbstractMatrix) = FiveArg()
else
    MulStyle(::AbstractMatrix) = ThreeArg()
end
MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...))

Base.isreal(A::LinearMap) = eltype(A) <: Real
LinearAlgebra.issymmetric(::LinearMap) = false # default assumptions
LinearAlgebra.ishermitian(A::LinearMap{<:Real}) = issymmetric(A)
LinearAlgebra.ishermitian(::LinearMap) = false # default assumptions
LinearAlgebra.isposdef(::LinearMap) = false # default assumptions

Base.ndims(::LinearMap) = 2
Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
Base.length(A::LinearMap) = size(A)[1] * size(A)[2]

# check dimension consistency for y = A*x and Y = A*X
function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector)
     m, n = size(A)
     (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!"))
     return nothing
 end
 function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)
     m, n = size(A)
     (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
     return nothing
 end

# conversion of AbstractMatrix to LinearMap
convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A)
convert_to_lmaps_(A::LinearMap) = A
convert_to_lmaps() = ()
convert_to_lmaps(A) = (convert_to_lmaps_(A),)
@inline convert_to_lmaps(A, B, Cs...) =
    (convert_to_lmaps_(A), convert_to_lmaps_(B), convert_to_lmaps(Cs...)...)

Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false)
    length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
    if isone(α)
        iszero(β) && (A_mul_B!(y, A, x); return y)
        isone(β) && (y .+= A * x; return y)
        # β != 0, 1
        rmul!(y, β)
        y .+= A * x
        return y
    elseif iszero(α)
        iszero(β) && (fill!(y, zero(eltype(y))); return y)
        isone(β) && return y
        # β != 0, 1
        rmul!(y, β)
        return y
    else # α != 0, 1
        iszero(β) && (A_mul_B!(y, A, x); rmul!(y, α); return y)
        isone(β) && (y .+= rmul!(A * x, α); return y)
        # β != 0, 1
        rmul!(y, β)
        y .+= rmul!(A * x, α)
        return y
    end
end
# the following is of interest in, e.g., subspace-iteration methods
function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false)
    (size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
    @inbounds @views for i = 1:size(X, 2)
        mul!(Y[:, i], A, X[:, i], α, β)
    end
    # starting from Julia v1.1, we could use the `eachcol` iterator
    # for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
    #     mul!(Yi, A, Xi, α, β)
    # end
    return Y
end

A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector)  = mul!(y, A, x)
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)

include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
include("transpose.jl") # transposing linear maps
include("linearcombination.jl") # defining linear combinations of linear maps
include("composition.jl") # composition of linear maps
include("functionmap.jl") # using a function as linear map
include("blockmap.jl") # block linear maps
include("kronecker.jl") # Kronecker product of linear maps
include("conversion.jl") # conversion of linear maps to matrices

"""
    LinearMap(A; kwargs...)
    LinearMap(J, M::Int)
    LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)

Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
with the purpose of redefining its properties via the keyword arguments `kwargs`;
a `UniformScaling` object `J` with specified (square) dimension `M`; or
from a function or callable object `f`. In the latter case, one also needs to specify
the size of the equivalent matrix representation `(M, N)`, i.e. for functions `f` acting
on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably,
also the `eltype` `T` of the corresponding matrix representation needs to be specified, i.e.
whether the action of `f` on a vector will be similar to e.g. multiplying by numbers of type `T`.
If not specified, the devault value `T=Float64` will be assumed. Optionally, a corresponding
function `fc` can be specified that implements the transpose/adjoint of `f`.

The keyword arguments and their default values for functions `f` are
*   issymmetric::Bool = false : whether `A` or `f` acts as a symmetric matrix
*   ishermitian::Bool = issymmetric & T<:Real : whether `A` or `f` acts as a Hermitian matrix
*   isposdef::Bool = false : whether `A` or `f` acts as a positive definite matrix.
For existing linear maps or matrices `A`, the default values will be taken by calling
`issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.

For functions `f`, there is one more keyword argument
*   ismutating::Bool : flags whether the function acts as a mutating matrix multiplication
    `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
    or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
    The default value is guessed by looking at the number of arguments of the first occurence
    of `f` in the method table.
"""
LinearMap(A::MapOrMatrix; kwargs...) = WrappedMap(A; kwargs...)
LinearMap(J::UniformScaling, M::Int) = UniformScalingMap(J.λ, M)
LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)
LinearMap(f, fc, M::Int; kwargs...) = LinearMap{Float64}(f, fc, M; kwargs...)
LinearMap(f, fc, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, fc, M, N; kwargs...)

LinearMap{T}(A::MapOrMatrix; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)

end # module


text:
export mutate

"Mutate: Change weights"
function mutate_weights(indiv::NEATInd, cfg::Dict)
    # TODO: check original weight mutation
    ind = NEATInd(indiv)
    for c in ind.connections
        if rand() < cfg["p_mut_weights"]
            c.weight = c.weight + randn()*cfg["weight_factor"]
        end
    end
    ind
end

"split a connection by adding a neuron"
function split_connection(c::Connection, neuron_nb::Int, cfg::Dict)
    # Create neuron
    n = Neuron(0.0, 0.0, false)

    # Create 2 new connections
    # TODO: neuron number does not correspond to order for recurrence checks
    c1 = Connection(c.in_node, neuron_nb, 1.0, true, cfg["innovation_max"]+1)
    c2 = Connection(neuron_nb, c.out_node, c.weight, true, cfg["innovation_max"]+2)
    cfg["innovation_max"] += 2

    c1, c2, n
end

"Mutation: Add a node in a connection"
function mutate_add_neuron(indiv::NEATInd, cfg::Dict)
    ind = NEATInd(indiv)

    ci = rand(1:length(ind.connections))
    c = ind.connections[ci]
    c1, c2, n = split_connection(c, length(ind.neurons), cfg)
    deleteat!(ind.connections, ci)
    append!(ind.connections, [c1, c2])
    push!(ind.neurons, n)

    ind
end

"Mutation: Add a connection between 2 random nodes"
function mutate_add_connection(indiv::NEATInd, cfg::Dict, allow_recurrent::Bool=false)
    ind = NEATInd(indiv)

    # Select 2 different neurons at random, can't create a connection that already exists
    nn = length(ind.neurons)
    valid = trues(nn, nn)
    for c in ind.connections
        valid[c.in_node, c.out_node] = false
    end
    conns = findall(valid)
    if length(conns) > 0
        shuffle!(conns)
        cfg["innovation_max"] += 1
        c = Connection(conns[1][1], conns[1][2], rand_weight(), true, cfg["innovation_max"])
        push!(ind.connections, c)
    end

    ind
end

"Mutation: Toggle Enabled"
function mutate_enabled(indiv::NEATInd; n_times=1)
    ind = NEATInd(indiv)
    ids = randperm(length(ind.connections))
    for i in 1:n_times
        ind.connections[ids[i]].enabled = !ind.connections[ids[i]].enabled
    end
    ind
end

"mutate(ind::NEATInd, cfg::Dict): return a new mutated individual"
function mutate(indiv::NEATInd, cfg::Dict)
    if rand() < cfg["p_mutate_add_neuron"]
        return mutate_add_neuron(indiv, cfg)
    elseif rand() < cfg["p_mutate_add_connection"]
        return mutate_add_connection(indiv, cfg)
    elseif rand() < cfg["p_mutate_weights"]
        return mutate_weights(indiv, cfg)
    elseif rand() < cfg["p_mutate_enabled"]
        return mutate_enabled(indiv, cfg)
    end
    # return clone if no mutation occurs
    NEATInd(indiv)
end


text:
using ITensors
using ITensorsVisualization
using Test

@testset "ITensorsVisualization.jl" begin
  N = 10
  s(n) = Index([QN("Sz", 0) => 1, QN("Sz", 1) => 1]; tags = "S=1/2,Site,n=$n")
  l(n) = Index([QN("Sz", 0) => 10, QN("Sz", 1) => 10]; tags = "Link,l=$n")
  h(n) = Index([QN("Sz", 0) => 5, QN("Sz", 1) => 5]; tags = "ham,Link,l=$n")
  s⃗ = [s(n) for n in 1:N]
  l⃗ = [l(n) for n in 1:N-1]
  h⃗ = [h(n) for n in 1:N-1]

  # Add some more indices between two of the tensors
  x = Index([QN("Sz", 0) => 2]; tags = "X")
  y = Index([QN("Sz", 0) => 2]; tags = "Y")

  n = 2
  ψn1n2 = randomITensor(l⃗[n-1], s⃗[n], s⃗[n+1], l⃗[n+1], dag(x), dag(y))
  hn1 = randomITensor(dag(h⃗[n-1]), s⃗[n]', dag(s⃗[n]), h⃗[n], x, y)
  hn2 = randomITensor(dag(h⃗[n]), s⃗[n+1]', dag(s⃗[n+1]), h⃗[n+1])
  ELn0 = randomITensor(l⃗[n-1]', h⃗[n-1], dag(l⃗[n-1]))
  ERn2 = randomITensor(l⃗[n+1]', dag(h⃗[n+1]), dag(l⃗[n+1]))

  R = @visualize ELn0 * ψn1n2 * hn1 * hn2 * ERn2
  @test R ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

  # Split it up into multiple contractions
  R1 = @visualize ELn0 * ψn1n2 * hn1
  R2 = @visualize R1 * hn2 * ERn2 labels = ["T1", "T2", "T3"]
  @test R2 ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

  @test_throws ErrorException @visualize R1 * hn2 * ERn2 labels = ["T1", "T2"]
end

max_stars_repo_name: mtfishman/ITensorsVisualization

max_issues_repo_name: mtfishman/ITensorsVisualization

max_forks_repo_name: mtfishman/ITensorsVisualization


text:
typealias LinearPipeline{T<:ImageTransformation} Vector{T}

function _transform(ops, N, types, img)
    ex = :(img_0 = img)
    for i in 1:N
        types[i] <: Any || throw(ArgumentError("Pipline must consist of ImageTransformation objects"))
        vin  = symbol("img_$(i-1)")
        vout = symbol("img_$i")
        ex = quote
            $ex
            $vout = transform(ops[$i], $vin)
        end
    end
    vlast = symbol("img_$N")
    ex = :($ex; $vlast)
    ex
end

@generated function transform{N,T}(ops::NTuple{N}, img::T)
    @inbounds res = _transform(ops, N, ops.parameters, img)
    res
end

function transform{I<:AbstractImage}(pl::LinearPipeline, imgs::AbstractArray{I})
    imgs_out = similar(imgs)
    pl_tuple = (pl...)
    for iter in eachindex(imgs)
        img = imgs[iter]
        T = typeof(img)
        imgs_out[iter] = transform(pl_tuple, img)::T
    end
    imgs_out
end

function transform{T}(pl::LinearPipeline, img::T)
    transform((pl...), img)::T
end


text:
mutable struct MinimaxObjective <: AbstractObjective
    name::String
    sense::String
    cores::Vector{MinimaxObjectiveCore}

    MinimaxObjective(cores) = new("minimax", "min", cores)
    MinimaxObjective() = new("minimax", "min", MinimaxObjectiveCore[])
end


text:
# TESTS FOR ALL 3D UNITCELLS
@testset "3D Unitcells" begin

    #----------------------
    #   CUBIC
    #----------------------
    @testset "Cubic" begin

        # test Default Wrapper
        @test_nowarn getUnitcellCubic();
        @test_nowarn getUnitcellSC();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellCubic()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}
        @test typeof(getUnitcellSC()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # test Default Wrapper getUnitcell(...)
        id = :cubic
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellCubic(String);
        @test_nowarn getUnitcellCubic(Symbol);
        @test_nowarn getUnitcellCubic(Int64, Symbol);
        @test_nowarn getUnitcellSC(String);
        @test_nowarn getUnitcellSC(Symbol);
        @test_nowarn getUnitcellSC(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellCubic(1);
        @test_nowarn getUnitcellCubic(String, 1);
        @test_nowarn getUnitcellCubic(String, Symbol, 1);
        @test_nowarn getUnitcellSC(1);
        @test_nowarn getUnitcellSC(String, 1);
        @test_nowarn getUnitcellSC(String, Symbol, 1);

    end

    #----------------------
    #   DIAMOND
    #----------------------
    @testset "Diamond" begin

        # test Default Wrapper
        @test_nowarn getUnitcellDiamond();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellDiamond()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # test Default Wrapper getUnitcell(...)
        id = :diamond
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellDiamond(String);
        @test_nowarn getUnitcellDiamond(Symbol);
        @test_nowarn getUnitcellDiamond(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellDiamond(1);
        @test_nowarn getUnitcellDiamond(String, 1);
        @test_nowarn getUnitcellDiamond(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellDiamond(2);
        @test_nowarn getUnitcellDiamond(String, 2);
        @test_nowarn getUnitcellDiamond(String, Symbol, 2);

    end

    #----------------------
    #   FCC
    #----------------------
    @testset "FCC" begin

        # test Default Wrapper
        @test_nowarn getUnitcellFCC();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellFCC()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # test Default Wrapper getUnitcell(...)
        id = :fcc
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellFCC(String);
        @test_nowarn getUnitcellFCC(Symbol);
        @test_nowarn getUnitcellFCC(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellFCC(1);
        @test_nowarn getUnitcellFCC(String, 1);
        @test_nowarn getUnitcellFCC(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellFCC(2);
        @test_nowarn getUnitcellFCC(String, 2);
        @test_nowarn getUnitcellFCC(String, Symbol, 2);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellFCC(4);
        @test_nowarn getUnitcellFCC(String, 4);
        @test_nowarn getUnitcellFCC(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcellFCC(4))))) == 3

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellFCC(5);
        @test_nowarn getUnitcellFCC(String, 5);
        @test_nowarn getUnitcellFCC(String, Symbol, 5);

    end

    #----------------------
    #   PYROCHLORE
    #----------------------
    @testset "Pyrochlore" begin

        # test Default Wrapper
        @test_nowarn getUnitcellPyrochlore();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellPyrochlore()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # test Default Wrapper getUnitcell(...)
        id = :pyrochlore
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellPyrochlore(String);
        @test_nowarn getUnitcellPyrochlore(Symbol);
        @test_nowarn getUnitcellPyrochlore(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellPyrochlore(1);
        @test_nowarn getUnitcellPyrochlore(String, 1);
        @test_nowarn getUnitcellPyrochlore(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellPyrochlore(2);
        @test_nowarn getUnitcellPyrochlore(String, 2);
        @test_nowarn getUnitcellPyrochlore(String, Symbol, 2);

    end

    #----------------------
    #   HYPERKAGOME
    #----------------------
    @testset "Hyperkagome" begin

        # test Default Wrapper
        @test_nowarn getUnitcellHyperkagome();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellHyperkagome()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # test Default Wrapper getUnitcell(...)
        id = :pyrochlore
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellHyperkagome(String);
        @test_nowarn getUnitcellHyperkagome(Symbol);
        @test_nowarn getUnitcellHyperkagome(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellHyperkagome(1);
        @test_nowarn getUnitcellHyperkagome(String, 1);
        @test_nowarn getUnitcellHyperkagome(String, Symbol, 1);

    end

    #----------------------
    #   (8,3)a
    #----------------------
    @testset "(8,3)a" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_8_3_a();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_8_3_a()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_8_3_a(String);
        @test_nowarn getUnitcell_8_3_a(Symbol);
        @test_nowarn getUnitcell_8_3_a(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_a(1);
        @test_nowarn getUnitcell_8_3_a(String, 1);
        @test_nowarn getUnitcell_8_3_a(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_a(4);
        @test_nowarn getUnitcell_8_3_a(String, 4);
        @test_nowarn getUnitcell_8_3_a(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_8_3_a(4))))) == 3

    end

    #----------------------
    #   (8,3)b
    #----------------------
    @testset "(8,3)b" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_8_3_b();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_8_3_b()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_8_3_b(String);
        @test_nowarn getUnitcell_8_3_b(Symbol);
        @test_nowarn getUnitcell_8_3_b(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_b(1);
        @test_nowarn getUnitcell_8_3_b(String, 1);
        @test_nowarn getUnitcell_8_3_b(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_b(4);
        @test_nowarn getUnitcell_8_3_b(String, 4);
        @test_nowarn getUnitcell_8_3_b(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_8_3_b(4))))) == 3

    end

    #----------------------
    #   (8,3)c
    #----------------------
    @testset "(8,3)c" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_8_3_c();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_8_3_c()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_8_3_c(String);
        @test_nowarn getUnitcell_8_3_c(Symbol);
        @test_nowarn getUnitcell_8_3_c(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_c(1);
        @test_nowarn getUnitcell_8_3_c(String, 1);
        @test_nowarn getUnitcell_8_3_c(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_c(4);
        @test_nowarn getUnitcell_8_3_c(String, 4);
        @test_nowarn getUnitcell_8_3_c(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_8_3_c(4))))) == 3

    end

    #----------------------
    #   (8,3)n
    #----------------------
    @testset "(8,3)n" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_8_3_n();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_8_3_n()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_8_3_n(String);
        @test_nowarn getUnitcell_8_3_n(Symbol);
        @test_nowarn getUnitcell_8_3_n(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_n(1);
        @test_nowarn getUnitcell_8_3_n(String, 1);
        @test_nowarn getUnitcell_8_3_n(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_8_3_n(4);
        @test_nowarn getUnitcell_8_3_n(String, 4);
        @test_nowarn getUnitcell_8_3_n(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_8_3_n(4))))) == 3

    end

    #----------------------
    #   (9,3)a
    #----------------------
    @testset "(9,3)a" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_9_3_a();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_9_3_a()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_9_3_a(String);
        @test_nowarn getUnitcell_9_3_a(Symbol);
        @test_nowarn getUnitcell_9_3_a(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_9_3_a(1);
        @test_nowarn getUnitcell_9_3_a(String, 1);
        @test_nowarn getUnitcell_9_3_a(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_9_3_a(2);
        @test_nowarn getUnitcell_9_3_a(String, 2);
        @test_nowarn getUnitcell_9_3_a(String, Symbol, 2);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_9_3_a(4);
        @test_nowarn getUnitcell_9_3_a(String, 4);
        @test_nowarn getUnitcell_9_3_a(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_9_3_a(4))))) == 3

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_9_3_a(5);
        @test_nowarn getUnitcell_9_3_a(String, 5);
        @test_nowarn getUnitcell_9_3_a(String, Symbol, 5);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_9_3_a(5))))) == 3

    end

    #----------------------
    #   (10,3)a
    #----------------------
    @testset "(10,3)a" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_10_3_a();
        @test_nowarn getUnitcellHyperoctagon();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_10_3_a()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}
        @test typeof(getUnitcellHyperoctagon()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_10_3_a(String);
        @test_nowarn getUnitcell_10_3_a(Symbol);
        @test_nowarn getUnitcell_10_3_a(Int64, Symbol);
        @test_nowarn getUnitcellHyperoctagon(String);
        @test_nowarn getUnitcellHyperoctagon(Symbol);
        @test_nowarn getUnitcellHyperoctagon(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_a(1);
        @test_nowarn getUnitcell_10_3_a(String, 1);
        @test_nowarn getUnitcell_10_3_a(String, Symbol, 1);
        @test_nowarn getUnitcellHyperoctagon(1);
        @test_nowarn getUnitcellHyperoctagon(String, 1);
        @test_nowarn getUnitcellHyperoctagon(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_a(4);
        @test_nowarn getUnitcell_10_3_a(String, 4);
        @test_nowarn getUnitcell_10_3_a(String, Symbol, 4);
        @test_nowarn getUnitcellHyperoctagon(4);
        @test_nowarn getUnitcellHyperoctagon(String, 4);
        @test_nowarn getUnitcellHyperoctagon(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_10_3_a(4))))) == 3
        @test length(unique(label.(bonds(getUnitcellHyperoctagon(4))))) == 3

    end

    #----------------------
    #   (10,3)b
    #----------------------
    @testset "(10,3)b" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_10_3_b();
        @test_nowarn getUnitcellHyperhoneycomb();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_10_3_b()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}
        @test typeof(getUnitcellHyperhoneycomb()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_10_3_b(String);
        @test_nowarn getUnitcell_10_3_b(Symbol);
        @test_nowarn getUnitcell_10_3_b(Int64, Symbol);
        @test_nowarn getUnitcellHyperhoneycomb(String);
        @test_nowarn getUnitcellHyperhoneycomb(Symbol);
        @test_nowarn getUnitcellHyperhoneycomb(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_b(1);
        @test_nowarn getUnitcell_10_3_b(String, 1);
        @test_nowarn getUnitcell_10_3_b(String, Symbol, 1);
        @test_nowarn getUnitcellHyperhoneycomb(1);
        @test_nowarn getUnitcellHyperhoneycomb(String, 1);
        @test_nowarn getUnitcellHyperhoneycomb(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_b(2);
        @test_nowarn getUnitcell_10_3_b(String, 2);
        @test_nowarn getUnitcell_10_3_b(String, Symbol, 2);
        @test_nowarn getUnitcellHyperhoneycomb(2);
        @test_nowarn getUnitcellHyperhoneycomb(String, 2);
        @test_nowarn getUnitcellHyperhoneycomb(String, Symbol, 2);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_b(4);
        @test_nowarn getUnitcell_10_3_b(String, 4);
        @test_nowarn getUnitcell_10_3_b(String, Symbol, 4);
        @test_nowarn getUnitcellHyperhoneycomb(4);
        @test_nowarn getUnitcellHyperhoneycomb(String, 4);
        @test_nowarn getUnitcellHyperhoneycomb(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_10_3_b(4))))) == 3
        @test length(unique(label.(bonds(getUnitcellHyperhoneycomb(4))))) == 3

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_b(5);
        @test_nowarn getUnitcell_10_3_b(String, 5);
        @test_nowarn getUnitcell_10_3_b(String, Symbol, 5);
        @test_nowarn getUnitcellHyperhoneycomb(5);
        @test_nowarn getUnitcellHyperhoneycomb(String, 5);
        @test_nowarn getUnitcellHyperhoneycomb(String, Symbol, 5);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_10_3_b(5))))) == 3
        @test length(unique(label.(bonds(getUnitcellHyperhoneycomb(5))))) == 3

    end

    #----------------------
    #   (10,3)c
    #----------------------
    @testset "(10,3)c" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_10_3_c();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_10_3_c()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_10_3_c(String);
        @test_nowarn getUnitcell_10_3_c(Symbol);
        @test_nowarn getUnitcell_10_3_c(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_c(1);
        @test_nowarn getUnitcell_10_3_c(String, 1);
        @test_nowarn getUnitcell_10_3_c(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_c(4);
        @test_nowarn getUnitcell_10_3_c(String, 4);
        @test_nowarn getUnitcell_10_3_c(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_10_3_c(4))))) == 3

    end

    #----------------------
    #   (10,3)d
    #----------------------
    @testset "(10,3)d" begin

        # test Default Wrapper
        @test_nowarn getUnitcell_10_3_d();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcell_10_3_d()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,3}, LB,B<:AbstractBond{LB,3}}

        # Test Wrapper with different label types
        @test_nowarn getUnitcell_10_3_d(String);
        @test_nowarn getUnitcell_10_3_d(Symbol);
        @test_nowarn getUnitcell_10_3_d(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_d(1);
        @test_nowarn getUnitcell_10_3_d(String, 1);
        @test_nowarn getUnitcell_10_3_d(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcell_10_3_d(4);
        @test_nowarn getUnitcell_10_3_d(String, 4);
        @test_nowarn getUnitcell_10_3_d(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcell_10_3_d(4))))) == 3

    end

end


text:
# SINGLE NODE/CORE VERSION

import Pkg; Pkg.instantiate()

include("/users/yh31/scratch/projects/gigaword_64k/src/gigaword_64k.jl")

const path_intermed_data = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/intermediate_files"

const output_dir = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/reductions"

reduce_jdfs_and_save(path_intermed_data, output_dir, ns_tuple = (70_000, 64_000))


text:
module SeparatingAxisTheorem2D

using LinearAlgebra
using StaticArrays
using CoordinateTransformations
using RecipesBase
import StaticArrays: SUnitRange

export Shape2D, AxisAlignedBoundingBox, AABB, LineSegment, Polygon, Triangle, Circle, CompoundShape
export intersecting, inflate, sweep, sweep_intersecting

# Utilities
   norm2(v)    = v⋅v
    perp(v)    = SVector(v[2], -v[1])
 project(v, w) = w*(v⋅w)/(w⋅w)
projectN(v, w) = w*(v⋅w)
 reflect(v, a) = v - 2*project(v, a)
reflectN(v, a) = v - 2*projectN(v, a)
  rotate(v, θ) = ((s, c) = sincos(θ); SVector(v[1]*c - v[2]*s, v[1]*s + v[2]*c))
  ccw(a, b, c) = (b - a)×(c - a) > 0

projectNextrema(pts, n)    = SVector{2}(extrema(p⋅n for p in pts))
voronoi_region(p, v)       = (d = p⋅v; ifelse(d < 0, -1, ifelse(d > v⋅v, 1, 0)))
overlapping(I1, I2)        = I1[1] <= I2[2] && I2[1] <= I1[2]
ininterval(x, I)           = I[1] <= x <= I[2]
enclosing_interval(I1, I2) = SVector(min(I1[1], I2[1]), max(I1[2], I2[2]))
minmaxV((x,y))             = SVector(minmax(x, y))
minmaxV(x, y)              = SVector(minmax(x, y))
wrap1(i, N)                = ifelse(i < 1, N, ifelse(i > N, 1, i))

  circshift1(v::StaticVector{N}) where {N} = pushfirst(v[SUnitRange(1,N-1)], v[end])
uncircshift1(v::StaticVector{N}) where {N} = push(v[SUnitRange(2,N)], v[1])
@generated function unrolled_any(f, ::Val{N}) where {N}
    N == 0 ? false : :(Base.@_inline_meta; $(foldr((x, y) -> :($x || $y), [:(f($i)) for i in 1:N])))
end
@generated function unrolled_all(f, ::Val{N}) where {N}
    N == 0 ? true  : :(Base.@_inline_meta; $(foldr((x, y) -> :($x && $y), [:(f($i)) for i in 1:N])))
end
@inline unrolled_any(f, tup::NTuple{N,Any}) where {N} = unrolled_any(i -> (Base.@_inline_meta; f(tup[i])), Val(N))
@inline unrolled_all(f, tup::NTuple{N,Any}) where {N} = unrolled_all(i -> (Base.@_inline_meta; f(tup[i])), Val(N))

# Shape Definitions
const Point = AbstractVector{<:Number}
abstract type Shape2D{T} end
abstract type AbstractPolygon{T} <: Shape2D{T} end
## BoundingBox
struct AxisAlignedBoundingBox{T} <: AbstractPolygon{T}
    bounds::SMatrix{2,2,T,4}
end
AxisAlignedBoundingBox((xl, xu), (yl, yu)) = AxisAlignedBoundingBox(SMatrix{2,2}(xl, yl, xu, yu))
AxisAlignedBoundingBox(x::Number, y::Number) = AxisAlignedBoundingBox(SMatrix{2,2}(-x/2, -y/2, x/2, y/2))
const AABB{T} = AxisAlignedBoundingBox{T}
Base.getindex(B::AABB, i, j=:) = B.bounds[i,j]
getAABB(aabb::AABB) = aabb
getAABB(X::Shape2D) = X.aabb
getAABB((x, y)::Point) = AABB((x, x), (y, y))

## Line Segments
abstract type AbstractLineSegment{T} <: AbstractPolygon{T} end
### SimpleLineSegment
struct SimpleLineSegment{T} <: AbstractLineSegment{T}
    v   ::SVector{2,T}
    w   ::SVector{2,T}
    edge::SVector{2,T}
    aabb::AABB{T}

    SimpleLineSegment{T}(v, w) where {T} = new(v, w, w - v, AABB(minmax(v[1], w[1]), minmax(v[2], w[2])))
end
SimpleLineSegment(v, w) = (T = eltype(v); SimpleLineSegment{T}(SVector{2,T}(v), SVector{2,T}(w)))

### LineSegment
struct LineSegment{T} <: AbstractLineSegment{T}
    v     ::SVector{2,T}
    w     ::SVector{2,T}
    edge  ::SVector{2,T}
    normal::SVector{2,T}
    ndotv ::T
    aabb  ::AABB{T}

    function LineSegment{T}(v, w; normalize_normals=false) where {T}
        edge   = w - v
        normal = normalize_normals ? normalize(perp(edge)) : perp(edge)
        ndotv  = normal⋅v
        aabb   = AABB(minmax(v[1], w[1]), minmax(v[2], w[2]))
        new(v, w, edge, normal, ndotv, aabb)
    end
end
LineSegment(v, w) = (T = eltype(v); LineSegment{T}(SVector{2,T}(v), SVector{2,T}(w)))

## Polygon
struct Polygon{N,T} <: AbstractPolygon{T}
    points  ::SVector{N,SVector{2,T}}
    edges   ::SVector{N,SVector{2,T}}
    normals ::SVector{N,SVector{2,T}}
    nextrema::SVector{N,SVector{2,T}}
    aabb    ::AABB{T}

    function Polygon{N,T}(points; normalize_normals=false) where {N,T}
        edges    = uncircshift1(points) .- points
        normals  = normalize_normals ? normalize.(perp.(edges)) : perp.(edges)
        nextrema = projectNextrema.(Ref(points), normals)
        @assert last.(nextrema) ≈ dot.(points, normals) "Polygon must be convex"
        aabb     = AABB(extrema(first.(points)), extrema(last.(points)))
        new(points, edges, normals, nextrema, aabb)
    end
end
Polygon(points) = (N = length(points); T = eltype(eltype(points)); Polygon{N,T}(SVector{N,SVector{2,T}}(points)))
Polygon(points...) = Polygon(points)
Polygon(B::AABB) = Polygon(SVector.(B[1][SVector(1,2,2,1)], B[2][SVector(1,1,2,2)]))

## Triangle
Triangle((P1, P2, P3)) = ccw(P1, P2, P3) ? Polygon(P1, P2, P3) : Polygon(P1, P3, P2)
Triangle(points...) = Triangle(points)

## Circle
struct Circle{T} <: Shape2D{T}
    c::SVector{2,T}
    r::T
    r2::T
    aabb::AABB{T}

    function Circle{T}(c, r) where {T}
        @assert r > 0 "Circle radius $r must be positive"
        aabb = AABB((c[1] - r, c[1] + r), (c[2] - r, c[2] + r))
        new(c, r, r*r, aabb)
    end
end
Circle(c, r) = (T = eltype(c); Circle{T}(SVector{2,T}(c), T(r)))
Circle(r) = (T = eltype(r); Circle{T}(zeros(SVector{2,T}), r))

## Compound Shapes
struct CompoundShape{PS,T} <: Shape2D{T}
    parts::PS
    aabb ::AABB{T}

    function CompoundShape{PS,T}(parts) where {PS,T}
        aabb = AABB((minimum(getAABB(P)[1,1] for P in parts), maximum(getAABB(P)[1,2] for P in parts)),
                    (minimum(getAABB(P)[2,1] for P in parts), maximum(getAABB(P)[2,2] for P in parts)))
        new(parts, aabb)
    end
end
CompoundShape(parts::Shape2D{T}...) where {T} = CompoundShape{typeof(parts),T}(parts)
CompoundShape(parts::Vector{S}) where {T,S<:Shape2D{T}} = CompoundShape{Vector{S},T}(parts)

# Projecting Shapes onto Axes
projectNextrema(B::AABB,                n) = minmaxV(n[1].*B[1]) + minmaxV(n[2].*B[2])
projectNextrema(L::AbstractLineSegment, n) = minmaxV(L.v⋅n, L.w⋅n)
projectNextrema(P::Polygon,             n) = projectNextrema(P.points, n)
projectNextrema(C::Circle,              n) = (d = C.c⋅n; r = C.r*norm(n); SVector(d - r, d + r))

# Separating Axis
## General Method (unused?)
separating_axis(S1::Shape2D, S2::Shape2D, a) = !overlapping(projectNextrema(S1, a), projectNextrema(S2, a))
## Specialized Methods (axis corresponding to first shape)
separating_axis(L::LineSegment, S::Shape2D) = !ininterval(L.ndotv, projectNextrema(S, L.normal))
separating_axis(P::Polygon, S::Shape2D, i::Integer) = !overlapping(P.nextrema[i], projectNextrema(S, P.normals[i]))

# Intersection Checking
intersecting(X, Y) = !AABBseparated(X, Y) && _intersecting(X, Y)
intersecting(X, Y, f) = intersecting(X, f(Y))
## Special Cases (Point/AABB)
intersecting(B1::AABB, B2::AABB)  = overlapping(B1[1], B2[1]) && overlapping(B1[2], B2[2])
intersecting( p::Point, B::AABB)  = ininterval(p[1], B[1]) && ininterval(p[2], B[2])
intersecting( B::AABB,  p::Point) = intersecting(p, B)

## Broadphase
AABBseparated(S1::Shape2D, S2::Shape2D) = !intersecting(getAABB(S1), getAABB(S2))
AABBseparated( p::Point,    S::Shape2D) = !intersecting(p, getAABB(S))
AABBseparated( S::Shape2D,  p::Point)   = AABBseparated(p, S)

## Narrowphase - Point
_intersecting(p::Point, L::LineSegment{T}) where {T} = (x = p - L.v; e = L.edge; abs(x×e) < eps(T) && 0 <= x⋅e <= e⋅e)
# _intersecting(p::Point, P::Polygon) = all(ininterval(p⋅n, I) for (n, I) in zip(P.normals, P.nextrema))
_intersecting(p::Point, P::Polygon{N}) where {N} = unrolled_all(i -> ininterval(p⋅P.normals[i], P.nextrema[i]), Val(N))
_intersecting(p::Point, C::Circle) = norm2(p - C.c) <= C.r2
_intersecting(S::Union{LineSegment,Polygon,Circle}, p::Point) = _intersecting(p, S)

## Narrowphase - AbstractPolygon (AABB/LineSegment/Polygon)
_intersecting(X::AbstractPolygon, Y::AbstractPolygon) = !_any_separating_axis(X, Y) && !_any_separating_axis(Y, X)
_any_separating_axis(B::AABB,        S::Shape2D) = false # established by !AABBseparated by this point
_any_separating_axis(L::LineSegment, S::Shape2D) = separating_axis(L, S)
# _any_separating_axis(P::Polygon{N},  S::Shape2D) where {N} = any(separating_axis(P, S, i) for i in 1:N)
_any_separating_axis(P::Polygon{N},  S::Shape2D) where {N} = unrolled_any(i -> separating_axis(P, S, i), Val(N))

## Narrowphase - AbstractLineSegment (SimpleLineSegment, LineSegment)
_intersecting(L1::AbstractLineSegment, L2::AbstractLineSegment) = ccw(L1.v, L2.v, L2.w) != ccw(L1.w, L2.v, L2.w) &&
                                                                  ccw(L1.v, L1.w, L2.v) != ccw(L1.v, L1.w, L2.w)

## Narrowphase - Circle
_intersecting(C::Circle, B::AABB) = _intersecting(clamp.(C.c, B[:,1], B[:,2]), C)
function _intersecting(C::Circle, L::AbstractLineSegment)
    _intersecting(L.v, C) && return true
    _intersecting(L.w, C) && return true
    vc = C.c - L.v
    d2 = norm2(L.edge)
    d2*C.r2 < (L.edge×vc)^2 && return false
    0 <= vc⋅L.edge <= d2
end
function _intersecting(C::Circle, P::Polygon{N}) where {N}
    for i in 1:N
        pi2c = C.c - P.points[i]
        vr = voronoi_region(pi2c, P.edges[i])
        if vr == 0
            pi2c⋅P.normals[i] > C.r*norm(P.normals[i]) && return false
        else
            j = wrap1(i + vr, N)
            pj2c = C.c - P.points[j]
            vr*voronoi_region(pj2c, P.edges[j]) < 0 && norm2(ifelse(vr > 0, pj2c, pi2c)) > C.r2 && return false
        end
    end
    true
end
_intersecting(X::AbstractPolygon, C::Circle) = _intersecting(C, X)
_intersecting(C1::Circle, C2::Circle) = norm2(C1.c - C2.c) <= (C1.r + C2.r)^2

## Compound Shape Narrowphase
_intersecting(X::Union{Point,Shape2D}, S::CompoundShape) = any(P -> intersecting(X, P), S.parts)
_intersecting(X::Union{Point,Shape2D}, S::CompoundShape{<:Tuple}) = unrolled_any(P -> intersecting(X, P), S.parts)
_intersecting(S::CompoundShape, X::Union{Point,AbstractPolygon,Circle}) = _intersecting(X, S)

# Transformations
## Inflation
inflate(p::Point, ε; round_corners=true) = Circle(p, ε)
inflate(B::AABB,  ε; round_corners=true) = inflate(Polygon(B), ε, round_corners=round_corners)
function inflate(L::AbstractLineSegment, ε; round_corners=true)
    n = normalize(L isa LineSegment ? L.normal : perp(L.edge))
    if round_corners
        CompoundShape(Polygon(SVector(L.v + ε*n, L.w + ε*n, L.w - ε*n, L.v - ε*n)), Circle(L.v, ε), Circle(L.w, ε))
    else
        m = perp(n)
        Polygon(SVector(L.v + ε*(n + m), L.w + ε*(n - m), L.w + ε*(-n - m), L.v + ε*(-n + m)))
    end
end
function inflate(P::Polygon, ε; round_corners=true)
    push_out_corner_vector(n0, n1) = (c = n0×n1; abs(c) < 1e-6 ? n0 : (perp(n1) - perp(n0))/c)
    normals = normalize.(P.normals)
    !round_corners && return Polygon(P.points .+ ε.*push_out_corner_vector.(circshift1(normals), normals))
    CompoundShape(
        Polygon(vec([P.points .+ ε.*circshift1(normals) P.points .+ ε.*normals]')),
        Circle.(P.points, ε)...
    )
end
inflate(C::Circle,        ε; round_corners=true) = Circle(C.c, C.r + ε)
inflate(S::CompoundShape, ε; round_corners=true) = CompoundShape(inflate.(S.parts, ε, round_corners=round_corners)...)

## Affine Transformation
(f::Translation)(S::Shape2D) = transform(f, S)    # julia#14919 --> AbstractAffineMap
(f::LinearMap)(S::Shape2D)   = transform(f, S)
(f::AffineMap)(S::Shape2D)   = transform(f, S)
transform(f::Translation,                B::AABB)          = AABB(B.bounds .+ f.translation)
transform(f::Union{LinearMap,AffineMap}, B::AABB)          = transform(f, Polygon(B))
transform(f::AbstractAffineMap,          L::LineSegment)   = LineSegment(f(L.v), f(L.w))
transform(f::AbstractAffineMap,          P::Polygon)       = Polygon(f.(P.points))
transform(f::Translation,                C::Circle)        = Circle(f(C.c), C.r)
transform(f::Union{LinearMap,AffineMap}, C::Circle)        = Circle(f(C.c), norm(f.linear[:,1]))    # isometries only
transform(f::AbstractAffineMap,          S::CompoundShape) = CompoundShape(f.(S.parts)...)

## Sweep
sweep(X, f1, f2) = sweep(f1(X), f2(X))
sweep(v::Point, w::Point) = LineSegment(v, w)
sweep(B1::AABB, B2::AABB) = sweep(Polygon(B1), Polygon(B2))
function sweep(L1::AbstractLineSegment, L2::AbstractLineSegment; check_degenerate=false)
    if intersecting(L1, L2)
        λ = [L1.edge L2.edge]\(L2.v - L1.v)
        p = L1.v + λ[1]*L1.edge
        if ccw(L1.v, L2.v, p)
            CompoundShape(Polygon(L1.v, L2.v, p), Polygon(L1.w, L2.w, p))
        else
            CompoundShape(Polygon(L2.v, L1.v, p), Polygon(L2.w, L1.w, p))
        end
    else
        if check_degenerate
            V, W = SimpleLineSegment(L1.v, L2.v), SimpleLineSegment(L1.w, L2.w)
            intersecting(V, W) && return sweep(V, W)
        end
        a = ccw(L1.v, L1.w, L2.w)
        b = ccw(L1.w, L2.w, L2.v)
        c = ccw(L2.w, L2.v, L1.v)
        d = ccw(L2.v, L1.v, L1.w)
        if a == c
            CompoundShape(Triangle(L1.v, L1.w, L2.w), Triangle(L2.w, L2.v, L1.v))
        else # b == d
            CompoundShape(Triangle(L1.w, L2.w, L2.v), Triangle(L2.v, L1.v, L1.w))
        end
    end
end
sweep(P1::Polygon, P2::Polygon) = CompoundShape(P1, sweep.(SimpleLineSegment.(P1.points, uncircshift1(P1.points)),
                                                           SimpleLineSegment.(P2.points, uncircshift1(P2.points)))...)
sweep(C1::Circle,  C2::Circle)  = inflate(SimpleLineSegment(C1.c, C2.c), C1.r, round_corners=true)
sweep(S1::CompoundShape, S2::CompoundShape) = CompoundShape(sweep.(S1.parts, S2.parts)...)

# Sweep Intersection Checking
combined_AABB(X, Y) = (BX = getAABB(X); BY = getAABB(Y); AABB([min.(BX[:,1], BY[:,1]) max.(BX[:,2], BY[:2])]))
## X Static, Y Dynamic
sweep_intersecting(X, Y, f1, f2) = sweep_intersecting(X, f1(Y), f2(Y))
sweep_intersecting(X, Y1, Y2) = !AABBseparated(X, combined_AABB(Y1, Y2)) && _sweep_intersecting(X, Y1, Y2)

_sweep_intersecting(X, Y1, Y2) = intersecting(X, sweep(Y1, Y2))

## X and Y Dynamic
sweep_intersecting(X, fX1, fX2, Y, fY1, fY2) = sweep_intersecting(X, Y, inv(fX1) ∘ fY1, inv(fX2) ∘ fY2)

# Distance Computation
## TODO: port over from old MotionPlanning.jl code

# Plot Recipes
@recipe function f(B::AABB; dims=(1,2))
    seriestype :=  :shape
    fillcolor  --> :match
    linecolor  --> :match
    label      --> ""
    x, y = dims
    coords = (B[1][SVector(1,2,2,1)], B[2][SVector(1,1,2,2)])
    coords[x], coords[y]
end

@recipe function f(L::AbstractLineSegment; dims=(1,2))
    seriestype :=  :shape
    fillcolor  --> :match
    linecolor  --> :match
    label      --> ""
    x, y = dims
    SVector(L.v[x], L.w[x]), SVector(L.v[y], L.w[y])
end

@recipe function f(P::Polygon; dims=(1,2))
    seriestype :=  :shape
    fillcolor  --> :match
    linecolor  --> :match
    label      --> ""
    x, y = dims
    getindex.(P.points, x), getindex.(P.points, y)
end

@recipe function f(C::Circle; dims=(1,2), n=64)
    seriestype :=  :shape
    fillcolor  --> :match
    linecolor  --> :match
    label      --> ""
    x, y = dims
    θ = range(-π, stop=π, length=n)
    C.c[x] .+ C.r*cos.(θ), C.c[y] .+ C.r*sin.(θ)
end

@recipe f(S::CompoundShape; dims=(1,2)) = (dims --> dims; foreach(P -> @series(begin P end), S.parts))

end # module


text:
using Nettle

function Base.trunc(s::AbstractString, n::Integer)
    n > 0 || throw(DomainError())
    l = length(s)
    l > n || return s
    n > 3 || return s[1:n]
    return s[1:n-3] * "..."
end

tests = [""    => "d41d8cd98f00b204e9800998ecf8427e",
         "a"   => "0cc175b9c0f1b6a831c399e269772661",
         "abc" => "900150983cd24fb0d6963f7d28e17f72",
         "message digest" => "f96b697d7cb7938d525a2f31aaf161d0",
         "abcdefghijklmnopqrstuvwxyz" => "c3fcd3d76192e4007dfb496cca67e13b",
         "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789" =>
         "d174ab98d277d9f5a5611c2c9f419d9f",
         "12345678901234567890123456789012345678901234567890123456789012345678901234567890" =>
         "57edf4a22be3c955ac49da2e2107b67a",
         "foobad" => "3858f62230ac3c915f300c664312c63f"]

println("Testing Julia's MD5 hash against RFC 1321.")
for (k, h) in sort(tests, by = length ∘ first)
    md5sum = hexdigest("md5", k)
    @printf("%20s → %s ", trunc(k, 15), md5sum)
    if md5sum == h
        println("MD5 OK")
    else
        println("MD5 Bad")
        println("* The sum should be  ", h)
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=

  Repeated measures of IQ in Turing.jl

  From WebPPL model
  https://mhtess.github.io/bdappl/chapters/03-simpleModels.html
  """
  The data are the measures xijx_{ij}xij
     for the i=1,...,ni = 1, . . . ,ni=1,...,n people
  and their j=1,...,mj = 1, . . . ,mj=1,...,m repeated test scores.

  We assume that the differences in repeated test scores are distributed as
  Gaussian error terms with zero mean and unknown precision. The mean of the
  Gaussian of a person’s test scores corresponds to their latent true IQ.
  This will be different for each person. The standard deviation of the
  Gaussians corresponds to the accuracy of the testing instruments in
  measuring the one underlying IQ value. We assume this is the same for
  every person, since it is conceived as a property of the tests themselves.
  """

  Note: The WebPPL model has sigma ~ Uniform(0,50) and get these values, with
  a quite low sigma (~ 7)
   """
   sigma: 6.996331408573459
   p1:94.81734864858849
   p2:109.78959481757205
   p3:154.67515393166423

   credibleInterval on sigma:
   [ '0.90', [ 3.3645266702554633, 10.402949581214315 ] ]
   credibleInterval on p1:
   [ '0.90', [ 88.48060304863435, 102.16143539836233 ] ]

   credibleInterval on p2:
   [ '0.90', [ 102.73276331774267, 116.84164858686201 ] ]

   credibleInterval on p3:
   [ '0.90', [ 149.12746700987051, 161.4365469824775 ] ]
   """

  This Turing.jl model with sigma ~ Uniform(0,50) give a much larger posterior sigma (~20)
  """
  parameters       mean       std   naive_se      mcse        ess      rhat   ess_per_sec
       sigma    19.9092   13.8640     0.1386    1.1993    26.2159    1.8145       54.6165
      mus[1]    89.2099   25.0480     0.2505    0.7799   881.3664    1.0013     1836.1801
      mus[2]   105.8519   23.4151     0.2342    0.9114   533.6730    1.0020     1111.8188
      mus[3]   146.5696   22.5392     0.2254    0.8805   442.2614    1.0029      921.3778
  ...
  Credible interval for sigma with mass 0.9: (7.636206..42.867467)

  """

  So I lower the interval range of sigma and get the following:
  """
  Summary Statistics
  parameters       mean       std   naive_se      mcse       ess      rhat   ess_per_sec
      Symbol    Float64   Float64    Float64   Float64   Float64   Float64       Float64

       sigma     6.0077    2.9381     0.0294    0.2733   29.9863    1.0457       13.0205
      mus[1]    91.2915    6.6268     0.0663    0.5185   63.0714    1.1298       27.3866
      mus[2]   103.7605    6.4221     0.0642    0.4227   89.3549    1.0412       38.7993
      mus[3]   148.7040    6.3352     0.0634    0.4894   98.1541    1.0009       42.6201

  Quantiles
  parameters       2.5%      25.0%      50.0%      75.0%      97.5%
      Symbol    Float64    Float64    Float64    Float64    Float64

       sigma     2.8435     3.4661     4.1825     8.8713    11.3996
      mus[1]    76.6420    89.0146    92.6321    94.7485   105.7608
      mus[2]    91.8652   100.3283   104.2820   106.5811   119.8080
      mus[3]   137.4726   146.9825   148.9408   148.9408   162.6477

  Credible interval for sigma with mass 0.9: (3.466133..11.095075)
  Credible interval for mus[1] with mass 0.9: (78.278793..97.360210)
  Credible interval for mus[2] with mass 0.9: (91.267542..109.106586)
  Credible interval for mus[3] with mass 0.9: (140.436663..160.814064)
  """

  And it's still larger credible intervals...

=#
using Turing
include("jl_utils.jl")

@model function repeated_iq_measurements(data)
    n,_ = size(data)

    # everyone shares same sigma (corresponding to measurement error)
    sigma ~ Uniform(0, 12)

    # each person has a separate latent IQ
    mus ~ filldist(Uniform(0,200),n)
    for i in 1:n
        for j in 1:length(data[i])
            data[i,j] ~ Normal(mus[i], sigma)
        end
    end

end

data = [
   90 95 100
   105 110 115
   150 155 160]
println("data:", data)
model = repeated_iq_measurements(data)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,10), 10_000)

display(chns)

credible_interval(chns,"sigma",0.90)
credible_interval(chns,"mus[1]",0.90)
credible_interval(chns,"mus[2]",0.90)
credible_interval(chns,"mus[3]",0.90)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#Semi-Supervised GAN
# Ref: https://machinelearningmastery.com/semi-supervised-generative-adversarial-network/
# The approach is based on the definition of the semi-supervised model in the 2016 paper by Tim Salimans,
# et al.from OpenAI titled “Improved Techniques for Training GANs.”

using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: binarycrossentropy, logitcrossentropy
using Images
using MLDatasets
using Statistics, StatsBase
using Parameters: @with_kw
using Random
using Printf
using CUDAapi
using Zygote
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 100
    nclasses::Int = 10
    sample_each::Int = 5000
    epochs::Int = 20
    verbose_freq::Int = 1000
    output_x::Int = 10
    output_y::Int = 10
    lr_dscr::Float64 = 0.0002
    lr_gen::Float64 = 0.0002
end

function generator(args)
    return Chain(Dense(args.latent_dim, 6272), x->leakyrelu.(x, 0.2f0), x-> reshape(x, 7,7,128, size(x,2)),
            ConvTranspose((4, 4), 128 => 128; stride = 2, pad = 1),
            BatchNorm(128, leakyrelu),
            Dropout(0.25),
            ConvTranspose((4, 4), 128 => 64; stride = 2, pad = 1),
            BatchNorm(64, leakyrelu),
            Conv((7, 7), 64 => 1, tanh; stride = 1, pad = 3)) |> gpu
end

function discriminator(args)
    return Chain(Conv((3,3), 1=>128, pad=(1,1), stride = (2,2)),
                  x->leakyrelu.(x, 0.2f0),
                  Dropout(0.4),
                  Conv((3,3), 128=>128, pad= (1,1), stride = (2,2), leakyrelu),
                  x->leakyrelu.(x, 0.2f0),
                  x -> reshape(x, :, size(x, 4)),
                  Dropout(0.4), Dense(6272, 10)) |> gpu
end

function load_data(hparams)
    # Load MNIST dataset
    images, labels = MLDatasets.MNIST.traindata(Float32)
    # Normalize to [-1, 1] and convert it to WHCN
    image_tensor = permutedims(reshape(@.(2f0 * images - 1f0), 28, 28, 1, :), (2, 1, 3, 4))
    img, y = [],[]
    for i in 0:hparams.nclasses-1
        x = findall(==(i), labels)
        img_for_class = image_tensor[:,:,:,x]
        labels_for_class = labels[x]
        idxs = sample(1:length(x), hparams.sample_each, replace = false)
        [push!(img, img_for_class[:,:,:,idx]) for idx in idxs]
        push!(y, (labels_for_class[idxs]...))
    end
    imgs = zeros(28,28,1,length(y))
    for i in 1:length(img)
        imgs[:,:,:,i] = img[i]
    end
    n = length(y)
    y = float.(Flux.onehotbatch(y, 0:hparams.nclasses-1))
    idxs = sample(1:n, n, replace=false)
    img, y = imgs[:,:,:,idxs], y[:,idxs]
    data = [(img[:,:,:,r],y[:,r]) |> gpu for r in partition(1:n, hparams.batch_size)]
    return data
end

# Loss functions
function discr_loss(real_output, fake_output)
    real_loss = mean(binarycrossentropy.(real_output, 1f0))
    fake_loss = mean(binarycrossentropy.(fake_output, 0f0))
    return real_loss + fake_loss
end

function aux_loss(output, labels)
   return Flux.crossentropy(output, labels)
end

generator_loss(fake_output) = mean(binarycrossentropy.(fake_output, 1f0))

function custom_activation(output)
    exp_sum = sum(exp.(output))
    return exp_sum / (exp_sum + 1.0f0)
end

function train_discr(discr, fake_data, original_data ,labels, opt_discr)
    ps = params(discr)
    loss = 0.0
    for i in 1:2
    loss, back = Zygote.pullback(ps) do
           original_common = discr(original_data) |> gpu
           fake_common =  discr(fake_data) |> gpu
           d1 = discr_loss(custom_activation(original_common), custom_activation(fake_common))
           x = softmax(original_common)
           a1 = aux_loss(x, labels)
           (d1 + a1)/2
           end
    grads = back(1f0)
    update!(opt_discr, ps, grads)
    end
    return loss
end

Zygote.@nograd train_discr

function train_gan(gen, discr, original_data, label, opt_gen, opt_discr, hparams)
    noise = randn!(similar(original_data, (hparams.latent_dim, hparams.batch_size))) |> gpu
    ps = params(gen)
    loss = Dict()
    loss["gen"], back = Zygote.pullback(ps) do
            fake = gen(noise) |> gpu
            loss["discr"] = train_discr(discr, fake, original_data, label, opt_discr)
            fake_common = discr(fake)
            generator_loss(custom_activation(fake_common))
    end
    grads = back(1f0)
    update!(opt_gen, ps, grads)
    return loss
end

function create_output_image(gen, fixed_noise, hparams)
    @eval Flux.istraining() = false
    fake_images = @. cpu(gen(fixed_noise))
    @eval Flux.istraining() = true
    image_array = dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function train()
    hparams = HyperParams()

    data = load_data(hparams)

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = discriminator(hparams)

    # Generator
    gen =  generator(hparams)

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr, (0.5,0.999))
    opt_gen = ADAM(hparams.lr_gen, (0.5,0.999))

    isdir("output")||mkdir("output")

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for (x,y) in data
            # Update discriminator and generator
            loss = train_gan(gen, dscr, x, y, opt_gen, opt_dscr, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss["discr"]), Generator loss = $(loss["gen"])")
                # Save generated fake image
                output_image = create_output_image(gen, fixed_noise, hparams)
                save(@sprintf("output/sgan_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(gen, fixed_noise, fixed_labels, hparams)
    save(@sprintf("output/sgan_steps_%06d.png", train_steps), output_image)
end

cd(@__DIR__)
train()


text:
#==============================================================================
    stuff to look up from the REPL
    usage uncertain so far

    TODO:   this needs to parse keyword tokens and search the array elements

    author:     mahiki@
 =============================================================================#

"""
    @remindme

Return helpful text based on keywords.

@remindme(keyword)

# Examples
```jldoctest
julia> @remindme perf;
avoid non-constant globals
focus on inplace operations
write for-loops
minimize type-instabilities
```
"""
macro remindme(keyword)
    x = string(keyword)

    for patt in keys(remindmelist)
        if occursin(x, patt)
            for item in remindmelist[patt]
                println(item)
            end
        end
    end
end

# need to match key and call those keys

remindmelist = Dict(
    "performance tips" =>
        [
            "avoid non-constant globals"
            "focus on inplace operations"
            "write for-loops"
            "minimize type-instabilities"
            "write code in functions"
        ]
    )

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
function RK4(::Type{T},
            N::Int,
            xyz::Array{Float64,1},
            σ::Float64,
            ρ::Float64,
            β::Float64,
            s::Float64,
            Δt::Float64) where {T<:AbstractFloat}

    # preallocate for storing results - store without scaling
    XYZout = Array{Float64,2}(undef,3,N+1)
    XYZout[:,1] = xyz

    # Runge Kutta 4th order coefficients including time step and sigma for x
    RKα = zeros(3,4)
    RKα[2,:] = [1/6.,1/3.,1/3.,1/6.]*Δt
    RKα[3,:] = RKα[2,:]
    RKα[1,:] = RKα[2,:]*σ

    RKβ = zeros(3,3)
    RKβ[2,:] = [1/2.,1/2.,1.]*Δt
    RKβ[3,:] = RKβ[2,:]
    RKβ[1,:] = RKβ[2,:]*σ

    # convert everything to the desired number system determined by T and scale
    xyz = T.(s*xyz)
    s_inv = T(1.0/s)
    ρ,β = T.([ρ,β])
    RKα = T.(RKα)
    RKβ = T.(RKβ)

    # preallocate memory for intermediate results
    xyz0 = deepcopy(xyz)
    xyz1 = deepcopy(xyz)
    dxyz = zero(xyz)       # tendencies

    for i = 1:N
        @simd for j in 1:3
            @inbounds xyz1[j] = xyz[j]
        end

        for rki = 1:4
            rhs!(dxyz,xyz1[1],xyz1[2],xyz1[3],ρ,β,s_inv)

            if rki < 4
                @simd for j in 1:3
                    @inbounds xyz1[j] = xyz[j] + dxyz[j] * RKβ[j,rki]
                end
            end

            # sum the RK steps on the go
            @simd for j in 1:3
                @inbounds xyz0[j] += dxyz[j] * RKα[j,rki]
            end
        end

        @simd for j in 1:3
            @inbounds xyz[j] = xyz0[j]
        end

        # store as 64bit, undo scaling
        XYZout[:,i+1] = Float64.(xyz)/s
    end

    return XYZout
end


text:
using ObjectDetector
ObjectDetector.YOLO.v2_tiny_416_COCO(silent=true)
ObjectDetector.YOLO.v3_tiny_416_COCO(silent=true)
ObjectDetector.YOLO.v3_320_COCO(silent=true)
ObjectDetector.YOLO.v3_416_COCO(silent=true)
ObjectDetector.YOLO.v3_608_COCO(silent=true)