text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module PyCall

export pycall, pyimport, pybuiltin, PyObject,
       PyPtr, pyincref, pydecref, pyversion, PyArray, PyArray_Info,
       pyerr_check, pyerr_clear, pytype_query, PyAny, @pyimport, PyDict,
       pyisinstance, pywrap, pytypeof, pyeval, PyVector, pystring,
       pyraise, pytype_mapping, pygui, pygui_start, pygui_stop,
       pygui_stop_all, @pylab, set!, PyTextIO, @pysym

import Base: size, ndims, similar, copy, getindex, setindex!, stride,
       convert, pointer, summary, convert, show, haskey, keys, values,
       eltype, get, delete!, empty!, length, isempty, start, done,
       next, filter!, hash, splice!, pop!, ==, isequal, push!,
       unshift!, shift!, append!, insert!, prepend!, writemime, mimewritable

# Python C API is not interrupt-save.  In principle, we should
# use sigatomic for every ccall to the Python library, but this
# should really be fixed in Julia (#2622).  However, we will
# use the sigatomic_begin/end functions to protect pycall and
# similar long-running (or potentially long-running) code.
import Base: sigatomic_begin, sigatomic_end

## Compatibility import for v0.3, v0.4
using Compat
if VERSION >= v"0.4.0-dev+3710"
    import Base.unsafe_convert
else
    const unsafe_convert = Base.convert
end

#########################################################################

const depfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(depfile) || error("PyCall not properly installed. Please run Pkg.build(\"PyCall\")")
include(depfile) # generated by Pkg.build("PyCall")

macro pysym(func)
    :(($func, libpython))
end
macro pyglobal(name)
    :(cglobal(($name, libpython)))
end
macro pyglobalobj(name)
    :(cglobal(($name, libpython), PyObject_struct))
end

#########################################################################

# Mirror of C PyObject struct (for non-debugging Python builds).
# We won't actually access these fields directly; we'll use the Python
# C API for everything.  However, we need to define a unique Ptr type
# for PyObject*, and we might as well define the actual struct layout
# while we're at it.
immutable PyObject_struct
    ob_refcnt::Int
    ob_type::Ptr{Void}
end

typealias PyPtr Ptr{PyObject_struct} # type for PythonObject* in ccall

#########################################################################
# Wrapper around Python's C PyObject* type, with hooks to Python reference
# counting and conversion routines to/from C and Julia types.

type PyObject
    o::PyPtr # the actual PyObject*
    function PyObject(o::PyPtr)
        po = new(o)
        finalizer(po, pydecref)
        return po
    end
end

PyNULL() = PyObject(convert(PyPtr, C_NULL))

function pydecref(o::PyObject)
    ccall(@pysym(:Py_DecRef), Void, (PyPtr,), o.o)
    o.o = convert(PyPtr, C_NULL)
    o
end

function pyincref(o::PyObject)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    o
end

# doing an incref *before* creating a PyObject may safer in the
# case of borrowed references, to ensure that no exception or interrupt
# induces a double decref.
function pyincref(o::PyPtr)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    PyObject(o)
end

pyisinstance(o::PyObject, t::PyObject) =
  t.o != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t.o) == 1

pyisinstance(o::PyObject, t::Union(Ptr{Void},PyPtr)) =
  t != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1

pyquery(q::Ptr{Void}, o::PyObject) =
  ccall(q, Cint, (PyPtr,), o) == 1

pytypeof(o::PyObject) = o.o == C_NULL ? throw(ArgumentError("NULL PyObjects have no Python type")) : pycall(TypeType, PyObject, o)

# conversion to pass PyObject as ccall arguments:
unsafe_convert(::Type{PyPtr}, po::PyObject) = po.o

# use constructor for generic conversions to PyObject
convert(::Type{PyObject}, o) = PyObject(o)
PyObject(o::PyObject) = o

#########################################################################

include("pyinit.jl")
include("exception.jl")
include("gui.jl")

#########################################################################

include("gc.jl")

# make a PyObject that embeds a reference to keep, to prevent Julia
# from garbage-collecting keep until o is finalized.
PyObject(o::PyPtr, keep::Any) = pyembed(PyObject(o), keep)

#########################################################################

include("pybuffer.jl")
include("conversions.jl")
include("pytype.jl")
include("callback.jl")
include("io.jl")

#########################################################################
# Pretty-printing PyObject

function pystring(o::PyObject)
    if o.o == C_NULL
        return "NULL"
    else
        s = ccall((@pysym :PyObject_Repr), PyPtr, (PyPtr,), o)
        if (s == C_NULL)
            pyerr_clear()
            s = ccall((@pysym :PyObject_Str), PyPtr, (PyPtr,), o)
            if (s == C_NULL)
                pyerr_clear()
                return string(o.o)
            end
        end
        return convert(AbstractString, PyObject(s))
    end
end

function show(io::IO, o::PyObject)
    print(io, "PyObject $(pystring(o))")
end

if VERSION < v"0.4.0-dev+6471"
    function Base.help(o::PyObject)
        pycall(pybuiltin("help"), PyAny, o)
    end
else
    function Base.Docs.doc(o::PyObject)
        Base.Docs.Text(haskey(o, "__doc__") ?
                       convert(AbstractString, o["__doc__"]) :
                       "Python object (no docstring found)")
    end
end

#########################################################################
# computing hashes of PyObjects

const pysalt = hash("PyCall.PyObject") # "salt" to mix in to PyObject hashes
hashsalt(x) = hash(x, pysalt)

function hash(o::PyObject)
    if o.o == C_NULL
        hashsalt(C_NULL)
    elseif is_pyjlwrap(o)
        # call native Julia hash directly on wrapped Julia objects,
        # since on 64-bit Windows the Python 2.x hash is only 32 bits
        hashsalt(unsafe_pyjlwrap_to_objref(o.o))
    else
        h = ccall((@pysym :PyObject_Hash), Py_hash_t, (PyPtr,), o)
        if h == -1 # error
            pyerr_clear()
            return hashsalt(o.o)
        end
        hashsalt(h)
    end
end

#########################################################################
# PyObject equality

const Py_EQ = convert(Cint, 2) # from Python's object.h

function ==(o1::PyObject, o2::PyObject)
    if o1.o == C_NULL || o2.o == C_NULL
        return o1.o == o2.o
    elseif is_pyjlwrap(o1)
        if is_pyjlwrap(o2)
            return unsafe_pyjlwrap_to_objref(o1.o) ==
                   unsafe_pyjlwrap_to_objref(o2.o)
        else
            return false
        end
    else
        val = ccall((@pysym :PyObject_RichCompareBool), Cint,
                    (PyPtr, PyPtr, Cint), o1, o2, Py_EQ)
        return val == -1 ? o1.o == o2.o : @compat Bool(val)
    end
end

isequal(o1::PyObject, o2::PyObject) = o1 == o2 # Julia 0.2 compatibility

#########################################################################
# For o::PyObject, make o["foo"] and o[:foo] equivalent to o.foo in Python,
# with the former returning an raw PyObject and the latter giving the PyAny
# conversion.

function getindex(o::PyObject, s::AbstractString)
    if (o.o == C_NULL)
        throw(ArgumentError("ref of NULL PyObject"))
    end
    p = ccall((@pysym :PyObject_GetAttrString), PyPtr,
              (PyPtr, Ptr{Uint8}), o, bytestring(s))
    if p == C_NULL
        pyerr_clear()
        throw(KeyError(s))
    end
    return PyObject(p)
end

getindex(o::PyObject, s::Symbol) = convert(PyAny, getindex(o, string(s)))

function setindex!(o::PyObject, v, s::AbstractString)
    if (o.o == C_NULL)
        throw(ArgumentError("assign of NULL PyObject"))
    end
    if -1 == ccall((@pysym :PyObject_SetAttrString), Cint,
                   (PyPtr, Ptr{Uint8}, PyPtr), o, bytestring(s), PyObject(v))
        pyerr_clear()
        throw(KeyError(s))
    end
    o
end

setindex!(o::PyObject, v, s::Symbol) = setindex!(o, v, string(s))

function haskey(o::PyObject, s::AbstractString)
    if (o.o == C_NULL)
        throw(ArgumentError("haskey of NULL PyObject"))
    end
    return 1 == ccall((@pysym :PyObject_HasAttrString), Cint,
                      (PyPtr, Ptr{Uint8}), o, bytestring(s))
end

haskey(o::PyObject, s::Symbol) = haskey(o, string(s))

#########################################################################
# keys(o) should return an iterator over the keys (members) of o

type PyObjectMembers
    members::PyVector{@compat Tuple{Symbol,PyObject}}
end

function show(io::IO, m::PyObjectMembers)
    print(io, "PyObjectMembers iterator (", length(m), " members)")
end

keys(o::PyObject) = PyObjectMembers(pycall(inspect["getmembers"],
                                           PyVector{(Symbol,PyObject)}, o))

start(m::PyObjectMembers) = 1
done(m::PyObjectMembers, i::Int) = i > length(m.members)
next(m::PyObjectMembers, i::Int) = (m.members[i][1], i+1)
length(m::PyObjectMembers) = length(m.members)

#########################################################################
# Create anonymous composite w = pywrap(o) wrapping the object o
# and providing access to o's members (converted to PyAny) as w.member.

# we skip wrapping Julia reserved words (which cannot be type members)
const reserved = Set{ASCIIString}(["while", "if", "for", "try", "return", "break", "continue", "function", "macro", "quote", "let", "local", "global", "const", "abstract", "typealias", "type", "bitstype", "immutable", "ccall", "do", "module", "baremodule", "using", "import", "export", "importall", "pymember", "false", "true", "Tuple"])

function pywrap(o::PyObject, mname::Symbol=:__anon__)
    members = convert(Vector{@compat Tuple{AbstractString,PyObject}},
                      pycall(inspect["getmembers"], PyObject, o))
    filter!(m -> !(m[1] in reserved), members)
    if VERSION >= v"0.4.0-dev+4448"
        m = Module(mname, false)
    else
        # Hack to create an anonymous bare module on older Julia
        m = Module(:pyimport)
        m = eval(m, Expr(:toplevel, Expr(:module, false, mname, Expr(:block)), mname))
    end
    consts = [Expr(:const, Expr(:(=), symbol(x[1]), convert(PyAny, x[2]))) for x in members]
    exports = try
                  convert(Vector{Symbol}, o["__all__"])
              catch
                  [symbol(x[1]) for x in filter(x -> x[1][1] != '_', members)]
              end
    eval(m, Expr(:toplevel, consts..., :(pymember(s) = $(getindex)($(o), s)),
                 Expr(:export, exports...)))
    m
end

#########################################################################

pyimport(name::AbstractString) =
    PyObject(@pycheckn ccall((@pysym :PyImport_ImportModule), PyPtr,
                             (Ptr{Uint8},), bytestring(name)))

pyimport(name::Symbol) = pyimport(string(name))

# convert expressions like :math or :(scipy.special) into module name strings
modulename(s::Symbol) = string(s)
function modulename(e::Expr)
    if e.head == :.
        string(modulename(e.args[1]), :., modulename(e.args[2]))
    elseif e.head == :quote
        modulename(e.args...)
    else
        throw(ArgumentError("invalid module"))
    end
