text:
# Copyright (C) 2019, by Bill Burdick, ZLIB licensed, https://github.com/zot/domdom

"""
    Domdom

A simple server for [Domdom](https://github.com/zot/domdom).

Remember that there will be some "soak time" at the start because of Julia's JIT behavior, it should be quite zippy after that.

To start it, use the start function like this:

start(JSON, DIRECTORY, HANDLER)

JSON is a JSON-compatible Julia object (arrays, dicts, strings, numbers, booleans, etc.) as detailed in [Domdom](https://github.com/zot/domdom).

DIRECTORY the server uses for your HTML files, etc.

HANDLER is an event handler. You can use patternhandler() for easy event handling (see documentation). A "raw" event handler (patternhandler creates one of these) is a function(EVENT, CON, LOCATION, ARG)

- EVENT is :key, :click, or :set
- CON is the Domdom connection object
- LOCATION the path in the JSON document (contained in CON.document) for the event
- ARG is the JSON property for the event
"""
module Domdom

using HTTP, Sockets, JSON

export start, useverbose, document, web2julia, julia2web, verbose
export Connection, Dom, DomArray, DomObject, DomValue
export connection, path, webpath, contents, parent, root, clone, domproperties, domproperties!
export getpath, isdirty, dirty!, clean!, cleanall!

usingverbose = false
error = nothing
useverbose(flag) = global usingverbose = flag

"A Domdom connection"
mutable struct Connection
    document
    queuing
    queue
    index
    handler::Function
    socket::HTTP.WebSockets.WebSocket
    host
    port
    properties
    dirty
    domproperties
    Connection(document, handler, socket, host, port) = new(document, false, [], Dict(), handler, socket, host, port, Dict(), WeakKeyDict(), WeakKeyDict())
end

struct DomProperties
    properties
    DomProperties() = new(Dict())
end

domid = 0

"""
    DomArray

An array within a Domdom document. This functions like an array but it also contains the path
of the array. Domdom documents are trees, placing the same array in more than one place in the
document is not allowed.
"""
mutable struct DomArray <: AbstractArray{Any, 1}
    connection::Connection
    path::Array
    contents::Array{Any, 1}
    id
end

"""
    DomObject

An object within a Domdom document. This functions like a dictionary but it also contains the path
of the array. Domdom documents are trees, placing the same object in more than one place in the
document is not allowed.
"""
mutable struct DomObject <: AbstractDict{String, Any}
    connection::Connection
    path::Array
    contents::Dict{Symbol, Any}
    id
end

"Either a DomObject or a DomArray"
const Dom = Union{DomObject, DomArray}
const DomValue = Union{Dom, String, Symbol, Number, Bool, Char, Nothing}
"Either a DomObject, a DomArray, or a Connection"
const DomCon = Union{Connection, Dom}

DomArray(con::Dom, values::DomValue...) = DomArray(connection(con), [], collect(DomValue, values))
DomArray(con::Connection, path::Array, values::Array) = DomArray(con, path, collect(DomValue, values), [global domid += 1])
DomArray(con::Connection, path::Array, values::Array{DomValue, 1}) = DomArray(con, path, values, [global domid += 1])
DomObject(con::DomCon, type::Symbol; props...) = DomObject(connection(con), [], Dict{Symbol, DomValue}([:type=>type, props...]))
DomObject(con::DomCon; props...) = DomObject(connection(con), [], Dict{Symbol, DomValue}([props...]))
DomObject(con::Connection, path::Array, contents::Dict{Symbol, DomValue}) = DomObject(con, path, contents, [global domid += 1])

id(el::Dom) = getfield(el, :id)
connection(el::Dom) = getfield(el, :connection)
connection(con::Connection) = con
webpath(el::Dom) = getfield(el, :path)
webpath!(el::Dom, path) = setfield!(el, :path, path)
path(el::Dom) = web2julia(webpath(el))
Base.parent(el::Dom) = getpath(connection(el), path(el)[1:end - 1])
function Base.replace(el::Dom, el2::Dom)
    verbose("REPLACING:\n  $(el)\nwith\n  $(el2)")
    con = connection(el)
    path = absolutepath(con, path(el))
    setpath(con, path, el2)
    queue(con, [:set, julia2web(path), juliaobj2webobj(value)])
end
clone(x) = x
clone(el::DomArray) = DomArray(connection(el), [], Any[clone(x) for x in el.contents])
clone(el::DomObject) = DomObject(connection(el), [], Dict(Any[k => clone(v) for (k, v) in contents(el)]))
contents(el::Dom) = getfield(el, :contents)
"Root of the document"
root(con::Connection) = con.document.contents
root(el::Dom) = root(connection(el))
domproperties(el::Dom) = get!(connection(el).domproperties, id(el), DomProperties())
Base.show(io::IO, el::DomArray) = print(io, "DomArray(",join(map(stringFor, path(el)), ", "),")[", join(map(repr, el), ", "), "]")
Base.show(io::IO, el::DomObject) = print(io, "DomObject(",join(map(stringFor, path(el)), ", "),")[", join(map(p->"$(String(p[1]))=>$(repr(p[2]))", collect(contents(el))), ", "), "]")
Base.getproperty(el::DomObject, name::Symbol) = el[name]
Base.getproperty(props::DomProperties, name::Symbol) = get(getfield(props, :properties), name, nothing)
Base.setproperty!(el::DomObject, name::Symbol, value::DomValue) = el[name] = value
Base.setproperty!(props::DomProperties, name::Symbol, value) = getfield(props, :properties)[name] = value
Base.getindex(el::DomObject, name) = get(contents(el), name, nothing)
Base.getindex(el::DomArray, index) = if index in 1:length(contents(el)) contents(el)[index] else nothing end
function Base.setindex!(el::DomObject, value::DomValue, name)
    subsetindex(el, value, get(contents(el), name, nothing), name)
    contents(el)[name] = value
end
function Base.setindex!(el::DomArray, value::DomValue, index::Integer)
    if length(el) + 1 == index
        push!(contents(el), value)
        subsetindex(el, value, nothing, index)
    else
        old = contents(el)[index]
        contents(el)[index] = value
        subsetindex(el, value, old, index)
    end
end
function subsetindex(el::Dom, value, oldValue, index)
    # verify that el is actually connected
    if connection(el) != nothing && getpath(connection(el), path(el)) === el
        newpath = checkpath([webpath(el)..., julia2web(index)], value)
        adjustindex(connection(el).index, newpath, oldValue, value)
        queue(connection(el), [:set, newpath, juliaobj2webobj(value)])
    end
end
Base.length(el::Dom) = length(contents(el))
Base.iterate(el::Dom, state...) = iterate(contents(el), state...)
Base.size(el::DomArray) = size(el.contents)
Base.IndexStyle(::Type{DomArray}) = IndexLinear()
Base.pairs(el::DomObject) = pairs(contents(el))
Base.keys(el::DomObject) = keys(contents(el))
Base.values(el::DomObject) = values(contents(el))
function Base.push!(el::DomArray, items::DomValue...)
    for item in items
        el[end + 1] = item
    end
end
function Base.pop!(dom::DomArray)
    adjustindex(connection(dom).index, [webpath(dom)..., length(dom) - 1], dom[end], nothing)
    pop!(dom.contents)
    #queue(connection(dom), [:deleteLast, webpath(dom)])
    queue(connection(dom), [:splice, webpath(dom), length(dom), 1]) # don't subtract one from length
end

Base.deleteat!(dom::DomArray, i::Integer) = deleteat!(dom, i:i)
function Base.deleteat!(dom::DomArray, r::UnitRange{<:Integer})
    for i in r
        adjustindex(connection(dom).index, [webpath(dom)..., i - 1], dom[i], nothing)
    end
    deleteat!(dom.contents, r)
    parentpath = path(dom)[1:end - 1]
    for i in r[1]:length(dom)
        changepath([parentpath..., i], dom[i])
    end
    queue(connection(dom), [:splice, webpath(dom), r[1] - 1, length(r)])
end

function changepath(elpath, value::Dom)
    webpath!(value, julia2web(elpath))
    for (key, v) in pairs(contents(value))
        changepath([elpath..., key], v)
    end
end
changepath(elpath, value) = nothing

function checkpath(elpath, value::Dom, set = false)
    if webpath(value) == [] || set
        for (k, v) in pairs(contents(value))
            if isa(v, Dom)
                checkpath([elpath..., julia2web(k)], v, true)
            end
        end
        setfield!(value, :path, elpath)
    elseif webpath(value) != elpath
        throw(ArgumentError("Attempt to move a tree value from [$(join(webpath(value), ", "))] to [$(join(elpath, ", "))]"))
    end
    elpath
end
checkpath(path, value) = path

function absolutepath(con::Connection, path)
    if isa(path, Symbol)
        con.index[path]
    elseif length(path) > 0 && isa(path[1], String) && startswith(path[1], "@")
        [con.index[Symbol(path[1][2:end])]..., path[2:end]...]
    else
        path
    end
end

"""
    getpath(CON::Connection, PATH) -> Any

Return the value at PATH in CON

CON is a domdom connection

PATH is a path in the conneciton
"""
getpath(con::Connection, path::AbstractArray{T, 1} where T) = walk(con, path)

walk(con::Connection, path::AbstractArray{T, 1} where T) = walk(con, con.document.contents, path)
function walk(con::Connection, obj::Dom, path::AbstractArray{T, 1} where T)
    local key, curindex

    path = absolutepath(con, path)
    try
        for index in eachindex(path)
            curindex = index
            key = path[index]
            obj = obj[key]
        end
    catch err
        if isa(err, ArgumentError)
            println(stderr, "ERROR: Could not find index $key at $curindex in path $(join([path[1:curindex]], ", "))")
            throw(err)
        end
    end
    obj
end

function setpath(con::Connection, path, value)
    path = absolutepath(con, path)
    adjustindex(con.index, path, getpath(con, path), value)
    local parent = getpath(con, path[1:end - 1])
    if isa(path[end], Integer) # append item if it's an array of len-1
        if length(parent) >= path[end]
            contents(parent)[path[end]] = value
        elseif length(parent) == path[end] - 1
            push!(parent, value)
        else
            throw(BoundsError(parent, path[end]))
        end
    else
        contents(parent)[path[end]] = value
    end
    dirty!(parent)
end

function adjustindex(index, path, oldjson, newjson)
    local oldids = findIds(path, oldjson)
    local newids = findIds(path, newjson)
    local oldkeys = Set(keys(oldids))

    for (k, v) in newids
        index[k] = v
        delete!(oldkeys, k)
    end
    for k in oldkeys
        delete!(index, k)
    end
end

function findIds(path, json, ids = Dict())
    if isa(json, Array)
        for k in 1:length(json)
            findIds([path..., k - 1], json[k], ids)
        end
    elseif isa(json, Dict)
        if haskey(json, :id)
            ids[Symbol(json[:id])] = path
        end
        for (k, v) in json
            findIds([path..., k], v, ids)
        end
    end
    ids
end

function insertpath(con::Connection, path, value)
    path = absolutepath(con, path)
    local endkey = path[end]
    local parent = getpath(con, path[1:end - 1])

    adjustindex(con.index, path, nothing, value)
    if isa(parent, Array)
        splice!(parent, endkey:0, [value])
    elseif isa(parent, Dict)
        parent[endkey] = value
    end
    #queue(con, [:insert, julia2web(path), juliaobj2webobj(value)])
    queue(con, [:splice, julia2web(parent), julia2web(endkey), 0, juliaobj2webobj(value)])
end

"""
    isdirty(dom)

Determine if an object is dirty
"""
isdirty(dom::Dom) = haskey(connection(dom).dirty, id(dom))
"""
    clean!(dom)

State that an object is clean
"""
clean!(dom::Dom) = delete!(connection(dom).dirty, id(dom))
"""
    cleanall(dom)
    cleanall(con)

State that all objects in a connection are clean
"""
cleanall!(dom::Dom) = cleanall!(connection(dom))
cleanall!(con::Connection) = con.dirty = WeakKeyDict()
"""
    dirty!(dom)
    dirty!(con, path)

Make a field dirty
"""
dirty!(con::Connection, path) = dirty!(getpath(con, path))
dirty!(dom::Dom) = connection(dom).dirty[id(dom)] = true

"print a message if usingverbose is true"
function verbose(args...)
    if usingverbose
        println(args...)
    end
