text:
function server(port, log_dir; tasks_file="tasks", finished_tasks_file="finished_tasks")
    # Load tasks
    tasks = open(tasks_file) do file
        map(readlines(file)) do line
            strip(line)
        end
    end

    # Load finished tasks and filter out those from the todo list
    finished_tasks = isfile(finished_tasks_file) ? readdlm(finished_tasks_file) : []
    todo = filter(1:length(tasks)) do i
        !(i in finished_tasks)
    end

    # Order is randomized to balance load in case of unequal task
    # execution time
    todo = todo[randperm(length(todo))]

    jobid = ENV["SLURM_JOB_ID"]

    open("$(log_dir)/server-$(jobid).log", "w+") do log_file
        server_logger = SimpleLogger(log_file, Logging.Debug)

        @debug "Starting server at $(gethostname()), $(getipaddr())"
        @debug "$(length(tasks)) tasks total, $(length(todo)) unfinished"

        @async begin
            server = listen(IPv4(0), port)
            clients_running = 0
            todo_id = 0
            while true
                try
                    sock = accept(server)
                    @debug "got connection"
                    cl_query = strip(readline(sock))
                    @debug "got query: $(cl_query)"
                    if cl_query == "-1"
                        clients_running += 1
                    else
                        run_id = parse(Int, cl_query)
                        @debug "client claims $(run_id) is finished"
                        push!(finished_tasks, run_id)
                        writedlm(finished_tasks_file, finished_tasks)
                        sleep(0.1)
                    end
                    todo_id += 1
                    if todo_id > length(todo)
                        write(sock, "-1\n")
                        clients_running -= 1
                        clients_running <= 0 && break
                    else
                        run_id = todo[todo_id]
                        write(sock, "$(run_id)\n")
                        cmd = tasks[run_id]
                        @debug "telling client to run $(cmd)"
                        write(sock, cmd)
                    end
                    close(sock)
                catch e
                    error(server_logger, "$(e)")
                    break
                end
            end
            close(server)
            @debug "Server finished"
        end
    end
end


text:
__precompile__()
module GadgetReader

using LinearAlgebra: rotate!

using Cosmology
using GadgetIO
using LinearAlgebra
using Statistics
using Unitful
using UnitfulAstro

include("gadget_config.jl")
include("snapshot.jl")
include("particles.jl")
include("read_particles.jl")
include("transformations.jl")
include("units.jl")
include("utils.jl")

# structs
export Snapshot,
       GadgetConfig,
       Particles,

       # read
       read_redshift,
       read_header_particle_mass,

       # transformations
       rotate,
       rotate!,
       rotate_edgeon,
       rotate_edgeon!,
       translate,
       translate!,
       center_of_mass_iterative,
       translate_to_center_of_mass_iterative,
       translate_to_center_of_mass_iterative!,
       center_of_velocity,
       translate_to_center_of_velocity,
       translate_to_center_of_velocity!

       # units

end

ext: jl


text:
module Orbits

using ConcreteStructs
using StaticArrays
using Printf

export SimpleOrbit, KeplerianOrbit

abstract type AbstractOrbit end

Base.broadcastable(o::AbstractOrbit) = Ref(o)

"""
    relative_position(::AbstractOrbit, t)

The relative position, `[x, y, z]`, of the companion compared to the host at time `t`. In
other words, this is the vector pointing from the host to the companion along the line of
sight. Nominally, the units of this distance are relative to the host's radius. For
example, a distance of 2 is 2 *stellar* radii.
"""
relative_position(::AbstractOrbit, t)

include("simple.jl")
include("keplerian.jl")

end # module


text:
" Implements a NonlinearSource operator "
type NonlinearForce <: Kernel
    u::Variable
end

import MOOSE.computeQpResidual

@inline function computeQpResidual(kernel::NonlinearForce, qp::Int64, i::Int64)
    u = kernel.u

    return -u.value[qp] * u.value[qp] * u.phi[qp][i]
end

import MOOSE.computeQpJacobian

@inline function computeQpJacobian(kernel::NonlinearForce, v::Variable, qp::Int64, i::Int64, j::Int64)::Float64
    u = kernel.u

    if u.id == v.id
        return -2 * u.value[qp] * v.phi[qp][j] * u.phi[qp][i]
    end

    return 0
end

import MOOSE.coupledVars

@inline function coupledVars(kernel::NonlinearForce)::Array{Variable}
    return [kernel.u]
end


text:
# estimate with MCMC (Turing) for conjugate prior
# compare with jlbayes_conjugate_regression.jl

using Turing
using Distributions
using Random
using LinearAlgebra
using StatsPlots
using Optim

# Generate moc data
n = 50
dim = 1
Random.seed!(99)
u = rand(Normal(0, 0.7), n) # true value σ² = 0.49
x = rand(Uniform(-sqrt(3.0), sqrt(3.0)), n, dim) # x generated from Uniform[-√3, √3]
X = [ones(n) x] # dependent variable matrix with ones
β = Vector(1:1:dim+1) # true coefficients
y = X*β + u # true data generating process

# set prior (conjugate prior for camparison with analytic result)
# p.78 eq (3.32)
k = size(X)[2]
β0 = zeros(k)
τ0 = 0.2
A0 = τ0 * I # Identity matrix
ν0 = 5.0
λ0 = 7.0
# H0 = (λ0 / ν0) * inv(A0) # scale matrix for marginal prior of multivariate coefficient (mv Tdist)
# h0 = diag(sqrt(Matrix(H0, k, k))) # scale params for marginal prior of single coefficient (T dist) H0の対角成分
sd0 = diag(sqrt(Matrix(A0, k, k))) # the parameter for conditional prior for β|σ² ~ N_k(β0, σ²*inv(A0)), (inv(A0) part)

@model function linear_regression(X, y)
    # conjugate priors
    # https://turing.ml/dev/tutorials/05-linear-regression/

    # Set variance priors
    σ² ~ InverseGamma(ν0/2, λ0/2)

    # Set the priors on our coefficients.
    nfeatures = size(X, 2)
    β ~ MvNormal(β0, σ² * inv(A0))

    # Write likelihood.
    y ~ MvNormal(X * β, σ² * I)
end

# Draw random number from posterior.
model = linear_regression(X, y)
n_draws = 5_000
n_chains = 4
n_tune = 1_000
# sample(model, sampler, parallel_type, n, n_chains)
# chain = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains, discard_initial=n_tune)
# How to implement burn-in / warmup? --> discard initial value automatically. ; discard_adapt=false to turn off
chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains)
chn = sample(model, NUTS(), MCMCThreads(), n_draws, n_chains; discard_adapt=false)

# Chain info.
typeof(chn)
size(chn)
describe(chn)
summarize(chn)
quantile(chn)
plot(chn)

# Indexing a Chains object.
β1_chn = chn["β[1]"]
β2_chn = chn[Symbol("β[2]")]
σ²_chn = chn[:σ²]
all_β_chn = group(chn, :β) # get subset of all parameters include "β"

typeof(β1_chn)
β1_chn[iter = 1000:5000, chain = 1]
mean(β1_chn)
std(β1_chn)

# Monte Carlo Standard Error
mcse(θ_chn) = sqrt((1 / (length(θ_chn) * (length(θ_chn) - 1))) * sum((θ_chn .- mean(θ_chn)).^2))

[β1_chn, β2_chn, σ²_chn]

mcse.([β1_chn, β2_chn, σ²_chn])

# without burn-in
mean(β1_chn[iter = 1:5000])

histogram(β1_chn[:, 1])
histogram!(β1_chn[:, 2])
histogram!(β1_chn[:, 3])
histogram!(β1_chn[:, 4])

# mode (MLE and MAP) estimates (Optimization with Optim.jl) as an option (continuous only)
# Note that loading Optim explicitly is required for mode estimation to function,
# as Turing does not load the opimization suite unless Optim is loaded as well.

# Generate a MLE estimate (default optimizer is LBFGS optimizer).
mle_estimate = optimize(model, MLE())
# Use Newton
# mle_estimate = optimize(model, MLE(), Newton())

# Generate a MAP estimate (default optimizer is LBFGS optimizer).
map_estimate = optimize(model, MAP())

# Analyze mode estimate
# Import StatsBase to use it's statistical methods.
# using StatsBase

# Print out the coefficient table.
# coeftable(mle_estimate)

# Analytical Solution
inv(X' * X) * X' * y


text:
# Implement the interface for breadth first in a dummy way
struct DummyBreadthFirstSearch <: AbstractBreadthFirstSearch{Symbol} end

@testset "Forced implementations" begin
    search = DummyBreadthFirstSearch()
    tree = BPTree("Dummy data")

    @test_throws MissingImplementationError BranchAndPrune.bisect(search, "Dummy element")
    @test_throws MissingImplementationError BranchAndPrune.process(search, tree)

    # TODO add test for get_leaf_id! and insert_leaf! for general searches
    # TODO add test for keyfunc for AbstractKeySearch
end

BranchAndPrune.process(::DummyBreadthFirstSearch, s::Symbol) = s, s
BranchAndPrune.bisect(::DummyBreadthFirstSearch, s::Symbol) = :store, :store

# TODO Add similar tests for depth first and key search
@testset "Breadth first search" begin
    #= Build the following tree for testing
                    (1)
                   /   \
                (2)     (5: to_bisect)
               /   \
    (3: to_store)  (4: to_discard)
    =#
    rt = BPNode(0, [2, 5])
    node = BPNode(1, [3, 4])
    to_store = BPLeaf(:store, 2, :working)
    to_discard = BPLeaf(:discard, 2, :working)
    to_bisect = BPLeaf(:bisect, 1, :working)

    tree = BPTree(Dict(1 => rt, 2 => node),
                  Dict(3 => to_store, 4 => to_discard, 5 => to_bisect),
                  [3, 4, 5])

    search = DummyBreadthFirstSearch()

    # First iteration processes the to_store leaf
    tree, _ = iterate(search, tree)

    @test nnodes(tree) == 5
    @test tree[3].status == :final
    @test length(tree.working_leaves) == 2

    # Second iteration processes the to_discard leaf
    tree, _ = iterate(search, tree)

    @test nnodes(tree) == 4
    @test length(tree.working_leaves) == 1

    # Second iteration processes the to_bisect leaf
    tree, _ = iterate(search, tree)

    @test nnodes(tree) == 6
    @test length(tree.working_leaves) == 2
    @test isa(tree[5], BPNode)
end


text:
include("./RPSAgents.jl")
import Pkg
Pkg.add("StatsBase")
Pkg.add("Distributions")
using StatsBase, Distributions
using .RPSAgents

function evaluate(agent1, agent2, steps)
  # Hàm đánh giá 2 agent đầu vào bằng cách cho chúng đấu với nhau
  score = 0 # Số lượt thắng của agent1
  tie = 0 # Số lượt hòa của agent1
  agent1_last_step, agent2_last_step = 0, 0
  for i in 0:steps-1
    agent1_observation = Dict(
      "step" => i,
      "lastOpponentAction" => agent2_last_step
    )
    agent2_observation = Dict(
      "step" => i,
      "lastOpponentAction" => agent1_last_step
    )
    # 2 agent lần lượt ra quyết định
    agent1_action = agent1.run(agent1_observation)
    agent2_action = agent2.run(agent2_observation)
    if(agent1_action == agent2_action)
      tie += 1
    else
      counter_agent1_action = (agent1_action + 1) % 3
      if(agent2_action == counter_agent1_action)
        score += 0
      else
        score += 1
      end
    end
    agent1_last_step = agent1_action
    agent2_last_step = agent2_action
  end
  return score, steps - score - tie
end

function compare()
  # Thực nghiệm so sánh ThompsonSamplingAgent với các agent khác
  # Danh sách các agent để so sánh với agent của chúng ta
  agents = Dict(
    "random_agent" => RandomAgent(),
    "only_scissor_agent" => OnlyScissorAgent(),
    "only_rock_agent" => OnlyRockAgent(),
    "only_paper_agent" => OnlyPaperAgent(),
    "copy_agent" => CopyAgent(),
    "freq_counting_agent" => FreqCountingAgent(),
    "counter_last_action_agent" => CounterLastActionAgent(),
    "markov_agent" => MarkovAgent(),
    "thompson_sampling_agent" = ThompsonSamplingAgent()
  )
  for (name, agent) in agents
    ourScores = []
    oppoScores = []
    usedAgents = Dict(
    "random_agent" => [],
    "only_scissor_agent" => [],
    "only_rock_agent" => [],
    "only_paper_agent" => [],
    "copy_agent" => [],
    "freq_counting_agent" => [],
    "counter_last_action_agent" => []
    )
    for i in (1:5)
      our_agent = ThompsonSamplingAgent()
      score1, score2 = evaluate(our_agent, agent, 1000)
      append!(ourScores, score1)
      append!(oppoScores, score2)
      for (agentName, agentCount) in our_agent.agent_count
        append!(usedAgents[agentName], agentCount)
      end
    end
    println("\n ----------------- \n")
    println(name)
    println("Our score: ", ourScores)
    println("Oppo score: ", oppoScores)
    println("Used agents: ")
    for (agentName, agentCounts) in usedAgents
      println("Mean used ", agentName, ": ", agentCounts)
    end
  end
end
compare()


text:
using LinearAlgebra
using Random
using Distributions

# plot features require Gadfly, etc
using ApproxManifoldProducts

# more packages required
using TransformUtils

const LinAlg = LinearAlgebra
const KDE = KernelDensityEstimate
const TU = TransformUtils
const AMP = ApproxManifoldProducts

# temporarily required for experimental versions of on-manifold products
KDE.setForceEvalDirect!(true)

## create two densities

# two densities on a cylinder
p = manikde!(randn(2,100), (:Euclid, :Circular) )

pts2a = 3.0*randn(1,100).+5.0
pts2b = TransformUtils.wrapRad.(0.5*randn(1,100).+pi)
q = manikde!([pts2a;pts2b], (:Euclid, :Circular) )

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q], (:Euclid, :Circular), Niter=2)