end

# separate this function in order to make it easier to write more
# pyimport-like functions
function pyimport_name(name, optional_varname)
    len = length(optional_varname)
    if len > 0 && (len != 2 || optional_varname[1] != :as)
        throw(ArgumentError("usage @pyimport module [as name]"))
    elseif len == 2
        optional_varname[2]
    elseif typeof(name) == Symbol
        name
    else
        mname = modulename(name)
        throw(ArgumentError("$mname is not a valid module variable name, use @pyimport $mname as <name>"))
    end
end

macro pyimport(name, optional_varname...)
    mname = modulename(name)
    Name = pyimport_name(name, optional_varname)
    quote
        if !isdefined($(Expr(:quote, Name)))
            const $(esc(Name)) = pywrap(pyimport($mname))
        elseif !isa($(esc(Name)), Module)
            error("@pyimport: ", $(Expr(:quote, Name)), " already defined")
        end
        nothing
    end
end

#########################################################################

# look up a global builtin
function pybuiltin(name)
    builtin[name]
end

#########################################################################

typealias TypeTuple{N} Union(Type,NTuple{N, Type})

function pycall(o::Union(PyObject,PyPtr), returntype::TypeTuple, args...; kwargs...)
    oargs = map(PyObject, args)
    nargs = length(args)
    sigatomic_begin()
    try
        arg = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,),
                                       nargs))
        for i = 1:nargs
            @pycheckz ccall((@pysym :PyTuple_SetItem), Cint,
                             (PyPtr,Int,PyPtr), arg, i-1, oargs[i])
            pyincref(oargs[i]) # PyTuple_SetItem steals the reference
        end
        if isempty(kwargs)
            ret = PyObject(@pycheckn ccall((@pysym :PyObject_Call), PyPtr,
                                          (PyPtr,PyPtr,PyPtr), o, arg, C_NULL))
        else
            kw = PyObject((AbstractString=>Any)[string(k) => v for (k, v) in kwargs])
            ret = PyObject(@pycheckn ccall((@pysym :PyObject_Call), PyPtr,
                                            (PyPtr,PyPtr,PyPtr), o, arg, kw))
        end
        jret = convert(returntype, ret)
        return jret
    finally
        sigatomic_end()
    end
end

#########################################################################
# Once Julia lets us overload ".", we will use [] to access items, but
# for now we can define "get"

function get(o::PyObject, returntype::TypeTuple, k, default)
    r = ccall((@pysym :PyObject_GetItem), PyPtr, (PyPtr,PyPtr), o,PyObject(k))
    if r == C_NULL
        pyerr_clear()
        default
    else
        convert(returntype, PyObject(r))
    end
end

get(o::PyObject, returntype::TypeTuple, k) =
    convert(returntype, PyObject(@pycheckn ccall((@pysym :PyObject_GetItem),
                                 PyPtr, (PyPtr,PyPtr), o, PyObject(k))))

get(o::PyObject, k, default) = get(o, PyAny, k, default)
get(o::PyObject, k) = get(o, PyAny, k)

function delete!(o::PyObject, k)
    e = ccall((@pysym :PyObject_DelItem), Cint, (PyPtr, PyPtr), o, PyObject(k))
    if e == -1
        pyerr_clear() # delete! ignores errors in Julia
    end
    return o
end

function set!(o::PyObject, k, v)
    @pycheckz ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     o, PyObject(k), PyObject(v))
    v
end

#########################################################################
# support IPython _repr_foo functions for writemime of PyObjects

for (mime, method) in ((MIME"text/html", "_repr_html_"),
                       (MIME"image/jpeg", "_repr_jpeg_"),
                       (MIME"image/png", "_repr_png_"),
                       (MIME"image/svg+xml", "_repr_svg_"),
                       (MIME"text/latex", "_repr_latex_"))
    T = istext(mime()) ? AbstractString : Vector{Uint8}
    @eval begin
        function writemime(io::IO, mime::$mime, o::PyObject)
            if o.o != C_NULL && haskey(o, $method)
                r = pycall(o[$method], PyObject)
                r.o != pynothing && return write(io, convert($T, r))
            end
            throw(MethodError(writemime, (io, mime, o)))
        end
        mimewritable(::$mime, o::PyObject) =
            o.o != C_NULL && haskey(o, $method) && let meth = o[$method]
                meth.o != pynothing &&
                pycall(meth, PyObject).o != pynothing
            end
    end
end

#########################################################################

const Py_single_input = 256  # from Python.h
const Py_file_input = 257
const Py_eval_input = 258

const pyeval_fname = bytestring("PyCall.jl") # filename for pyeval

# evaluate a python string, returning PyObject, given a dictionary
# (string/symbol => value) of local variables to use in the expression
function pyeval_(s::AbstractString, locals::PyDict, input_type)
    sb = bytestring(s) # use temp var to prevent gc before we are done with o
    sigatomic_begin()
    try
        o = PyObject(@pycheckn ccall((@pysym :Py_CompileString), PyPtr,
                                     (Ptr{Uint8}, Ptr{Uint8}, Cint),
                                     sb, pyeval_fname, input_type))
        main = @pycheckn ccall((@pysym :PyImport_AddModule),
                                PyPtr, (Ptr{Uint8},),
                                bytestring("__main__"))
        maindict = @pycheckn ccall((@pysym :PyModule_GetDict), PyPtr,
                                    (PyPtr,), main)
        return PyObject(@pycheckn ccall((@pysym :PyEval_EvalCode),
                                         PyPtr, (PyPtr, PyPtr, PyPtr),
                                         o, maindict, locals))
    finally
        sigatomic_end()
    end
end

function pyeval(s::AbstractString, returntype::TypeTuple=PyAny,
                locals=PyDict{AbstractString, PyObject}(),
                input_type=Py_eval_input; kwargs...)
    for (k, v) in kwargs
        locals[string(k)] = v
    end
    return convert(returntype, pyeval_(s, locals, input_type))
end

#########################################################################
# Precompilation: just an optimization to speed up initialization.
# Here, we precompile functions that are passed to cfunction by __init__,
# for the reasons described in JuliaLang/julia#12256.

precompile(jl_Function_call, (PyPtr,PyPtr,PyPtr))
precompile(pyio_repr, (PyPtr,))
precompile(pyjlwrap_dealloc, (PyPtr,))
precompile(pyjlwrap_repr, (PyPtr,))
precompile(pyjlwrap_hash, (PyPtr,))
precompile(pyjlwrap_hash32, (PyPtr,))

for f in (jl_IO_close, jl_IO_fileno, jl_IO_flush, jl_IO_isatty, jl_IO_readable, jl_IO_writable, jl_IO_readline, jl_IO_readlines, jl_IO_seek, jl_IO_seekable, jl_IO_tell, jl_IO_writelines, jl_IO_read_text, jl_IO_read, jl_IO_readall_text, jl_IO_readall, jl_IO_readinto, jl_IO_write)
    precompile(f, (PyPtr,PyPtr))
end

end # module PyCall


text:
"""
$(README)
Depends on:
$(IMPORTS)
"""
module Vulkan

using Reexport
using DocStringExtensions
using AutoHashEquals: @auto_hash_equals

using VulkanCore: vk
using .vk
const core = vk
export vk, core

using Base: RefArray
import Base: convert, cconvert, unsafe_convert, &, |, xor, isless, ==, typemax, in, parent
using MLStyle

@reexport using ResultTypes: unwrap, unwrap_error, iserror
using ResultTypes: ResultTypes

@template FUNCTIONS =
    """
    $(DOCSTRING)
    $(METHODLIST)
    """

@template (METHODS, MACROS) =
    """
    $(DOCSTRING)
    $(TYPEDSIGNATURES)
    """

@template TYPES =
    """
    $(DOCSTRING)
    $(TYPEDEF)
    $(TYPEDFIELDS)
    """

include("preferences.jl")

# generated wrapper
include("prewrap.jl")

include("CEnum/CEnum.jl")
using .CEnum

@static if Sys.islinux()
    include("../generated/linux.jl")
elseif Sys.isapple()
    include("../generated/macos.jl")
elseif Sys.isbsd()
    include("../generated/bsd.jl")
elseif Sys.iswindows()
    include("../generated/windows.jl")
end

include("utils.jl")
include("driver.jl")
include("validation.jl")
include("instance.jl")
include("device.jl")
include("dispatch.jl")
include("print.jl")

const global_dispatcher = APIDispatcher()

export
        # Wrapper
        VulkanStruct,
        ReturnedOnly,
        Handle,
        to_vk,
        from_vk,

        # Driver
        set_driver,
        @set_driver,

        # Printing
        print_app_info,
        print_available_devices,
        print_device_info,
        print_devices,

        # Device
        find_queue_family,

        # Debugging
        default_debug_callback,

        # Pointer utilities
        function_pointer,
        pointer_length,

        # Bitmask manipulation utilities
        BitMask,
        @bitmask_flag,

        # Error handling
        VulkanError,
        @check,
        iserror

end # module Vulkan

ext: jl


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler multicomponent equations
gamma1              = 1.4
gamma2              = 1.648
gas_constant1       = 0.287
gas_constant2       = 1.578
equations           = CompressibleEulerMulticomponentEquations2D(gamma1, gamma2, gas_constant1, gas_constant2)

initial_condition   = initial_condition_shock_bubble

surface_flux        = flux_lax_friedrichs
volume_flux         = flux_chandrashekar
basis               = LobattoLegendreBasis(3)
indicator_sc        = IndicatorHennemannGassner(equations, basis,
                                                alpha_max=0.5,
                                                alpha_min=0.001,
                                                alpha_smooth=true,
                                                variable=density_pressure)
volume_integral     = VolumeIntegralShockCapturingHG(indicator_sc;
                                                     volume_flux_dg=volume_flux,
                                                     volume_flux_fv=surface_flux)
solver              = DGSEM(basis, surface_flux, volume_integral)

coordinates_min     = (-2.25, -2.225)
coordinates_max     = ( 2.20,  2.225)
mesh                = TreeMesh(coordinates_min, coordinates_max,
                               initial_refinement_level=7,
                               n_cells_max=1_000_000)

semi                = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan               = (0.0, 0.012) # 0.0021..
ode                 = semidiscretize(semi, tspan)

summary_callback    = SummaryCallback()

analysis_interval   = 300 # change to 200+ for whole calc.
analysis_callback   = AnalysisCallback(semi, interval=analysis_interval)

alive_callback      = AliveCallback(analysis_interval=analysis_interval)

save_solution       = SaveSolutionCallback(interval=300,     # 40 or change to 200+ for whole calc.
                                           save_initial_solution=true,
                                           save_final_solution=true,
                                           solution_variables=:primitive)

stepsize_callback   = StepsizeCallback(cfl=0.3)

callbacks           = CallbackSet(summary_callback,
                                  analysis_callback,
                                  alive_callback,
                                  save_solution,
                                  stepsize_callback)

###############################################################################
# run the simulation
sol                 = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
                                                      dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                                                      save_everystep=false,
                                                      callback=callbacks,
                                                      maxiters=1e5);
summary_callback() # print the timer summary


text:
"""
Apply the rule for each cell in the grid, using optimisations
allowed for the supertype of the rule.
"""
maprule!(data::AbstractSimData, rule::Rule) = blockrun!(data, rule)