end

"low-level event handler invoked by a socket message"
function handleset(con, location, value)
    verbose("SET $(repr(location)) = $(repr(value))")
    location = web2julia(location)
    adjustindex(con.index, location, getpath(con, location), value)
    setpath(con, location, value)
    con.handler(:set, con, location, value)
end

const cmds = Dict([
    :click => (con::Connection, name, location)-> con.handler(:click, con, web2julia(location), name)
    :set => handleset
    :key => (con::Connection, name, location)-> con.handler(:key, con, web2julia(location), name)
])

const nullhandler(args...) = nothing

const emptyDict = Dict()

web2julia(path::Number) = path + 1
web2julia(path::String) = Symbol(path)
web2julia(path::AbstractArray) = web2julia.(path)

julia2web(path::Number) = path - 1
julia2web(path::Symbol) = String(path)
julia2web(path::AbstractArray) = julia2web.(path)

stringFor(item) = repr(item)
stringFor(item::Symbol) = String(item)

juliaobj2webobj(x) = x
juliaobj2webobj(array::DomArray) = juliaobj2webobj.(array.contents)
juliaobj2webobj(obj::DomObject) = Dict(map(p->(p[1], juliaobj2webobj(p[2])), collect(contents(obj))))

call(con::Connection, cmd, args...) = cmds[Symbol(cmd)](con, args...)

function queue(con::Connection, item)
    if con.queuing > 0
        verbose("QUEUING ", JSON.json(item))
        push!(con.queue, item)
    else
        verbose("SENDING ", JSON.json(item))
        send(con, item)
    end
end

startqueue(con) = con.queuing += 1

function flushqueue(con::Connection)
    if con.queuing > 0
        con.queuing -= 1
        if con.queuing == 0 # send
            if length(con.queue) == 1
                send(con, con.queue[1])
            elseif length(con.queue) > 1
                send(con, [:batch, con.queue])
            end
        end
        empty!(con.queue)
    end
end

send(con::Connection, item) = write(con.socket, JSON.json(item))

function document(con::Connection, dom = con.document)
    con.document = dom
    local ids = findIds([], con.document[:contents])

    checkpath([], con.document[:contents])
    verbose("set document to ", repr(dom))
    if !isempty(ids)
        push!(con.index, ids...)
    end
    queue(con, ["document", juliaobj2webobj(con.document)])
end