## better wrap around example

pts1a = 6.0*rand(1,100).-3.0
pts1b = TransformUtils.wrapRad.(0.3*randn(1,100).-0.6*pi)
p = manikde!([pts1a;pts1b], (:Euclid, :Circular) )

pts2a = 6.0*rand(1,100).-3.0
pts2b = TransformUtils.wrapRad.(0.3*randn(1,100).+0.6*pi)
q = manikde!([pts2a;pts2b], (:Euclid, :Circular) )

# approximate the product between hybrid manifold densities
pq = manifoldProduct([p;q], (:Euclid, :Circular), Niter=3)

## import plotting

using Cairo
using Gadfly, Colors
using KernelDensityEstimatePlotting

## Different plots

plotKDE([p;q; pq], levels=3, c=["red";"green";"magenta"])

plotKDE([p;q; pq], dims=[1],c=["red";"green";"magenta"])

pl = plotKDECircular( [marginal(p, [2;]);marginal(q, [2;]); marginal(pq, [2;])] )

pl = Gadfly.plot(y=pGM[1,:],x=pGM[2,:], Geom.histogram2d(xbincount=50,ybincount=30))

## export the plot to image file

pl |> PDF("/tmp/test.pdf",20cm,15cm)
@async run(`evince /tmp/test.pdf`)


text:
# syntax: proto3
const Syntax = (;[
    Symbol("SYNTAX_PROTO2") => Int32(0),
    Symbol("SYNTAX_PROTO3") => Int32(1),
]...)