blockrun!(data, context, args...) = begin
    nrows, ncols = framesize(data)
    r = radius(data)
    if r > 0
        blocksize = 2r
        status = sourcestatus(data)

        @inbounds for bj in 1:size(status, 2) - 1, bi in 1:size(status, 1) - 1
            status[bi, bj] || continue
            # Convert from padded block to init dimensions
            istart = blocktoind(bi, blocksize) - r
            jstart = blocktoind(bj, blocksize) - r
            # Stop at the init row/column size, not the padding or block multiple
            istop = min(istart + blocksize - 1, nrows)
            jstop = min(jstart + blocksize - 1, ncols)
            # Skip the padding
            istart = max(istart, 1)
            jstart = max(jstart, 1)

            for j in jstart:jstop, i in istart:istop
                ismasked(data, i, j) && continue
                blockdo!(data, context, (i, j), args...)
            end
        end
    else
        for j in 1:ncols, i in 1:nrows
            ismasked(data, i, j) && continue
            blockdo!(data, context, (i, j), args...)
        end
    end
end

@inline blockdo!(data, rule::Rule, I) = begin
    @inbounds state = source(data)[I...]
    @inbounds dest(data)[I...] = applyrule(rule, data, state, I)
    nothing
end

maprule!(data::AbstractSimData, rule::PartialRule) = begin
    data = WritableSimData(data)
    # Update active blocks in the dest array
    @inbounds parent(dest(data)) .= parent(source(data))
    # Run the rule for active blocks
    blockrun!(data, rule)
    updatestatus!(sourcestatus(data), deststatus(data))
end

@inline blockdo!(data::WritableSimData, rule::PartialRule, I) = begin
    state = source(data)[I...]
    state == zero(state) && return
    applyrule!(rule, data, state, I)
end

# struct UpdateDest end
# @inline blockdo!(data::WritableSimData, ::UpdateDest, I) = begin
    # out = source(data)[I...]
    # dest(data)[I...] = out
# end

"""
Run the rule for all cells, writing the result to the dest array
The neighborhood is copied to the rules neighborhood buffer array for performance
# TODO test 1d
"""
# maprule!(data::AbstractSimData{T,1}, rule::Union{NeighborhoodRule, Tuple{NeighborhoodRule,Vararg}},
#           args...)  where T = begin
#     # The rule provides the neighborfdood buffer
#     r = radius(data)
#     sze = hoodsize(r)
#     buf = similar(init(data), sze, sze)
#     src, dst = source(data), dest(data)
#     nrows = framesize(data)
#     handleoverflow!(data, r)
#     # Setup buffer array between rows
#     # Ignore the first column, it wil be copied over in the main loop
#     for i in 2:sze
#         @inbounds buf[i] = src[i-1-r]
#     end
#     # Run rule for a row
#     @inbounds for i in 1:nrows
#         copyto!(buf, 1, buf, 2)
#         buf[sze] = src[i+r]
#         state = buf[r+1]
#         dst[i] = applyrule(rule, data, state, (i,))
#     end
# end

maprule!(data::SingleSimData{T,2}, rule::Union{NeighborhoodRule,Chain{<:Tuple{NeighborhoodRule,Vararg}}},
         args...) where T = begin
    # The rule provides the neighborhood buffer
    r = radius(rule)
    blocksize = 2r
    hoodsize = 2r + 1
    src, dst = parent(source(data)), parent(dest(data))
    srcstatus, dststatus = sourcestatus(data), deststatus(data)
    sumstatus = localstatus(data)
    bufs = buffers(data)
    nrows, ncols = framesize(data)

    handleoverflow!(data, r)
    # Center of the buffer for both axes
    center = r + 1
    deststatus(data) .= false

    # Run the rule row by row. When we move along a row by one cell, we access only
    # a single new column of data the same hight of the nighborhood, and move the existing
    # data in the neighborhood buffer array accross by one column. This saves on reads
    # from the main array, and focusses reads and writes in the small buffer array that
    # should be in fast local memory.
    @inbounds for bi = 1:size(srcstatus, 1) - 1
        sumstatus .= false
        i = blocktoind(bi, blocksize)
        rowsinblock = min(blocksize, nrows - blocksize * (bi - 1))
        skippedlastblock = true
        freshbuffer = true

        b11, b12 = srcstatus[bi,     1], srcstatus[bi,     2]
        b21, b22 = srcstatus[bi + 1, 1], srcstatus[bi + 1, 2]
        sumstatus[1, 2] = false
        sumstatus[2, 2] = false
        for bj = 1:size(srcstatus, 2) - 1
            sumstatus[1, 1] = sumstatus[1, 2]
            sumstatus[2, 1] = sumstatus[2, 2]
            sumstatus[1, 2] = false
            sumstatus[2, 2] = false

            b11, b21 = b12, b22
            b12, b22 = srcstatus[bi, bj + 1], srcstatus[bi + 1, bj + 1]

            jstart = blocktoind(bj, blocksize)
            jstop = min(jstart + blocksize - 1, ncols)

            # Use this block unless it or its neighbors are active
            if !(b11 | b12 | b21 | b22)
                # Skip this block
                skippedlastblock = true
                continue
            end

            if skippedlastblock
                for y = 1:hoodsize
                    for b in 1:rowsinblock
                        for x = 1:hoodsize
                            val = src[i + b + x - 2, jstart + y - 1]
                            bufs[b][x, y] = val
                        end
                    end
                end
                skippedlastblock = false
                freshbuffer = true
            end

            for j in jstart:jstop
                centerbj = j - jstart < r ? 1 : 2
                if freshbuffer
                    freshbuffer = false
                else
                    # Move the neighborhood buffers accross one column
                    for b in 1:rowsinblock
                        buf = bufs[b]
                        # copyto! uses linear indexing, so 2d dims are transformed manually
                        copyto!(buf, 1, buf, hoodsize + 1, (hoodsize - 1) * hoodsize)
                    end
                    # Copy a new column to each neighborhood buffer
                    for b in 1:rowsinblock
                        buf = bufs[b]
                        for x in 1:hoodsize
                            buf[x, hoodsize] = src[i + b + x - 2, j + 2r]
                        end
                    end
                end
                for b in 1:rowsinblock
                    ii = i + b - 1
                    ismasked(data, ii, j) && continue

                    centerbi = b <= r ? 1 : 2
                    # Run the rule using buffer b
                    buf = bufs[b]
                    state = buf[center, center]
                    # @assert state == src[ii + r, j + r]
                    newstate = applyrule(rule, data, state, (ii, j), buf)
                    sumstatus[centerbi, centerbj] |= newstate != zero(newstate)
                    dst[ii + r, j + r] = newstate
                end
                # Combine blocks with the previous rows / cols
                # TODO only write the first column
                dststatus[bi, bj] |= sumstatus[1, 1]
                dststatus[bi, bj+1] |= sumstatus[1, 2]
                dststatus[bi+1, bj] |= sumstatus[2, 1]
                # Start new block fresh to remove old status
                dststatus[bi+1, bj+1] = sumstatus[2, 2]
            end
        end
    end
    updatestatus!(sourcestatus(data), deststatus(data))
end

"""
Wrap overflow where required. This optimisation allows us to ignore
bounds checks on neighborhoods and still use a wraparound grid.
"""
handleoverflow!(data::SingleSimData, r::Integer) = handleoverflow!(data, overflow(data), r)
handleoverflow!(data::SingleSimData{T,1}, overflow::WrapOverflow, r::Integer) where T = begin
    # Copy two sides
    @inbounds copyto!(source, 1-r:0, source, nrows+1-r:nrows)
    @inbounds copyto!(source, nrows+1:nrows+r, source, 1:r)
end
handleoverflow!(data::SingleSimData{T,2}, overflow::WrapOverflow, r::Integer) where T = begin
    # TODO optimise this. Its mostly a placeholder so wrapping still works in GOL tests.
    src = source(data)
    nrows, ncols = framesize(data)
    # Left
    @inbounds copyto!(src, CartesianIndices((1:nrows, 1-r:0)),
                      src, CartesianIndices((1:nrows, ncols+1-r:ncols)))
    # Right
    @inbounds copyto!(src, CartesianIndices((1:nrows, ncols+1:ncols+r)),
                      src, CartesianIndices((1:nrows, 1:r)))
    # Top
    @inbounds copyto!(src, CartesianIndices((1-r:0, 1:ncols)),
                      src, CartesianIndices((ncols+1-r:ncols, 1:ncols)))
    # Bottom
    @inbounds copyto!(src, CartesianIndices((ncols+1:ncols+r, 1:ncols)),
                      src, CartesianIndices((1:r, 1:ncols)))

    # Copy four corners
    @inbounds copyto!(src, CartesianIndices((1-r:0, 1-r:0)),
                      src, CartesianIndices((nrows+1-r:nrows, ncols+1-r:ncols)))
    @inbounds copyto!(src, CartesianIndices((1-r:0, ncols+1:ncols+r)),
                      src, CartesianIndices((nrows+1-r:nrows, 1:r)))
    @inbounds copyto!(src, CartesianIndices((nrows+1:nrows+r, ncols+1:ncols+r)),
                      src, CartesianIndices((1:r, 1:r)))
    @inbounds copyto!(src, CartesianIndices((nrows+1:nrows+r, 1-r:0)),
                      src, CartesianIndices((1:r, ncols+1-r:ncols)))

    # Wrap status
    status = sourcestatus(data)
    # status[:, 1] .|= status[:, end-1] .| status[:, end-2]
    # status[1, :] .|= status[end-1, :] .| status[end-2, :]
    # status[end-1, :] .|= status[1, :]
    # status[:, end-1] .|= status[:, 1]
    # status[end-2, :] .|= status[1, :]
    # status[:, end-2] .|= status[:, 1]
    # FIXME: Buggy currently, just running all in Wrap mode
    status .= true
end
handleoverflow!(data, overflow::RemoveOverflow, r) = nothing

combinestatus(x::Number, y::Number) = x + y
combinestatus(x::Integer, y::Integer) = x | y

updatestatus!(copyto::AbstractArray, copyfrom::AbstractArray) =
    @inbounds copyto .= copyfrom
updatestatus!(copyto, copyfrom) = nothing

"""
    applyrule(rules::Chain, data, state, (i, j))

Chained rules. If a `Chain` of rules is passed to applyrule, run them sequentially for each
cell.  This can have much beter performance as no writes occur between rules, and they are
essentially compiled together into compound rules. This gives correct results only for
CellRule, or for a single NeighborhoodRule followed by CellRule.
"""
@inline applyrule(rules::Chain{<:Tuple{<:NeighborhoodRule,Vararg}}, data, state, index, buf) = begin
    state = applyrule(rules[1], data, state, index, buf)
    applyrule(tail(rules), data, state, index)
end
@inline applyrule(rules::Chain, data, state, index) = begin
    state == zero(state) && return state
    newstate = applyrule(rules[1], data, state, index)
    applyrule(tail(rules), data, newstate, index)
end
@inline applyrule(rules::Chain{Tuple{}}, data, state, index) = state


text:
import Base.show
export show

# A Constant is simply a wrapper around a native Julia constant
# Hence, we simply display its value
function show(io::IO, x::Constant)
  print(io, "Constant($(x.value))")