# open_file copied from Gadfly.jl, MIT license: https://github.com/GiovineItalia/Gadfly.jl
#
function open_file(filename)
    if Sys.isapple()
        run(`open $(filename)`)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $(filename)`)
    elseif Sys.iswindows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    else
        @warn "Opening browseers is not supported on OS $(string(Sys.KERNEL))"
    end
end

"""
    start(DIR::String, PORT = 8085; [browse = false]) do CONNECTION, EVENT
        event.onset(TYPE, PROPERTY) do dom, key, arg, obj, event
            ...
        end
        event.onclick(TYPE, PROPERTY) do dom, key, arg, obj, event
            ...
        end
        event.onchanged(TYPE) do dom, key, arg, obj, event
            ...
        end
        ...
    end

Start a Domdom server, given a function that configures the events and the new connection

The provided function must return the top-level document
"""
function start(eventfunc, dir::String, port = 8085; browse = false, config = (args...)->())
    set = Dict()
    changed = Dict()
    click = Dict()
    start(dir, patternhandler(Dict([:set=>set :changed=>changed :click=>click])), port, browse=browse, config=config) do con
        eventfunc(con, (
            onset = (hnd, type, attr)-> set[(type, attr)] = hnd,
            onchanged = (hnd, type)-> changed[type] = hnd,
            onclick = (hnd, type, attr)-> click[(type, attr)] = hnd,
        ))
    end
end
function start(conFunc, dir::String, handler::Function, port = 8085; browse = false, config = (args...)->())
    verbose("STARTING HTTP ON PORT $port, DIR $dir")
    #host = Sockets.localhost
    local host = ip"0.0.0.0"

    if browse
        @async begin
            #sleep(3)
            open_file("http://localhost:8085")
        end
    end
    config(dir, host, port)
    HTTP.listen(host, port) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do ws
                try
                    connection = Connection(Dict(), handler, ws, host == ip"0.0.0.0" ? Sockets.localhost : host, port)
                    startqueue(connection)
                    document(connection, conFunc(connection))
                    flushqueue(connection)
                    while !eof(ws)
                        frame = HTTP.WebSockets.readframe(ws)
                        startqueue(connection)
                        verbose("RECEIVED: ", String(frame))
                        call(connection, JSON.Parser.parse(String(frame))...)
                        flushqueue(connection)
                    end
                catch err
                    global error
                    error = err
                    if isa(err, HTTP.WebSockets.WebSocketError)
                        println("WebSocket Error $(repr(err.status))")
                    else
                        #showerror(stderr, err, catch_backtrace())
                        Base.display_error(err, catch_backtrace())
                    end
                end
            end
        else
            local req = http.message
            local file = dir * replace(http.message.target, r"\?.*$" => "")

            req.body = read(http)
            closeread(http)
            if isdir(file) file = file * "/index.html" end
            body = read(file)
            req.response = HTTP.Response(body)
            HTTP.setheader(req.response, "Content-Type" => sniff(file, body))
            HTTP.setheader(req.response, "Access-Control-Allow-Origin" => "*")
            req.response.request = req
            startwrite(http)
            write(http, req.response.body)
        end
    end
    (dir, host, port)
end

filetypes = Dict([
    "css" => "text/css"
])

function sniff(file, body)
    local m = match(r"\.(.*)$", file)

    if m == nothing
        HTTP.sniff(body)
    else
        get(filetypes, m.captures[1]) do
            HTTP.sniff(body)
        end
    end
end

"""
    patternhandler(DICT; clickhandler = nullhandler, sethandler = nullhandler, keyhandler = nullhandler)

DICT is a dictionary with optional :set, :click, and :key entries.

Each entry is (TYPE, FIELD)=> function(DOMITEM, KEY, ARG, OBJ, EVENT)

- TYPE is the type of the event's JSON object
- FIELD is the field in the event's JSON object that is being changed or clicked
- DOMITEM is the item in the document
- KEY is the field name
- OBJ is the JSON object
- EVENT is the
"""
function patternhandler(dict; clickhandler = nullhandler, sethandler = nullhandler, keyhandler = nullhandler)
    function (event, con, location, arg)
        local obj = event == :key ? con.document : getpath(con, event == :click ? location : location[1:end - 1])
        local objtype = obj[:type]
        local key = if event in [:click :key]
            location = [location..., arg]
            isa(arg, String) ? Symbol(arg) : arg
        elseif event == :set
            location[end]
        end
        local handlers = get(dict, event, emptyDict)
        local defaulthandler = event == :click ? clickhandler : event == :set ? sethandler : keyhandler

        verbose("EVENT $event LOCATION $(repr(location)) TYPE $(typeof(objtype)) $objtype KEY $(typeof(key)) $(repr(key)) ARG $(repr(arg)) VALUE $(repr(getpath(con, location)))")
        get(handlers, (objtype, key), get(handlers, objtype, defaulthandler))(getpath(con, location[1:end - 1]), key, arg, obj, event)
        if event == :set
            changedhandlers = get(dict, :changed, emptyDict)
            get(changedhandlers, objtype, nullhandler)(getpath(con, location[1:end - 1]), key, arg, obj, event)
        end
    end
end

end #module

max_stars_repo_name: zot/domdom

max_issues_repo_name: zot/domdom

max_forks_repo_name: zot/domdom


text:
using ModelingToolkit, ParameterizedFunctions

### Capture MATLABDiffEq.jl type issues
@parameters t a b c d
@variables x(t) y(t)
D = Differential(t)
eqs = [D(x) ~ a*x - b*x*y
       D(y) ~ -c*y + d*x*y]
sys = ODESystem(eqs)
equations(sys)
matstr = Symbolics.build_function(map(x->x.rhs,equations(sys)),states(sys),
                                        parameters(sys),independent_variable(sys),
                                        target = ModelingToolkit.MATLABTarget())
@test matstr == "diffeqf = @(t,internal_var___u) [
  internal_var___p(1) * internal_var___u(1) + -1 * internal_var___p(2) * internal_var___u(1) * internal_var___u(2);
  -1 * internal_var___p(3) * internal_var___u(2) + internal_var___p(4) * internal_var___u(1) * internal_var___u(2);
];
"

f = @ode_def_bare LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
p = [1.5,1,3,1]
tspan = (0.0,10.0)
u0 = [1.0,1.0]
prob = ODEProblem(f,u0,tspan,p)

sys = modelingtoolkitize(prob)
[eq.rhs for eq ∈ equations(sys)]
matstr = Symbolics.build_function(map(x->x.rhs,equations(sys)),states(sys),
                                        parameters(sys),independent_variable(sys),
                                        target = ModelingToolkit.MATLABTarget())
@test matstr == "diffeqf = @(t,internal_var___u) [
  internal_var___p(1) * internal_var___u(1) - internal_var___p(2) * internal_var___u(1) * internal_var___u(2);
  internal_var___u(2) * -(internal_var___p(3)) + internal_var___p(4) * internal_var___u(1) * internal_var___u(2);
];
"

ext: jl


text:
using InteractionNetworkModels, Distributions, BenchmarkTools, Plots

# The Model(s)
model_mode = Hollywood(-3.0, Poisson(7), 10)
S = sample(model_mode, 10)
S = [[1,1,1,1], [2,2,2,2], [3,3,3,3]]
V = 1:20

d_lsp = MatchingDist(FastLSP(100))
# d_f = FastMatchingDist(FastLCS(100), 51)

K_inner = DimensionRange(2, 50)
K_outer = DimensionRange(1, 50)

model = SIM(S, 3.2, d_lsp, V, K_inner, K_outer)

# model_f = SIM(S, 4.0, d_f, V, 50, 50)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=5, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.1, model.K_inner.l, model.K_inner.u),
    lag=1,
    K=200
)
mcmc_sampler_sp = SimMcmcInsertDeleteSubpath(
    ν_ed=4, β=0.6, ν_td=3,
    len_dist=TrGeometric(0.1, model.K_inner.l, model.K_inner.u),
    lag=1,
    K=200
)

@time out=mcmc_sampler(
    model,
    lag=1,
    init=model.mode,
    desired_samples=1000,
    burn_in=0
)
plot(out)

@time out_sp=mcmc_sampler_sp(
    model,
    lag=1,
    init=model.mode,
    desired_samples=1000,
    burn_in=0
)
plot!(out_sp)

summaryplot(out)
summaryplot!(out_sp)


text:
# All functions specific to the Zosterops experiments.
# Initialise predefined species and redefine dispersal and reproduction

# IMPORTANT NOTES (when using `zosterops` mode)
#
# 1. the second environmental niche is now `above-ground carbon` instead of
#    `precipitation` (unfortunately, we can't actually rename it everywhere)
# 2. the `cellsize` setting now determines the patch carrying capacity in
#    individuals, not grams
# 3. the `fertility` setting is the absolute number of offspring per breeding pair
#    (instead of a metabolic coefficient)
# 4. the `tolerance` setting now determines the probability that a mate of
#    another species is accepted
# 5. `degpleiotropy` must be set to 0, otherwise the species initialisation will fail

let zosterops = Individual[]
    """
        initzosteropsspecies(settings)

    Initialise the predefined Zosterops species archetypes (silvanus/highland, jubaensis/lowland).
    """
    function initzosteropsspecies(settings::Dict{String, Any})
        # ensure that "one gene, one trait" is true and that we have species definitions
        (settings["degpleiotropy"] != 0) && simlog("degpleiotropy must be 0", 'e')
        (isnothing(settings["species"])) && simlog("no species defined", 'e')
        # load per-species AGC optimum and tolerance values from settings
        for (s, v) in pairs(settings["species"])
            push!(zosterops, initzosteropsspecies(s, Float64(v[1]), Float64(v[2]), settings))
        end
    end

    """
        initzosteropsspecies(name, precopt, prectol, settings)

    Create a new individual then modify its traits to create a Zosterops archetype.
    """
    function initzosteropsspecies(name::String, precopt::Float64, prectol::Float64, settings::Dict{String, Any})
        archetype = createind(settings)
        archetype.lineage = name
        # Find the gene that codes for the relevant trait and change that trait's value
        for chromosome in archetype.genome
            for gene in chromosome.genes
                isempty(gene.codes) && continue
                genetrait = settings["traitnames"][gene.codes[1].nameindex]
                #XXX do we also need to set tempopt/temptol?
                if genetrait == "precopt"
                    gene.codes[1].value = precopt
                elseif genetrait == "prectol"
                    gene.codes[1].value = prectol
                end
            end
        end
        # then recalculate the individual's trait dict and return the finished archetype
        archetype.traits = gettraitdict(archetype.genome, settings["traitnames"])
        return archetype
    end

    """
        getzosteropsspecies(name, sex, settings)

    Return a new individual of the named species
    """
    global function getzosteropsspecies(name::String, sex::Sex, settings::Dict{String, Any})
        # This function has to be global to escape the let-block of the species list
        isempty(zosterops) && initzosteropsspecies(settings)
        bird = nothing
        for z in zosterops
            if z.lineage == name
                bird = deepcopy(z)
                break
            end
        end
        (isnothing(bird)) && simlog("Unknown species name: $name", 'e')
        bird.id = rand(UInt32)
        varyalleles!(bird.genome, settings, rand()) #XXX do we want mutations?
        bird.traits = gettraitdict(bird.genome, settings["traitnames"])
        bird.size = bird.traits["repsize"] #if we want mutations, we need to fix the size afterwards
        bird.sex = sex
        return bird
    end

    """
        getzosteropsnames(settings)

    Return the names of all defined Zosterops species.
    """
    #XXX utility function - but do I need it?
    global function getzosteropsnames(settings::Dict{String, Any})
        isempty(zosterops) && initzosteropsspecies(settings)
        return map(s->s.lineage, zosterops)
    end
end

"""
    zgenesis(patch, settings)

Create a new community of Zosterops breeding pairs (possibly of multiple species).
Returns an array of individuals.
"""
function zgenesis(patch::Patch, settings::Dict{String, Any})
    community = Array{Individual, 1}()
    (isnothing(settings["species"])) && simlog("no species defined", 'e')
    # check which species can inhabit this patch
    species = Array{String, 1}()
    for s in settings["species"]
        # s[1] = species name, s[2][1] = AGC opt, s[2][2] = AGC tol
        if abs(s[2][1]-patch.prec) <= s[2][2]
            push!(species, s[1])
        end
    end
    (isempty(species)) && return community
    # calculate the number of initial breeding pairs and add a male and a female for each
    npairs = Integer(rand(0:round(settings["cellsize"]/2)))
    for i in 1:npairs
        sp = rand(species)
        m = getzosteropsspecies(sp, male, settings)
        f = getzosteropsspecies(sp, female, settings)
        f.partner = m.id
        m.partner = f.id
        push!(community, m)
        push!(community, f)
        simlog("Adding a pair of Z. $sp", 'd')
    end
    community
end

"""
    zdisperse!(world, settings)

Dispersal of bird individuals within the world. Males disperse first, looking
for suitable habitats within their dispersal range to establish territories.
Females disperse second, looking for available mates. (Cf. Aben et al. 2016)
"""
function zdisperse!(world::Array{Patch,1}, settings::Dict{String, Any}, sex::Sex=male)
    for patch in world
        newseedbank = Array{Individual,1}()
        while length(patch.seedbank) > 0
            #XXX This is probably rather inefficient (creating a new list each time)
            juvenile = pop!(patch.seedbank)
            if juvenile.sex == sex
                zdisperse!(juvenile, world, patch.location,
                           Integer(settings["cellsize"]), settings["tolerance"])
            else
                push!(newseedbank, juvenile)
            end
        end
        patch.seedbank = newseedbank
    end
    (sex == male) && zdisperse!(world, settings, female)
end

"""
    zdisperse!(bird, world, location, cellsize)

Dispersal of a single bird.
"""
function zdisperse!(bird::Individual, world::Array{Patch,1}, location::Tuple{Int, Int},
                    cellsize::Int, tolerance::Float64)
    # keep track of where we've been and calculate the max dispersal distance
    x, y = location
    route = [location]
    #XXX disperse!() adds `/sqrt(2)`??
    #FIXME we're probably going to need another distribution
    maxdist = rand(Logistic(bird.traits["dispmean"], bird.traits["dispshape"]))
    while maxdist > 0
        # calculate the best habitat patch in the surroundings (i.e. the closest to AGC optimum)
        target = [(x-1, y-1), (x, y-1), (x+1, y-1),
                  (x-1, y),             (x+1, y),
                  (x-1, y+1), (x, y+1), (x+1, y+1)]
        filter!(c -> !in(c, route), target)
        possdest = map(p -> coordinate(p[1], p[2], world), target)
        filter!(p->!isnothing(p), possdest)
        if iszero(length(possdest))
            simlog("A Z.$(bird.lineage) died after failed dispersal.", 'd')
            return
        end
        bestdest = possdest[1]
        bestfit = abs(bestdest.prec - bird.traits["precopt"])
        for patch in possdest[2:end]
            patchfit = abs(patch.prec - bird.traits["precopt"])
            if patchfit < bestfit
                bestdest = patch
                bestfit = patchfit
            end
        end
        # check if the patch is full, within the bird's AGC range, and (for females) has a mate
        if bird.sex == female #XXX not ideal to do this here
            partner = findfirst(b -> ziscompatible(bird, b, tolerance), bestdest.community)
        end
        if (bestfit <= bird.traits["prectol"] &&
            length(bestdest.community) < cellsize &&
            (bird.sex == male || !isnothing(partner)))
            if bird.sex == female
                partner = bestdest.community[partner]
                bird.partner = partner.id
                partner.partner = bird.id
            end
            bird.marked = true
            push!(bestdest.community, bird)
            simlog("A Z.$(bird.lineage) moved to $(bestdest.location[1])/$(bestdest.location[2]).", 'd')
            return #if we've found a spot, we're done
        end
        x, y = bestdest.location
        push!(route, bestdest.location)
        maxdist -= 1
    end #if the max dispersal distance is reached, the individual simply dies
    simlog("A Z.$(bird.lineage) died after failed dispersal.", 'd')
end

"""
    ziscompatible(female, male, tolerance)

Check to see whether two birds are reproductively compatible.
"""
function ziscompatible(f::Individual, m::Individual, tolerance::Float64)
    !(m.sex == male && f.sex == female) && return false
    !(m.partner == 0 && f.partner == 0) && return false
    #FIXME When birds are looking for a mate, they're still juvenile, but they only mate as adults
    #!(m.size >= m.traits["repsize"] && f.size >= f.traits["repsize"]) && return false
    (m.lineage != f.lineage && rand(Float64) > tolerance) && return false
    #TODO how about seqsimilarity and genetic compatibility?
    simlog("Found a partner: $(f.id) and $(m.id).", 'd')
    return true
end

"""
    zreproduce!(patch, settings)

Reproduction of Zosterops breeding pairs in a patch.
"""
function zreproduce!(patch::Patch, settings::Dict{String, Any})
    #FIXME reproduction happens too seldom?
    noffs = Integer(settings["fertility"])
    for bird in patch.community
        if bird.sex == female && bird.partner != 0
            pt = findfirst(b -> b.id == bird.partner, patch.community)
            if isnothing(pt)
                simlog("Z.$(bird.lineage) no. $(bird.id) has no partner.", 'd')
                continue
            end
            partner = patch.community[pt]
            simlog("A Z.$(bird.lineage) mated with a Z.$(partner.lineage), $noffs offspring.", 'd')
            #TODO Offspring are assigned the lineage of their mother. Is that what we want?
            append!(patch.seedbank, createoffspring(noffs, bird, partner,
                                                    settings["traitnames"], true))
        end
    end
end

let width = 0, height = 0
    """
        coordinate(x, y, world)

    A utility function to perform a fast look-up for the patch at coordinate x/y.
    Important: this assumes a rectangular world with coordinates in row-major order!
    """
    global function coordinate(x::Int, y::Int, world::Array{Patch,1})
        if iszero(width)
            width = maximum(p -> p.location[1], world)
            height = maximum(p -> p.location[2], world)
        end
        (x <= 0 || y <= 0 || x > width || y > height) && return
        i = ((y-1) * width) + x
        return world[i]
    end
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
module PopGen

##   o O       o O       o ############ O       o O       o O
## o | | O   o | | O   o |              | O   o | | O   o | | O
## | | | | O | | | | O | | Dependencies | | O | | | | O | | | | O
## O | | o   O | | o   O |              | o   O | | o   O | | o
##   O o       O o       O ############ o       O o       O o

using DataFrames, PlotlyJS, GeneticVariation, Distributions, MultipleTesting

export PopObj,
    summary,
    nancycats,
    gulfsharks,
    csv,
    genepop,
    vcf,
    sample_names,
    isolate_genotypes,
    locations, locations!,
    population, populations, population!, populations!,
    remove_inds!,
    remove_loci!,
    missing,
    heterozygosity, het, He,
    hwe_test,
    plot_missing,
    plot_locations

##   o O       o O       o ############ O       o O       o O
## o | | O   o | | O   o |              | O   o | | O   o | | O
## | | | | O | | | | O | |  Load files  | | O | | | | O | | | | O
## O | | o   O | | o   O |              | o   O | | o   O | | o
##   O o       O o       O ############ o       O o       O o

include("PopObj.jl")
include("Read.jl")
include("Manipulate.jl")
include("Plotting.jl")
include("Datasets.jl")
include("AlleleFreq.jl")
include("HardyWeinberg.jl")

end # module PopGen


text:
#!/usr/bin/env julia

using Luxor

# some of these tests produce errors on Linux versions of Cairo
# but no errors on MacOS X Cairo.

function matrix_tests(fname)
    # matrix tests

    # translate(dx,dy) =	transform([1,  0, 0,  1, dx, dy])                shift by
    # scale(fx, fy)    =    transform([fx, 0, 0, fy,  0, 0])                 scale by
    # rotate(A)        =    transform([c, s, -c, c,   0, 0])                 rotate to A radians
    # x-skew(a)        =    transform([1,  0, tan(a), 1,   0, 0])            xskew by A
    # y-skew(a)        =    transform([1, tan(a), 0, 1, 0, 0])               yskew by A
    # flip HV          =    transform([fx, 0, 0, fy, cx(*1-fx), cy* (fy-1)]) flip

    Drawing(1000,1000, fname)

    # absolute position 200,250 relative to top left origin (0, 0)
    setmatrix([1, 0, 0, 1, 200, 250.0])
    sethue("red")
    fontsize(20)
    text("hello world")

    sethue("green")
    # absolute position 400,500 from top left origin (0, 0)
    setmatrix([1, 0, 0, 1, 400, 500.0])
    text("hello world")

    sethue("white")
    # absolute position 230 230
    setmatrix([1, 0, 0, 1, 230, 230.0])
    text("hello world")

    # shift relative by 10 x (right) and 20 y (down) after moving to absolute 300 300
    setmatrix([1, 0, 0, 1, 300, 300])
    transform([1, 0, 0, 1, 10, 20])
    sethue("yellow")
    text("hello world")

    setmatrix([1, 0, 0, 1, 300, 300])
    transform([1, 0, 0, 1, 0, 0])
    sethue("gray")
    text("hello world")

    # scale by 2 x (right) and 3 y after moving to 300 300
    setmatrix([1, 0, 0, 1, 300, 300.0])
    transform([2, 0, 0, 3, 0, 0.0])
    sethue("purple")
    text("hello world")

    # using rotate by angle, not using matrix
    # absolute move to 100, 100
    setmatrix([1, 0, 0, 1, 100.0, 100])
    fontsize(20)
    sethue("brown")
    for i in 1:4
        # rotate by pi/2 radians
        rotate(pi/2)
        text("hello world")
    end

    # rotate by angle, around current position
    # absolute move to 500 across 500 down
    setmatrix([1, 0, 0, 1, 300.0, 600])
    fontsize(20)
    sethue("cyan")
    n = 10
    angle = 2pi/n
    for i in 1:n
        # rotate another 'angle' radians
        transform([cos(angle), sin(angle), -sin(angle), cos(angle), 0, 0])
        text("$i hello world")
    end

    fontsize(20)
    # skew X
    setopacity(0.5)
    sethue("purple")
    # move to 800 across 300 down
    setmatrix([1, 0, 0, 1, 800, 300])
    # Cairo doesn't like tan(90) skews (or the matrix that results)
    for i in 0:5:40
        transform([1, 0, tand(i), 1, 0, 20])    # xskew by i° and move down by 20 for each one
        text("hello world")
    end

    # skew y
    setopacity(0.75)
    sethue("orange")
    setmatrix([1, 0, 0, 1, 800, 500])
    # Cairo doesn't like tan(90) skews
    for i in 0:5:40
        transform([1, tand(i), 0, 1, 0, 20])    # yskew by i° and move down by 20 for each one
        text("hello world")
    end

    # flip hv
    # Flip H/V with center shifted by cx:cy

    sethue("magenta")
    setmatrix([1, 0, 0, 1, 800, 200.0])
    fx = 1
    fy = -1
    cx = cy = 50
    transform([fx, 0, 0, fy, cx * (1-fx), cy * (fy-1)])
    text("hello world")

    finish()
    println("finished test: output in $(fname)")
end

matrix_tests("/tmp/matrix-tests.pdf")


text:
using BinDeps

@BinDeps.setup

libta_lib = library_dependency("libta_lib", aliases=["libta_lib", "libta_lib.so"])

if is_apple()
    #const libta_lib = "/usr/local/lib/libta_lib.0.0.0.dylib"

    if Pkg.installed("Homebrew") === nothing
        error("Homebrew package not installed, please run Pkg.add(\"Homebrew\")")
    end
    using Homebrew
    provides(Homebrew.HB, "ta-lib", libta_lib, os = :Darwin)

elseif is_linux()
    #const libta_lib = "/usr/lib/libta_lib.so.0.0.0"

    provides(Sources, URI("http://prdownloads.sourceforge.net/ta-lib/ta-lib-0.4.0-src.tar.gz"), libta_lib, unpacked_dir="ta-lib")
    provides(BuildProcess, Autotools(libtarget = "libta_lib", configure_options=["--prefix=/usr"]), libta_lib, os = :Unix)

else
    error("TALib.jl doesn't support this OS")
end

@BinDeps.install Dict(:libta_lib => :libta_lib)


text:
using XRayQuadrics
using Documenter

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

makedocs(;
    modules=[XRayQuadrics],
    authors="Thanasi Pantazides <thanasipantazides@gmail.com> and contributors",
    repo="https://github.com/thanasipantazides/XRayQuadrics.jl/blob/{commit}{path}#{line}",
    sitename="XRayQuadrics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Theory" => [
            "Introduction" => "introduction.md",
            "Quadric Surfaces" => "quadrics.md",
            "Intersections and Collisions" => "collisions.md"
        ]
    ],
)


text:
function sma(ta::TimeArray, n::Integer)
  tstamps = timestamp(ta)[n:end]

  vals = zeros(size(values(ta),1) - (n-1), size(values(ta),2))
  for i in 1:size(values(ta),1) - (n-1)
    for j in 1:size(values(ta),2)
      vals[i,j] = nanmean(values(ta)[i:i+(n-1),j])[1]
    end
  end

  cname = Symbol[]
  cols  = colnames(ta)
  for c in 1:length(cols)
    push!(cname, Symbol(string(cols[c], "_sma_", n)))
  end

  TimeArray(tstamps, vals, cname, meta(ta))
end

function ema(ta::TimeArray, n::Integer; wilder::Bool = false)
  k = (1 + !wilder) / (n + !wilder)

  tstamps = timestamp(ta)[n:end]

  vals    =  zeros(size(values(ta),1), size(values(ta),2))
  # seed with first value with an sma value
  vals[n,:] = values(sma(ta, n))[1,:]

  for i in n+1:size(values(ta),1)
    for j in 1:size(values(ta),2)
      _v = values(ta)[i,j]
      vals[i,j] = isnan(_v) ? vals[i-1, j] : _v * k + vals[i-1, j] * (1-k)
    end
  end

  cname   = Symbol[]
  cols    = colnames(ta)
  for c in 1:length(cols)
    push!(cname, Symbol(string(cols[c], "_ema_", n)))
  end

  TimeArray(tstamps, vals[n:length(ta),:], cname, meta(ta))
end

function kama(ta::TimeArray, n::Integer = 10, fn::Integer = 2, sn::Integer = 30)
  vola = moving(nansum, abs.(ta .- lag(ta)), n)
  change = abs.(ta .- lag(ta, n))
  er = safediv.(change, vola)  # Efficiency Ratio

  # Smooth Constant
  sc = (er .* (2 / (fn + 1) - 2 / (sn + 1)) .+ 2 / (sn + 1)).^2

  cl = ta[n+1:end]

  vals = similar(Array{Float64}, axes(values(cl)))
  # using simple moving average as initial kama
  pri_kama = nanmean(values(ta[1:n]))

  @assert length(cl) == length(sc)

  for idx ∈ 1:length(cl)
    _p_k = pri_kama .+ values(sc[idx]) .* (values(cl[idx]) .- pri_kama)

    #Added check for NaN and adequately backfill with the correct last available pri_kama
    if typeof(pri_kama) == Float64
      pri_kama = isnan.(_p_k)[1] ? pri_kama : _p_k
    else
      _p_k[isnan.(_p_k)] .= pri_kama[isnan.(_p_k)]
      pri_kama = _p_k
    end

    vals[idx, :] = pri_kama

  end

  cols =
  if length(colnames(ta)) == 1
    [:kama]
  else
    [Symbol("$(c)_kama") for c in colnames(ta)]
  end

  TimeArray(timestamp(cl), vals, cols)
end

function env(ta::TimeArray, n::Integer; e::AbstractFloat = 0.1)
  tstamps = timestamp(ta)[n:end]

  s = sma(ta, n)

  upper = values(s) .* (1 + e)
  lower = values(s) .* (1 - e)

  up_cnames = Symbol.(["$(string.(c_name))_env_$(n)_up" for c_name in colnames(ta)])
  lw_cnames = Symbol.(["$(string.(c_name))_env_$(n)_low" for c_name in colnames(ta)])

  u = TimeArray(tstamps, upper, up_cnames, meta(ta))
  l = TimeArray(tstamps, lower, lw_cnames, meta(ta))

  merge(l, u; method = :inner)
end

# Array dispatch for use by other algorithms

function sma(a::Array, n::Integer)
  vals = zeros(size(a,1) - (n-1), size(a,2))

  for i in 1:size(a,1) - (n-1)
    for j in 1:size(a,2)
      vals[i,j] = nanmean(a[i:i+(n-1),j])
    end
  end

  vals
end

function ema(a::Array, n::Integer; wilder = false)
  k = (1 + !wilder) / (n + !wilder)

  vals = zeros(size(a,1), size(a,2))
  # seed with first value with an sma value
  vals[n,:] = sma(a, n)[1,:]

  for i in n+1:size(a,1)
    for j in 1:size(a,2)
      vals[i,j] = a[i,j] * k + vals[i-1, j] * (1-k)
    end
  end

  vals[n:end, :]
end

function env(a::AbstractArray, n::Integer; e::AbstractFloat = 0.1)
  s = sma(a, n)

  upper = @. s * (1 + e)
  lower = @. s * (1 - e)

  [lower upper]
end

@doc raw"""
    sma(ta::TimeArray, n)
    sma(A::Array, n)

Simple Moving Average

# Formula

```math
\text{SMA} = \frac{\sum_i^n{P_i}}{n}
```
"""
sma

@doc raw"""
    ema(ta::TimeArray, n, wilder = false)
    ema(A::Array, n, wilder = false)

Exponemtial Moving Average

A.k.a. exponentially weighted moving average (EWMA)

# Formula

Let ``k`` denote the degree of weighting decrease.

If parameter `wilder` is `true`, ``k = \frac{1}{n}``,
else ``k = \frac{2}{n + 1}``.

```math
\text{EMA}_t = k \times P_t + (1 - k) \times \text{EMA}_{t - 1}
```
"""
ema

@doc raw"""
    kama(ta::TimeArray, n = 10, fn = 2, sn = 30)

Kaufman's Adaptive Moving Average

# Arguments

- `n`: period

- `fn`: the fastest EMA constant

- `sn`: the slowest EMA constant

# Formula

```math
\begin{align*}
  \text{KAMA}_t     & = \text{KAMA}_{t-1} + \text{SC} \times (\text{Price} - \text{KAMA}_{t-1}) \\
  \text{SC}         & = (\text{ER} \times (\frac{2}{fn + 1} - \frac{2}{sn + 1}) + \frac{2}{sn + 1})^2 \\
  \text{ER}         & = \frac{Change}{Volatility} \\
  \text{Change}     & = | \text{Price} - \text{Price}_{t-n} | \\
  \text{Volatility} & = \sum_{i}^{n} | \text{Price}_i - \text{Price}_{i - 1} |