mutable struct Option <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Option(; kwargs...)
        obj = new(meta(Option), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Option
const __meta_Option = Ref{ProtoMeta}()
function meta(::Type{Option})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Option)
            __meta_Option[] = target = ProtoMeta(Option)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :value => _Any]
            meta(target, Option, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Option[]
    end
end
function Base.getproperty(obj::Option, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::_Any
    else
        getfield(obj, name)
    end
end

mutable struct EnumValue <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function EnumValue(; kwargs...)
        obj = new(meta(EnumValue), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct EnumValue
const __meta_EnumValue = Ref{ProtoMeta}()
function meta(::Type{EnumValue})
    ProtoBuf.metalock() do
        if !isassigned(__meta_EnumValue)
            __meta_EnumValue[] = target = ProtoMeta(EnumValue)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :number => Int32, :options => Base.Vector{Option}]
            meta(target, EnumValue, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_EnumValue[]
    end
end
function Base.getproperty(obj::EnumValue, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :number
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    else
        getfield(obj, name)
    end
end

const Field_Kind = (;[
    Symbol("TYPE_UNKNOWN") => Int32(0),
    Symbol("TYPE_DOUBLE") => Int32(1),
    Symbol("TYPE_FLOAT") => Int32(2),
    Symbol("TYPE_INT64") => Int32(3),
    Symbol("TYPE_UINT64") => Int32(4),
    Symbol("TYPE_INT32") => Int32(5),
    Symbol("TYPE_FIXED64") => Int32(6),
    Symbol("TYPE_FIXED32") => Int32(7),
    Symbol("TYPE_BOOL") => Int32(8),
    Symbol("TYPE_STRING") => Int32(9),
    Symbol("TYPE_GROUP") => Int32(10),
    Symbol("TYPE_MESSAGE") => Int32(11),
    Symbol("TYPE_BYTES") => Int32(12),
    Symbol("TYPE_UINT32") => Int32(13),
    Symbol("TYPE_ENUM") => Int32(14),
    Symbol("TYPE_SFIXED32") => Int32(15),
    Symbol("TYPE_SFIXED64") => Int32(16),
    Symbol("TYPE_SINT32") => Int32(17),
    Symbol("TYPE_SINT64") => Int32(18),
]...)

const Field_Cardinality = (;[
    Symbol("CARDINALITY_UNKNOWN") => Int32(0),
    Symbol("CARDINALITY_OPTIONAL") => Int32(1),
    Symbol("CARDINALITY_REQUIRED") => Int32(2),
    Symbol("CARDINALITY_REPEATED") => Int32(3),
]...)

mutable struct Field <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Field(; kwargs...)
        obj = new(meta(Field), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Field
const __meta_Field = Ref{ProtoMeta}()
function meta(::Type{Field})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Field)
            __meta_Field[] = target = ProtoMeta(Field)
            fnum = Int[1,2,3,4,6,7,8,9,10,11]
            allflds = Pair{Symbol,Union{Type,String}}[:kind => Int32, :cardinality => Int32, :number => Int32, :name => AbstractString, :type_url => AbstractString, :oneof_index => Int32, :packed => Bool, :options => Base.Vector{Option}, :json_name => AbstractString, :default_value => AbstractString]
            meta(target, Field, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Field[]
    end
end
function Base.getproperty(obj::Field, name::Symbol)
    if name === :kind
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :cardinality
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :number
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :type_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :oneof_index
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :packed
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    elseif name === :json_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :default_value
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct _Enum <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function _Enum(; kwargs...)
        obj = new(meta(_Enum), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct _Enum
const __meta__Enum = Ref{ProtoMeta}()
function meta(::Type{_Enum})
    ProtoBuf.metalock() do
        if !isassigned(__meta__Enum)
            __meta__Enum[] = target = ProtoMeta(_Enum)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :enumvalue => Base.Vector{EnumValue}, :options => Base.Vector{Option}, :source_context => SourceContext, :syntax => Int32]
            meta(target, _Enum, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta__Enum[]
    end
end
function Base.getproperty(obj::_Enum, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :enumvalue
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{EnumValue}
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    elseif name === :source_context
        return (obj.__protobuf_jl_internal_values[name])::SourceContext
    elseif name === :syntax
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct _Type <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function _Type(; kwargs...)
        obj = new(meta(_Type), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct _Type
const __meta__Type = Ref{ProtoMeta}()
function meta(::Type{_Type})
    ProtoBuf.metalock() do
        if !isassigned(__meta__Type)
            __meta__Type[] = target = ProtoMeta(_Type)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :fields => Base.Vector{Field}, :oneofs => Base.Vector{AbstractString}, :options => Base.Vector{Option}, :source_context => SourceContext, :syntax => Int32]
            meta(target, _Type, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta__Type[]
    end
end
function Base.getproperty(obj::_Type, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :fields
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Field}
    elseif name === :oneofs
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :options
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Option}
    elseif name === :source_context
        return (obj.__protobuf_jl_internal_values[name])::SourceContext
    elseif name === :syntax
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export Syntax, _Type, Field_Kind, Field_Cardinality, Field, _Enum, EnumValue, Option


text:
module julia_bot

using Base.Threads
using Printf
using Random
using SQLite
using Serialization
using Sockets
using Tables
using UUIDs

include("../../jjson/jjson.jl")
include("bot_state.jl")
include("parse_config.jl")
include("tcp.jl")
include("mcmc.jl")
include("irc.jl")
include("repl.jl")

export run_bot

function run_bot(config_filename::String)
    tcp_sock = TCPSocket(;delay=true)
    bot = Bot(tcp_sock)

    config::Config = parseconfig!(config_filename, bot)
    oauth::String = read_oauth_file(config.token_file)

    status = irc_connect(bot.sock, config)
    println(stdout, "Status after irc_connect: ", status)

    status = irc_auth(bot.sock, config.nickname, oauth)
    println(stdout, "Status after irc_auth: ", status)

    if !isempty(config.channel)
        irc_join(bot.sock, config.channel)
    end # if
    Threads.@spawn local_repl(bot.sock)

    println(stdout, "Initialise:")
    status, msg = irc_readlines(bot.sock, config)
end

end # module

max_stars_repo_name: Iaggelis/bald-bot

max_issues_repo_name: Iaggelis/bald-bot

max_forks_repo_name: Iaggelis/bald-bot


text:
"""
    twoD(X)

Project `X` to two dmensions using PCA
"""
function twoD(X)
    X = X .- mean(X,dims=1)
    X = X ./ std(X,dims=1)
    s = svd(X)
    Y = s.U[:,1:2]#.*s.S[1:2]'
    Y[:,1], Y[:,2]
end

"""
    threeD(X)

Project `X` to three dimensions using PCA
"""
function threeD(X)
    X = X .- mean(X,dims=1)
    X = X ./ std(X,dims=1)
    s = svd(X)
    Y = s.U[:,1:3]#.*s.S[1:2]'
    Y[:,1], Y[:,2], Y[:,3]
end

"""
    s1(x, dims=:)

normalize x sums to 1
"""
function s1(x, dims=:)
    if isderiving()
        m = minimum(real(x), dims=dims)
        if any(<(0), m)
            x = x .- m
        end
        return x ./ sum(x, dims=dims)
    else
        m = minimum(x, dims=dims)
        x = float.(x)
        if any(<(0), m)
            x .-= m
        end
        x ./= sum(x, dims=dims)
    end
end

"""
    n1(x)

normalize x norm 1
"""
n1(x) = x./norm(x)
n1(x::AbstractMatrix, dims) = mapslices(n1, x, dims=dims)
"""
    v1(x, dims=:)

normalize x var 1
"""
function v1(x, dims=:)
    x = x .- mean(x, dims=dims)
    x .= x./std(x, dims=dims)
end

function v1!(x, dims=:)
    x .-= mean(x, dims=dims)
    x .= x./std(x, dims=dims)
end

"""
    median1(x, dims=:)

normalize x median 0, median absolute deviation = 1
"""
function m1(x, dims=:)
    x = x .- median(x, dims=dims)
    # x .= x./mapslices(mad, x, dims=dims)
end
function m1(x::AbstractVector)
    x = x .- median(x)
    x .= x./mad(x, normalize=true)
end

"""
    bp_filter(x, passband)

Band-pass filter, passband is tuple, `fs` assumed = 1
"""
function bp_filter(x, passband)
    responsetype = Bandpass(passband..., fs=1)
    designmethod = Butterworth(2)
    filt(digitalfilter(responsetype, designmethod), x)
end

"""
    lp_filter(x, cutoff)

Low-pass filter, `fs` assumed = 1
"""
function lp_filter(x, cutoff)
    responsetype = Lowpass(cutoff, fs=1)
    designmethod = Butterworth(2)
    filt(digitalfilter(responsetype, designmethod), x)
end


text:
⊗(A,B) = kron(A, B)

function transforms_e(Lw, r̂, l)

    A = (Lw - Lw*r̂ - Lw)/(2*tanh((r̂-1)/l) + tanh(-2/l)*(r̂ - 1))
    b = (A*tanh(-2/l) + Lw)/2
    c = Lw - b
    xt = (r,s) -> (A .* tanh.((r .- 1) ./ l) .+ b .* r .+ c,
                   ((A .* sech.((r .- 1) ./ l).^2) ./ l) .+ b,
                   zeros(size(s)))
    yt = (r,s) -> (A .* tanh.((s .- 1) ./ l) .+ b.*s .+ c,
                   zeros(size(r)),
                   ((A .* sech.((s .- 1) ./ l).^2) ./ l) .+ b)

    #=
    a = (Lw*r̂ + Lw - r̂ - 1)/(2 * (r̂*tanh(1/l) - tanh(r̂/l)))
    b = (Lw - 2*a*tanh(1/l))/2
    c = b - a*tanh(-1/l)
    xt = (r,s) -> (a .* tanh.((r) ./ l) .+ b .* r .+ c,
                   ((a .* sech.((r) ./ l).^2) ./ l) .+ b,
                   zeros(size(s)))
    yt = (r,s) -> (a .* tanh.((s) ./ l) .+ b.* s .+ c,
                   zeros(size(r)),
                   ((a .* sech.((s) ./ l).^2) ./ l) .+ b)
    =#

    #=
    r = -1:.01:1
    s = -1:.01:1
    scatter(r, xt(r,s), legend=false)
    gui()
    =#

    return xt, yt

end

function transforms_ne(Lw, el_x, el_y)

    xt = (r,s) -> (el_x .* tan.(atan((Lw)/el_x).* (0.5*r .+ 0.5)),
                   el_x .* sec.(atan((Lw)/el_x).* (0.5*r .+ 0.5)).^2 * atan((Lw)/el_x) * 0.5 ,
                   zeros(size(s)))

    yt = (r,s) -> (el_y .* tan.(atan((Lw)/el_y).* (0.5*s .+ 0.5)) ,
                   zeros(size(r)),
                   el_y .* sec.(atan((Lw)/el_y).*(0.5*s .+ 0.5)) .^2 * atan((Lw)/el_y) * 0.5 )

    return xt, yt
end

function transforms_n(Lw)

    xt = (r,s) -> (Lw/2 .* (r .+ 1),
                   fill(Lw/2, size(r)),
                   zeros(size(s)))

    yt = (r,s) -> (Lw/2 .* (s .+ 1),
                   zeros(size(r)),
                   fill(Lw/2, size(s)))
    return xt, yt
end

function create_metrics(Nr, Ns, B_p, μ, ρ,
                        xf=(r,s)->(r, ones(size(r)), zeros(size(r))),
                        yf=(r,s)->(s, zeros(size(s)), ones(size(s))))

    Nrp = Nr + 1
    Nsp = Ns + 1
    Np = Nrp * Nsp

    r = range(-1, stop=1, length=Nrp)
    s = range(-1, stop=1, length=Nsp)

    # Create the mesh
    r = ones(1, Nsp) ⊗ r
    s = s' ⊗ ones(Nrp)

    (x, xr, xs) = xf(r, s)
    (y, yr, ys) = yf(r, s)

    J = xr .* ys - xs .* yr

    @assert minimum(J) > 0

    JI = 1 ./ J

    rx =  ys ./ J
    sx = -yr ./ J
    ry = -xs ./ J
    sy =  xr ./ J

    μx = μy = μ(x, y, B_p)

    # variable coefficient matrix components
    crr = J .* (rx .* μx .* rx + ry .* μy .* ry)
    crs = J .* (sx .* μx .* rx + sy .* μy .* ry)
    css = J .* (sx .* μx .* sx + sy .* μy .* sy)

    # surface matrices
    (xf1, yf1) = (view(x, 1, :), view(y, 1, :))
    nx1 = -ys[1, :]
    ny1 =  xs[1, :]
    sJ1 = hypot.(nx1, ny1)
    nx1 = nx1 ./ sJ1
    ny1 = ny1 ./ sJ1
    μf1 = μ(xf1, yf1, B_p)
    η = μf1 ./ (2 * sqrt.(μf1 ./ ρ(xf1, yf1, B_p)))

    (xf2, yf2) = (view(x, Nrp, :), view(y, Nrp, :))
    nx2 =  ys[end, :]
    ny2 = -xs[end, :]
    sJ2 = hypot.(nx2, ny2)
    nx2 = nx2 ./ sJ2
    ny2 = ny2 ./ sJ2
    μf2 = μ(xf2, yf2, B_p)

    (xf3, yf3) = (view(x, :, 1), view(y, :, 1))
    nx3 =  yr[:, 1]
    ny3 = -xr[:, 1]
    sJ3 = hypot.(nx3, ny3)
    nx3 = nx3 ./ sJ3
    ny3 = ny3 ./ sJ3

    (xf4, yf4) = (view(x, :, Nsp), view(y, :, Nsp))
    nx4 = -yr[:, end]
    ny4 =  xr[:, end]
    sJ4 = hypot.(nx4, ny4)
    nx4 = nx4 ./ sJ4
    ny4 = ny4 ./ sJ4

    (coord = (x,y),
     facecoord = ((xf1, xf2, xf3, xf4), (yf1, yf2, yf3, yf4)),
     μf2 = μf2,
     η = η,
     crr = crr, css = css, crs = crs,
     J=J,
     JI = JI,
     sJ = (sJ1, sJ2, sJ3, sJ4),
     nx = (nx1, nx2, nx3, nx4),
     ny = (ny1, ny2, ny3, ny4),
     rx = rx, ry = ry, sx = sx, sy = sy)
end

function transfinite_blend(α1, α2, α3, α4, α1s, α2s, α3r, α4r, r, s)
  # +---4---+
  # |       |
  # 1       2
  # |       |
  # +---3---+
  @assert [α1(-1) α2(-1) α1( 1) α2( 1)] ≈ [α3(-1) α3( 1) α4(-1) α4( 1)]

  x = (1 .+ r) .* α2(s)/2 + (1 .- r) .* α1(s)/2 +
      (1 .+ s) .* α4(r)/2 + (1 .- s) .* α3(r)/2 -
     ((1 .+ r) .* (1 .+ s) .* α2( 1) +
      (1 .- r) .* (1 .+ s) .* α1( 1) +
      (1 .+ r) .* (1 .- s) .* α2(-1) +
      (1 .- r) .* (1 .- s) .* α1(-1)) / 4

  xr =  α2(s)/2 - α1(s)/2 +
        (1 .+ s) .* α4r(r)/2 + (1 .- s) .* α3r(r)/2 -
      (+(1 .+ s) .* α2( 1) +
       -(1 .+ s) .* α1( 1) +
       +(1 .- s) .* α2(-1) +
       -(1 .- s) .* α1(-1)) / 4

  xs = (1 .+ r) .* α2s(s)/2 + (1 .- r) .* α1s(s)/2 +
       α4(r)/2 - α3(r)/2 -
      (+(1 .+ r) .* α2( 1) +
       +(1 .- r) .* α1( 1) +
       -(1 .+ r) .* α2(-1) +
       -(1 .- r) .* α1(-1)) / 4

  return (x, xr, xs)
end

function transfinite(x1, x2, x3, x4,
                     y1, y2, y3, y4)

        ex = [(α) -> x1 * (1 .- α) / 2 + x3 * (1 .+ α) / 2,
          (α) -> x2 * (1 .- α) / 2 + x4 * (1 .+ α) / 2,
          (α) -> x1 * (1 .- α) / 2 + x2 * (1 .+ α) / 2,
          (α) -> x3 * (1 .- α) / 2 + x4 * (1 .+ α) / 2]
    exα = [(α) -> -x1 / 2 + x3 / 2,
           (α) -> -x2 / 2 + x4 / 2,
           (α) -> -x1 / 2 + x2 / 2,
           (α) -> -x3 / 2 + x4 / 2]
    ey = [(α) -> y1 * (1 .- α) / 2 + y3 * (1 .+ α) / 2,
          (α) -> y2 * (1 .- α) / 2 + y4 * (1 .+ α) / 2,
          (α) -> y1 * (1 .- α) / 2 + y2 * (1 .+ α) / 2,
          (α) -> y3 * (1 .- α) / 2 + y4 * (1 .+ α) / 2]
    eyα = [(α) -> -y1 / 2 + y3 / 2,
           (α) -> -y2 / 2 + y4 / 2,
           (α) -> -y1 / 2 + y2 / 2,
           (α) -> -y3 / 2 + y4 / 2]

    xt(x,y) = transfinite_blend(ex[1], ex[2], ex[3], ex[4],
                                exα[1], exα[2], exα[3], exα[4],
                                x, y)
    yt(x,y) = transfinite_blend(ey[1], ey[2], ey[3], ey[4],
                                eyα[1], eyα[2], eyα[3], eyα[4],
                                x, y)

    return xt, yt

end


text:
#=

Data structures for exponents and polynomials. Orderings. Arithmetic operations. All experimental.

For the moment, lex order is assumed and almost enforced.

=#

import Base.+, Base.*, Base.-, Base.^, Base.repr

immutable Indeterminate
    degree::Int
    Indeterminate(degree::Int) = degree > 0 ? new(degree) :
        error("Degree must be > 0")
end

type Monomial
    indeterminate::Indeterminate
    exponent::Array{Int,1}
    coefficient::Number

    function Monomial(x::Indeterminate)
        new(x, fill(1,x.degree), 1)
    end

    function Monomial(x::Indeterminate, c::Number)
        new(x, fill(1, x.degree), c)
    end

    function Monomial(x::Indeterminate, ex::Array{Int,1})
        Monomial(x, ex, 1)
    end

    function Monomial(x::Indeterminate, ex::Array{Int,1}, c::Number)
        length(ex) == x.degree || error("Length of exponent, $ex, must match degree of indeterminate, $x.degree")
        new(x,ex,c)
    end

end

typealias Multinomial Array{Monomial,1}

function repr(x::Monomial)
    ex = x.exponent
    c = x.coefficient
    "$c*x^$ex"
end

function repr(A::Multinomial)
    str = ""
    for a in A
        str = string(str, " + ",repr(a))
    end
    str
end

function ltlex(a::Monomial, b::Monomial)
    b.indeterminate == a.indeterminate || error("Monomials $a and $b have different indeterminates.")
    return lexcmp(b.exponent,a.exponent) == 1
end

""" collect_terms(a)

Collect terms Monomial array in place, modifying a, and return the length of the valid entries. Entries will be sorted in lexicographically increasing order.

"""
function collect_terms!(A::Multinomial)
    if length(A) <= 1 return length(A) end
    sort!(A,lt=ltlex)
    mt = 0
    current = 1
    for i in 2:length(A)
        if A[i].exponent == A[current].exponent
            # If A[i] has the same exponent as a[current], add the coefficient
            # of A[i] to that of A[current]. If mt=0 set mt=i.
            A[current].coefficient += A[i].coefficient
            if mt == 0 mt = i end
        elseif mt > 0
            # Else, if mt > 0, set A[mt] = A[i], set current = mt,
            # and set mt = mt+1. Note mt < i always.
            A[mt] = A[i]
            current = mt
            mt += 1
        else
            current = i
        end
    end
    # eliminate 0 terms
    mt = 0
    n = 0
    for i = 1:current
        if isapprox(A[i].coefficient, 0)
            if mt == 0
                mt = i
            end
        elseif mt > 0
            A[mt] = A[i]
            n = mt
            mt += 1
        else
            n = i
        end
    end
    return max(1,n) # return at least 1 term, even if 0.
end

function +(A::Multinomial, B::Multinomial)
    C = deepcopy(vcat(A,B))
    n = collect_terms!(C)
    return C[1:n]
end

function *(A::Multinomial, B::Multinomial)
    na = length(A)
    nb = length(B)
    if na==0 || nb==0 return Array(Monomial,0) end
    x = A[1].indeterminate
    C = Array(Monomial, na*nb)
    for i in 1:na
        Ae = A[i].exponent
        Ac = A[i].coefficient
        offset = (i-1)*nb
        for j in 1:nb
            C[j+offset] = Monomial(x, Ae+B[j].exponent, Ac*B[j].coefficient)
        end
    end
    n = collect_terms!(C)
    if !(typeof(n) == Int)
        @show C
        @show n
    end
    return C[1:n]
end

# syntactic sugar

function -(a::Monomial)  Monomial(a.indeterminate, a.exponent, -a.coefficient) end

function -(A::Multinomial)
    B = deepcopy(A)
    for b in B
        b.coefficient = -b.coefficient
    end
    return B
end

function -(A::Multinomial, B::Multinomial)
    C = vcat(deepcopy(A), -(B)) # -() deepcopies B
    n = collect_terms!(C)
    return C[1:n]
end

+(A::Monomial, B::Multinomial) = [A] + B
+(A::Multinomial, B::Monomial) = A + [B]
+(A::Monomial, B::Monomial)    = [A] + [B]
-(A::Monomial, B::Multinomial) = [A] - B
-(A::Multinomial, B::Monomial) = A - [B]
-(A::Monomial, B::Monomial)    = [A] - [B]
*(A::Monomial, B::Multinomial) = [A]*B
*(A::Multinomial, B::Monomial) = A*[B]
*(A::Monomial, B::Monomial)    = [A]*[B]

+(x::Number, a::Monomial) = Monomial(a.indeterminate,
                                     zeros(Int,length(a.exponent)), x) + a
+(a::Monomial, x::Number) = x+a
+(x::Number, A::Multinomial) = Monomial(A[1].indeterminate,
                                        zeros(Int,length(A[1].exponent)), x)+A
+(A::Multinomial, x::Number) = x + A
-(x::Number, a::Monomial) = (-a)+x
-(a::Monomial, x::Number) = a+(-x)
-(x::Number, A::Multinomial) = x + (-A)
-(A::Multinomial, x::Number) = A + (-x)
*(x::Number, a::Monomial) = Monomial(a.indeterminate,
                                     deepcopy(a.exponent), x*a.coefficient)
*(a::Monomial, x::Number) = x*a
*(x::Number, A::Multinomial) = [Monomial(a.indeterminate,
                                         deepcopy(a.exponent),
                                         x*a.coefficient) for a in A]
*(A::Multinomial, x::Number) = x*A

^(s::Indeterminate, xp::Array{Int,1}) = Monomial(s,xp,1)

function demo_arithmetic()
    @show x = Indeterminate(2)
    println()
    # A has an int and a real coefficient
    # B has an int and an imaginary coefficient
    # Both A and B have constant terms with coefficient 1
    @show A = 1 + 3.4*x^[1,2]
    println()
    @show B = 1 + 6im*x^[2,1]
    println()
    @show A+B
    println()
    # A-B has a 0 term which is eliminated during term collection
    @show A-B
    println()
    @show A*B
    nothing
end


text:
using StringDistances
using Clustering

# Functions for preprocessing raw names data in preparation for SOM algorithm

function generate_basis_set(name_list,set_size=1000)
   names_data = distance_matrix(name_list)
   dgram = hclust(Matrix(names_data),linkage= :average)
   cluster_indicies = Clustering.cutree(dgram, k=set_size)
   get_clustroids(name_list,cluster_indicies,set_size)
end

function generate_distance_df(name_list)
   dm = distance_matrix(name_list,my_compare=long_compare)
   df = dm |> DataFrame
   names!(df,Symbol.(name_list))
   copy(df)
end

#--------

function short_compare(a,b)
    1- compare(a,b, DamerauLevenshtein())
end

function long_compare(a,b)
    (
    (1-compare(a,b, DamerauLevenshtein())) *
    (1-compare(a,b, RatcliffObershelp())) *
          dictionaryDistance(a,b)
     )^(1/2)
end

function dictionaryDistance(a,b)
   min,max = minmax(length(a),length(b))
   for i in 1:min
      if a[i] != b[i]
         return 1 - (i-1)/max
      end
   end
   return (max-min)/max
end

#---------------

function distance_matrix(names;my_compare=short_compare)
    l = length(names)
    m = zeros(l,l)
    for i in 1:l
        for j in 1:l
            name1 = names[i]
            name2 = names[j]
            m[i,j] = my_compare(name1,name2)
        end
    end
    m
end

function get_cluster(j,names_list,cluster_indicies)
   out = []
   k=0
   for i in cluster_indicies
      k += 1
      if i == j
        push!(out,names_list[k])
      end
   end
   out
end

function get_clustroids(names_list,cluster_indicies,n=1000)
   out = []
   for i in 1:n
      k = findfirst(x->x==i,cluster_indicies)
      push!(out,names_list[k])
   end
   out
end

function generate_labels(basis_set)
   permutedims(basis_set[:,:],(2,1))
end


text:
############################################################################################
# Define the basic type system
############################################################################################
"Most abstract supertype for all SEMs"
abstract type AbstractSem end

"Supertype for all single SEMs, e.g. SEMs that have at least the fields `observed`, `imply`, `loss` and `optimizer`"
abstract type AbstractSemSingle{O, I, L, D} <: AbstractSem end

"Supertype for all collections of multiple SEMs"
abstract type AbstractSemCollection <: AbstractSem end

"Supertype for all loss functions of SEMs. If you want to implement a custom loss function, it should be a subtype of `SemLossFunction`."
abstract type SemLossFunction end

"""
    SemLoss(args...; loss_weights = nothing, ...)

Constructs the loss field of a SEM. Can contain multiple `SemLossFunction`s, the model is optimized over their sum.
See also [`SemLossFunction`](@ref).

# Arguments
- `args...`: Multiple `SemLossFunction`s.
- `loss_weights::Vector`: Weights for each loss function. Defaults to unweighted optimization.

# Examples
```julia
my_ml_loss = SemML(...)
my_ridge_loss = SemRidge(...)
my_loss = SemLoss(SemML, SemRidge; loss_weights = [1.0, 2.0])
```
"""
mutable struct SemLoss{F <: Tuple, T}
    functions::F
    weights::T
end

function SemLoss(functions...; loss_weights = nothing, kwargs...)

    if !isnothing(loss_weights)
        loss_weights = SemWeight.(loss_weights)
    else
        loss_weights = Tuple(SemWeight(nothing) for _ in 1:length(functions))
    end

    return SemLoss(
        functions,
        loss_weights
        )
end

# weights for loss functions or models. If the weight is nothing, multiplication returs second argument
struct SemWeight{T}
    w::T
end

Base.:*(x::SemWeight{Nothing}, y) = y
Base.:*(x::SemWeight, y) = x.w*y

"""
Supertype of all objects that can serve as the `optimizer` field of a SEM.
Connects the SEM to its optimization backend and controls options like the optimization algorithm.
If you want to connect the SEM package to a new optimization backend, you should implement a subtype of SemOptimizer.
"""
abstract type SemOptimizer end

"""
Supertype of all objects that can serve as the observed field of a SEM.
Pre-processes data and computes sufficient statistics for example.
If you have a special kind of data, e.g. ordinal data, you should implement a subtype of SemObserved.
"""
abstract type SemObserved end

"""
Supertype of all objects that can serve as the imply field of a SEM.
Computed model-implied values that should be compared with the observed data to find parameter estimates,
e. g. the model implied covariance or mean.
If you would like to implement a different notation, e.g. LISREL, you should implement a subtype of SemImply.
"""
abstract type SemImply end

"Subtype of SemImply for all objects that can serve as the imply field of a SEM and use some form of symbolic precomputation."
abstract type SemImplySymbolic <: SemImply end

"""
    Sem(;observed = SemObservedData, imply = RAM, loss = SemML, optimizer = SemOptimizerOptim, kwargs...)

Constructor for the basic `Sem` type.
All additional kwargs are passed down to the constructors for the observed, imply, loss and optimizer fields.

# Arguments
- `observed`: object of subtype `SemObserved` or a constructor.
- `imply`: object of subtype `SemImply` or a constructor.
- `loss`: object of subtype `SemLossFunction`s or constructor; or a tuple of such.
- `optimizer`: object of subtype `SemOptimizer` or a constructor.

Returns a Sem with fields
- `observed::SemObserved`: Stores observed data, sample statistics, etc. See also [`SemObserved`](@ref).
- `imply::SemImply`: Computes model implied statistics, like Σ, μ, etc. See also [`SemImply`](@ref).
- `loss::SemLoss`: Computes the objective and gradient of a sum of loss functions. See also [`SemLoss`](@ref).
- `optimizer::SemOptimizer`: Connects the model to the optimizer. See also [`SemOptimizer`](@ref).
"""
mutable struct Sem{O <: SemObserved, I <: SemImply, L <: SemLoss, D <: SemOptimizer} <: AbstractSemSingle{O, I, L, D}
    observed::O
    imply::I
    loss::L
    optimizer::D
end

############################################################################################
# automatic differentiation
############################################################################################
"""
    SemFiniteDiff(;observed = SemObservedData, imply = RAM, loss = SemML, optimizer = SemOptimizerOptim, kwargs...)

Constructor for `SemFiniteDiff`.
All additional kwargs are passed down to the constructors for the observed, imply, loss and optimizer fields.

# Arguments
- `observed`: object of subtype `SemObserved` or a constructor.
- `imply`: object of subtype `SemImply` or a constructor.
- `loss`: object of subtype `SemLossFunction`s or constructor; or a tuple of such.
- `optimizer`: object of subtype `SemOptimizer` or a constructor.

Returns a Sem with fields
- `observed::SemObserved`: Stores observed data, sample statistics, etc. See also [`SemObserved`](@ref).
- `imply::SemImply`: Computes model implied statistics, like Σ, μ, etc. See also [`SemImply`](@ref).
- `loss::SemLoss`: Computes the objective and gradient of a sum of loss functions. See also [`SemLoss`](@ref).
- `optimizer::SemOptimizer`: Connects the model to the optimizer. See also [`SemOptimizer`](@ref).
"""
struct SemFiniteDiff{O <: SemObserved, I <: SemImply, L <: SemLoss, D <: SemOptimizer} <: AbstractSemSingle{O, I, L, D}
    observed::O
    imply::I
    loss::L
    optimizer::D
end

############################################################################################
# ensemble models
############################################################################################
"""
    SemEnsemble(models..., optimizer = SemOptimizerOptim, weights = nothing, kwargs...)

Constructor for ensemble models.

# Arguments
- `models...`: `AbstractSem`s.
- `optimizer`: object of subtype `SemOptimizer` or a constructor.
- `weights::Vector`:  Weights for each model. Defaults to the number of observed data points.

All additional kwargs are passed down to the constructor for the optimizer field.

Returns a SemEnsemble with fields
- `n::Int`: Number of models.
- `sems::Tuple`: `AbstractSem`s.
- `weights::Vector`: Weights for each model.
- `optimizer::SemOptimizer`: Connects the model to the optimizer. See also [`SemOptimizer`](@ref).
- `identifier::Dict`: Stores parameter labels and their position.
"""
struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, D, I} <: AbstractSemCollection
    n::N
    sems::T
    weights::V
    optimizer::D
    identifier::I
end

function SemEnsemble(models...; optimizer = SemOptimizerOptim, weights = nothing, kwargs...)
    n = length(models)
    npar = n_par(models[1])

    # default weights

    if isnothing(weights)
        nobs_total = sum(n_obs.(models))
        weights = [n_obs(model)/nobs_total for model in models]
    end

    # check identifier equality
    id = identifier(models[1])
    for model in models
        if id != identifier(model)
            throw(ErrorException("The identifier of your models do not match. \n
            Maybe you tried to specify models of an ensemble via ParameterTables. \n
            In that case, you may use RAMMatrices instead."))
        end
    end

    # optimizer
    if !isa(optimizer, SemOptimizer)
        optimizer = optimizer(;kwargs...)
    end

    return SemEnsemble(
        n,
        models,
        weights,
        optimizer,
        id
        )
end

"""
    n_models(ensemble::SemEnsemble) -> Integer

Returns the number of models in an ensemble model.
"""
n_models(ensemble::SemEnsemble) = ensemble.n
"""
    models(ensemble::SemEnsemble) -> Tuple{AbstractSem}

Returns the models in an ensemble model.
"""
models(ensemble::SemEnsemble) = ensemble.sems
"""
    weights(ensemble::SemEnsemble) -> Vector

Returns the weights of an ensemble model.
"""
weights(ensemble::SemEnsemble) = ensemble.weights
"""
    optimizer(ensemble::SemEnsemble) -> SemOptimizer

Returns the optimizer part of an ensemble model.
"""
optimizer(ensemble::SemEnsemble) = ensemble.optimizer

############################################################################################
# additional methods
############################################################################################
"""
    observed(model::AbstractSemSingle) -> SemObserved

Returns the observed part of a model.
"""
observed(model::AbstractSemSingle) = model.observed

"""
    imply(model::AbstractSemSingle) -> SemImply

Returns the imply part of a model.
"""
imply(model::AbstractSemSingle) = model.imply

"""
    loss(model::AbstractSemSingle) -> SemLoss

Returns the loss part of a model.
"""
loss(model::AbstractSemSingle) = model.loss

"""
    optimizer(model::AbstractSemSingle) -> SemOptimizer

Returns the optimizer part of a model.
"""
optimizer(model::AbstractSemSingle) = model.optimizer


text:
module JsonBuilder

import JSON: print
import Base: push!, getindex, done, next, eof

export @json

macro json(s)
    quote
        io = IOBuffer()
        $(json(s)...)
        takebuf_string(io)
    end
end

macro json(io, s)
    quote
        io = $io
        $(json(s)...)
        nothing
    end
end

function json(s)
    if !isa(s, Expr) || s.head != :string
        isa(s, String) ? (s = Expr(:string, s)) :
        error("invalid invocation of json macro")
    end

    if !isa(s.args[1], String) # quick fix for invocations like `@json "$x"`
        unshift!(s.args, " ")
    end

    x = Parser([], s.args, 1, 1)
    parse_value!(x)
    code_gen(x)
end

abstract Token
type ObjectMixin <: Token end
type ArrayMixin  <: Token end
type EOF <: Token end
type Var <: Token content end
type Str <: Token content end
type Raw <: Token content end

type Parser
    result::Vector{Token}; s; i::Int; j::Int
end

push!(p::Parser, x) = push!(p.result, x)
getindex(p::Parser, x) = p.s[x]
getindex(p::Parser, x, y) = p.s[x][y]
done(p::Parser) = p.j > length(p[p.i])
eof(p::Parser)  = p.i >= length(p.s)
next(p::Parser) = !done(p) ? p[p.i][p.j] :
                  !eof(p)  ? error("unexpected interpolation") :
                             error("unexpected EOF")

function parse_space!(p::Parser)
    while !done(p)
        if !isspace(next(p))
            if next(p) == '#'
                while !done(p) && next(p) != '\n'
                    p.j += 1
                end
            else
                return
            end
        else
            p.j += 1
        end
    end
end

function parse_value!(p::Parser)
    parse_space!(p)
    if done(p)
        if eof(p)
            error("malformed json")
        else
            push!(p, Var(p[p.i+1]))
            p.i, p.j = p.i + 2, 1
        end
    elseif next(p) == '{'
        parse_object!(p)
    elseif next(p) == '['
        parse_array!(p)
    elseif next(p) in ('"', '\'')
        parse_string!(p)
    else
        parse_literal!(p)
    end
end

function parse_literal!(p::Parser)
    start = p.j
    while !done(p) && (next(p) in "+-.0123456789Eaeflnrstu") # true, false, null, numbers; just hope the user won't misspell them :)
        p.j += 1
    end
    push!(p, Raw(p[p.i, start:p.j-1]))
end

function parse_object!(p::Parser)
    parse_char!('{', p)
    parse_space!(p)
    if !done(p) && next(p) == '}'
        parse_char!('}', p)
        return
    end

    while true
        parse_key!(p)
        if p.j+2 <= length(p[p.i]) && p[p.i, p.j:p.j+2] == "..."
            p.j += 3
            push!(p, ObjectMixin())
        else
            parse_space!(p)
            parse_char!(':', p)
            parse_space!(p)
            parse_value!(p)
        end
        parse_space!(p)
        if next(p) == '}'
            parse_char!('}', p)
            return
        else
            parse_char!(',', p)
            parse_space!(p)
        end
    end
end

function parse_array!(p::Parser)
    parse_char!('[', p)
    parse_space!(p)
    if !done(p) && next(p) == ']'
        parse_char!(']', p)
        return
    end

    while true
        parse_value!(p)
        if p.j+2 <= length(p[p.i]) && p[p.i, p.j:p.j+2] == "..."
            p.j += 3
            push!(p, ArrayMixin())
        end
        parse_space!(p)
        if next(p) == ']'
            parse_char!(']', p)
            return
        else
            parse_char!(',', p)
            parse_space!(p)
        end
    end
end

function parse_string!(p::Parser)
    start, q = p.j, next(p)
    p.j += 1
    while !done(p)
        if next(p) == q
            push!(p, Str(p[p.i, start+1:p.j-1])) # JSON.print will add quotes around them, so no need to insert quotes here
            p.j += 1
            return
        end
        p.j += 1
    end
    error("non-terminated string literal. hint: if you need interpolations within string literals, interpolate the whole string like `\$(\"some \$var\")`")
end

function parse_key!(p::Parser)
    if done(p)
        push!(p, Var(p[p.i+1]))
        p.i, p.j = p.i + 2, 1
    elseif next(p) in ('"', '\'')
        parse_string!(p)
    else # simple form, without quote marks
        start = p.j
        while !done(p) && any(next(p) in x for x in ('a':'z', 'A':'Z', '0':'9', "\$_"))
            p.j += 1
        end
        push!(p, Raw('"'), Raw(p[p.i, start:p.j-1]), Raw('"'))
    end
end

function parse_char!(x::Char, p::Parser)
    if next(p) == x
        p.j += 1
        push!(p, Raw(x))
    else
        error("unexpected $(next(p))")
    end
end

macro gen(ex)
    :( push!(result, $(Expr(:quote, ex))) )
end

function code_gen(x)
    i, result, x = 1, Expr[], x.result

    push!(x, EOF()) # ensure that x[i+1] always valid

    while x[i] != EOF()
        if isa(x[i], Str)
            @gen print(io, $(x[i].content))
        elseif isa(x[i], Raw)
            if isa(x[i+1], Raw)
                x[i+1] = Raw(string(x[i].content, x[i+1].content))
            else
                @gen write(io, $(x[i].content))
            end
        elseif isa(x[i], Var)
            if isa(x[i+1], ObjectMixin)
                @gen join(io, $(esc(x[i].content)), ',') do io, x
                    k, v = x
                    print(io, string(k))
                    write(io, ':')
                    print(io, v)
                end
                i += 1 # skip the Mixin Token
            elseif isa(x[i+1], ArrayMixin)
                @gen join(print, io, $(esc(x[i].content)), ',')
                i += 1 # skip the Mixin Token
            else
                @gen print(io, $(esc(x[i].content)))
            end
        else
            error("BUG, Please fire an issue with your json template string.")
        end

        i += 1
    end

    result
end

function join(f, io::IO, iter, delim)
    i = start(iter)

    if !done(iter, i)
        str, i  = next(iter, i)
        f(io, str)
    end

    while !done(iter, i)
        write(io, delim)
        str, i  = next(iter,i)
        f(io, str)
    end
end

end

using JsonBuilder


text:
using DataStructures
using TypedPolynomials

# Utility functions
LT(f) = leadingterm(f)
LTS(F) = [LT(f) for f in F]
LM(f) = leadingmonomial(f)
LMS(F) = unique([LM(f) for f in F])
LC(f) = leadingcoefficient(f)
LCS(F) = [LC(f) for f in F]
LCMLM(f1, f2) = lcm(LM(f1), LM(f2))
LCMLM(ftuple) = lcm(LM(ftuple[1]), LM(ftuple[2]))
LCMLT(f1, f2) = lcm(LT(f1), LT(f2))
LCMLT(ftuple) = lcm(LT(ftuple[1]), LT(ftuple[2]))
Lpart(f1, f2) = div(LCMLM(f1, f2), LT(f1)) * f1
Mpart(f1, f2) = div(LCMLM(f1, f2), LT(f1))
spoly(f1, f2) = Lpart(f1, f2) - Lpart(f2, f1)
allpairs(s) = [(i, j) for i in 1:s for j in i+1:s if i ≠ j]
Mon(L) = unique(collect(Iterators.flatten([monomials(l) for l in L])))
row2poly(mat, row, monoms) = polynomial(mat[row, :] .* Mon(monoms))
rows(N, monoms) = [row2poly(N, i, monoms) for i in 1:size(N)[1]]

"Return sorted minimum degree pairs by LCMLM"
function mindegpairs(F::Array{T}) where {T <: AbstractPolynomialLike}
    map(x->x[1], sort([((f1, f2), degree(LCMLM(f1, f2)))
                       for (i, f1) in enumerate(F) for f2 in F[i+1:end] if f1 ≠ f2],
                      lt=(mdp1, mdp2) -> mdp1[2] < mdp2[2] ? true : false,
                      rev=true))
end

"Normal selection strategy"
function Select(B, G)
    min = Inf
    for pair in B
        d = degree(LCMLM(G[pair[1]], G[pair[2]]))
        if d < min
            min = d
        end
    end

    minpairs = []
    polys = []
    for pair in B
        d = degree(LCMLM(G[pair[1]], G[pair[2]]))
        if d == min
            push!(minpairs, pair)
            push!(polys, (G[pair[1]], G[pair[2]]))
        end
    end
    (minpairs, polys)
end

"Basic F4 from CLO"
function F4(F::Array{T}) where {T <: AbstractPolynomialLike}
    G = copy(F)
    t = length(G)
    B = allpairs(t)
    count = 0
    while !isempty(B)
        Bprime, polys = Select(B, G)
        B = setdiff(B, Bprime)
        L₁ = [Lpart(pair[1], pair[2]) for pair in polys]
        L₂ = [Lpart(pair[2], pair[1]) for pair in polys]
        L = unique(cat(L₁, L₂, dims=1))
        M, monoms = ComputeM(L, G)
        N = RREF(M)
        Nplus = [p for p in rows(N, monoms) if LM(p) ∉ LMS(rows(M, monoms)) && !isapproxzero(p)]
        for n in Nplus
            G = union(G, [n])
            t = length(G)
            B = union(B, [(i, t) for i in 1:t-1])
        end
    end
    G
end

function ComputeM(L, G)
    H = copy(L)
    done = Set(LMS(H))
    MonH = Set(Mon(H))

    while MonH ≠ done
        xβ = maximum(setdiff(Mon(H), done))

        done = union(done, [xβ])
        for f in G
            if divides(LM(f), xβ)
                H = union(H, [div(xβ, LM(f)) * f])
                MonH = union(MonH, Mon(H))
            end
        end
    end
    t = typeof(coefficients(G[1])[1])
    M = zeros(t, length(H), length(MonH))
    for (i, poly) in enumerate([zip(coefficients(h), monomials(h)) for h in H])
        for (coeff, monom) in poly
            for (j, m) in enumerate(MonH)
                if monom == m
                    M[i, j] = coeff
                end
            end
        end
    end
    return M, MonH
end

function REF(M)
    A = copy(M)
    r = 0
    rows, cols = size(A)
    for i in 1:cols
        piv_found = false
        for j in r+1:rows
            if A[j, i] ≠ 0
                r += 1
                A[j, :], A[r, :] = A[r, :], A[j, :]
                A[r, :] *= A[r, i]^-1
                piv_found = true
                break
            end
        end
        if piv_found
            for j in r+1:rows
                A[j, :] -= A[j, i] * A[r, :]
            end
        end
    end
    A
end

function RREF(M)
    A = copy(M)
    rows, cols = size(A)
    r = 0
    for i in 1:cols
        piv_found = false
        for j in r+1:rows
            if A[j, i] ≠ 0
                r += 1
                A[j, :], A[r, :] = A[r, :], A[j, :]
                A[r, :] *= A[r,i]^-1
                piv_found = true
                break
            end
        end
        if piv_found
            for j in [n for n in 1:rows if n != r]
                A[j, :] -= A[j, i] * A[r, :]
            end
        end
    end
    A
end

"Reduce a polynomial modulo P"
function REDPOL(f, P)
    qp = [0 for p in P]
    g = copy(f)
    while any([divides(LT(p), LT(g)) for p in P])
        for (i, p) in enumerate(P)
            if divides(LT(p), LT(g))
                m = div(g, p)
                g -= m*p
                qp[i] += m
            end
        end
    end
    @assert maximum(LTS(qp .* P)) ≤ LT(f)
    (qp, g)
end

"Reduce a system of Polynomials"
function REDUCTION(P)
    Q = copy(P)
    for p in Q
        if any([divides(LT(q),LT(p)) for q in setdiff(Q, [p])])
            setdiff!(Q, [p])
            h = rem(p, Q)
            if ~isapproxzero(h)
                union!(Q, [h])
            end
        end
    end
    [LC(q)^(-1) * q for q in Q]
end

"Test if a system of polynomials is a Groebner basis"
function GROEBNERTEST(G)
    B = mindegpairs(G)
    while !isempty(B)
        g1, g2 = pop!(B)
        h = spoly(g1, g2)
        h = rem(h, G)
        if !isapproxzero(h)
            println("h reduced by G: $h")
            return false
        end
    end
    true
end

"Basic Buchberger algorithm"
function GROEBNER(F)
    G = copy(F)
    B = mindegpairs(G)
    while !isempty(B)
        g1, g2 = pop!(B)
        h = spoly(g1, g2)
        h = rem(h, G)
        if !isapproxzero(h)
            for g in G
                union!(B, Set([(g, h)]))
            end
            union!(G, [h])
        end
    end
    G
end

"Produce unique Reduced Groebner Basis"
function REDGROEBNER(G)
    F = copy(G)
    H = Vector{eltype(F)}()
    while ~isempty(F)
        f0 = pop!(F)
        if all([!divides(LM(f), LM(f0)) for f in F]) &&
            all([!divides(LM(h), LM(f0)) for h in H])
            union!(H, [f0])
        end
    end
    REDUCTION(H)
end

function UPDATE(Gold, Bold, h::Polynomial)
    C = [(h, g) for g in Gold]
    D = Vector{Tuple{eltype(h), eltype(h)}}()
    while !isempty(C)
        h, g1 = pop!(C)
        if isdisjoint(variables(LM(h)), variables(LM(g1))) ||
            (all([!divides(LCMLM(p), LCMLM(h, g1)) for p in C]) &&
            all([!divides(LCMLM(p), LCMLM(h, g1)) for p in D]))
            push!(D, (h, g1))
        end
    end

    E = Vector{eltype(D)}()
    while !isempty(D)
        h, g = pop!(D)
        if !isdisjoint(variables(LM(h)), variables(LM(g)))
            push!(E, (h, g))
        end
    end

    Bnew = Vector{eltype(E)}()
    while !isempty(Bold)
        g1, g2 = pop!(Bold)
        if !divides(LM(h), LCMLM(g1, g2)) ||
            LCMLM(g1, h) == LCMLM(g1, g2) ||
            LCMLM(h, g2) == LCMLM(g1, g2)
            push!(Bnew, (g1, g2))
        end
    end

    union!(Bnew, E)
    Gnew = Vector{eltype(Gold)}()
    while !isempty(Gold)
        g = pop!(Gold)
        if !divides(LM(h), LM(g))
            union!(Gnew, [g])
        end
    end
    union!(Gnew, [h])
    (Gnew, Bnew)
end

function GROEBNERNEW2(F::Array{T}) where {T <: AbstractPolynomialLike}
    K = copy(F)
    G = Vector{eltype(F)}()
    B = Set{Tuple{eltype(F), eltype(F)}}()
    while !isempty(K)
        h = pop!(K)
        (G, B) = UPDATE(G, B, h)
    end
    while !isempty(B)
        g1, g2 = pop!(B)
        h = spoly(g1, g2)
        h = rem(h, G)
        if !isapproxzero(h)
            (G, B) = UPDATE(G, B, h)
        end
    end
    G
end

function GEBAUERMOLLER(F::Array{T}) where {T <: AbstractPolynomialLike}
    REDGROEBNER(GROEBNERNEW2(F))
end

function EXTGROEBNER(F::Array{T}) where {T <: AbstractPolynomialLike}
    G = copy(F)
    #𝒢 = [(f, REDPOL(f, G)) for f in F]
    B = mindegpairs(G)

    while !isempty(B)
        g1, g2 = pop!(B)
        g = spoly(g1, g2)
        (qp, h) = REDPOL(g, G)
        if !isapproxzero(h)
            union!(B, [(p, h) for p in G])
     #       union!(𝒢, (h, REDPOL(h, F)))
            union!(G, [h])
        end
    end

    ℱ = [(f, REDPOL(f, G)) for f in F]
    (G, ℱ)
end

function EXTREDGROEBNER(F::Array{T}) where {T <: AbstractPolynomialLike}
    G = copy(F)
    𝒢 = [(f, f == g ? 1 : 0) for f in F for g in G]
    B = mindegpairs(G)

    while !isempty(B)
        g1, g2 = pop!(B)
        m1 = Mpart(g1, g2)
        m2 = Mpart(g2, g1)
        g = m1*g1 - m2*g2
        (qp, h) = REDPOL(g, setdiff(G, [g]))
        if !isapproxzero(h)
            union!(B, [(p, h) for p in setdiff(G, [g])])
            union!(𝒢, [(h, REDPOL(f, setdiff(G, [g]))[1]) for f in F])
            union!(G, [h])
        end
    end

end

@polyvar x[1:4]

fz = 4x[1]*x[2]^2*x[3] + 4x[3]^2 - 5x[1]^3 + 7x[1]^2*x[3]^2

f1 = x[1]^2 + x[1]*x[2] - 1.0
f2 = x[1]^2 - x[3]^2
f3 = x[1]*x[2] + 1.0
I = [f1, f2, f3]

g1 = (1//1)*x[1]^2 + x[2]^2 + x[3]^2 - (1//1)
g2 = (1//1)*x[1]^2 - x[2] + x[3]^2
g3 = (1//1)*x[1] - x[3]

G = [g1, g2, g3]

h1 = 3.0*x[1]^2 + 2.0*x[2]*x[3] - 2.0*x[1]*x[4]
h2 = 2.0*x[1]*x[3] - 2.0*x[2]*x[4]
h3 = 2.0*x[1]*x[2] - 2.0*x[3] - 2.0*x[3]*x[4]
h4 = x[1]^2 + x[2]^2 + x[3]^2 - 1.0

H = [h1, h2, h3, h4]

z1 = x[1]*x[2] + x[2]^2
z2 = x[3] + 1//1
z3 = x[2]^3*x[3] + x[1]^2
z4 = x[3]^5*x[2]^2 + 4//1

Z = [z1, z2, z3, z4]

c1 = x[1] + x[2] + x[3] + x[4]
c2 = x[1]*x[2] + x[2]*x[3] + x[3]*x[4] + x[1]*x[4]
c3 = x[1]*x[2]*x[3] + x[2]*x[3]*x[4] + x[3]*x[4]*x[1] + x[4]*x[1]*x[2]
c4 = x[1]*x[2]*x[3]*x[4] - 1//1
Cyclic4 = [c1, c2, c3, c4]


text:
using GeoMakie, Test

@testset "GeoMakie" begin
    lons = -180:180
    lats = -90:90
    field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

    fig = Figure()
    ax = GeoAxis(fig[1,1], coastlines=true)
    el = surface!(ax, lons, lats, field; shading = false)
    @test true
    # display(fig)
end


text:
import Base: rand,
             ∈,
             isempty,
             convert

export HalfSpace, LinearConstraint,
       an_element,
       constrained_dimensions,
       halfspace_left, halfspace_right

"""
    HalfSpace{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a (closed) half-space of the form ``a⋅x ≤ b``.

### Fields

- `a` -- normal direction (non-zero)
- `b` -- constraint

### Examples

The set ``y ≥ 0`` in the plane:

```jldoctest
julia> HalfSpace([0, -1.], 0.)
HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0], 0.0)
```
"""
struct HalfSpace{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
    a::VN
    b::N

    function HalfSpace{N, VN}(a::VN, b::N
                             ) where {N<:Real, VN<:AbstractVector{N}}
        @assert !iszero(a) "a half-space needs a non-zero normal vector"
        return new{N, VN}(a, b)
    end
end

isoperationtype(::Type{<:HalfSpace}) = false

# convenience constructor without type parameter
HalfSpace(a::VN, b::N) where {N<:Real, VN<:AbstractVector{N}} =
    HalfSpace{N, VN}(a, b)

"""
    LinearConstraint

Alias for `HalfSpace`
"""
const LinearConstraint = HalfSpace

# --- LazySet interface functions ---

"""
    dim(hs::HalfSpace)::Int

Return the dimension of a half-space.

### Input

- `hs` -- half-space

### Output

The ambient dimension of the half-space.
"""
function dim(hs::HalfSpace)::Int
    return length(hs.a)
end

"""
    ρ(d::AbstractVector{N}, hs::HalfSpace{N})::N where {N<:Real}

Evaluate the support function of a half-space in a given direction.

### Input

- `d`  -- direction
- `hs` -- half-space

### Output

The support function of the half-space.
If the set is unbounded in the given direction, the result is `Inf`.
"""
function ρ(d::AbstractVector{N}, hs::HalfSpace{N})::N where {N<:Real}
    v, unbounded = σ_helper(d, Hyperplane(hs.a, hs.b); error_unbounded=false,
                            halfspace=true)
    if unbounded
        return N(Inf)
    end
    return dot(d, v)
end

"""
    σ(d::AbstractVector{N}, hs::HalfSpace{N}) where {N<:Real}

Return the support vector of a half-space.

### Input

- `d`  -- direction
- `hs` -- half-space

### Output

The support vector in the given direction, which is only defined in the
following two cases:
1. The direction has norm zero.
2. The direction is the half-space's normal direction.
In both cases the result is any point on the boundary (the defining hyperplane).
Otherwise this function throws an error.
"""
function σ(d::AbstractVector{N}, hs::HalfSpace{N}) where {N<:Real}
    v, unbounded = σ_helper(d, Hyperplane(hs.a, hs.b); error_unbounded=true,
                            halfspace=true)
    return v
end

"""
    isbounded(hs::HalfSpace)::Bool

Determine whether a half-space is bounded.

### Input

- `hs` -- half-space

### Output

`false`.
"""
function isbounded(::HalfSpace)::Bool
    return false
end

"""
    isuniversal(hs::HalfSpace{N}, [witness]::Bool=false
               )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a half-space is universal.

### Input

- `P`       -- half-space
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `false`
* If `witness` option is activated: `(false, v)` where ``v ∉ P``

### Algorithm

Witness production falls back to `isuniversal(::Hyperplane)`.
"""
function isuniversal(hs::HalfSpace{N}, witness::Bool=false
                    )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    if witness
        return isuniversal(Hyperplane(hs.a, hs.b), true)
    else
        return false
    end
end

"""
    an_element(hs::HalfSpace{N})::Vector{N} where {N<:Real}

Return some element of a half-space.

### Input

- `hs` -- half-space

### Output

An element on the defining hyperplane.
"""
function an_element(hs::HalfSpace{N})::Vector{N} where {N<:Real}
    return an_element_helper(Hyperplane(hs.a, hs.b))
end

"""
    ∈(x::AbstractVector{N}, hs::HalfSpace{N})::Bool where {N<:Real}

Check whether a given point is contained in a half-space.

### Input

- `x` -- point/vector
- `hs` -- half-space

### Output

`true` iff ``x ∈ hs``.

### Algorithm

We just check if ``x`` satisfies ``a⋅x ≤ b``.
"""
function ∈(x::AbstractVector{N}, hs::HalfSpace{N})::Bool where {N<:Real}
    return dot(x, hs.a) <= hs.b
end

"""
    rand(::Type{HalfSpace}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::HalfSpace{N}

Create a random half-space.

### Input

- `HalfSpace` -- type for dispatch
- `N`         -- (optional, default: `Float64`) numeric type
- `dim`       -- (optional, default: 2) dimension
- `rng`       -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`      -- (optional, default: `nothing`) seed for reseeding

### Output

A random half-space.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
Additionally, the constraint `a` is nonzero.
"""
function rand(::Type{HalfSpace};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing
             )::HalfSpace{N}
    rng = reseed(rng, seed)
    a = randn(rng, N, dim)
    while iszero(a)
        a = randn(rng, N, dim)
    end
    b = randn(rng, N)
    return HalfSpace(a, b)
end

"""
    isempty(hs::HalfSpace)::Bool

Return if a half-space is empty or not.

### Input

- `hs` -- half-space

### Output

`false`.
"""
function isempty(hs::HalfSpace)::Bool
    return false
end

"""
    constraints_list(hs::HalfSpace{N}) where {N<:Real}

Return the list of constraints of a half-space.

### Input

- `hs` -- half-space

### Output

A singleton list containing the half-space.
"""
function constraints_list(hs::HalfSpace{N}) where {N<:Real}
    return [hs]
end

"""
    constraints_list(A::AbstractMatrix{N}, b::AbstractVector{N}) where {N<:Real}

Convert a matrix-vector representation to a linear-constraint representation.

### Input

- `A` -- matrix
- `b` -- vector

### Output

A list of linear constraints.
"""
function constraints_list(A::AbstractMatrix{N}, b::AbstractVector{N}
                         ) where {N<:Real}
    m = size(A, 1)
    @assert m == length(b) "a matrix with $m rows is incompatible with a " *
                           "vector of length $(length(b))"

    return [LinearConstraint(A[i, :], b[i]) for i in 1:m]
end

"""
    constrained_dimensions(hs::HalfSpace{N})::Vector{Int} where {N<:Real}

Return the indices in which a half-space is constrained.

### Input

- `hs` -- half-space

### Output

A vector of ascending indices `i` such that the half-space is constrained in
dimension `i`.

### Examples

A 2D half-space with constraint ``x1 ≥ 0`` is constrained in dimension 1 only.
"""
function constrained_dimensions(hs::HalfSpace{N})::Vector{Int} where {N<:Real}
    return nonzero_indices(hs.a)
end

"""
    halfspace_left(p::AbstractVector{N},
                   q::AbstractVector{N})::HalfSpace{N} where {N<:Real}

Return a half-space describing the 'left' of a two-dimensional line segment
through two points.

### Input

- `p` -- first point
- `q` -- second point

### Output

The half-space whose boundary goes through the two points `p` and `q` and which
describes the left-hand side of the directed line segment `pq`.

### Algorithm

The implementation is simple: the half-space ``a⋅x ≤ b`` is calculated as
`a = [dy, -dx]`, where ``d = (dx, dy)`` denotes the line segment
`pq`, that is, ``\\vec{d} = \\vec{p} - \\vec{q}``, and `b = dot(p, a)`.

### Examples

The left half-space of the "east" and "west" directions in two-dimensions are
the upper and lower half-spaces:

```jldoctest halfspace_left
julia> using LazySets: halfspace_left

julia> halfspace_left([0.0, 0.0], [1.0, 0.0])
HalfSpace{Float64,Array{Float64,1}}([0.0, -1.0], 0.0)

julia> halfspace_left([0.0, 0.0], [-1.0, 0.0])
HalfSpace{Float64,Array{Float64,1}}([0.0, 1.0], 0.0)
```

We create a box from the sequence of line segments that describe its edges:

```jldoctest halfspace_left
julia> H1 = halfspace_left([-1.0, -1.0], [1.0, -1.0]);

julia> H2 = halfspace_left([1.0, -1.0], [1.0, 1.0]);

julia> H3 = halfspace_left([1.0, 1.0], [-1.0, 1.0]);

julia> H4 = halfspace_left([-1.0, 1.0], [-1.0, -1.0]);

julia> H = HPolygon([H1, H2, H3, H4]);

julia> B = BallInf([0.0, 0.0], 1.0);

julia> B ⊆ H && H ⊆ B
true
```
"""
function halfspace_left(p::AbstractVector{N},
                        q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
    @assert length(p) == length(q) == 2 "the points must be two-dimensional"
    @assert p != q "the points must not be equal"
    a = [q[2] - p[2], p[1] - q[1]]
    return HalfSpace(a, dot(p, a))
end

"""
    halfspace_right(p::AbstractVector{N},
                    q::AbstractVector{N})::HalfSpace{N} where {N<:Real}

Return a half-space describing the 'right' of a two-dimensional line segment
through two points.

### Input

- `p` -- first point
- `q` -- second point

### Output

The half-space whose boundary goes through the two points `p` and `q` and which
describes the right-hand side of the directed line segment `pq`.

### Algorithm

See the documentation of `halfspace_left`.
"""
function halfspace_right(p::AbstractVector{N},
                         q::AbstractVector{N})::HalfSpace{N} where {N<:Real}
    return halfspace_left(q, p)
end

"""
    is_tighter_same_dir_2D(c1::LinearConstraint{N},
                           c2::LinearConstraint{N}) where {N<:Real}

Check if the first of two two-dimensional constraints with equivalent normal
direction is tighter.

### Input

- `c1`     -- first linear constraint
- `c2`     -- second linear constraint
- `strict` -- (optional; default: `false`) check for strictly tighter
              constraints?

### Output

`true` iff the first constraint is tighter.
"""
function is_tighter_same_dir_2D(c1::LinearConstraint{N},
                                c2::LinearConstraint{N};
                                strict::Bool=false) where {N<:Real}
    @assert dim(c1) == dim(c2) == 2 "this method requires 2D constraints"
    @assert c1.a <= c2.a <= c1.a "the constraints must have the same " *
        "normal direction"

    lt = strict ? (<) : (<=)
    if c1.a[1] == zero(N)
        @assert c2.a[1] == zero(N)
        return lt(c1.b, c1.a[2] / c2.a[2] * c2.b)
    end
    return lt(c1.b, c1.a[1] / c2.a[1] * c2.b)
end

function _linear_map_hrep(M::AbstractMatrix{N}, P::HalfSpace{N}, use_inv::Bool;
                          inverse::Union{Nothing, AbstractMatrix{N}}=nothing
                         ) where {N<:Real}
    constraint = _linear_map_hrep_helper(M, P, use_inv; inverse=inverse)[1]
    return HalfSpace(constraint.a, constraint.b)
end

"""
    translate(hs::HalfSpace{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:Real}

Translate (i.e., shift) a half-space by a given vector.

### Input

- `hs`    -- half-space
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated half-space.

### Notes

The normal vectors of the halfspace (vector `a` in `a⋅x ≤ b`) is shared with the
original halfspace if `share == true`.

### Algorithm

A half-space ``a⋅x ≤ b`` is transformed to the half-space ``a⋅x ≤ b + a⋅v``.
In other words, we add the dot product ``a⋅v`` to ``b``.
"""
function translate(hs::HalfSpace{N}, v::AbstractVector{N}; share::Bool=false
                  ) where {N<:Real}
    @assert length(v) == dim(hs) "cannot translate a $(dim(hs))-dimensional " *
                                 "set by a $(length(v))-dimensional vector"
    a = share ? hs.a : copy(hs.a)
    b = hs.b + dot(hs.a, v)
    return HalfSpace(a, b)
end


text:
module MPI

using Libdl, Serialization
using Requires
using DocStringExtensions

function serialize(x)
    s = IOBuffer()
    Serialization.serialize(s, x)
    take!(s)
end

function deserialize(x)
    s = IOBuffer(x)
    Serialization.deserialize(s)
end

primitive type SentinelPtr Sys.WORD_SIZE
end

primitive type MPIPtr Sys.WORD_SIZE
end
Base.cconvert(::Type{MPIPtr}, x::SentinelPtr) = x
Base.unsafe_convert(::Type{MPIPtr}, x::SentinelPtr) = reinterpret(MPIPtr, x)

function _doc_external(fname)
"""
- `$fname` man page: [OpenMPI](https://www.open-mpi.org/doc/current/man3/$fname.3.php), [MPICH](https://www.mpich.org/static/docs/latest/www3/$fname.html)
"""
end

include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
include("error.jl")
include("handle.jl")
include("info.jl")
include("comm.jl")
include("environment.jl")
include("datatypes.jl")
include("buffers.jl")
include("operators.jl")
include("pointtopoint.jl")
include("collective.jl")
include("topology.jl")
include("onesided.jl")
include("io.jl")

include("deprecated.jl")

function __init__()
    @static if Sys.isunix()
        # need to open libmpi with RTLD_GLOBAL flag for Linux, before
        # any ccall cannot use RTLD_DEEPBIND; this leads to segfaults
        # at least on Ubuntu 15.10
        Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
    end

    # disable UCX memory hooks since it can mess up dlopen
    # https://github.com/openucx/ucx/issues/4001
    ENV["UCX_MEM_MMAP_RELOC"] = "no"
    ENV["UCX_MEM_MALLOC_HOOKS"] = "no"
    ENV["UCX_MEM_MALLOC_RELOC"] = "no"
    ENV["UCX_MEM_EVENTS"] = "no"

    if filesize(dlpath(libmpi)) != libmpi_size
        error("MPI library has changed, re-run Pkg.build(\"MPI\")")
    end

    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" include("cuda.jl")
end

end


text:
#
# PseudoArrays.jl -
#
# Types and methods to facilitate the definition of custom array-like types.
#

module PseudoArrays

export
    CartesianArray,
    PseudoArray,
    LinearArray

using Base: axes1, elsize, tail, OneTo, throw_boundserror, @propagate_inbounds
import Base: getindex, setindex!, checkbounds

"""

Abstract type `PseudoArray{T,N,S}` is to be derived by types that want to
provide an array-like interface.  Parameter `T` is the element type, parameter
`N` is the number of dimensions and parameter `S` is the index style:
`IndexCartesian` or `IndexLinear`.

!!! note
    The indexing style must be part of the signature because it must be
    possible to call `IndexStyle()` on the data type not the instance.  Another
    possibility would have been to have the type of the embedded array be part
    of the signature but this is more restrictive.

Alias `LinearArray{T,N}` is an abstract type that can be derived by types that
want to provide an array-like interface with array values stored in an array
whose index style is linear.

Usage can be as simple as:

    struct CustomArray{T,N,...} <: LinearArray{T,N}
        arr::Array{T,N} # can be any array type with linear index style
        ...             # anything else
    end

    @inline Base.parent(A::CustomArray) = A.arr

As a result, instances of `CustomArray{T,N}` will be seen as instances of
`AbstractArray{T,N}` and behave as if they implement linear indexing.  Apart
from the needs to extend the `Base.parent` method, the interface to
`LinearArray{T,N}` should provide any necessary methods for indexation, getting
the dimensions, the element type, *etc.* for the derived custom type.  You may
however override these definitions by more optimized or more suitable methods
specialized for your custom array-like type.

Similarly, alias `CartesianArray{T,N}` is an abstract type that can be derived
by types that want to provide an array-like interface with array values stored
in an array whose index style is Cartesian.  For such array-like object, index
checking requires an efficient implementation of the `Base.axes()` method which
you may have to specialize.  The default implementation is:

    @inline Base.axes(A::PseudoArray) = axes(parent(A))

"""
abstract type PseudoArray{T,N,S<:IndexStyle} <: AbstractArray{T,N} end
const CartesianArray{T,N} = PseudoArray{T,N,IndexCartesian}
const LinearArray{T,N} = PseudoArray{T,N,IndexLinear}

@doc @doc(PseudoArray) CartesianArray
@doc @doc(PseudoArray) LinearArray

# Make PseudoArray instances behave like arrays (indexing is considered later).
#Base.eltype(::PseudoArray{T,N}) where {T,N} = T # FIXME: not needed
#Base.ndims(::PseudoArray{T,N}) where {T,N} = N # FIXME: not needed
@inline Base.length(A::PseudoArray) = length(parent(A))
@inline Base.size(A::PseudoArray) = size(parent(A))
Base.size(A::PseudoArray, d) = size(parent(A), d)
@inline Base.axes(A::PseudoArray) = axes(parent(A))
Base.axes(A::PseudoArray, d) = axes(parent(A), d)
@inline Base.axes1(A::PseudoArray) = axes1(parent(A))
@inline Base.IndexStyle(::Type{<:PseudoArray{T,N,S}}) where {T,N,S} = S()
Base.parent(A::T) where {T<:PseudoArray} =
    error(string("method parent() must be extended for instances of ", T))
Base.elsize(::Type{<:PseudoArray{T,N}}) where {T,N} = elsize(Array{T,N})
Base.sizeof(A::PseudoArray) = sizeof(parent(A))
Base.pairs(S::IndexCartesian, A::PseudoArray) = pairs(S, parent(A))
Base.pairs(S::IndexLinear, A::PseudoArray) = pairs(S, parent(A))

# Make LinearArray instances efficient iterators.
@inline Base.iterate(A::LinearArray, i=1) =
    ((i % UInt) - 1 < length(A) ? (@inbounds A[i], i + 1) : nothing)

@inline function getindex(A::LinearArray{T,N}, i::Int) where {T,N}
    @boundscheck checkbounds(A, i)
    @inbounds r = parent(A)[i]
    r
end

@inline function getindex(A::CartesianArray{T,N},
                          I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    @inbounds r = parent(A)[I...]
    r
end

@inline function setindex!(A::LinearArray{T,N}, val, i::Int) where {T,N}
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[i] = val
    A
end

@inline function setindex!(A::CartesianArray{T,N}, val,
                           I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    @inbounds parent(A)[I...] = val
    A
end

end # module


text:
import Random: seed!
import LinearAlgebra: norm, dot

"""
    genprob(num_pts, u; sd=abs(rand(Int64)), K=5)

Generate a list of `num_pts + 1` functions to build function `min_{i =
1:num_pts+1} f_i`. The functions are defined inside the box `[0, u]` in
R^n. The function is created so that there are several (random)
quadratic functions with increasing minimum values. The `num_pts +
1`-th function is a plane that achieves the global solution at (0,
0). The values of the local minima are powers of `K`. To generate
always the same functions, one should provide the seed parameter `sd`.

Returns the list of functions and the bounds `l` and `u`.

"""
function genprob(num_pts, u; sd::Int64=abs(rand(Int64)), K=5)

    fmin_list = Vector{Function}(undef, num_pts + 1)

    n = length(u)

    l = zeros(Float64, n)

    sd = 38904589 + sd

    seed!(sd)

    K = 5.0

    for i = 1:num_pts

        fmin_list[i] = let

            local _n = n
            local k = K^i
            local d = 1000 * rand(n)
            local x0 = rand(n) .* u

            (x) -> k + 0.5 * sum((d[j] * (x[j] - x0[j])^2 for j = 1:_n))

        end

    end

    fmin_list[num_pts + 1] = let

        local c = rand(n)

        @. c = c / norm(c, 2)

        @views c[1] = (K^(num_pts + 2) - dot(c[2:end], u[2:end])) / u[1]

        (x) -> dot(c, x)

    end

    return fmin_list, l, u

end

"""
    gen_fmin(flist)

This function simply returns the `f_min` function, defined by `min
flist[i]`, which is usefull only for drawing purposes and intuition.

"""
function gen_fmin(flist)

    fmin = let

        fl = copy(flist)

        m = length(fl)

        (x) -> minimum((fl[i](x) for i = 1:m))

    end

    return fmin

end


text:
function cai(f,c)
    return x -> isnan(x) ? c : complex(f(x))
end

function zai(f)
    return x -> isnan(x) ? 0.0im : complex(f(x))
end

function horner(c::AbstractVector,kr::AbstractRange{Int},x)
    T = promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end
    ret = zero(T)
    @inbounds for k in reverse(kr)
        ret = muladd(x,ret,c[k])
    end
    ret
end

function horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractVector)
    n,T = length(x),promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zeros(T,n)
    @inbounds for k in reverse(kr)
        ck = c[k]
        @simd for i = 1:n
            ret[i] = muladd(x[i],ret[i],ck)
        end
    end
    ret
end

horner(c::AbstractVector,x) = horner(c,1:length(c),x)
horner(c::AbstractVector,x::AbstractArray) = horner(c,1:length(c),x)
horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractArray) = reshape(horner(c,kr,vec(x)),size(x))

function evaluate(f::AbstractVector,S::OscLaurent{D,R},z) where {D,R}
    #z = tocanonical(domain(S),z)
    zz = mappoint(domain(S),Circle(),z)
    invz = 1/zz
    (horner(f,1:2:length(f),zz) + horner(f,2:2:length(f),invz).*invz)*exp.(1im*z*S.exp)
end

function mobiusdiff(v::AbstractVector{T}) where T<:Number #need to fix
    n = length(v)
    if n == 1
        w = zeros(T,1)
    else
        w = Array{T}(undef, n+2)
        w[1] = zero(T)
        n = length(v)
        k = 1 #just for clarity
        @inbounds w[3] = 1.0im*k*v[3] + 0.5im*(k+1)*v[5]
        @inbounds w[2] = -1.0im*k*v[2] - 0.5im*(k+1)*v[4]
        l = (isodd(n) ? n÷2 : n÷2-1)
        for k=2:l-1
            @inbounds w[2k] = -1.0im*k*v[2k] - 0.5im*(k+1)*v[2k+2] - 0.5im*(k-1)*v[2k-2]
            @inbounds w[2k+1] = 1.0im*k*v[2k+1] + 0.5im*(k+1)*v[2k+3] + 0.5im*(k-1)*v[2k-1]
        end
        #### Deal with large n case
        if iseven(n) # if n is even, then l = n/2 -1, n = 2l + 2
            # positive modes
            @inbounds w[2l] = -1.0im*l*v[2l] - 0.5im*(l+1)*v[2l+2] - 0.5im*(l-1)*v[2l-2]
            @inbounds w[2l+2] = -1.0im*(l+1)*v[2l+2] - 0.5im*(l)*v[2l]
            @inbounds w[2l+4] = -0.5im*(l+1)*v[2l+2] # index -(n/2+2)
            # negative modes
            @inbounds w[2l+1] = 1.0im*l*v[2l+1] + 0.5im*(l-1)*v[2l-1]
            @inbounds w[2l+3] = 0.5im*(l)*v[2l+1] # index n/2
        else # if n is odd, then things are symmetric
            @inbounds w[2l+1] = 1.0im*l*v[2l+1] + 0.5im*(l-1)*v[2l-1]
            @inbounds w[2l+3] = 0.5im*(l)*v[2l+1]
            @inbounds w[2l] = -1.0im*l*v[2l] - 0.5im*(l-1)*v[2l-2]
            @inbounds w[2l+2] = -0.5im*(l)*v[2l]
        end
        w[1] = -fsum(w) # needed because fsum(w) = 0 is needed to enforce decay at infinity
        w
    end
end

function mob_conj(x::AbstractVector{T}) where T <: Number
    y = x
    if isodd(length(x))
        y[2:2:end] = conj(x[3:2:end])
        y[3:2:end] = conj(x[2:2:end])
    else  # This is not the exact conjugate!
        y[2:2:end-1] = conj(x[3:2:end-1])
        y[3:2:end-1] = conj(x[2:2:end-1])
    end

end


text:
# functions that allow physics to register their own options checking

"""
  Dictionary to hold the physics-specific options checking functions.

  Do not put functions into this dictionary directly, use the accessor
  function [`registerOptionsChecker`](@ref)
"""
global const PhysicsOptionsFuncs = Dict{String, Function}()

"""
  API function for registering an options checking function for a given
  physics.

  Note that users generally do not have to call this function directly,
  it is handled by [`register_physics`](@ref)

  This function is run after the non-physics specific default options
  have been supplied.  The function must have a a method with signature

  ```
    func(opts::Dict)
  ```

  which modifies the dictionary in-place.

  It is illegal to register more than one options checker per physics
  or a single options checker for more than one physics.

  **Inputs**

   * physics_name: the name defined by the physics module
   * func: the options checking function
"""
function registerOptionsChecker(physics_name::String, func::Function)

  if haskey(PhysicsOptionsFuncs, physics_name)
    error("An options checker has already been registered for physics $(physics_name)")
  end

  for (key, val) in PhysicsOptionsFuncs
    if val == func
      error("Options checker for physics $(physics_name) is already registered for physics $(key)")
    end
  end

  if length(methods(func, (Dict,))) == 0
    error("The options checker for physics $(physics_name) does not have a method matching func(opts::Dict)")
  end

  # function appears acceptable, add it to the list
  PhysicsOptionsFuncs[physics_name] = func

  return nothing
end

"""
  Prints the name of all currently registered physics modules and the module
  name itself
"""
function printPhysicsModules(f::IO=STDOUT)

  n = length(keys(PhysicsOptionsFuncs))
  println(f, "Currently registered physics modules: ", n)
  for (name, val) in PhysicsOptionsFuncs
    println(f, "  ", name)
  end

  return nothing
end


text:
# ADD TESTS SET FOR INVERSE
#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic
using IntervalContractors

#Preamble
setprecision(53)
setprecision(Interval, Float64)

# Using approximate checks for validaty update later?
import Base.isapprox
isapprox(x::Interval,y::Interval) = isapprox(x.lo,y.lo,atol=1E-4) && isapprox(x.hi,y.hi,atol=1E-4)

@testset "inv_rev_test" begin
    @test inv_rev(∅, entireinterval(Float64))[2] == ∅
    @test inv_rev(Interval(0.0, 1.0), ∅)[2] == ∅
    @test inv_rev(∅, Interval(0.0, 1.0))[2] == ∅
    @test inv_rev(Interval(-2.0, -1.0), entireinterval(Float64))[2] == Interval(-1.0, -0.5)
    @test inv_rev(Interval(1.0, 1.0), entireinterval(Float64))[2] == Interval(1.0, 1.0)
    @test inv_rev(entireinterval(Float64), entireinterval(Float64))[2] == Interval(-∞, ∞)
    @test inv_rev(Interval(-Inf, 0.0), entireinterval(Float64))[2] == Interval(-∞, 0.0)
end