end

# A variable and affine expression will display its size
function show(io::IO, x::AffineExpr)
  if x.head == :variable
    print(io, "Variable($(x.size[1]), $(x.size[2]))")
  else
    print(io, "Affine Expression($(x.size[1]), $(x.size[2]))")
  end
end

# A sum of squares expression simply prints its name
function show(io::IO, x::SumSquaresExpr)
  print(io, "Sum of Squares Expression")
end

# A linear equality constraint will print
function show(io::IO, c::EqConstraint)
  print(io, "Equality Constraint: $(c.lhs) == $(c.rhs)")
end


text:
const DEFAULT_SUBSAMPLE_NUMBER = 5 # default number of samples used in sub-pixel smoothing (total number is square of this)

################################################################################
const TRANSVERSE_MODE_RESOLUTION_MULTIPLIER = 3
const PROPAGATION_CONSTANT_IMAGINARY_TOLERANCE = 1e-2
const ANALYSIS_QUADRATURE_POINTS = 31
const PROGRESS_UPDATE_TIME = .05 # in seconds


text:
function Gradient(arg0::DerivativeStructure)
    return Gradient((DerivativeStructure,), arg0)
end

function Gradient(arg0::jdouble, arg1::Vector{jdouble})
    return Gradient((jdouble, Vector{jdouble}), arg0, arg1)
end

function abs(obj::Gradient)
    return jcall(obj, "abs", Gradient, ())
end

function acos(obj::Gradient)
    return jcall(obj, "acos", Gradient, ())
end

function acosh(obj::Gradient)
    return jcall(obj, "acosh", Gradient, ())
end

function add(obj::Gradient, arg0::Gradient)
    return jcall(obj, "add", Gradient, (Gradient,), arg0)
end

function add(obj::Gradient, arg0::jdouble)
    return jcall(obj, "add", Gradient, (jdouble,), arg0)
end

function asin(obj::Gradient)
    return jcall(obj, "asin", Gradient, ())
end

function asinh(obj::Gradient)
    return jcall(obj, "asinh", Gradient, ())
end

function atan(obj::Gradient)
    return jcall(obj, "atan", Gradient, ())
end

function atan2(obj::Gradient, arg0::Gradient)
    return jcall(obj, "atan2", Gradient, (Gradient,), arg0)
end

function atanh(obj::Gradient)
    return jcall(obj, "atanh", Gradient, ())
end

function cbrt(obj::Gradient)
    return jcall(obj, "cbrt", Gradient, ())
end

function ceil(obj::Gradient)
    return jcall(obj, "ceil", Gradient, ())
end

function compose(obj::Gradient, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "compose", Gradient, (jdouble, jdouble), arg0, arg1)
end

function constant(::Type{Gradient}, arg0::jint, arg1::jdouble)
    return jcall(Gradient, "constant", Gradient, (jint, jdouble), arg0, arg1)
end

function copy_sign(obj::Gradient, arg0::Gradient)
    return jcall(obj, "copySign", Gradient, (Gradient,), arg0)
end

function copy_sign(obj::Gradient, arg0::jdouble)
    return jcall(obj, "copySign", Gradient, (jdouble,), arg0)
end

function cos(obj::Gradient)
    return jcall(obj, "cos", Gradient, ())
end

function cosh(obj::Gradient)
    return jcall(obj, "cosh", Gradient, ())
end

function divide(obj::Gradient, arg0::Gradient)
    return jcall(obj, "divide", Gradient, (Gradient,), arg0)
end

function divide(obj::Gradient, arg0::jdouble)
    return jcall(obj, "divide", Gradient, (jdouble,), arg0)
end

function equals(obj::Gradient, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::Gradient)
    return jcall(obj, "exp", Gradient, ())
end

function expm1(obj::Gradient)
    return jcall(obj, "expm1", Gradient, ())
end

function floor(obj::Gradient)
    return jcall(obj, "floor", Gradient, ())
end

function get_exponent(obj::Gradient)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::Gradient)
    return jcall(obj, "getField", GradientField, ())
end

function get_free_parameters(obj::Gradient)
    return jcall(obj, "getFreeParameters", jint, ())
end

function get_gradient(obj::Gradient)
    return jcall(obj, "getGradient", Vector{jdouble}, ())
end

function get_order(obj::Gradient)
    return jcall(obj, "getOrder", jint, ())
end

function get_partial_derivative(obj::Gradient, arg0::Vector{jint})
    return jcall(obj, "getPartialDerivative", jdouble, (Vector{jint},), arg0)
end

function get_partial_derivative(obj::Gradient, arg0::jint)
    return jcall(obj, "getPartialDerivative", jdouble, (jint,), arg0)
end

function get_real(obj::Gradient)
    return jcall(obj, "getReal", jdouble, ())
end

function get_value(obj::Gradient)
    return jcall(obj, "getValue", jdouble, ())
end

function hash_code(obj::Gradient)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::Gradient, arg0::Gradient)
    return jcall(obj, "hypot", Gradient, (Gradient,), arg0)
end

function linear_combination(obj::Gradient, arg0::Gradient, arg1::Gradient, arg2::Gradient, arg3::Gradient)
    return jcall(obj, "linearCombination", Gradient, (Gradient, Gradient, Gradient, Gradient), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Gradient, arg0::Gradient, arg1::Gradient, arg2::Gradient, arg3::Gradient, arg4::Gradient, arg5::Gradient)
    return jcall(obj, "linearCombination", Gradient, (Gradient, Gradient, Gradient, Gradient, Gradient, Gradient), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Gradient, arg0::Gradient, arg1::Gradient, arg2::Gradient, arg3::Gradient, arg4::Gradient, arg5::Gradient, arg6::Gradient, arg7::Gradient)
    return jcall(obj, "linearCombination", Gradient, (Gradient, Gradient, Gradient, Gradient, Gradient, Gradient, Gradient, Gradient), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::Gradient, arg0::Vector{Gradient}, arg1::Vector{Gradient})
    return jcall(obj, "linearCombination", Gradient, (Vector{Gradient}, Vector{Gradient}), arg0, arg1)
end

function linear_combination(obj::Gradient, arg0::Vector{jdouble}, arg1::Vector{Gradient})
    return jcall(obj, "linearCombination", Gradient, (Vector{jdouble}, Vector{Gradient}), arg0, arg1)
end

function linear_combination(obj::Gradient, arg0::jdouble, arg1::Gradient, arg2::jdouble, arg3::Gradient)
    return jcall(obj, "linearCombination", Gradient, (jdouble, Gradient, jdouble, Gradient), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Gradient, arg0::jdouble, arg1::Gradient, arg2::jdouble, arg3::Gradient, arg4::jdouble, arg5::Gradient)
    return jcall(obj, "linearCombination", Gradient, (jdouble, Gradient, jdouble, Gradient, jdouble, Gradient), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Gradient, arg0::jdouble, arg1::Gradient, arg2::jdouble, arg3::Gradient, arg4::jdouble, arg5::Gradient, arg6::jdouble, arg7::Gradient)
    return jcall(obj, "linearCombination", Gradient, (jdouble, Gradient, jdouble, Gradient, jdouble, Gradient, jdouble, Gradient), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::Gradient)
    return jcall(obj, "log", Gradient, ())
end

function log10(obj::Gradient)
    return jcall(obj, "log10", Gradient, ())
end

function log1p(obj::Gradient)
    return jcall(obj, "log1p", Gradient, ())
end

function multiply(obj::Gradient, arg0::Gradient)
    return jcall(obj, "multiply", Gradient, (Gradient,), arg0)
end

function multiply(obj::Gradient, arg0::jdouble)
    return jcall(obj, "multiply", Gradient, (jdouble,), arg0)
end

function multiply(obj::Gradient, arg0::jint)
    return jcall(obj, "multiply", Gradient, (jint,), arg0)
end

function negate(obj::Gradient)
    return jcall(obj, "negate", Gradient, ())
end

function new_instance(obj::Gradient, arg0::jdouble)
    return jcall(obj, "newInstance", Gradient, (jdouble,), arg0)
end

function pow(::Type{Gradient}, arg0::jdouble, arg1::Gradient)
    return jcall(Gradient, "pow", Gradient, (jdouble, Gradient), arg0, arg1)
end

function pow(obj::Gradient, arg0::Gradient)
    return jcall(obj, "pow", Gradient, (Gradient,), arg0)
end

function pow(obj::Gradient, arg0::jdouble)
    return jcall(obj, "pow", Gradient, (jdouble,), arg0)
end

function pow(obj::Gradient, arg0::jint)
    return jcall(obj, "pow", Gradient, (jint,), arg0)
end

function reciprocal(obj::Gradient)
    return jcall(obj, "reciprocal", Gradient, ())
end

function remainder(obj::Gradient, arg0::Gradient)
    return jcall(obj, "remainder", Gradient, (Gradient,), arg0)
end

function remainder(obj::Gradient, arg0::jdouble)
    return jcall(obj, "remainder", Gradient, (jdouble,), arg0)
end

function rint(obj::Gradient)
    return jcall(obj, "rint", Gradient, ())
end

function root_n(obj::Gradient, arg0::jint)
    return jcall(obj, "rootN", Gradient, (jint,), arg0)
end

function scalb(obj::Gradient, arg0::jint)
    return jcall(obj, "scalb", Gradient, (jint,), arg0)
end

function signum(obj::Gradient)
    return jcall(obj, "signum", Gradient, ())
end

function sin(obj::Gradient)
    return jcall(obj, "sin", Gradient, ())
end

function sin_cos(obj::Gradient)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::Gradient)
    return jcall(obj, "sinh", Gradient, ())
end

function sqrt(obj::Gradient)
    return jcall(obj, "sqrt", Gradient, ())
end

function subtract(obj::Gradient, arg0::Gradient)
    return jcall(obj, "subtract", Gradient, (Gradient,), arg0)
end

function subtract(obj::Gradient, arg0::jdouble)
    return jcall(obj, "subtract", Gradient, (jdouble,), arg0)
end

function tan(obj::Gradient)
    return jcall(obj, "tan", Gradient, ())
end

function tanh(obj::Gradient)
    return jcall(obj, "tanh", Gradient, ())
end

function taylor(obj::Gradient, arg0::Vector{jdouble})
    return jcall(obj, "taylor", jdouble, (Vector{jdouble},), arg0)
end

function to_degrees(obj::Gradient)
    return jcall(obj, "toDegrees", Gradient, ())
end

function to_derivative_structure(obj::Gradient)
    return jcall(obj, "toDerivativeStructure", DerivativeStructure, ())
end

function to_radians(obj::Gradient)
    return jcall(obj, "toRadians", Gradient, ())
end