\end{align*}
```
"""
kama

@doc raw"""
    env(ta::TimeArray, n; e = 0.1)
    env(A::AbstractArray, n; e = 0.1)

Moving Average Envelope

```math
\begin{align*}
  \text{Upper Envelope} & = \text{SMA}(n) \times (1 + e) \\
  \text{Lower Envelope} & = \text{SMA}(n) \times (1 - e)
\end{align*}
```

# Arguments

- `e`: the envelope, `0.1` implies the `10%` envelope.

# References

- [TradingView](https://www.tradingview.com/wiki/Envelope_(ENV))
"""
env


text:
function lab_13(r::Robot)
    for i in ((Nord, Ost), (Nord, West), (Sud, Ost), (Sud, West))
        go_special!(r, i)
        return_special!(r, i)
    end
    putmarker!(r)
end

function return_special!(r,side)
    while (ismarker(r))
        move!(r, inverse(side[1]))
        move!(r, inverse(side[2]))
    end
end

function go_special!(r, side::NTuple{2,HorizonSide})
    while isborder(r,side[1]) == false  &&   isborder(r,side[2]) == false
        move!(r,side[1])
        move!(r,side[2])
        putmarker!(r)
    end
end

inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))


text:
using Tmatrix
import Zygote
import FiniteDifferences

finite_difference_method = FiniteDifferences.central_fdm(5, 1; factor=1e-6)

# inputs ====================================================
n_θ_points = 5
n_ϕ_points = 2
m, n, m_, n_ = 2, 3, 2, 5
n_max = 1
k1_r, k1_i, k2_r, k2_i = 1e5, 1e3, 2e5, 3e3
rotationally_symmetric = true
symmetric_about_plane_perpendicular_z = false
rx, rz = 1e-6, 1.2e-6
kind = "regular"
J_superscript = 11
Q_superscript = J_superscript

θ_1D_array = LinRange(1e-16, π, n_θ_points);
ϕ_1D_array = LinRange(1e-16, 2π, n_ϕ_points);
if rotationally_symmetric
    θ_array = collect(θ_1D_array)
    ϕ_array = zeros(size(θ_array))
else
    # eventually, this should be removed. I just keep it for sanity checks.
    θ_array, ϕ_array = Tmatrix.meshgrid(θ_1D_array, ϕ_1D_array);
end
r_array, n̂_array = Tmatrix.ellipsoid(rx, rz, θ_array);

epss = 1e-6

# test "J_mn_m_n__integrand_SeparateRealImag" ====================================================================================
J = Tmatrix.J_mn_m_n__integrand_SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript)
J_ = Tmatrix.J_mn_m_n__integrand(m, n, m_, n_, complex(k1_r, k1_i) .* r_array, complex(k2_r, k2_i) .* r_array, r_array, θ_array, ϕ_array, n̂_array; kind=kind, J_superscript=J_superscript)

println("==================================================================")
println("testing <<J_mn_m_n__integrand_SeparateRealImag>>")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(J))
println("Displaying result from complex function ------------------------------")
display(J_)

println("testing autodiff for <<J_mn_m_n__integrand_SeparateRealImag>> -----------------------------------------")
@time J_mn_m_n__integrand_SeparateRealImag_jacobian = Zygote.jacobian(Tmatrix.J_mn_m_n__integrand_SeparateRealImag, m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array,kind,J_superscript)
function J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)
    Tmatrix.J_mn_m_n__integrand_SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript)
end
@time J_mn_m_n__integrand_SeparateRealImag_lean_jacobian = Zygote.jacobian(J_mn_m_n__integrand_SeparateRealImag_lean, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)
@time J_jacobian_finite_diff = FiniteDifferences.jacobian(finite_difference_method, J_mn_m_n__integrand_SeparateRealImag_lean, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)
@time J_jacobian_manual_numerical_diff = (
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r+epss, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r-epss, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)) / (2*epss),
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i+epss, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i-epss, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)) / (2*epss),
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r+epss, k2_i, r_array, θ_array, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r-epss, k2_i, r_array, θ_array, ϕ_array, n̂_array)) / (2*epss),
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i+epss, r_array, θ_array, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i-epss, r_array, θ_array, ϕ_array, n̂_array)) / (2*epss),
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i, r_array.+epss, θ_array, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i, r_array.-epss, θ_array, ϕ_array, n̂_array)) ./ (2*epss),
    (J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i, r_array, θ_array.+epss, ϕ_array, n̂_array) - J_mn_m_n__integrand_SeparateRealImag_lean(k1_r, k1_i, k2_r, k2_i, r_array, θ_array.-epss, ϕ_array, n̂_array)) ./ (2*epss),
)
finite_difference_method = FiniteDifferences.central_fdm(5, 1; factor=1e-10)
@time J_jacobian_finite_diff = FiniteDifferences.jacobian(finite_difference_method, J_mn_m_n__integrand_SeparateRealImag_lean, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)
hcat(J_mn_m_n__integrand_SeparateRealImag_lean_jacobian[1], J_jacobian_finite_diff[1])

# test "J_mn_m_n__SeparateRealImag" ================================================================================================
J = Tmatrix.J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript, rotationally_symmetric)
J_ = Tmatrix.J_mn_m_n_(m, n, m_, n_, complex(k1_r, k1_i) .* r_array, complex(k2_r, k2_i) .* r_array, r_array, θ_array, ϕ_array, n̂_array; kind=kind, J_superscript=J_superscript, rotationally_symmetric=rotationally_symmetric)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(J))
println("Displaying result from complex function ------------------------------")
display(J_)

# the problem that Zygote is trying to differentiate "rotationally_symmetric".
# how can I prevent Zygote from differentiating with respect to a given argument?
Zygote.jacobian(
    (m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) -> Tmatrix.J_mn_m_n__SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, J_superscript, rotationally_symmetric),
    m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array
)

Zygote.jacobian(J_mn_m_n__integrand_SeparateRealImag_lean, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array)[5]

epss = big(1e-6);
vcat(
    (
        J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array.+epss), big.(θ_array), big.(ϕ_array), [big.(x) for x in n̂_array]) -
        J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array.-epss), big.(θ_array), big.(ϕ_array), [big.(x) for x in n̂_array])) ./
        (2*epss)...
)

single_one_array = zeros(size(r_array)); single_one_array[1] = 1

for epss in [big(1e-50)]
    display(vcat(
        (
            J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array+epss*single_one_array), big.(θ_array), big.(ϕ_array), [big.(x) for x in n̂_array]) -
            J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array-epss*single_one_array), big.(θ_array), big.(ϕ_array), [big.(x) for x in n̂_array])) ./
            (2*epss)...
    ))
    println()

    display(vcat(
        (
            J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array), big.(θ_array+epss*single_one_array), big.(ϕ_array), [big.(x) for x in n̂_array]) -
            J_mn_m_n__integrand_SeparateRealImag_lean(big(k1_r), big(k1_i), big(k2_r), big(k2_i), big.(r_array), big.(θ_array-epss*single_one_array), big.(ϕ_array), [big.(x) for x in n̂_array])) ./
            (2*epss)...
    ))
    println()
end

# test "Q_mn_m_n_SeparateRealImag"
@time Q = Tmatrix.Q_mn_m_n_SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, Q_superscript, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
@time Q_ = Tmatrix.Q_mn_m_n_(m, n, m_, n_, complex(k1_r, k1_i), complex(k2_r, k2_i), complex(k1_r, k1_i) .* r_array, complex(k2_r, k2_i) .* r_array, r_array, θ_array, ϕ_array, n̂_array; kind=kind, Q_superscript=Q_superscript, rotationally_symmetric=rotationally_symmetric)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(Q))
println("Displaying result from complex function ------------------------------")
display(Q_)

Zygote.jacobian(
    (m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) -> Tmatrix.Q_mn_m_n_SeparateRealImag(m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, Q_superscript, rotationally_symmetric, symmetric_about_plane_perpendicular_z),
    m, n, m_, n_, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array
)

# test "Q_matrix_SeparateRealImag"
@time Q = Tmatrix.Q_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
@time Q_ = Tmatrix.Q_matrix(n_max, complex(k1_r, k1_i), complex(k2_r, k2_i), complex(k1_r, k1_i) .* r_array, complex(k2_r, k2_i) .* r_array, r_array, θ_array, ϕ_array, n̂_array, kind=kind, rotationally_symmetric=rotationally_symmetric, symmetric_about_plane_perpendicular_z=false, verbose=false)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(Q))
println("Displaying result from complex function ------------------------------")
display(Q_)

# The jacobian calculation is so slow. For n_max = 1, it takes 52.409485 seconds. This is slow compared to 0.015419 seconds required to evaluate the function
Zygote.jacobian(
    (n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array) -> Tmatrix.Q_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, kind, rotationally_symmetric, symmetric_about_plane_perpendicular_z),
     n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array
)

# test "T_matrix_SeparateRealImag"
@time T = Tmatrix.T_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
@time T_ = Tmatrix.T_matrix(n_max, complex(k1_r, k1_i), complex(k2_r, k2_i), complex(k1_r, k1_i) .* r_array, complex(k2_r, k2_i) .* r_array, r_array, θ_array, ϕ_array, n̂_array; rotationally_symmetric=rotationally_symmetric, symmetric_about_plane_perpendicular_z=false)

println("==================================================================")
println("Displaying results from _SeparateRealImag and the complex function")
println("Displaying result from _SeparateRealImag function --------------------")
display(Tmatrix.get_complex_matrix_from_concatenated_real_imag(T))
println("Displaying result from complex function ------------------------------")
display(T_)

@time T = Tmatrix.T_matrix_SeparateRealImag(n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, rotationally_symmetric, symmetric_about_plane_perpendicular_z)
@time ∂T = Zygote.jacobian(Tmatrix.T_matrix_SeparateRealImag, n_max, k1_r, k1_i, k2_r, k2_i, r_array, θ_array, ϕ_array, n̂_array, rotationally_symmetric, symmetric_about_plane_perpendicular_z)


text:
# constants in cgs units
const c = 2.997824582e10     # cm s^-1
const σ = 5.6704e-5          # erg cm^-2 s^-1 K^-4
const kB = 1.3806485e-16     # erg K^-1
const mH = 1.6737236e-24     # gm

# constants in other units
const c_ms = 2.99782458e8    # m s^-1
const c_kms = 2.99782458e5   # km s^-1

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
# Returns index corresponding to the name in NamedTupleDist.
# If the name corresponds to a multivariate distribution,
# an array containig all indices of this distribution is returned.
function asindex(vs::NamedTupleShape, name::Symbol)
    name = string(name)
    if name in (rn = repeatednames(vs))
        indxs = findall(x-> x == string(name), rn)
        length(indxs) > 1 ? (return indxs) : (return indxs[1])
    else
        throw(ArgumentError("Symbol :$name found not in $vs."))
    end
end

# Returns index corresponding to the name in NamedTupleDist.
function asindex(vs::NamedTupleShape, name::Expr)
    name = string(name)
    if name in (an = all_active_names(vs))
        return findfirst(x-> x == name, an)[1]
    else
        throw(ArgumentError("Symbol :$name found not in $vs."))
    end
end

# for samples
function isshaped(samples::BAT.DensitySampleVector)
    isa(varshape(samples), NamedTupleShape) ? (return true) : (return false)
end

function asindex(samples::BAT.DensitySampleVector, name::Union{Expr, Symbol})
    if isshaped(samples)
        return asindex(varshape(samples), name)
    else
        throw(ArgumentError("Samples are unshaped. Key :$name cannot be matched. Use index instead."))
    end
end

# for prior
function asindex(ntd::NamedTupleDist, name::Union{Expr, Symbol})
    return asindex(varshape(ntd), name)
end

function asindex(
    x::Union{DensitySampleVector, NamedTupleDist, MarginalDist},
    key::Integer
)
    return key
end

#for MarginalDist
function asindex(marg::MarginalDist, name::Union{Expr, Symbol})
    idx = asindex(marg.origvalshape, name)
    if idx in marg.dims
        return idx
    else
        throw(ArgumentError("Key :$name not in MarginalDist"))
    end
end

# Return the name corresponding to the index as Symbol (for univariate)
# or Expr for (multivariate) distributions
function getname(vs::NamedTupleShape, idx::Integer)
    names = all_active_names(vs)(vs)
    return Meta.parse(names[idx])
end

# Return the name corresponding to the index as String
function getstring(vs::NamedTupleShape, idx::Integer)
    names = all_active_names(vs)
    return names[idx]
end

function getstring(prior::NamedTupleDist, idx::Integer)
    vs = varshape(prior)
    names = all_active_names(vs)
    return names[idx]
end

function getstring(samples::BAT.DensitySampleVector, idx::Integer)
    if isshaped(samples)
        vs = varshape(samples)
        names = all_active_names(vs)
        return names[idx]
    else
        throw(ArgumentError("Samples are unshaped. Key :$name cannot be matched. Use index instead."))
    end
end

function getstring(marg::MarginalDist, idx::Integer)
    if idx in marg.dims
        vs = marg.origvalshape
        names = all_active_names(vs)
        return names[idx]
    else
        throw(ArgumentError("Index $idx not in MarginalDist"))
    end
end

# Return array of strings with the names of all indices.
# For a multivariate distribution, names for each dimension are created by appending "[i]" to the name.
function allnames(vs::NamedTupleShape)
    subindxs = [collect(1:getproperty(vs, k).len) for k in keys(vs)]
    names = [length(subindxs[i]) > 1 ? subname.(string(keys(vs)[i]), subindxs[i]) : [string(keys(vs)[i])] for i in 1:length(vs)]

    return reduce(vcat, names)
end

function all_active_names(vs::NamedTupleShape)
    subindxs = [collect(1:getproperty(vs, k).len) for k in keys(vs)]
    names = [length(subindxs[i]) > 1 ? subname.(string(keys(vs)[i]), subindxs[i]) : [string(keys(vs)[i])] for i in 1:length(vs) if getproperty(vs, keys(vs)[i]).len>0]

    return reduce(vcat, names)
end

# Return array of strings with the names of all indices.
# For a multivariate distribution, the name is repeated for all its dimensions.
function repeatednames(vs::NamedTupleShape)
    names = [[string(k) for i in 1:getproperty(vs, k).len] for k in keys(vs) if getproperty(vs, k).len>0]
    return reduce(vcat, names)
end

# Generate sub-name for one dimension of a multivariate distribution
function subname(name::String, indx::Integer)
    return name*"["*string(indx)*"]"
end

function get_fixed_names(vs::NamedTupleShape)
    active_names = all_active_names(vs)
    all_names = allnames(vs)
    fixed_names = [n for n in all_names if !in(n, active_names)]
end


text:
using ControlBenchmarks
using Test
using SafeTestsets

@testset "ControlBenchmarks.jl" begin

@testset "Benchmarks" begin
    # Run tests for all the benchmarks
    include( "benchmarks/runtests.jl" )
end

end


text:
using CollisionDetection

"""
    weld(mesh1, mesh2, ...;boundary=false) -> welded_mesh

Build a mesh by welding or pasting together the inputs. Vertices from different
meshes that coincide up to the tolerance will be merged into one. The order
cells appear in the output mesh is equal to the order in the inputs.
boundary = true; will merge only vertices of different meshes if one of
them exists on the 'open' boundary of the mesh.
"""
function weld end

function weld(Γ₁, Γ₂; boundary=false)

    # T = eltype(eltype(Γ₁.vertices))
    T = coordtype(Γ₁)
    tol = sqrt(eps(T))

    verts1 = skeleton(Γ₁,0) # TODO: check this, it might break for unused verts
    radii = zeros(numcells(verts1))
    indcs = [v[1] for v in cells(verts1)]
    cntrs = [cartesian(center(chart(verts1,v))) for v in cells(verts1)]
    tree = Octree(cntrs, radii)

    nv1 = numvertices(Γ₁)
    nv2 = numvertices(Γ₂)

    idmap = collect(nv1 .+ (1:nv2))
    num_equal_vertices = 0

    V1 = vertices(Γ₁)
    V2 = similar(V1,0)

    in_interior1 = interior_vpredicate(Γ₁)
    in_interior2 = interior_vpredicate(Γ₂)
    cpred(u,v,i,j) = boundary ? norm(u-v) < tol && (!in_interior1(indcs[i])
                || !in_interior2(j)) : norm(u-v) < tol

    for (j,v) in enumerate(Γ₂.vertices)
        found = false
        pred(c,s) = fitsinbox(Array(v), 0.0, c, s+tol)
        for box in boxes(tree, pred)
            for i in box
                u = cntrs[i]
                if cpred(u,v,i,j) #norm(u-v) < tol
                    idmap[j] = indcs[i]
                    num_equal_vertices += 1
                    found = true
                    # we implicitly assume that the first copy of a point
                    # appearing in the vertex buffer is the one actually
                    # referred to in the index buffer. DANGEROUS!
                    break
                end
            end
            found && break
        end
        if !found
            idmap[j] -= num_equal_vertices
            push!(V2, v)
        end
    end

    V = [V1; V2]
    F = [Γ₁.faces; Γ₂.faces]

    nc1 = numcells(Γ₁)
    # for (i,c) in enumerate(Γ₂.faces)
    for (i,c) in enumerate(cells(Γ₂))
        F[nc1+i] = map_ids(c, idmap)
    end

    return Mesh(V,unique(F))

end

function weld(G1::SComplex2D, G2::SComplex2D; seam)

    T = coordtype(G1)
    P = vertextype(G1)
    tol = sqrt(eps(T))

    nodes1 = skeleton(G1,0)
    edges1 = skeleton(G1,1)

    node_ctrs1 = [cartesian(center(chart(nodes1,c))) for c in cells(nodes1)]
    edge_ctrs1 = [cartesian(center(chart(edges1,c))) for c in cells(edges1)]

    node_radii1 = zeros(T, length(node_ctrs1))
    edge_radii1 = [volume(chart(edges1,c)) for c in cells(edges1)] / 2

    node_tree1 = Octree(node_ctrs1, node_radii1)
    edge_tree1 = Octree(edge_ctrs1, edge_radii1)

    nodes2 = skeleton(G2,0)
    edges2 = skeleton(G2,1)

    node_is_on_seam = inclosure_gpredicate(seam)
    edge_is_on_seam = overlap_gpredicate(seam)
    # is_interior_edge = interior_tpredicate(G1)
    # is_interior_node = interior_vpredicate(G1)

    Nodes2 = similar(G2.nodes, 0)
    Edges2 = similar(G2.edges, 0)

    num_equal_edges = 0
    edges_map = collect(numcells(edges1) .+ (1:numcells(edges2)))
    for (j,edge2) in enumerate(cells(edges2))
        found = false
        ctr2 = cartesian(center(chart(edges2, edge2)))
        rad2 = volume(chart(edges2,edge2)) / 2
        for box in boxes(edge_tree1, (c,s)->fitsinbox(ctr2, rad2, c, s+tol))
            for i in box
                ctr1 = edge_ctrs1[i]
                # For now we assume edges are the same if their centroids coincide
                if norm(ctr1-ctr2) < tol && edge_is_on_seam(chart(edges1,cells(edges1)[i]))
                # if norm(ctr1-ctr2) < tol && !is_interior_edge(cells(edges1)[i])
                    edges_map[j] = i
                    num_equal_edges += 1
                    found = true
                    break
                end
                found && break
            end
        end
        if !found
            push!(Edges2, edge2)
            edges_map[j] -= num_equal_edges
        end
    end

    Faces2 = similar(G2.faces, 0)
    for face in cells(G2)
        push!(Faces2,map(p->sign(p)*getindex(edges_map,abs(p)),face))
    end

    num_equal_nodes = 0
    nodes_map = collect(numcells(nodes1) .+ (1:numcells(nodes2)))
    for (j,node2) in enumerate(cells(nodes2))
        found = false
        ctr2 = cartesian(center(chart(nodes2, node2)))
        rad2 = 0.0
        for box in boxes(node_tree1, (c,s)->fitsinbox(ctr2, rad2, c, s+tol))
            for i in box
                ctr1 = node_ctrs1[i]
                # For now we assume edges are the same if their centroids coincide
                if norm(ctr1-ctr2) < tol && node_is_on_seam(ctr1)
                # if norm(ctr1-ctr2) < tol && !is_interior_node(cells(nodes1)[i][1])
                    nodes_map[j] = i
                    num_equal_nodes += 1
                    found = true
                    break
                end
                found && break
            end
        end
        if !found
            push!(Nodes2, node2)
            nodes_map[j] -= num_equal_nodes
        end
    end

    # Edges2 = similar(G2.edges, 0)
    for (i,edge) in enumerate(Edges2)
        edge_mapped = map(p->sign(p)*getindex(nodes_map,abs(p)),edge)
        Edges2[i] = edge_mapped
        # push!(Edges2, edge_mapped)
    end

    SComplex2D{T,P}(
        vcat(G1.vertices, G2.vertices),
        vcat(G1.nodes, Nodes2),
        vcat(G1.edges, Edges2),
        vcat(G1.faces, Faces2))
end

@generated function map_ids(c::SVector{N,T}, idmap) where {N,T}
    xp = :(SVector{$N,$T}())
    for i in 1:N
        push!(xp.args, :(idmap[c[$i]]))
    end
    return xp
end

function weld(meshes...; boundary=false)
    @assert length(meshes) > 1
    G = weld(meshes[1], meshes[2], boundary=boundary)
    for i in 3:length(meshes)
        G = weld(G, meshes[i], boundary=boundary)
    end
    G
end


text:
module PowerVega

import InfrastructureModels
import PowerModels
using Statistics, LightGraphs

using VegaLite, VegaDatasets
using DataFrames
using JSON
using FilePaths

import PyCall
const nx = PyCall.PyNULL()
const scipy = PyCall.PyNULL()

function __init__()
    copy!(nx, PyCall.pyimport_conda("networkx", "networkx"))
    copy!(scipy, PyCall.pyimport_conda("scipy", "scipy"))
end

include("core/data.jl")
include("core/layout.jl")
include("core/plot.jl")

include("core/export.jl")  # must be last to properly export all functions

end # module


text:
abstract type Grid end

function _is_monotonic_increasing{T}(x::Array{T, 1})
    if length(x) < 2
        return true
    end

    all(elem > 0 for elem in x[2:end] .- x[1:end-1])
end

struct RegularGrid <: Grid
    left :: Float64
    right :: Float64
    points :: Int
    vertices :: Array{Float64, 1}

    function RegularGrid(left=0, right=1, points=50)
        return new(left, right, points, linspace(left, right, points))
    end
end

struct ArbitraryGrid <: Grid
    vertices :: Array{Float64, 1}
    monotonic :: Bool

    ArbitraryGrid(data) = new(data, _is_monotonic_increasing(data))
end

grid_step(g::RegularGrid) = (g.right - g.left) / g.points

vertices(g::RegularGrid) = g.vertices
vertices(g::ArbitraryGrid) = g.vertices

is_monotonic(g::RegularGrid) = true
is_monotonic(g::ArbitraryGrid) = g.monotonic


text:
mutable struct HSDEStatus <: AbstractStatus
    m::Int64
    n::Int64
    i::Int64
    model::FOSMathProgModel
    status::Symbol
    checki::Int64
    eps::Float64
    verbose::Int64
    checked::Bool
    direct::Bool
    init_time::UInt64     #As reported by time_ns()
    init_duration::UInt64 #In nano-seconds
    debug::Int64
end

""" checkstatus(stat::HSDEStatus, x)
Returns `false` if no check was made
If convergence check is done, returns `true` and sets stat.status to one of:
    :Continue, :Optimal, :Unbounded, :Infeasible
"""
function checkstatus(stat::HSDEStatus, z; override = false)
    #TODO fix m, n
    if stat.i % stat.checki == 0 || override
        t = time_ns() - stat.init_time
        m, n, i, model, verbose, debug, ϵ = stat.m, stat.n, stat.i, stat.model, stat.verbose, stat.debug, stat.eps
        x, y, s, r, τ, κ = getvalues(model, m, n, z, i, verbose, debug)
        ϵpri = ϵdual = ϵgap = ϵinfeas = ϵunbound = ϵ
        p   = norm(model.A*x/τ + s/τ - model.b)/norm(1+norm(model.b))
        d   = norm(model.A'*y/τ + model.c - r/τ)/norm(1+norm(model.c))
        ctx = (model.c'*x)[1]
        bty = (model.b'*y)[1]
        g   = norm(ctx/τ + bty/τ)/(1+abs(ctx/τ)+abs(bty/τ))
        if debug > 0
            savedata(i, p, d, g, ctx, bty, κ, τ, x, y, s, t, model, debug)
        end
        # TODO test u^T v
        if verbose > 0
            if !stat.direct
                cgiter = getcgiter(stat.model.data)
                push!(model.history, :cgiter, i, cgiter)
                printstatusiter(i, p, d, g, ctx, bty, κ/τ, cgiter, t)
            else
                printstatusiter(i, p, d, g, ctx, bty, κ/τ, t)
            end
        end

        status = :Continue
        if p <= ϵpri*(1+norm(model.b)) && d <= ϵdual*(1+norm(model.c)) && g <= ϵgap*(1+norm(ctx/τ)+norm(bty/τ))
            if verbose > 0
                println("Found solution i=$i")
            end
            status =  :Optimal
        elseif norm(model.A*x + s) <= ϵunbound*(-ctx/norm(model.c))
            status = :Unbounded
        elseif norm(model.A'*y) <= ϵinfeas*(-bty/norm(model.b))
            status = :Infeasible
        end
        stat.status = status
        stat.checked = true
        return stat.checked
    else
        stat.checked = false
        return stat.checked
    end
end

function printstatusheader(stat::HSDEStatus)
    if stat.verbose > 0
        println("Time to initialize: $(stat.init_duration/1e9)s")
        width = 76 + (stat.direct ? 0 : 5)
        println("-"^width)
        print(" Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau")
        !stat.direct && print(" | cg ")
        println(" | time")
        println("-"^width)
    end
end

function printstatusiter(i, p, d, g, stx, bty, κτ, t)
    @printf("%6d|% 9.2e % 9.2e % 9.2e % 9.2e % 9.2e % 9.2e % .1es\n",i,p,d,g,stx,-bty,κτ,t/1e9)
end

function printstatusiter(i, p, d, g, stx, bty, κτ, cgiter, t)
    @printf("%6d|% 9.2e % 9.2e % 9.2e % 9.2e % 9.2e % 9.2e % 4d % .1es\n",i,p,d,g,stx,-bty,κτ,cgiter,t/1e9)
end

function getvalues(model, m, n, z, i, verbose, debug)
    nu = n+m+1
    x = view(z, 1:n)
    y = view(z, (n+1):(n+m))
    r = view(z, (nu+1):(nu+n))
    s = view(z, (nu+n+1):(nu+n+m))
    τ = z[nu]
    κ = z[2nu]
    x, y, s, r, τ, κ
end

# """
# Checks primal and dual feasibility assuming that x,y,s are satisfying the cone constraints
# """
# function printstatus(model, x, y, s, r, τ, κ, i, verbose, debug)
#     p   = norm(model.A*x/τ + s/τ - model.b)/norm(1+norm(model.b))
#     d   = norm(model.A'*y/τ + model.c - r/τ)/norm(1+norm(model.c))
#     ctx = (model.c'*x)[1]
#     bty = (model.b'*y)[1]
#     g   = norm(ctx/τ + bty/τ)/(1+abs(ctx/τ)+abs(bty/τ))
#     if debug > 0
#         savedata(i, p, d, g, ctx, bty, κ, τ, x, y, s, model, debug)
#     end
#     # TODO test u^T v
#     if verbose > 0
#         printstatusiter(i, p, d, g, ctx, bty, κ/τ)
#     end
#     return
# end

@generated function savedata(i, p, d, g, ctx, bty, κ, τ, x, y, s, t, model, debug)
    ex = :(history = model.history)
    for (v,vs) in [(:p, "p"), (:d, "d"), (:g, "g"),
                   (:ctx, "ctx"), (:bty, "bty"), (:κ, "κ"),
                   (:τ, "τ"), (:t, "t")]
        ex = :($ex ; push!(history, Symbol($vs), i, $v))
    end
    ex2 = :()
    for (v,vs) in [(:x, "x"), (:y,"y"), (:s,"s")]
        ex2 = :($ex2 ; push!(history, Symbol($vs), i, $v))
    end
    ex = :($ex; if debug > 1; $ex2; end)
    ex = :($ex; return)
    return ex
end


text:
@inline function U(chain::EAPChain)
  return (sum(chain.us) + U_interaction(chain)
          - dot(end_to_end(chain), [chain.Fx; 0.0; chain.Fz]));
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


text:
mutable struct GridRoomsDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(GridRoomsDirected)
@generate_setters(GridRoomsDirected)

mutable struct GridRoomsUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Goal}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    goal_pos::CartesianIndex{2}
    done::Bool
end

@generate_getters(GridRoomsUndirected)
@generate_setters(GridRoomsUndirected)

#####
# Directed
#####

function GridRoomsDirected(; T = Float32, grid_size = (2, 2), room_size = (5, 5), rng = Random.GLOBAL_RNG)
    grid_height = grid_size[1]
    grid_width = grid_size[2]
    room_height = room_size[1]
    room_width = room_size[2]

    height = room_height * grid_height - grid_height + 1
    width = room_width * grid_width - grid_width + 1
    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    origins = get_room_origins(grid_size, room_size)

    for origin in origins
        room = Room(origin, room_height, room_width)
        place_room!(world, room)

        world[WALL, room.region[(end + 1) ÷ 2, 1]] = false
        world[WALL, room.region[(end + 1) ÷ 2, end]] = false
        world[WALL, room.region[1, (end + 1) ÷ 2]] = false
        world[WALL, room.region[end, (end + 1) ÷ 2]] = false
    end

    world[WALL, [1, height], :] .= true
    world[WALL, :, [1, width]] .= true

    agent_pos = CartesianIndex(2, 2)
    agent_dir = RIGHT
    reward = zero(T)
    terminal_reward = one(T)
    goal_pos = CartesianIndex(height - 1, width - 1)
    done = false

    env = GridRoomsDirected(world, agent_pos, agent_dir, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::GridRoomsDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const GRID_ROOMS_DIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::GridRoomsDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), GRID_ROOMS_DIRECTED_LAYERS)

RLBase.state_space(env::GridRoomsDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::GridRoomsDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::GridRoomsDirected, ::RLBase.DefaultPlayer) = DIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::GridRoomsDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::GridRoomsDirected) = get_done(env)

function RLBase.reset!(env::GridRoomsDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::GridRoomsDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::GridRoomsDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)

    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

function GridRoomsUndirected(; T = Float32, grid_size = (2, 2), room_size = (5, 5), rng = Random.GLOBAL_RNG)
    grid_height = grid_size[1]
    grid_width = grid_size[2]
    room_height = room_size[1]
    room_width = room_size[2]

    height = room_height * grid_height - grid_height + 1
    width = room_width * grid_width - grid_width + 1
    objects = (AGENT, WALL, GOAL)
    world = GridWorldBase(objects, height, width)

    origins = get_room_origins(grid_size, room_size)

    for origin in origins
        room = Room(origin, room_height, room_width)
        place_room!(world, room)

        world[WALL, room.region[(end + 1) ÷ 2, 1]] = false
        world[WALL, room.region[(end + 1) ÷ 2, end]] = false
        world[WALL, room.region[1, (end + 1) ÷ 2]] = false
        world[WALL, room.region[end, (end + 1) ÷ 2]] = false
    end

    world[WALL, [1, height], :] .= true
    world[WALL, :, [1, width]] .= true

    agent_pos = CartesianIndex(2, 2)
    reward = zero(T)
    terminal_reward = one(T)
    goal_pos = CartesianIndex(height - 1, width - 1)
    done = false

    env = GridRoomsUndirected(world, agent_pos, reward, rng, terminal_reward, goal_pos, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::GridRoomsUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const GRID_ROOMS_UNDIRECTED_LAYERS = SA.SVector(2, 3)
RLBase.state(env::GridRoomsUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), GRID_ROOMS_UNDIRECTED_LAYERS)

RLBase.state_space(env::GridRoomsUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::GridRoomsUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::GridRoomsUndirected, player::RLBase.DefaultPlayer) = UNDIRECTED_NAVIGATION_ACTIONS
RLBase.reward(env::GridRoomsUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::GridRoomsUndirected) = get_done(env)

function RLBase.reset!(env::GridRoomsUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[GOAL, get_goal_pos(env)] = false

    new_goal_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_goal_pos!(env, new_goal_pos)
    world[GOAL, new_goal_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_reward!(env, zero(T))
    set_done!(env, false)

    return nothing
end

function (env::GridRoomsUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)

    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[GOAL, get_agent_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Common
#####

function get_room_origins(grid_size, room_size)
    grid_height = grid_size[1]
    grid_width = grid_size[2]
    room_height = room_size[1]
    room_width = room_size[2]

    vertical_steps = 1 : room_height - 1 : 1 + (grid_height - 1) * (room_height - 1)
    horizontal_steps = 1 : room_width - 1 : 1 + (grid_width - 1) * (room_width - 1)

    return [CartesianIndex(i, j) for i in vertical_steps for j in horizontal_steps]
end


text:
#@reexport module BlazePolynomialModel

#export fit_blaze_1d, make_flat_blaze, make_blaze, calc_blaze, eval_blaze
#export calc_normalizations, renormalize

#import ..AbstractBlazeType, ..AbstractBlaze1DType, ..AbstractBlaze2DType
import ..AbstractEchelleObservation, ..AbstractEchelleObservationSet
import ..default_pixel_buffer_lo, ..default_pixel_buffer_hi
import ..max_num_pixels_per_order
import ..num_pixels_per_order, ..num_orders, ..num_obs
import ..default_blaze_degree
#import ..calc_blaze
#using ..BlazeCommon

# Constants to module
const max_blaze_degree = 6
export max_blaze_degree
#const default_blaze_degree = 4 # TODO: move to each instrument?
#export default_blaze_degree

using Polynomials     # For blaze function
using Statistics

""" Type for storing a polynomial model for one order's blaze function. """
struct BlazePolynomial{T1,T2} <: AbstractBlaze1DType where {T1<:Integer, T2<:Real}
   order::T1
   log_lambda_lo::T2
   log_lambda_hi::T2
   pixel_range::UnitRange{T1}
   model::Poly{T2}
end
export BlazePolynomial

struct BlazePolynomialList{T1,T2} <: AbstractBlaze2DType where {T1<:Integer, T2<:Real}
   blaze_list::Vector{BlazePolynomial{T1,T2}}
end
export BlazePolynomialList

function fit_blaze_1d(x::AbstractRange{T1}, flux::AbstractArray{T2,1}; degree::Int = default_blaze_degree
      ) where {T1<:Real, T2<:Real}
   polyfit(x,flux,degree) # replaced x.-mean(x)
end
export fit_blaze_1d

function make_flat_blaze_polynomial_1d(;num_pixels_per_order::Integer = max_num_pixels_per_order,
      pixel_range::AbstractRange = 1:num_pixels_per_order,
      order::Integer = 0, degree::Integer = 0, log_lambda_lo::Real = 0.0, log_lambda_hi::Real = 0.0)
      @assert 0 <= degree <= max_blaze_degree
      coeff = vcat([1.0],zeros(degree))
      return BlazePolynomial(order, log_lambda_lo, log_lambda_hi, pixel_range, Poly(coeff) )
end

function make_flat_blaze_polynomial_1d(data::EOST, order::Integer;
      pixel_range::AbstractRange = 1:size(data.lambda,1), degree::Integer = 0
      ) where { EOST <: Union{AbstractEchelleObservationSet,AbstractEchelleObservation} }
   @assert 1 <= order <= max_num_orders
   @assert 0 <= degree <= max_blaze_degree
   pixel_lo = max(1,pixel_range[1] - 1)
   pixel_hi = min(size(data.lambda,1), pixel_range[end]+1)
   log_lambda_lo = log(minimum(data.lambda[pixel_lo,order,:]))
   log_lambda_hi = log(maximum(data.lambda[pixel_hi,order,:]))
   make_flat_blaze_polynomial_1d(order=order, log_lambda_lo=log_lambda_lo, log_lambda_hi=log_lambda_hi,
                  pixel_range=pixel_range, degree=degree)
end
export make_flat_blaze_polynomial_1d

function make_flat_blaze_polynomial_2d(data::EOST;
      pixel_range::AbstractRange = 1:size(data.lambda,1), degree::Integer = 0,
      ) where { EOST <: Union{AbstractEchelleObservationSet,AbstractEchelleObservation} }
   map(order-> make_flat_blaze_polynomial_1d(pixel_range=pixel_range, order=order, degree=degree), 1:num_orders(data) )
end
export make_flat_blaze_polynomial_2d

function make_blaze_polynomial_1d(data::EOST, order::Integer;
      pixel_range::AbstractRange = 1:size(data.lambda,1),
      pixel_buffer_lo::Integer = default_pixel_buffer_lo,
      pixel_buffer_hi::Integer = default_pixel_buffer_hi,
      degree::Integer = default_blaze_degree,
      verbose::Integer = 0
      ) where { EOST <: AbstractEchelleObservationSet }
      @assert 0 <= degree <= max_blaze_degree
      @assert 1 <= order <= num_orders(data)
      max_itterations_fit_blaze = 20
      # TODO: Insert code to truncate pixel_range to enforce positive blaze
      # blaze = fit_blaze_1d(pixel_range, vec(mean(data.flux[pixel_range,order,:],dims=2)), degree = degree)
      #println("# Order = ", order, " pixel_range = ", pixel_range)
      new_pixel_range = pixel_range
      last_pixel_range = 0:0
      local blaze
      i = 0
      while !((2*abs(new_pixel_range[1] - last_pixel_range[1]) <= pixel_buffer_lo) && (2*abs(new_pixel_range[end] - last_pixel_range[end]) <= pixel_buffer_hi )) &&
                  (i<max_itterations_fit_blaze)
            last_pixel_range = new_pixel_range
            blaze = fit_blaze_1d(new_pixel_range, vec(mean(data.flux[new_pixel_range,order,:],dims=2)), degree = degree)
            pixel_lo = max(pixel_range[1],  findfirst(x->x>0,polyval(blaze,pixel_range)) ) +pixel_buffer_lo
            pixel_hi = min(pixel_range[end],findlast( x->x>0,polyval(blaze,pixel_range)) ) -pixel_buffer_hi
            if verbose > 1
                  println("# order: ", order, " it=", i," range: ", pixel_lo, ":", pixel_hi, " blaze = ", polyval(blaze,max(pixel_lo-1,1)), " ", polyval(blaze,pixel_lo), " ... ", polyval(blaze, pixel_hi) )
            end
            new_pixel_range = pixel_lo:pixel_hi
            i += 1
      end
      if verbose > -1 && !((2*abs(new_pixel_range[1] - last_pixel_range[1]) <= pixel_buffer_lo) && (2*abs(new_pixel_range[end] - last_pixel_range[end]) <= pixel_buffer_hi ))
            println("# Warning: Reached maximum number of itterations fitting blaze for order ", order)
            println("# last_pixel_range = ", last_pixel_range)
            println("# new_pixel_range = ", new_pixel_range)
      end
      # blaze = blaze / mean(polyval(blaze,new_pixel_range))
      pixel_lo = max(1,new_pixel_range[1] - 1)
      pixel_hi = min(size(data.lambda,1), new_pixel_range[end]+1)
      log_lambda_lo = log(minimum(data.lambda[pixel_lo,order,:]))
      log_lambda_hi = log(maximum(data.lambda[pixel_hi,order,:]))
      return BlazePolynomial(order, log_lambda_lo, log_lambda_hi, new_pixel_range, blaze )
end
export make_blaze_polynomial_1d

function make_blaze_polynomial_2d(data::EOST;
      pixel_range::AbstractRange = 1:size(data.lambda,1), degree::Integer = default_blaze_degree
      ) where { EOST <: AbstractEchelleObservationSet}
      map(order-> make_blaze_polynomial_1d(data, order, pixel_range=pixel_range, degree=degree), 1:num_orders(data) )
end
export make_blaze_polynomial_2d

# TODO: Change calc_blaze to eval_blaze
function calc_blaze(b::BlazePolynomial, pixel::Real)
      @assert b.pixel_range[1] <= pixel <= b.pixel_range[end]
      polyval(b.model, pixel)
end

function calc_blaze(b::BlazePolynomial, pixel_range::AbstractRange{T2}) where {T2<:Real}
      @assert b.pixel_range[1]-0.5 <= pixel_range[1] <= b.pixel_range[end]+0.5
      @assert b.pixel_range[1]-0.5 <= pixel_range[end] <= b.pixel_range[end]+0.5
      blaze = polyval(b.model, pixel_range)
end

function calc_blaze(b::BlazePolynomial) where {T2<:Real}
      blaze = polyval(b.model, b.pixel_range)
end
export calc_blaze

#=
function calc_blaze(blaze_list::BlazeLT, data::EOT;
      output::AbstractArray{T2,2} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOT <: AbstractEchelleObservation  }
      @assert length(blaze_list) == num_orders(data)
      for blaze in blaze_list
            output[blaze.pixel_range,blaze.order] .= calc_blaze(blaze)
      end
      return output
end

function calc_blaze(blaze_list::BlazeLT, data::EOST;
      output::AbstractArray{T2,3} = zeros(eltype(data.flux),size(data.lambda))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      #BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet  }
      @assert length(blaze_list) == num_orders(data)
      for blaze in blaze_list
            output[blaze.pixel_range,blaze.order,:] .= calc_blaze(blaze)
      end
      return output
end
export calc_blaze
=#

#=
function calc_normalizations(blaze_list::BlazeLT, data::EOST;
      output::OutputT = zeros(eltype(data.flux),num_obs(data))
      ) where { T2<:Real,
      BlazeLT <: AbstractBlaze2DType,
      # BlazeT <: AbstractBlazeType, BlazeLT <: AbstractArray{BlazeT,1},
      EOST <: AbstractEchelleObservationSet, OutputT <: AbstractArray{T2,1}  }
   @assert length(blaze_list) == num_orders(data)
   @assert length(output) == num_obs(data)
   output .= map(obsid->
      mean(map(order->
            mean(data.flux[blaze_list[order].pixel_range,order,obsid] ./
                  calc_blaze(blaze_list[order]) )
            , 1:num_orders(data) ))
      , 1:num_obs(data))
end
export calc_normalizations
=#

function renormalize(blaze::BlazePolynomial, norm::Real)
      return BlazePolynomial(blaze.order, blaze.log_lambda_lo, blaze.log_lambda_hi,
                  blaze.pixel_range, blaze.model*norm )
end
export renormalize

#include("common.jl")

#end # module BlazePolynomialModel


text:
# # How to make a Balance law

# ```@meta
# CurrentModule = ClimateMachine.BalanceLaws
# ```

# Defining the set of solved PDEs in ClimateMachine revolve around defining a
# [`BalanceLaw`](@ref). A balance law solves equations of the form:

# ```
# ∂Y
# -- = - ∇ • ( F_{first_order}(Y) + F_{second_order}(Y, Σ) ) + S_{non_conservative}(Y, Σ)
# ∂t
# ```

# Here, `Y`, `Σ`, `F_{first_order}`, `F_{second_order}`, and
# `S_{non_conservative}` can be thought of column vectors[^1] expressing:

#  - `Y` the prognostic state variables, or unknowns of the PDEs to be solved

#  - `Σ` the gradients of functions of the prognostic state variables

#  - `F_{first-order}` contains all first-order fluxes (e.g. **not** functions
#  of gradients of any variables)

#  - `F_{second-order}` contains all second-order and higher-order fluxes
#  (e.g. functions of gradients of any variables)

#  - `S_{non_conservative}` non-conservative sources[^2]

# In order to alleviate users from being concerned with the burden of
# spatial discretization, users must provide their own implementations of
# the following methods, which are computed locally at each nodal point:

# ## Variable name specification methods

# The `vars_state` method is used to specify the names of the variables for the following state variable types:

# | **Type** | Purpose |
# |:-----|:----|
# | [`Prognostic`](@ref) | the variables in the prognostic state vector, typically mass, momentum, and various tracers. |
# | [`Auxiliary`](@ref) | any variables required for the balance law that aren't related to derivatives of the state variables (e.g. spatial coordinates or various integrals) or those needed to solve expensive auxiliary equations (e.g., temperature via a non-linear equation solve)     |
# | [`Gradient`](@ref) | the gradients of functions of the prognostic state variables. used to represent values before **and** after differentiation |
# | [`GradientFlux`](@ref) | the gradient fluxes necessary to impose Neumann boundary conditions. typically the product of a diffusivity tensor with a gradient state variable, potentially equivalent to the second-order flux for a prognostic state variable |
# | [`UpwardIntegrals`](@ref) | any one-dimensional vertical integrals from **bottom to top** of the domain required for the balance law. used to represent both the integrand **and** the resulting indefinite integral |
# | [`DownwardIntegrals`](@ref) | any one-dimensional vertical integral from **top to bottom** of the domain required for the balance law. each variable here must also exist in `vars_state` since the reverse integral kernels use subtraction to reverse the integral instead of performing a new integral. use to represent the value before **and** after reversing direction |