function variable(::Type{Gradient}, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(Gradient, "variable", Gradient, (jint, jint, jdouble), arg0, arg1, arg2)
end


text:
"""
    MattssonAlmquistVanDerWeide2018Accurate()

Coefficients of the optimized SBP operators with nonuniform grid given in
- Mattsson, Almquist, van der Weide (2018)
  Boundary optimized diagonal-norm SBP operators.
  Journal of Computational Physics 374, pp. 1261-1266.
"""
struct MattssonAlmquistVanDerWeide2018Accurate <: SourceOfCoefficients end

function Base.show(io::IO, source::MattssonAlmquistVanDerWeide2018Accurate)
  if get(io, :compact, false)
    summary(io, source)
  else
    print(io,
        "Mattsson, Almquist, van der Weide (2018) \n",
        "  Boundary optimized diagonal-norm SBP operators ('Accurate'). \n",
        "  Journal of Computational Physics 374, pp. 1261-1266.")
  end
end

function construct_grid(::MattssonAlmquistVanDerWeide2018Accurate, accuracy_order, xmin, xmax, N)
  @argcheck N > 6
  T = promote_type(typeof(xmin), typeof(xmax))

  if accuracy_order == 4
    xstart = SVector(
      T(0.0000000000000e+00),
      T(6.8764546205559e-01),
      T(1.8022115125776e+00),
      T(2.8022115125776e+00),
      T(3.8022115125776e+00),
    )
  elseif accuracy_order == 6
    xstart = SVector(
      T(0.0000000000000e+00),
      T(4.4090263368623e-01),
      T(1.2855984345073e+00),
      T(2.2638953951239e+00),
      T(3.2638953951239e+00),
      T(4.2638953951239e+00),
      T(5.2638953951239e+00),
    )
  elseif accuracy_order == 8
    xstart = SVector(
      T(0.0000000000000e+00),
      T(3.8118550247622e-01),
      T(1.1899550868338e+00),
      T(2.2476300175641e+00),
      T(3.3192851303204e+00),
      T(4.3192851303204e+00),
      T(5.3192851303204e+00),
      T(6.3192851303204e+00),
      T(7.3192851303204e+00),
    )
  else
    throw(ArgumentError("Accuracy order $accuracy_order not implemented/derived."))
  end

  BoundaryAdaptedGrid(xmin, xmax, xstart, N)
end

function first_derivative_coefficients(source::MattssonAlmquistVanDerWeide2018Accurate,
                                       order::Int, T=Float64, mode=FastMode())
  if order == 4
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,4}(SVector(T(-2.3518634632279443),
                                                T(3.0858932129093573),
                                                T(-0.9349061508864488),
                                                T(0.20087640120503583) )),
        # d2
        DerivativeCoefficientRow{T,1,4}(SVector(T(-0.6394095807755721),
                                                T(0),
                                                T(0.7917608885016885),
                                                T(-0.15235130772611635) )),
        # d3
        DerivativeCoefficientRow{T,1,5}(SVector(T(0.18446061657605967),
                                                T(-0.7539307161467195),
                                                T(0),
                                                T(0.6468087329936288),
                                                T(-0.07733863342297213) )),
        # d4
        DerivativeCoefficientRow{T,1,6}(SVector(T(-0.043308971097943676),
                                                T(0.1585246807786971),
                                                T(-0.7067880256335334),
                                                T(0),
                                                T(0.676082646803181),
                                                T(-0.08451033085039762) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(2//3), T(-1//12))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(0.21259737557798),
                            T(1.0260290400758),
                            T(1.0775123588954),
                            T(0.98607273802835) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  elseif order == 6
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,6}(SVector(T(-3.837232862086773),
                                                T(5.068376175016296),
                                                T(-1.6963756420892036),
                                                T(0.5851295073874213),
                                                T(-0.13275449453716695),
                                                T(0.012857316309426029) )),
        # d2
        DerivativeCoefficientRow{T,1,6}(SVector(T(-0.9591958116075092),
                                                T(0),
                                                T(1.2687130518352143),
                                                T(-0.3860731957925361),
                                                T(0.08345276954506707),
                                                T(-0.006896813980236088) )),
        # d3
        DerivativeCoefficientRow{T,1,6}(SVector(T(0.23226893958753303),
                                                T(-0.9178972823223088),
                                                T(0),
                                                T(0.8585728949642006),
                                                T(-0.19840158580452047),
                                                T(0.0254570335750956) )),
        # d4
        DerivativeCoefficientRow{T,1,7}(SVector(T(-0.07693308450770664),
                                                T(0.2682207441654631),
                                                T(-0.8244592785032373),
                                                T(0),
                                                T(0.7749305933808132),
                                                T(-0.15857634324473752),
                                                T(0.016817368709408833) )),
        # d5
        DerivativeCoefficientRow{T,1,8}(SVector(T(0.017248604775379643),
                                                T(-0.057293724956320936),
                                                T(0.1882698460378496),
                                                T(-0.7657841952075574),
                                                T(0),
                                                T(0.7505104744873311),
                                                T(-0.14956988077877112),
                                                T(0.01661887564208568) )),
        # d6
        DerivativeCoefficientRow{T,1,9}(SVector(T(-0.0016761725437960326),
                                                T(0.004750928272020083),
                                                T(-0.02423857479897864),
                                                T(0.15723370478311835),
                                                T(-0.753044107168037),
                                                T(0),
                                                T(0.7503740531217644),
                                                T(-0.15007481062435288),
                                                T(0.016674978958261432) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(0.13030223027124),
                            T(0.68851501587715),
                            T(0.95166202564389),
                            T(0.99103890475697),
                            T(1.0028757074552),
                            T(0.99950151111941) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  elseif order == 8
    left_boundary = (
        # d1
        DerivativeCoefficientRow{T,1,8}(SVector(T(-4.647545021331372),
                                                T(6.254178662563343),
                                                T(-2.4139100510438665),
                                                T(1.256639509571857),
                                                T(-0.647667696659897),
                                                T(0.24570663393330786),
                                                T(-0.05204535767697368),
                                                T(0.004643320643692943) )),
        # d2
        DerivativeCoefficientRow{T,1,8}(SVector(T(-1.0868211647678514),
                                                T(0),
                                                T(1.5195364183857292),
                                                T(-0.6543667982790691),
                                                T(0.3128620695707949),
                                                T(-0.11086808101811908),
                                                T(0.021234895991755504),
                                                T(-0.0015773398832359302) )),
        # d3
        DerivativeCoefficientRow{T,1,8}(SVector(T(0.26780878369504263),
                                                T(-0.9701235413358851),
                                                T(0),
                                                T(0.9726229663976735),
                                                T(-0.3684397832422662),
                                                T(0.1161859354443521),
                                                T(-0.018907619882971544),
                                                T(0.0008532589240538837) )),
        # d4
        DerivativeCoefficientRow{T,1,8}(SVector(T(-0.12264219453174283),
                                                T(0.3675044946198741),
                                                T(-0.8555984502334516),
                                                T(0),
                                                T(0.7955264157053945),
                                                T(-0.22405034810099697),
                                                T(0.04290097854980221),
                                                T(-0.0036408960088885555) )),
        # d5
        DerivativeCoefficientRow{T,1,9}(SVector(T(0.0680142412945559),
                                                T(-0.18906569319897437),
                                                T(0.34874693977149884),
                                                T(-0.855998538940948),
                                                T(0),
                                                T(0.7918634438230676),
                                                T(-0.19783080393205543),
                                                T(0.03775653807587011),
                                                T(-0.0034861268930073425)  )),
        # d6
        DerivativeCoefficientRow{T,1,10}(SVector(T(-0.026557903380662235),
                                                T(0.06895974248171861),
                                                T(-0.11319486586973267),
                                                T(0.24813802673084387),
                                                T(-0.8150412273910146),
                                                T(0),
                                                T(0.8048396208993126),
                                                T(-0.20182899393713088),
                                                T(0.038273766032179475),
                                                T(-0.0035881655655168633) )),
        # d7
        DerivativeCoefficientRow{T,1,11}(SVector(T(0.005594622888681003),
                                                T(-0.013135637494343391),
                                                T(0.0183198514986364),
                                                T(-0.047252720307139756),
                                                T(0.20250470348441063),
                                                T(-0.8004261411357496),
                                                T(0),
                                                T(0.7997353219580023),
                                                T(-0.19983539676215634),
                                                T(0.03806388509755359),
                                                T(-0.0035684892278956487) )),
        # d8
        DerivativeCoefficientRow{T,1,12}(SVector(T(-0.000499585189812272),
                                                T(0.0009766036909553027),
                                                T(-0.000827480914165871),
                                                T(0.0040138390353474905),
                                                T(-0.038683470124513085),
                                                T(0.20090351330033354),
                                                T(-0.8004576105033867),
                                                T(0),
                                                T(0.8000635199885898),
                                                T(-0.20001587999714746),
                                                T(0.038098262856599514),
                                                T(-0.0035717121428062043) )),
    )
    right_boundary = .- left_boundary
    upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
    central_coef = zero(T)
    lower_coef = -upper_coef
    left_weights = SVector( T(0.1075836807831),
                            T(0.61909685107891),
                            T(0.96971176519117),
                            T(1.1023441350947),
                            T(1.0244688965833),
                            T(0.99533550116831),
                            T(1.0008236941028),
                            T(0.99992060631812) )
    right_weights = left_weights
    left_boundary_derivatives = Tuple{}()
    right_boundary_derivatives = left_boundary_derivatives

    DerivativeCoefficients(left_boundary, right_boundary,
                            left_boundary_derivatives, right_boundary_derivatives,
                            lower_coef, central_coef, upper_coef,
                            left_weights, right_weights, mode, 1, order, source)
  # elseif order == 8
  #   left_boundary = (
  #       # d1
  #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d2
  #       DerivativeCoefficientRow{T,1,6}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d3
  #       DerivativeCoefficientRow{T,1,7}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d4
  #       DerivativeCoefficientRow{T,1,8}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d5
  #       DerivativeCoefficientRow{T,1,9}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #       # d6
  #       DerivativeCoefficientRow{T,1,10}(SVector(T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T(),
  #                                               T() )),
  #   )
  #   right_boundary = .- left_boundary
  #   upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
  #   central_coef = zero(T)
  #   lower_coef = -upper_coef
  #   left_weights = SVector( T(),
  #                           T(),
  #                           T(),
  #                           T(),
  #                           T(),
  #                           T() )
  #   right_weights = left_weights
  #   left_boundary_derivatives = Tuple{}()
  #   right_boundary_derivatives = left_boundary_derivatives

  #   DerivativeCoefficients(left_boundary, right_boundary,
  #                           left_boundary_derivatives, right_boundary_derivatives,
  #                           lower_coef, central_coef, upper_coef,
  #                           left_weights, right_weights, mode, 1, order, source)
  else
    throw(ArgumentError("Order $order not implemented/derived."))
  end
end


text:
using BitConverter
@testset "Script" begin
    @testset "Parse" begin
        script_pubkey = IOBuffer(hex2bytes("6a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937"))
        script = Script(script_pubkey)
        want = hex2bytes("304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a71601")
        @test script.data[1] == want
        want = hex2bytes("035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937")
        @test script.data[2] == want
        script_pubkey = IOBuffer(hex2bytes("fdfe0000483045022100c679944ff8f20373685e1122b581f64752c1d22c67f6f3ae26333aa9c3f43d730220793233401f87f640f9c39207349ffef42d0e27046755263c0a69c436ab07febc01483045022100eadc1c6e72f241c3e076a7109b8053db53987f3fcc99e3f88fc4e52dbfd5f3a202201f02cbff194c41e6f8da762e024a7ab85c1b1616b74720f13283043e9e99dab8014c69522102b0c7be446b92624112f3c7d4ffc214921c74c1cb891bf945c49fbe5981ee026b21039021c9391e328e0cb3b61ba05dcc5e122ab234e55d1502e59b10d8f588aea4632102f3bd8f64363066f35968bd82ed9c6e8afecbd6136311bb51e91204f614144e9b53aeffffffff05a08601000000000017a914081fbb6ec9d83104367eb1a6a59e2a92417d79298700350c00000000001976a914677345c7376dfda2c52ad9b6a153b643b6409a3788acc7f341160000000017a914234c15756b9599314c9299340eaabab7f1810d8287c02709000000000017a91469be3ca6195efcab5194e1530164ec47637d44308740420f00000000001976a91487fadba66b9e48c0c8082f33107fdb01970eb80388ac00000000"))
        want = hex2bytes("522102b0c7be446b92624112f3c7d4ffc214921c74c1cb891bf945c49fbe5981ee026b21039021c9391e328e0cb3b61ba05dcc5e122ab234e55d1502e59b10d8f588aea4632102f3bd8f64363066f35968bd82ed9c6e8afecbd6136311bb51e91204f614144e9b53ae")
        script = Script(script_pubkey)
        @test script.data[4] == want
    end
    @testset "Serialize" begin
        want = "6a47304402207899531a52d59a6de200179928ca900254a36b8dff8bb75f5f5d71b1cdc26125022008b422690b8461cb52c3cc30330b23d574351872b7c361e9aae3649071c1a7160121035d5c93d9ac96881f19ba1f686f15f009ded7c62efe85a872e6a19b43c15a2937"
        script_pubkey = IOBuffer(hex2bytes(want))
        script = Script(script_pubkey)
        @test bytes2hex(serialize(script)) == want
    end
    @testset "P2WPKH" begin
        h160 = hex2bytes("74d691da1574e6b3c192ecfb52cc8984ee7b6c56")
        script = BitcoinPrimitives.script(h160, type=:P2WPKH)
        @test BitcoinPrimitives.type(script) == :P2WPKH
    end
end


text:
using Distributed
import NODAL

addprocs()

@everywhere begin
    using NODAL
    function rosenbrock(x::Configuration, parameters::Dict{Symbol, Any})
        return (1.0 - x["i0"].value)^2 + 100.0 * (x["i1"].value - x["i0"].value^2)^2
    end
end

function optimize_rosenbrock()
    configuration = Configuration([FloatParameter(-2.0, 2.0, 0.0,"i0"),
                                   FloatParameter(-2.0, 2.0, 0.0,"i1")],
                                   "rosenbrock_config")

    tuning_run = Run(cost                = rosenbrock,
                     starting_point      = configuration,
                     stopping_criterion  = elapsed_time_criterion,
                     report_after        = 10,
                     reporting_criterion = elapsed_time_reporting_criterion,
                     duration            = 60,
                     methods             = [[:simulated_annealing 1];
                                            [:iterative_first_improvement 1];
                                            [:iterated_local_search 1];
                                            [:randomized_first_improvement 1];
                                            [:iterative_greedy_construction 1];
                                            [:iterative_probabilistic_improvement 1];])

    @spawn optimize(tuning_run)
    result = take!(tuning_run.channel)

    print(result)
    while !result.is_final
        result = take!(tuning_run.channel)
        print(result)
    end
end

optimize_rosenbrock()


text:
## Preloads
module DHC_Utils

    using Statistics
    using FFTW
    using Plots
    using LinearAlgebra
    using StaticArrays
    using HybridArrays

    export finklet
    export fink_filter_bank_slow
    export fink_filter_bank
    export fink_filter_list
    export speedy_DHC

    ## Define wavelet in Fourier space
    function finklet(j, l)
        # -------- set filters
        jrad = 7-j
        dθ = π/8        # 8 angular bins hardwired
        θ_l = dθ*l
        # -------- define coordinates
        nx = 256
        xbox = LinRange(-nx/2, nx/2-1 , nx)
        # make a 256x256 grid of X
        sx = xbox' .* ones(nx)
        sy = ones(nx)' .* xbox
        r  = sqrt.((sx).^2 + (sy).^2)
        θ  = mod.(atan.(sy, sx).+π .-θ_l,2*π)
        nozeros = r .> 0
        logr = log2.(r[nozeros])
        r[nozeros] = logr
        # -------- in Fourier plane, envelope of psi_j,l
        mask = (abs.(θ.-π).<= dθ) .& (abs.(r.-jrad) .<= 1)
        # -------- angular part
        ang = cos.((θ.-π).*4)
        # -------- radial part
        rad = cos.((r.-jrad).*π./2)
        psi = mask.*ang.*rad             #mask times angular part times radial part
        return psi
    end

    ## Compute the whole filter bank.  Legacy code for comparison
    function fink_filter_bank_slow(J,L)
        fink_filter = Array{Float64, 4}(undef, 256, 256, J, L)
        for l = 1:L
            for j = 1:J
                @inbounds fink_filter[:,:,j,l]=fftshift(finklet(j-1,l-1))
            end
        end
        return fink_filter
    end

    ## Faster filter bank generation.  Only square filters allowed.
    function fink_filter_bank(J::Integer, L::Integer; nx::Integer=256, wid::Integer=1)

        # -------- set parameters
        dθ   = π/8        # 8 angular bins hardwired
        dx   = nx÷2-1

        # -------- allocate output array of zeros
        filt = zeros(Float64, nx, nx, J, L)

        # -------- allocate theta and logr arrays
        logr = zeros(Float64, nx, nx)
        θ    = zeros(Float64, nx, nx)

        for l = 0:L-1
            θ_l = dθ*l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx, nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx, nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= (dθ*wid)

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r = sqrt(sx^2 + sy^2)
                        logr[y, x] = log2(max(1,r))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        #          the angular factor is the same for all j
            F_angular = cos.((θ[angmask].-π).*(4/wid))

        # -------- loop over j for the radial part
            for j = 0:J-1
                jrad  = 7-j
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1)

        # -------- radial part
                F_radial = cos.(Δj[rmask] .* (π/2))
                ind      = angmask[rmask]
                filt[ind,j+1,l+1] = F_radial .* F_angular[rmask]
            end
        end
        return filt
    end

    ## Make a list of non-zero pixels for the (Fourier plane) filters
    function fink_filter_list(filt)
        (ny,nx,J,L) = size(filt)

        # Allocate output arrays
        filtind = fill(CartesianIndex{2}[], J, L)
        filtval = fill(Float64[], J, L)

        # Loop over J,L and record non-zero values
        for l=1:L
            for j=1:J
                f = @view filt[:,:,j,l]
                ind = findall(f .> 1E-13)
                val = f[ind]
                filtind[j,l] = ind
                filtval[j,l] = val
            end
        end
        return [filtind, filtval]
    end

    ## Todo list
    # Check if 2 threads really help FFT when computer is busy
    function speedy_DHC(image::Array{Float64,2}, filter_list)
        # Use 2 threads for FFT
        FFTW.set_num_threads(2)

        # array sizes
        (Nx, Ny)  = size(image)
        (J,L)     = size(filter_list[1])

        # allocate coeff arrays
        out_coeff = []
        S0  = zeros(Float64, 2)
        S1  = zeros(Float64, J, L)
        S20 = zeros(Float64, J, L, J, L)
        S12 = zeros(Float64, J, L, J, L)

        # allocate image arrays for internal use
        im_fdf_0_1 = zeros(Float64,           Nx, Ny, J, L)   # this must be zeroed!
        im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, J, L)

        ## 0th Order
        S0[1]   = mean(image)
        norm_im = image.-S0[1]
        S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
        norm_im ./= sqrt(Nx*Ny*S0[2])

        append!(out_coeff,S0[:])

        ## 1st Order
        im_fd_0 = fft(norm_im)

        # unpack filter_list
        f_ind   = filter_list[1]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_list[2]

        zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

        # make a FFTW "plan" for an array of the given size and type
        P   = plan_ifft(im_fd_0)   # P is an operator, P*im is ifft(im)

        ## Main 1st Order and Precompute 2nd Order
        for l = 1:L
            for j = 1:J
                S1tot = 0.0
                f_i = f_ind[j,l]  # CartesianIndex list for filter
                f_v = f_val[j,l]  # Values for f_i
                # for (ind, val) in zip(f_i, f_v)   # this is slower!
                if length(f_i) > 0
                    for i = 1:length(f_i)
                        ind       = f_i[i]
                        zval      = f_v[i] * im_fd_0[ind]
                        S1tot    += abs2(zval)
                        zarr[ind] = zval
                        im_fdf_0_1[ind,j,l] = abs(zval)
                    end
                    S1[j,l] = S1tot
                    im_rd_0_1[:,:,j,l] .= abs2.(P*zarr)
                    zarr[f_i] .= 0
                end
            end
        end
        append!(out_coeff, S1[:])

        # we stored the abs()^2, so take sqrt (this is faster to do all at once)
        im_rd_0_1 .= sqrt.(im_rd_0_1)

        ## 2nd Order
        Amat = reshape(im_fdf_0_1, Nx*Ny, J*L)
        S12  = reshape(Amat' * Amat, J, L, J, L)
        Amat = reshape(im_rd_0_1, Nx*Ny, J*L)
        S20  = reshape(Amat' * Amat, J, L, J, L)

        append!(out_coeff, S20)
        append!(out_coeff, S12)

        return out_coeff
    end

end # of module

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
abstract type AbstractDomainDistribution end

"""
    DomainDistribution{D<:Distribution, O<:AbstractFisheryDomain}
        dist::D
        domain::O

A distribution over a fishery domain. Methods for `rand` (without size) are
defined for univariate, multivariate, and `MatrixNormal` distributions, and
return realizations reshaped to the size of the fishery domain.
"""
struct DomainDistribution{D, O<:AbstractFisheryDomain} <: AbstractDomainDistribution
    dist::D
    domain::O
end

## Construct a log-normal domain distribution where the μ vector provided
## defaults to the mean, but can be changed via the `s` argument (following the
## conventions of the `location` function from Distributions.jl).
function DomainDistribution(t::Type{<:Distributions.AbstractMvLogNormal},
                            μ::AbstractVector, Σ::AbstractMatrix,
                            domain::AbstractFisheryDomain,
                            s::Symbol = :mean)
    loc = location(t, s, μ, Σ)
    DomainDistribution(MvLogNormal(loc, Σ), domain)
end

function DomainDistribution(t::Type{<:Distributions.AbstractMvLogNormal},
                            μ::AbstractVector, Σ::AbstractPDMat,
                            domain::AbstractFisheryDomain,
                            s::Symbol = :mean)
    loc = location(t, s, μ, Σ.mat)
    DomainDistribution(MvLogNormal(loc, Σ), domain)
end

function DomainDistribution(ddist::DomainDistribution, domain::AbstractFisheryDomain)
    DomainDistribution(ddist.dist, domain)
end

## Pass an instance, not a type (exceptions above)
function DomainDistribution(t::Type, domain::AbstractFisheryDomain)
    error("Passed $t as a type, need an instance e.g. $t()")
end

domain(dd::DomainDistribution) = dd.domain

function rand(DD::DomainDistribution{D, O}) where {D<:UnivariateDistribution, O}
    rand(DD.dist, size(DD.domain)...)
end

function rand(DD::DomainDistribution{D, O}) where {D<:MultivariateDistribution, O}
    x = rand(DD.dist)
    reshape(x, size(DD.domain))
end

function rand(dd::DomainDistribution{D, O}) where {D<:NeutralLandscapeMaker, O}
    rand(dd.dist, size(dd.domain))
end

# Multiple domain distribution
struct MultiDomainDistribution <: AbstractDomainDistribution
    ddists::Vector{AbstractDomainDistribution}

    function MultiDomainDistribution(ddists::Vector{AbstractDomainDistribution})
        all(==(domain(ddists[1])), domain.(ddists)) || error("All components must have the same domain")
        new(ddists)
    end
end

function MultiDomainDistribution(dists::Vector, domain::AbstractFisheryDomain)
    ddist_vec = AbstractDomainDistribution[]
    for d in dists
        if d isa AbstractDomainDistribution
            push!(ddist_vec, d)
        else
            push!(ddist_vec, DomainDistribution(d, domain))
        end
    end

    MultiDomainDistribution(ddist_vec)
end

# If vector of domain distributions is provided without a domain, extract the
# domain from distributions that have it, check for equality, then use above
# constructor.
function MultiDomainDistribution(ddists::Vector)
    dom_vec = AbstractFisheryDomain[]
    for dist in ddists
        if dist isa AbstractDomainDistribution
            push!(dom_vec, domain(dist))
        end
    end
    length(dom_vec) > 0 || error("Must provide a domain")
    all(==(dom_vec[1]), dom_vec) || error("Each distribution must use the same domain")
    MultiDomainDistribution(ddists, dom_vec[1])
end

length(mdd::MultiDomainDistribution) = length(mdd.ddists)
domain(mdd::MultiDomainDistribution) = domain(mdd.ddists[1])
rand(mdd::MultiDomainDistribution) = rand.(mdd.ddists)
getindex(mdd::MultiDomainDistribution, i) = getindex(mdd.ddists, i)

# Blended distributions - combine multiple, weighted distributions over the domain
struct BlendedDomainDistribution{T} <: AbstractDomainDistribution
    ddists::MultiDomainDistribution
    wts::Vector{T}

    function BlendedDomainDistribution(ddists::MultiDomainDistribution, wts::Vector{T} = ones(length(ddists))) where {D<:AbstractDomainDistribution, T<:Real}
        new{T}(ddists, wts)
    end
end

function BlendedDomainDistribution(dists::Vector, domain::AbstractFisheryDomain, wts::Vector{<:Real} = ones(length(dists)))
    mdd = MultiDomainDistribution(dists, domain)
    BlendedDomainDistribution(mdd, wts)
end

domain(bdd::BlendedDomainDistribution) = domain(bdd.ddists)

function rand(bddist::BlendedDomainDistribution)
    rvec = rand(bddist.ddists)
    blend(rvec, bddist.wts)
end

# Classified domain distributions - classify
struct ClassifiedDomainDistribution{D, T, M, N} <: AbstractDomainDistribution
    ddist::D
    wts::Vector{T}
    mask::M

    function ClassifiedDomainDistribution(ddist::D, wts::Vector{T}, mask::M = nothing) where {D<:AbstractDomainDistribution, T<:Real, M}
        N = length(wts)
        new{D, T, M, N}(ddist, wts, mask)
    end
end

function ClassifiedDomainDistribution(ddist::AbstractDomainDistribution, nclass::Integer, mask = nothing)
    wts = ones(nclass)
    ClassifiedDomainDistribution(ddist, wts, mask)
end

function ClassifiedDomainDistribution(dist, domain::AbstractFisheryDomain, classes, mask = nothing)
    ClassifiedDomainDistribution(DomainDistribution(dist, domain), classes, mask)
end

domain(cdd::ClassifiedDomainDistribution) = cdd.ddist.domain

function rand(cddist::ClassifiedDomainDistribution{D, T, M, 2}) where {D, T, M}
    res = rand(cddist.ddist)
    classify!(res, cddist.wts, cddist.mask)
    res .== 2
end

function rand(cddist::ClassifiedDomainDistribution)
    res = rand(cddist.ddist)
    classify!(res, cddist.wts, cddist.mask)
    trunc.(Int, res)
end


text:
#
# Run examples defined in `Examples.jl` module.
#

include("./Examples.jl")

const ex = Examples

print(
    # s(2),
    # s(2.0),
    # s([2, 2])
    )

print(
    ex.factorial(1),
    ex.factorial(2),
    ex.factorial(3)
    )

max_stars_repo_name: groundf/course-julia

max_issues_repo_name: groundf/course-julia

max_forks_repo_name: groundf/course-julia


text:
using Test

_modulenames = quote
    SIMD
    MappedArrays
    SIMDIndexing
end
modulenames = [x for x in _modulenames.args if x isa Symbol]
@debug "Loading dependencies:" modulenames
for m in modulenames
    @eval import $m
end

modules = eval.(modulenames)
ambiguities = detect_ambiguities(modules..., imported=true, recursive=true)
@test ambiguities == []


text:
# A disjoint set implementation adapted from
# https://github.com/JuliaCollections/DataStructures.jl/blob/f57330a3b46f779b261e6c07f199c88936f28839/src/disjoint_set.jl
# under the MIT license: https://github.com/JuliaCollections/DataStructures.jl/blob/master/License.md

# imports
import ._TOP_MOD:
    length,
    eltype,
    union!,
    push!
# usings
import ._TOP_MOD:
    OneTo, collect, zero, zeros, one, typemax

# Disjoint-Set

############################################################
#
#   A forest of disjoint sets of integers
#
#   Since each element is an integer, we can use arrays
#   instead of dictionary (for efficiency)
#
#   Disjoint sets over other key types can be implemented
#   based on an IntDisjointSet through a map from the key
#   to an integer index
#
############################################################

_intdisjointset_bounds_err_msg(T) = "the maximum number of elements in IntDisjointSet{$T} is $(typemax(T))"

"""
    IntDisjointSet{T<:Integer}(n::Integer)

A forest of disjoint sets of integers, which is a data structure
(also called a union–find data structure or merge–find set)
that tracks a set of elements partitioned
into a number of disjoint (non-overlapping) subsets.
"""
mutable struct IntDisjointSet{T<:Integer}
    parents::Vector{T}
    ranks::Vector{T}
    ngroups::T
end

IntDisjointSet(n::T) where {T<:Integer} = IntDisjointSet{T}(collect(OneTo(n)), zeros(T, n), n)
IntDisjointSet{T}(n::Integer) where {T<:Integer} = IntDisjointSet{T}(collect(OneTo(T(n))), zeros(T, T(n)), T(n))
length(s::IntDisjointSet) = length(s.parents)

"""
    num_groups(s::IntDisjointSet)

Get a number of groups.
"""
num_groups(s::IntDisjointSet) = s.ngroups
eltype(::Type{IntDisjointSet{T}}) where {T<:Integer} = T

# find the root element of the subset that contains x
# path compression is implemented here
function find_root_impl!(parents::Vector{T}, x::Integer) where {T<:Integer}
    p = parents[x]
    @inbounds if parents[p] != p
        parents[x] = p = _find_root_impl!(parents, p)
    end
    return p
end

# unsafe version of the above
function _find_root_impl!(parents::Vector{T}, x::Integer) where {T<:Integer}
    @inbounds p = parents[x]
    @inbounds if parents[p] != p
        parents[x] = p = _find_root_impl!(parents, p)
    end
    return p
end

"""
    find_root!(s::IntDisjointSet{T}, x::T)

Find the root element of the subset that contains an member `x`.
Path compression happens here.
"""
find_root!(s::IntDisjointSet{T}, x::T) where {T<:Integer} = find_root_impl!(s.parents, x)

"""
    in_same_set(s::IntDisjointSet{T}, x::T, y::T)

Returns `true` if `x` and `y` belong to the same subset in `s`, and `false` otherwise.
"""
in_same_set(s::IntDisjointSet{T}, x::T, y::T) where {T<:Integer} = find_root!(s, x) == find_root!(s, y)

"""
    union!(s::IntDisjointSet{T}, x::T, y::T)

Merge the subset containing `x` and that containing `y` into one
and return the root of the new set.
"""
function union!(s::IntDisjointSet{T}, x::T, y::T) where {T<:Integer}
    parents = s.parents
    xroot = find_root_impl!(parents, x)
    yroot = find_root_impl!(parents, y)
    return xroot != yroot ? root_union!(s, xroot, yroot) : xroot
end

"""
    root_union!(s::IntDisjointSet{T}, x::T, y::T)

Form a new set that is the union of the two sets whose root elements are
`x` and `y` and return the root of the new set.
Assume `x ≠ y` (unsafe).
"""
function root_union!(s::IntDisjointSet{T}, x::T, y::T) where {T<:Integer}
    parents = s.parents
    rks = s.ranks
    @inbounds xrank = rks[x]
    @inbounds yrank = rks[y]

    if xrank < yrank
        x, y = y, x
    elseif xrank == yrank
        rks[x] += one(T)
    end
    @inbounds parents[y] = x
    s.ngroups -= one(T)
    return x
end

"""
    push!(s::IntDisjointSet{T})

Make a new subset with an automatically chosen new element `x`.
Returns the new element. Throw an `ArgumentError` if the
capacity of the set would be exceeded.
"""
function push!(s::IntDisjointSet{T}) where {T<:Integer}
    l = length(s)
    l < typemax(T) || throw(ArgumentError(_intdisjointset_bounds_err_msg(T)))
    x = l + one(T)
    push!(s.parents, x)
    push!(s.ranks, zero(T))
    s.ngroups += one(T)
    return x
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: September 2020

using InvertibleNetworks, LinearAlgebra, Test, Statistics, Random

# Random seed
Random.seed!(1);

###############################################################################
# Initialization

# Input
nx = 28
ny = 28*2
nc = 4
n_hidden = 64
batchsize = 5
X = rand(Float32, nx, ny, nc, batchsize)

# Layers and initialization
logdet = true
N1 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N2 = ActNorm(nc; logdet=logdet)
N3 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N = Composition(N1, N2, N3)

###############################################################################
# Test coherency of composition and manual composition

Y = copy(X)
Y, l3 = N3.forward(Y)
Y, l2 = N2.forward(Y)
Y, l1 = N1.forward(Y)
l = l1+l2+l3

Y_, l_ = N.forward(X)

@test isapprox(Y, Y_; rtol=1f-3)
@test isapprox(l, l_; rtol=1f-3)

###############################################################################
# Test coherency of composition and ∘

# Composition with ∘
N_ = N1∘N2∘N3

Y, l = N.forward(X)
Y_, l_ = N_.forward(X)

@test isapprox(Y, Y_; rtol=1f-3)
@test isapprox(l, l_; rtol=1f-3)

###############################################################################
# Test invertibility

@test isapprox(X, N.inverse(N.forward(X)[1]); rtol=1f-3)
@test isapprox(X, N.forward(N.inverse(X))[1]; rtol=1f-3)

###############################################################################
# Gradient Test

function loss(N, X, Y)

    # Forward pass
    Y_, lgdet = N.forward(X)

    # Residual and function value
    ΔY = Y_ - Y
    f = 0.5f0/batchsize*norm(ΔY)^2f0-lgdet

    # Back propagation
    ΔX, _ = N.backward(ΔY./batchsize, Y_)

    return f, ΔX, deepcopy(get_grads(N))

end

# Initializing nets
N1 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N2 = ActNorm(nc; logdet=logdet)
N3 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N0 = Composition(N1, N2, N3); N0.forward(rand(Float32, nx, ny, nc, batchsize))
θ0 = deepcopy(get_params(N0))
N1 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N2 = ActNorm(nc; logdet=logdet)
N3 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N = Composition(N1, N2, N3); N.forward(rand(Float32, nx, ny, nc, batchsize))
θ = deepcopy(get_params(N))

# Inputs
X = rand(Float32, nx, ny, nc, batchsize)
X0 = rand(Float32, nx, ny, nc, batchsize)
dX = rand(Float32, nx, ny, nc, batchsize)
Y = rand(Float32, nx, ny, nc, batchsize)
dθ = θ-θ0
for i = 1:length(dθ)
    (norm(θ0[i].data) != 0f0) && (dθ[i].data .*= norm(θ0[i].data)/norm(dθ[i].data))
end

# Gradient test (input)
set_params!(N0, θ0)
f0, ΔX, Δθ = loss(N0, X, Y)
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)
print("\nGradient test sequential network: input\n")
for j=1:maxiter
    f = loss(N0, X+h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h /= 2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test (parameters)
h = 0.1f0
maxiter = 5
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)
print("\nGradient test sequential network: parameters\n")
for j=1:maxiter
    set_params!(N0, θ0+h*dθ)
    f = loss(N0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(dθ, Δθ))
    print(err3[j], "; ", err4[j], "\n")
    global h /= 2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# Jacobian-related tests

# Gradient test

# Initializing nets
N1 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N2 = ActNorm(nc; logdet=logdet)
N3 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N0 = Composition(N1, N2, N3); N0.forward(rand(Float32, nx, ny, nc, batchsize))
θ0 = deepcopy(get_params(N0))
N1 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N2 = ActNorm(nc; logdet=logdet)
N3 = CouplingLayerHINT(nc, n_hidden; permute="full", logdet=logdet)
N = Composition(N1, N2, N3); N.forward(rand(Float32, nx, ny, nc, batchsize))
θ = deepcopy(get_params(N))
X = randn(Float32, nx, ny, nc, batchsize)

# Perturbation (normalized)
dθ = θ-θ0; dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)
dX = randn(Float32, nx, ny, nc, batchsize); dX *= norm(X)/norm(dX)

# Jacobian eval
dY, Y, _, _ = N.jacobian(dX, dθ, X)

# Test
print("\nJacobian test\n")
h = 0.1f0
maxiter = 5
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(N, θ+h*dθ)
    Y_loc, _ = N.forward(X+h*dX)
    err5[j] = norm(Y_loc - Y)
    err6[j] = norm(Y_loc - Y - h*dY)
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Adjoint test

set_params!(N, θ)
dY, Y, _, _ = N.jacobian(dX, dθ, X)
dY_ = randn(Float32, size(dY))
dX_, dθ_, _, _ = N.adjointJacobian(dY_, Y)
a = dot(dY, dY_)
b = dot(dX, dX_)+dot(dθ, dθ_)
@test isapprox(a, b; rtol=1f-3)


text:
## gradient function defintions

## GRAD RESULT

struct GradResult
    tape::Tape
    gvars::Dict{Int, Any}  # gradient vars: argid -> gradient var
end

function GradResult(tape::Tape)
    gvars = Dict{Int,Any}()
    # struct fields
    for (argid, dct) in tape.sfields
        gvars[argid] = Dict(field_path => tape.derivs[var_id]
                            for (field_path, var_id) in dct
                            if haskey(tape.derivs, var_id))  # not all fields may have derivatives
    end
    # other arguments
    struct_arg_ids = Set(keys(tape.sfields))
    for op in tape
        if op isa Input && !in(op.argid, struct_arg_ids)
            gvars[op.argid] = tape.derivs[op.var.id]
        end
    end
    return GradResult(tape, gvars)
end

Base.show(io::IO, g::GradResult) = print(io, "GradResult($(length(g.gvars)))")

function getindex(g::GradResult, argid::Int)
    tape = g.tape
    gvar = g.gvars[argid]
    if isa(gvar, Dict)
        return Dict(f => tape[id].var.val for (f, id) in gvar)
    else
        return tape[gvar].var.val
    end
end

## GRAD

getderiv(tape::Tape, id::Int) = tape[tape.derivs[id]].var
getderiv(tape::Tape, var::TAny) = getderiv(tape, var.id)
setderiv!(tape::Tape, var_id::Int, grad_var_id::Int) = (tape.derivs[var_id] = grad_var_id)
setderiv!(tape::Tape, var::TAny, grad_var::TAny) = (tape.derivs[var.id] = grad_var.id)

function rev_step!(op::Union{Call, Bcast}, i::Int)
    tape = op.var.tape
    y = op.var
    x = op.args[i]
    dy = getderiv(tape, y)
    dx = grad!(dy, Val(i), op)
    if !haskey(tape.derivs, x.id)
        setderiv!(tape, x, dx)
    else
        old_dx = getderiv(tape, x)
        new_dx = record!(tape, Call, +, (dx, old_dx))
        setderiv!(tape, x, new_dx)
    end
end

function back!(tape::Tape)
    # z - final variable, y - resulting variable of current op, x - dependencies of y
    # dy - derivative of z w.r.t. y
    z = tape[end].var
    # using Float32 for seed since for 64-bit args it will be expanded anyway
    dy = record!(tape, Constant, 1.0f0)
    # set initial derivative value
    tape.derivs[z.id] = dy.id
    for op in reverse(tape.ops[1:end-1])
        if op isa Call || op isa Bcast
            for i=1:length(op.args)
                # println("op = $op; i = $i")
                # backpropagate only non-constant tracked vars
                arg_op = tape[getid(op.args[i])]
                if op.args[i] isa TAny && !isa(arg_op, Constant)
                    rev_step!(op, i)
                end
            end
        end
    end
end

function make_tracked_args(tape::Tape, args...)
    targs = []
    for (argid, arg) in enumerate(args)
        if isstruct(arg)
            targ = record_struct!(tape, arg, argid)
        else
            targ = record!(tape, Input, arg; argid=argid)
        end
        push!(targs, targ)
    end
    return targs
end

"""
For each input that has a derivative on this tape check if the derivative
has the same size as the input.
"""
function check_deriv_sizes(tape::Tape)
    for (var_id, grad_var_id) in tape.derivs
        var_size = size(getvalue(tape[var_id]))
        grad_var_size = size(getvalue(tape[grad_var_id]))
        if  var_size != grad_var_size
            @warn "Gradient %$grad_var_id has size $grad_var_size, " *
                "but original variable %$var_id has size $var_size"
        end
    end
end

function _grad(f::Function, args...)
    tape = Tape(guess_device(args))
    # wrap args into tracked data
    targs = make_tracked_args(tape, args...)
    # execute function to fill in the tape
    tres = f(targs...)
    tape.resultid = getid(tres)
    # backpropagate gradients
    back!(tape)
    # consistency check
    check_deriv_sizes(tape)
    # construct GradResult object that wraps tape and provide accessors for computed derivatives
    return tres.val, GradResult(tape)
end

const GRAD_CACHE = Dict{Any, Tape}()

"""
Find gradient of `f` w.r.t. its arguments.
Example:

    val, g = grad(sum, rand(3))

where:
  - val is the value of `f` at this point
  - g::GradResult is a collection of gradients

GradResult is indexed by argument index and contains gradients
in a format most suitable for that argument, namely:

  - for arrays: arrays of the same type and size
  - for reals: reals
  - for mutable structs: dictionary of {(:field, :path) => value} pairs.

All gradients can be applied to original variables using `update!()` function.
"""
function grad(f::Function, args...; static=true)
    if static
        # key conists of function type and type of argument (for structs) or its size
        cache_key = (f, ([isstruct(arg) ? typeof(arg) : size(arg) for arg in args]...,))
        if haskey(GRAD_CACHE, cache_key)
            tape = GRAD_CACHE[cache_key]
            play!(tape, args...)
            return getvalue(tape[tape.resultid]), GradResult(tape)
        else
            val, g = _grad(f, args...)
            compile!(g.tape)
            GRAD_CACHE[cache_key] = g.tape
            return val, g
        end
    else
        return _grad(f, args...)
    end
end


text:
global warnings_enabled = true;

function warnings(val)
    global warnings_enabled = val
end

function warn(msg)
    if warnings_enabled
        @warn(msg)
    end
end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
"""
Author: Alex Wise and collaborators (see header of each file)
Created: October 2021
Contact: https://github.com/RvSpectML/RvLineList.jl
"""

module RvLineList

using JLD2, FileIO #objects used: save()

using CSV #objects used: CSV

using DataFrames #objects used: DataFrame

using Query #objects used: @filter

using RvSpectMLBase #objects used: read_data_paths

using EchelleInstruments #objects used: NEID, EXPRES, HARPSN
export NEID, EXPRES, HARPS, HARPSN

import EchelleCCFs:λ_air_to_vac #objects used: EchelleCCFs.λ_air_to_vac
#λ_air_to_vac = EchelleCCFs.λ_air_to_vac
#export λ_air_to_vac

using RvSpectML #objects used: extract_orders

using Dates #objects used: Dates.date
#filter datetime objects to get datetime list to use
#solar noon should use shift by longitude/15 to get from UTC to the thing we want

import Printf #objects used: @sprintf

using Statistics #objects used: quantile

using PyCall #makes it possible to call python functions from make_VALD_line_list.py in a julia environment
import Pandas.DataFrame as pd_df #used to convert julia DataFrame to python pandas.DataFrame

#mask manipulation
include("masks.jl")
export read_mask_air, read_mask_vacuum, binMask
export get_param_range, get_lambda_range, get_depth_range

#normaliation of EXPRES spectra - this file will not be needed once these functions are replaced with working versions from EchelleInstruments
include("expres_norm.jl")

#empirical mask generation
include("empirical_line_lists.jl")
export generateEmpiricalMask

# Packages we'll use in many places, but try to keep these minimal.
#using LinearAlgebra
#using NaNMath

end

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
function support_function_test(p::Rep)
    @test support_function([1, 1], p) ≈ 1
    @test support_function([1, 0], p) ≈ 1
    @testset "Scale" begin
        @test support_function([1, 1], 2p) ≈ 2
        @test support_function([1, 0], 2p) ≈ 2
    end
end
function support_function_test(lib::Polyhedra.Library)
    v = convexhull([1, 0], [0, 1], [-1, -1])
    @testset "V-representation" begin
        support_function_test(v)
    end
    p = polyhedron(v, lib)
    @testset "Polyhedron with V-representation" begin
        support_function_test(p)
    end
    @test !hrepiscomputed(p)
    h = hrep(polyhedron(v, lib))
    # We don't test just with the H-rep since it needs a solver
    p = polyhedron(h, lib)
    @testset "Polyhedron with H-representation" begin
        support_function_test(p)
    end
    @test !vrepiscomputed(p)
end