# ## Methods to compute gradients and integrals
# | **Method** |  Purpose |
# |:-----|:-----|
# | [`compute_gradient_argument!`](@ref) | specify how to compute the arguments to the gradients. can be functions of prognostic state  and auxiliary variables. |
# | [`compute_gradient_flux!`](@ref) | specify how to compute gradient fluxes. can be a functions of the gradient state, the prognostic state, and auxiliary variables.|
# | [`integral_load_auxiliary_state!`](@ref) | specify how to compute integrands. can be functions of the prognostic state and auxiliary variables. |
# | [`integral_set_auxiliary_state!`](@ref) | specify which auxiliary variables are used to store the output of the integrals. |
# | [`reverse_integral_load_auxiliary_state!`](@ref) | specify auxiliary variables need their integrals reversed. |
# | [`reverse_integral_set_auxiliary_state!`](@ref) | specify which auxiliary variables are used to store the output of the reversed integrals. |
# | [`update_auxiliary_state!`](@ref) | perform any updates to the auxiliary variables needed at the beginning of each time-step. Can be used to solve non-linear equations, calculate integrals, and apply filters. |
# | [`update_auxiliary_state_gradient!`](@ref) | same as above, but after computing gradients and gradient fluxes in case these variables are needed during the update. |

# ## Methods to compute fluxes and sources
# | **Method** | Purpose |
# |:-----|:-----|
# | [`flux_first_order!`](@ref) |  specify `F_{first_order}` for each prognostic state variable. can be functions of the prognostic state and auxiliary variables. |
# | [`flux_second_order!`](@ref)    |  specify `F_{second_order}` for each prognostic state variable. can be functions of the prognostic state, gradient flux state, and auxiliary variables. |
# | [`source!`](@ref)            |  specify `S_{non_conservative}`  for each prognostic state variable. can be functions of the prognostic state, gradient flux state, and auxiliary variables. |

# ## Methods to compute numerical fluxes
# | **Method** | Purpose |
# |:-----|:-----|
# | [`wavespeed`](@ref) | specify how to compute the local wavespeed if using the `RusanovNumericalFlux`. |
# | [`boundary_state!`](@ref) | define exterior nodal values of the prognostic state and gradient flux state used to compute the numerical boundary fluxes. |

# ## Methods to set initial conditions
# | **Method** | Purpose |
# |:-----|:-----|
# | [`init_state_prognostic!`](@ref) | provide initial values for the prognostic state as a function of time and space. |
# | [`init_state_auxiliary!`](@ref) | provide initial values for the auxiliary variables as a function of the geometry. |

# ## General Remarks

# While `Y` can be thought of a column vector (each _row_ of which corresponds
# to each state variable and its prognostic equation), the _second_ function
# argument inside these methods behave as dictionaries, for example:

# ```julia
# struct MyModel <: BalanceLaw end

# function vars_state(m::MyModel, ::Prognostic, FT)
#     @vars begin
#         ρ::FT
#         T::FT
#     end
# end

# function source!(m::MyModel, source::Vars, args...)
#     source.ρ = 1 # adds a source of 1 to RHS of ρ equation
#     source.T = 1 # adds a source of 1 to RHS of T equation
# end
# ```

# All equations are marched simultaneously in time.

# # Reference links

# [^1]: [Column Vectors](https://en.wikipedia.org/wiki/Row_and_column_vectors)

# [^2]: Note that using non-conservative sources should be a final resort,
# as this can leak conservation of the unknowns and lead to numerical
# instabilities. It is recommended to use either `F_{diffusive}` or
# `F_{non_diffusive}`, as these fluxes are communicated across elements[^3]

# [^3]: MPI communication occurs only across elements, not within each element,
# where there may be many [Gauss-Lobatto][^4] points

# [^4]: [Gauss-Lobatto](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Lobatto_rules)

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
export rule

@rule NormalMixture((:m, k), Marginalisation) (q_out::Any, q_switch::Any, q_p::GammaDistributionsFamily) = begin
    z_bar = clamp.(probvec(q_switch), tiny, 1.0 - tiny)
    return NormalMeanVariance(mean(q_out), inv(z_bar[k] * mean(q_p)))
end

@rule NormalMixture((:m, k), Marginalisation) (q_out::Any, q_switch::Any, q_p::Wishart) = begin
    z_bar = clamp.(probvec(q_switch), tiny, 1.0 - tiny)
    return MvNormalMeanCovariance(mean(q_out), cholinv(z_bar[k] * mean(q_p)))
end


text:
#To Do: Make it possible to load from file

#Create AbstractSpacecraft type
abstract type AbstractSpacecraft end

#original Code definition of RigidBodyVehicle with rotational inertia
struct RigidBodyVehicle <: AbstractSpacecraft
  Name::String
  DryMass::Float64
  Fuel::Float64
  RefArea::Float64
  RefLengthLon::Float64
  RefLengthLat::Float64
  ReflecCoeff::Float64
  SRParea::Float64
  InertiaTensor::Matrix{Float64}
  InvInertiaTensor::Matrix{Float64}
  VehicleAeroDataBase::Type{<:AbstractVehicleDataBase}
  Actuators::Vector{Type{<:AbstractActuator}}
  MissionSystem::Type{<:AbstractMissionSystem}
  GuidanceSystem::Type{<:AbstractGuidanceSystem}
  ControlSystem::Type{<:AbstractControlSystem}
  NavigationSystem::Type{<:AbstractNavigationSystem}
  InitialState::SimState
  Trajectory::Vector{SimState}

  #Constructor for vehicle type
  function RigidBodyVehicle(;Name::String = "MySpacecraft", Drymass::Float64 = 1000.,
      Fuel::Float64 = 0., ReferenceArea::Float64 = 100., ReferenceLengthLongitude::Float64 = 23.,
      ReferenceLengthLatitude::Float64 = 12., ReflectionCoefficient::Float64 = 0.5,
      SRP_Area::Float64 = 20., InertiaTensor::Matrix{Float64} = diagm([10000., 10000., 10000.]),
      InverseInertiaTensor::Matrix{Float64} = inv(InertiaTensor),
      VehicleAeroDatabase::Type{<:AbstractVehicleDataBase} = VehicleAeroDataBase(),
      Actuators::Vector{Type{<:AbstractActuator}} = Vector{AbstractActuator}(),
      MissionSystem::Type{<:AbstractMissionSystem} = MissionSystem(),
      GuidanceSystem::Type{<:AbstractGuidanceSystem} = GuidanceSystem(),
      ControlSystem::Type{<:AbstractControlSystem} = ControlSystem(),
      NavigationSystem::Type{<:AbstractNavigationSystem} = NavigationSystem(),
      InitialState::SimState = SimState(Vector{Float64}(6), Vector{Float64}(3)),
      Trajectory::Vector{SimState} = Vector{SimState}())

      #Create Vehicle
      new(Name, Drymass, Fuel, ReferenceArea, ReferenceLengthLongitude, ReferenceLengthLatitude,
            ReflectionCoefficient, SRP_Area, InertiaTensor, InverseInertiaTensor,
            VehicleAeroDataBase, Actuators, MissionSystem, GuidanceSystem, ControlSystem,
            NavigationSystem, InitialState, Trajectory)
  end #End Constructor
end