text:
"""
    AbstractBridgeOptimizer

A bridge optimizer applies given constraint bridges to a given optimizer thus
extending the types of supported constraints. The attributes of the inner
optimizer are automatically transformed to make the bridges transparent, e.g.
the variables and constraints created by the bridges are hidden.

By convention, the inner optimizer should be stored in a `model` field and
the dictionary mapping constraint indices to bridges should be stored in a
`bridges` field. If a bridge optimizer deviates from these conventions, it
should implement the functions `MOI.optimize!` and `bridge` respectively.
"""
abstract type AbstractBridgeOptimizer <: MOI.AbstractOptimizer end

# AbstractBridgeOptimizer interface

"""
    is_bridged(b::AbstractBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
              S::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether `b` tries to bridge `F`-in-`S` constraints
instead of passing it as is to its internal model.
"""
function is_bridged end
# Syntactic sugar
function is_bridged(b::AbstractBridgeOptimizer, ::Type{CI{F, S}}) where {F, S}
    return is_bridged(b, F, S)
end
# We don't bridge variables.
is_bridged(b::AbstractBridgeOptimizer, ::Type{VI}) = false

"""
    supports_bridging_constraint(b::AbstractBridgeOptimizer,
                               F::Type{<:MOI.AbstractFunction},
                               S::Type{<:MOI.AbstractSet})::Bool

Return a `Bool` indicating whether `b` supports bridging `F`-in-`S` constraints.
"""
function supports_bridging_constraint(::AbstractBridgeOptimizer,
                                    ::Type{<:MOI.AbstractFunction},
                                    ::Type{<:MOI.AbstractSet})
    return false
end

"""
    bridge_type(b::AbstractBridgeOptimizer,
                F::Type{<:MOI.AbstractFunction},
                S::Type{<:MOI.AbstractSet})

Return the `AbstractBridge` type to be used to bridge `F`-in-`S` constraints.
This function should only be called if `is_bridged(b, F, S)`.
"""
function bridge_type end

function concrete_bridge_type(b::AbstractBridgeOptimizer,
                              F::Type{<:MOI.AbstractFunction},
                              S::Type{<:MOI.AbstractSet})
    return concrete_bridge_type(bridge_type(b, F, S), F, S)
end

"""
    bridge(b::AbstractBridgeOptimizer, ci::CI)

Return the `AbstractBridge` used to bridge the constraint with index `ci`.
"""
bridge(b::AbstractBridgeOptimizer, ci::CI) = b.bridges[ci]
# By convention, they should be stored in a `bridges` field using a
# dictionary-like object.

# Implementation of the MOI interface for AbstractBridgeOptimizer

MOI.optimize!(b::AbstractBridgeOptimizer) = MOI.optimize!(b.model)
# By convention, the model should be stored in a `model` field

function MOI.is_empty(b::AbstractBridgeOptimizer)
    return isempty(b.bridges) && MOI.is_empty(b.model)
end
function MOI.empty!(b::AbstractBridgeOptimizer)
    MOI.empty!(b.model)
    MOI.empty!(b.bridged)
    empty!(b.bridges)
end
function MOI.supports(b::AbstractBridgeOptimizer,
                      attr::Union{MOI.AbstractModelAttribute,
                                  MOI.AbstractOptimizerAttribute})
    return MOI.supports(b.model, attr)
end
function MOI.copy_to(b::AbstractBridgeOptimizer, src::MOI.ModelLike;
                   copy_names = true)
    return MOIU.default_copy_to(b, src, copy_names)
end

# References
MOI.is_valid(b::AbstractBridgeOptimizer, vi::VI) = MOI.is_valid(b.model, vi)
function MOI.is_valid(b::AbstractBridgeOptimizer, ci::CI)
    if is_bridged(b, typeof(ci))
        return MOI.is_valid(b.bridged, ci)
    else
        return MOI.is_valid(b.model, ci)
    end
end
MOI.delete(b::AbstractBridgeOptimizer, vi::VI) = MOI.delete(b.model, vi)
function MOI.delete(b::AbstractBridgeOptimizer, ci::CI)
    if is_bridged(b, typeof(ci))
        if !MOI.is_valid(b, ci)
            throw(MOI.InvalidIndex(ci))
        end
        MOI.delete(b, bridge(b, ci))
        delete!(b.bridges, ci)
        MOI.delete(b.bridged, ci)
    else
        MOI.delete(b.model, ci)
    end
end

# Attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    if is_bridged(b, F, S)
        list = MOI.get(b.bridged, attr)
    else
        list = MOI.get(b.model, attr)
    end
    for bridge in values(b.bridges)
        for c in MOI.get(bridge, attr)
            i = something(findfirst(isequal(c), list), 0)
            if !iszero(i)
                MOI.deleteat!(list, i)
            end
        end
    end
    return list
end
function _numberof(b::AbstractBridgeOptimizer, model::MOI.ModelLike,
                   attr::Union{MOI.NumberOfConstraints, MOI.NumberOfVariables})
    s = MOI.get(model, attr)
    for v in values(b.bridges)
        s -= MOI.get(v, attr)
    end
    return s
end
function MOI.get(b::AbstractBridgeOptimizer, attr::MOI.NumberOfVariables)
    return _numberof(b, b.model, attr)
end
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.NumberOfConstraints{F, S}) where {F, S}
    if is_bridged(b, F, S)
        # The constraints contained in `b.bridged` may have been added by
        # bridges
        return _numberof(b, b.bridged, attr)
    else
        return _numberof(b, b.model, attr)
    end
end
function MOI.get(b::AbstractBridgeOptimizer, attr::MOI.ListOfConstraints)
    list_of_types = [MOI.get(b.model, attr); MOI.get(b.bridged, attr)]
    # Some constraint types show up in `list_of_types` even when all the
    # constraints of that type have been created by bridges and not by the user.
    # The code in `NumberOfConstraints` takes care of removing these constraints
    # from the counter so we can rely on it to remove these constraint types.
    types_to_remove = findall(iszero.(
        map(FS -> MOI.get(b, MOI.NumberOfConstraints{FS...}()), list_of_types)))
    deleteat!(list_of_types, types_to_remove)
    return list_of_types
end

# Model an optimizer attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::Union{MOI.AbstractModelAttribute,
                             MOI.AbstractOptimizerAttribute})
    return MOI.get(b.model, attr)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::Union{MOI.AbstractModelAttribute,
                              MOI.AbstractOptimizerAttribute},
                  value)
    return MOI.set(b.model, attr, value)
end

# Variable attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractVariableAttribute,
                 index::VI)
    return MOI.get(b.model, attr, index)
end
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractVariableAttribute,
                 indices::Vector{VI})
    return MOI.get(b.model, attr, indices)
end
function MOI.supports(b::AbstractBridgeOptimizer,
                      attr::MOI.AbstractVariableAttribute,
                      IndexType::Type{<:MOI.Index})
    return MOI.supports(b.model, attr, IndexType)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::MOI.AbstractVariableAttribute,
                  index::MOI.Index, value)
    return MOI.set(b.model, attr, index, value)
end
function MOI.set(b::AbstractBridgeOptimizer,
                  attr::MOI.AbstractVariableAttribute,
                  indices::Vector{<:MOI.Index}, values::Vector)
    return MOI.set(b.model, attr, indices, values)
end

# Constraint attributes
function MOI.get(b::AbstractBridgeOptimizer,
                 attr::MOI.AbstractConstraintAttribute,
                 ci::CI)
    if is_bridged(b, typeof(ci))
        if MOI.is_set_by_optimize(attr)
            MOI.get(b, attr, bridge(b, ci))
        else
            MOI.get(b.bridged, attr, ci)
        end
    else
        MOI.get(b.model, attr, ci)
    end
end
## Setting names
function MOI.supports(b::AbstractBridgeOptimizer, attr::MOI.ConstraintName,
                      Index::Type{<:CI})
    if is_bridged(b, Index)
        return MOI.supports(b.bridged, attr, Index)
    else
        return MOI.supports(b.model, attr, Index)
    end
end
function MOI.set(b::AbstractBridgeOptimizer, attr::MOI.ConstraintName,
                  constraint_index::CI, name::String)
    if is_bridged(b, typeof(constraint_index))
        MOI.set(b.bridged, attr, constraint_index, name)
    else
        MOI.set(b.model, attr, constraint_index, name)
    end
end
## Setting functions and sets
function MOI.set(b::AbstractBridgeOptimizer, ::MOI.ConstraintSet,
                  constraint_index::CI{F, S}, set::S) where {F, S}
    if is_bridged(b, typeof(constraint_index))
        MOI.set(b, MOI.ConstraintSet(), bridge(b, constraint_index), set)
        MOI.set(b.bridged, MOI.ConstraintSet(), constraint_index, set)
    else
        MOI.set(b.model, MOI.ConstraintSet(), constraint_index, set)
    end
end
function MOI.set(b::AbstractBridgeOptimizer, ::MOI.ConstraintFunction,
                  constraint_index::CI{F, S}, func::F) where {F, S}
    if is_bridged(b, typeof(constraint_index))
        MOI.set(b, MOI.ConstraintFunction(), bridge(b, constraint_index), func)
        MOI.set(b.bridged, MOI.ConstraintFunction(), constraint_index, func)
    else
        MOI.set(b.model, MOI.ConstraintFunction(), constraint_index, func)
    end
end

# Name
function MOI.get(b::AbstractBridgeOptimizer, IdxT::Type{<:MOI.Index},
                 name::String)
    if is_bridged(b, IdxT)
        return MOI.get(b.bridged, IdxT, name)
    else
        return MOI.get(b.model, IdxT, name)
    end
end

# We have no information as to whether the constraint is in the bridge or the
# model. Therefore, we try the model first, and then the bridge if that fails.
function MOI.get(b::AbstractBridgeOptimizer, IdxT::Type{CI},
                 name::String)
    ci = MOI.get(b.model, IdxT, name)
    if ci === nothing
        return MOI.get(b.bridged, IdxT, name)
    else
        return ci
    end
end

# Constraints
function MOI.supports_constraint(b::AbstractBridgeOptimizer,
                                F::Type{<:MOI.AbstractFunction},
                                S::Type{<:MOI.AbstractSet})
    if is_bridged(b, F, S)
        return supports_bridging_constraint(b, F, S) &&
            MOI.supports_constraint(b.bridged, F, S)
    else
        return MOI.supports_constraint(b.model, F, S)
    end
end
function MOI.add_constraint(b::AbstractBridgeOptimizer, f::MOI.AbstractFunction,
                            s::MOI.AbstractSet)
    if is_bridged(b, typeof(f), typeof(s))
        ci = MOI.add_constraint(b.bridged, f, s)
        @assert !haskey(b.bridges, ci)
        b.bridges[ci] = concrete_bridge_type(b, typeof(f), typeof(s))(b, f, s)
        return ci
    else
        return MOI.add_constraint(b.model, f, s)
    end
end
function MOI.modify(b::AbstractBridgeOptimizer, ci::CI,
                     change::MOI.AbstractFunctionModification)
    if is_bridged(b, typeof(ci))
        MOI.modify(b, bridge(b, ci), change)
        MOI.modify(b.bridged, ci, change)
    else
        MOI.modify(b.model, ci, change)
    end
end

# Objective
function MOI.modify(b::AbstractBridgeOptimizer, obj::MOI.ObjectiveFunction,
                     change::MOI.AbstractFunctionModification)
    MOI.modify(b.model, obj, change)
end

# Variables
MOI.add_variable(b::AbstractBridgeOptimizer) = MOI.add_variable(b.model)
MOI.add_variables(b::AbstractBridgeOptimizer, n) = MOI.add_variables(b.model, n)


text:
struct PharmClass
    system::String
    value::String
end

@inline function Base.isless(x::PharmClass, y::PharmClass)
    x_system = x.system::String
    y_system = y.system::String
    if x_system == y_system
        x_value = x.value::String
        y_value = y.value::String
        return Base.isless(x_value, y_value)
    end
    return Base.isless(x_system, y_system)
end

struct PharmGraph{I, F}
    g::MetaGraphs.MetaDiGraph{I, F}
end

Base.@kwdef struct Config
    rxnrel::String
    rxnsat::String
    showprogress::Bool = true
end

ext: jl


text:
#!/usr/bin/env julia

"""
This program generates batch sricpt to run IGV
Loads BAM files from localhost tunnel
Reads Variant positions from WIG file for each chromosome
"""

fh = ARGS[1] #      File with list of .wig file names
file = open(fh)
wig = readlines(file)

fP = ARGS[2] # GBM-pairs.txt
fil2 = open(fP)
bamList = readlines(fil2)

here = pwd() # presently working directory

for line in bamList#[1:5]
    out = []
    target = chomp(line)
    bam = split(target, ',')
    tumor, normal = bam[1], bam[2]
    #rm("$here/PNG/$target", recursive=true)
    mkdir("$here/PNG/$target")  # Creating snapshot directory
    push!(out, "new")
    push!(out, "genome hg19")
    push!(out, "load http://localhost:8000/bam/$tumor.bam,$here/sm10p.wig,$here/germline.wig.tdf,http://localhost:8000/bam/$normal.bam")
    push!(out, "maxPanelHeight 1000")
    push!(out, "snapshotDirectory $here/PNG/$target")
    for file in wig
        file = chomp(file)
        ln = open(file)
        lines = readlines(ln)
        chr = split(file, '.')
        for po in lines[3:end]
            pos = split(po, '\t')
            position = "goto "*chr[1]*":"*pos[1]*"-"*pos[1]
            push!(out, position)
            push!(out, "sort position")
            push!(out, "snapshot")
        end
    end
    push!(out, "exit")
    writedlm("$target.igv.txt", out)
end


text:
module CostHeatmapVisualizer
const project_root_dir = realpath(joinpath(@__DIR__, "../.."))
unique!(push!(LOAD_PATH, realpath(joinpath(project_root_dir, "test/utils"))))

import LinearAlgebra
import FileIO
import TestDynamics
using VegaLite: @vlplot

export cost_viz

"Returns a cost mapping that can be evaluated at (x, y) to get the running cost for that position."
function position_cost_mapping(;
    x_sequence,
    player_idx,
    control_system,
    positional_cost_weights,
    x_lane_center = nothing,
    y_lane_center = nothing,
    prox_min_regularization = 0.01,
    fix_cost_weights = (; state_lane = 0.1),
)

    opponent_positions = let
        opponent_indices = filter(!=(player_idx), eachindex(control_system.subsystems))
        opponent_position_indices = map(opponent_indices) do jj
            TestDynamics.state_indices(control_system, jj)[1:2]
        end
        map(opponent_position_indices) do opp_position_idxs
            x_sequence[opp_position_idxs, :]
        end
    end

    function (px, py)
        pos_ego = [px, py]

        J̃_position = (;
            state_lane = let
                # Note the switched indices. This is on purpose. The lane cost for the lane along y acts
                # on the x position and vice versa.
                x_lane_cost = isnothing(x_lane_center) ? 0 : (py - x_lane_center)^2
                y_lane_cost = isnothing(y_lane_center) ? 0 : (px - y_lane_center)^2
                x_lane_cost + y_lane_cost
            end,
            state_proximity = sum(opponent_positions) do pos_opponent_sequence
                sum(eachcol(pos_opponent_sequence)) do pos_opponent
                    -log(LinearAlgebra.norm_sqr(pos_ego - pos_opponent) + prox_min_regularization)
                end
            end,
        )
        sum(w * J̃_position[k] for (k, w) in pairs(positional_cost_weights)) +
        sum(w * J̃_position[k] for (k, w) in pairs(fix_cost_weights))
    end

end

"Generates a visualiztion of the inferred positional stage for a player."
function cost_viz(
    player_idx,
    player_config,
    ;
    x_sequence,
    control_system = control_system,
    show_y_label = false,
)
    cost_map = position_cost_mapping(;
        x_sequence,
        player_idx,
        control_system,
        y_lane_center = player_config.target_lane,
        positional_cost_weights = (; state_proximity = player_config.prox_cost),
    )

    x_domain = -1:0.09:2
    y_domain = -4:0.09:14

    max_size = 500
    x_position_domain = x_domain
    y_position_domain = y_domain
    x_range = only(diff(extrema(x_position_domain) |> collect))
    y_range = only(diff(extrema(y_position_domain) |> collect))
    max_range = max(x_range, y_range)
    canvas =
        @vlplot(width = max_size * x_range / max_range, height = max_size * y_range / max_range)

    normalized_cost = let
        data = [(; x, y, cost = cost_map(x, y)) for x in x_domain, y in y_domain]
        cmin, cmax = extrema(d.cost for d in data)
        [(; d.x, d.y, cost = (d.cost - cmin) / (cmax - cmin)) for d in data]
    end

    normalized_cost |>
    canvas + @vlplot(
        mark = {"point", shape = "square", filled = true, strokeWidth = 0},
        x = {"x:q", title = "Position x [m]"},
        y = {"y:q", title = "Position y [m]", axis = show_y_label ? true : nothing},
        color = {
            "cost:q",
            scale = {scheme = "viridis"},
            legend = {gradientLength = 200},
            title = "Normalized Cost",
        },
    )
end

end


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

module ProtoBuf

import Base.show, Base.copy!

export writeproto, readproto, ProtoMeta, ProtoMetaAttribs, meta, protobuild
export filled, isfilled, isfilled_default, which_oneof, fillset, fillset_default, fillunset
export show, copy!, set_field, set_field!, get_field, clear, add_field, add_field!, has_field, isinitialized
export ProtoEnum, lookup, enumstr
export ProtoServiceException, ProtoRpcChannel, ProtoRpcController, MethodDescriptor, ServiceDescriptor, ProtoService,
       AbstractProtoServiceStub, GenericProtoServiceStub, ProtoServiceStub, ProtoServiceBlockingStub,
       find_method, get_request_type, get_response_type, get_descriptor_for_type, call_method

using Compat

if isless(Base.VERSION, v"0.4.0-")
import Base.rsplit
rsplit{T<:AbstractString}(str::T, splitter; limit::Integer=0, keep::Bool=true) = rsplit(str, splitter, limit, keep)
end

if isless(Base.VERSION, v"0.4.0-")
fld_type(o, fld) = fieldtype(o, fld)
fld_names(x) = x.names
else
fld_type{T}(o::T, fld) = fieldtype(T, fld)
fld_names(x) = x.name.names
end

if isless(Base.VERSION, v"0.5.0-")
byte2str(x) = bytestring(x)
else
byte2str(x) = String(x)
end

# enable logging only during debugging
macro logmsg(s)
end
#macro logmsg(s)
#    quote
#        open("/tmp/protobuf.log", "a") do f
#            println(f, $(esc(s)))
#        end
#    end
#end

include("codec.jl")
include("svc.jl")

include("google/google.jl")

include("gen.jl")
include("utils.jl")

end # module

# Include Google ProtoBuf well known types (https://developers.google.com/protocol-buffers/docs/reference/google.protobuf).
# These are part of the `google.protobuf` module and are included automatically by the code generator.
# For hand coded modules, include them with: `using ProtoBuf; using google.protobuf`.
include("google/wellknown.jl")


text:
# Interface
const SerialZmqRelay = Base.Process

struct SerialRelayError <: Exception
    msg::String
end

SERIAL_PORTS = Dict{String,SerialZmqRelay}()

function launch_relay(
    port_name::String,
    baudrate::Int64,
    sub_endpoint::String,
    pub_endpoint::String,
)::SerialZmqRelay

    if port_name in keys(SERIAL_PORTS) && process_running(SERIAL_PORTS[port_name])
        @info "Serial port: $(port_name), already has associate ZMQ relay running!"
        return SERIAL_PORTS[port_name]
    end

    @info "Launching serial relay for: $(port_name), with write port: $(sub_endpoint) and read port: $(pub_endpoint)"
    relay_exe = joinpath(dirname(pathof(Mercury)), "..", "deps", "build", "relay_launch")
    cmd = `$relay_exe $port_name $baudrate $sub_endpoint $pub_endpoint`

    err = Pipe()
    serial_relay = run(pipeline(cmd, stderr = err), wait = false)
    close(err.in)
    sleep(1.0) # Wait for c constructor in relay_launch to either run successfully or fail

    if !process_running(serial_relay)
        err_msg = String(read(err))
        throw(SerialRelayError(err_msg))
    end

    SERIAL_PORTS[port_name] = serial_relay
    return serial_relay
end

function check_relays()
    if port_name in keys(SERIAL_PORTS)
        if !process_running(SERIAL_PORTS[port_name])
            @warn "Serial-ZMQ relay running on serialport: $port_name has stopped running!"
            delete!(SERIAL_PORTS, port_name)
        end
    end
end

function Base.close(serial_relay::SerialZmqRelay)
    if process_running(serial_relay) # true
        @info "Closing down Serial-ZMQ relay."
        kill(serial_relay)
    end
end

# TODO: Check if there is a way to capture STDOUT/STDERR from process without waiting
function check_relay_running(serial_relay::SerialZmqRelay)
    if !process_running(serial_relay) # true
        @warn "Serial-ZMQ relay shut down!."
        throw(ProcessFailedException)
    end
end

function closeall_relays()
    for port_name in keys(SERIAL_PORTS)
        @info "Closing down Serial-ZMQ relay running on serialport: $port_name"

        close(SERIAL_PORTS[port_name])
        delete!(SERIAL_PORTS, port_name)
    end
end

# Make sure all processes are cleaned up
Base.atexit(closeall_relays)


text:
## Defines all cipher functionality
## As usual, check out
#http://www.lysator.liu.se/~nisse/nettle/nettle.html#Cipher-functions

import Base: show
export CipherType, get_cipher_types, Encryptor, Decryptor, decrypt, decrypt!, encrypt, encrypt!

# This is a mirror of the nettle-meta.h:nettle_cipher struct
immutable NettleCipher
    name::Ptr{UInt8}
    context_size::Cuint
    block_size::Cuint
    key_size::Cuint
    set_encrypt_key::Ptr{Void}
    set_decrypt_key::Ptr{Void}
    encrypt::Ptr{Void}
    decrypt::Ptr{Void}
end

# For much the same reasons as in hash_common.jl, we define a separate, more "Julia friendly" type
immutable CipherType
    name::AbstractString
    context_size::Cuint
    block_size::Cuint
    key_size::Cuint
    set_encrypt_key::Ptr{Void}
    set_decrypt_key::Ptr{Void}
    encrypt::Ptr{Void}
    decrypt::Ptr{Void}
end

# These are the user-facing types that are used to actually {en,de}cipher stuff
immutable Encryptor
    cipher_type::CipherType
    state::Array{UInt8,1}
end
immutable Decryptor
    cipher_type::CipherType
    state::Array{UInt8,1}
end

# The function that maps from a NettleCipher to a CipherType
function CipherType(nc::NettleCipher)
    CipherType( uppercase(bytestring(nc.name)),
                nc.context_size, nc.block_size, nc.key_size,
                nc.set_encrypt_key, nc.set_decrypt_key, nc.encrypt, nc.decrypt)
end

# The global dictionary of hash types we know how to construct
const _cipher_types = Dict{AbstractString,CipherType}()

# We're going to load in each NettleCipher struct individually, deriving
# HashAlgorithm types off of the names we find, and storing the output
# and context size from the data members in the C structures
function get_cipher_types()
    # If we have already gotten the hash types from libnettle, don't query again
    if isempty(_cipher_types)
        cipher_idx = 1
        # nettle_ciphers is an array of pointers ended by a NULL pointer, continue reading hash types until we hit it
        while( true )
            ncptr = unsafe_load(cglobal(("nettle_ciphers",nettle),Ptr{Ptr{Void}}),cipher_idx)
            if ncptr == C_NULL
                break
            end
            cipher_idx += 1
            nc = unsafe_load(convert(Ptr{NettleCipher}, ncptr))
            cipher_type = CipherType(nc)
            _cipher_types[cipher_type.name] = cipher_type
        end
    end
    return _cipher_types
end

function Encryptor(name::AbstractString, key)
    cipher_types = get_cipher_types()
    name = uppercase(name)
    if !haskey(cipher_types, name)
        throw(ArgumentError("Invalid cipher type $name: call Nettle.get_cipher_types() to see available list"))
    end
    cipher_type = cipher_types[name]

    if length(key) != cipher_type.key_size
        throw(ArgumentError("Key must be $(cipher_type.key_size) bytes long"))
    end

    state = Array(UInt8, cipher_type.context_size)
    if nettle_major_version >= 3
        ccall( cipher_type.set_encrypt_key, Void, (Ptr{Void}, Ptr{UInt8}), state, pointer(key))
    else
        ccall( cipher_type.set_encrypt_key, Void, (Ptr{Void}, Cuint, Ptr{UInt8}), state, length(key), pointer(key))
    end

    return Encryptor(cipher_type, state)
end

function Decryptor(name::AbstractString, key)
    cipher_types = get_cipher_types()
    name = uppercase(name)
    if !haskey(cipher_types, name)
        throw(ArgumentError("Invalid cipher type $name: call Nettle.get_cipher_types() to see available list"))
    end
    cipher_type = cipher_types[name]

    if length(key) != cipher_type.key_size
        throw(ArgumentError("Key must be $(cipher_type.key_size) bytes long"))
    end

    state = Array(UInt8, cipher_type.context_size)
    if nettle_major_version >= 3
        ccall( cipher_type.set_decrypt_key, Void, (Ptr{Void}, Ptr{UInt8}), state, pointer(key))
    else
        ccall( cipher_type.set_decrypt_key, Void, (Ptr{Void}, Cuint, Ptr{UInt8}), state, length(key), pointer(key))
    end

    return Decryptor(cipher_type, state)
end

function decrypt!(state::Decryptor, result, data)
    if length(result) < length(data)
        throw(ArgumentError("Output array of length $(length(result)) insufficient for input data length ($(length(data)))"))
    end
    ccall(state.cipher_type.decrypt, Void, (Ptr{Void},Csize_t,Ptr{UInt8},Ptr{UInt8}),
        state.state, sizeof(data), pointer(result), pointer(data))
    return result
end

function decrypt(state::Decryptor, data)
    result = Array(UInt8, length(data))
    decrypt!(state, result, data)
    return result
end

function encrypt!(state::Encryptor, result, data)
    if length(result) < length(data)
        throw(ArgumentError("Output array of length $(length(result)) insufficient for input data length ($(length(data)))"))
    end
    ccall(state.cipher_type.encrypt, Void, (Ptr{Void},Csize_t,Ptr{UInt8},Ptr{UInt8}),
        state.state, sizeof(data), pointer(result), pointer(data))
    return result
end

function encrypt(state::Encryptor, data)
    result = Array(UInt8, length(data))
    encrypt!(state, result, data)
    return result
end

# The one-shot functions that make this whole thing so easy
decrypt(name::AbstractString, key, data) = decrypt(Decryptor(name, key), data)
encrypt(name::AbstractString, key, data) = encrypt(Encryptor(name, key), data)

# Custom show overrides make this package have a little more pizzaz!
function show(io::IO, x::CipherType)
    write(io, "$(x.name) Cipher\n")
    write(io, "  Context size: $(x.context_size) bytes\n")
    write(io, "  Block size: $(x.block_size) bytes\n")
    write(io, "  Key size: $(x.key_size) bytes")
end
show(io::IO, x::Encryptor) = write(io, "$(x.cipher_type.name) Encryption state")
show(io::IO, x::Decryptor) = write(io, "$(x.cipher_type.name) Decryption state")


text:
function Base.findmax(in::AFArray, dim::Int)
  dim > 0 || throw("Invalid dimension: $(dim).")
  out = af_array[0]
  idx = af_array[0]
  af_imax(out, idx, in, dim-1)
  AFArray(out[1]), AFArray(idx[1])
end

function Base.maximum(in::AFArray, dim::Int)
  dim > 0 || throw("Invalid dimension: $(dim).")
  out = af_array[0]
  af_max(out, in, dim-1)
  AFArray(out[1])
end

function Base.sum(in::AFArray, dim::Int)
  #dim > 0 || throw("Invalid dimension: $(dim).")
  out = af_array[0]
  af_sum(out, in, dim-1)
  AFArray(out[1])
end


text:
Base.Array{T,N}(dims::Base.OneTo...) where {T,N} = Array{T,N}(map(x -> x.stop, dims))

domain(x::AbstractArray, n) = indices(x, n)

function λopen(v::IVertex)
  prewalk(v) do v
    isa(value(v), Lambda) ? DataFlow.λopen(v) : v
  end
end

function λclose(v::IVertex)
  postwalk(v) do v
    isa(value(v), OLambda) ? DataFlow.λclose(v) : v
  end
end

function applybody(f, v::IVertex)
  @assert v.value isa Lambda
  vertex(Lambda(v.value.args, f(v.value.body)), v.inputs...)
end

function prewalkλ(f, v::IVertex)
  prewalk(v) do v
    v = f(v)
    v.value isa Lambda ? applybody(v -> prewalkλ(f, v), v) : v
  end
end

withopen(f, v) = v |> λopen |> f |> λclose

function dependents(v::IVertex)
  deps = ObjectIdDict()
  deps[v] = []
  DataFlow.prefor(v) do v
    foreach(w -> push!(get!(deps, w, []), v), v.inputs)
  end
  return deps
end

vcall(args...) = vertex(Call(), constant.(args)...)

# After creating a lambda, make sure the types of closed-over variables are
# propagated inside the new scope.
function closure_types(v)
  vars = [withtype(inputnode(i), vtype(v[i])) for i = 1:length(v.inputs)]
  vertex(Lambda(v.value.args, spliceinputs(v.value.body, vars...)),
         v.inputs...)
end

function tolambda(v::IVertex, args::IVertex...)
  λ = OLambda(length(args))
  is = [vertex(DataFlow.Split(i), constant(DataFlow.LooseEnd(λ.id)))
        for i = 1:length(args)]
  vertex(λ, DataFlow.postwalk!(v) do v
    get(is, findfirst(args, v), v)
  end) |> DataFlow.λclose |> closure_types
end

tolambda(v::IVertex, n::Integer) = vertex(Lambda(n, v))

max_stars_repo_name: JuliaTagBot/Tokamak

max_issues_repo_name: JuliaTagBot/Tokamak

max_forks_repo_name: FluxML/Tokamak


text:
module Millboard

__precompile__(true)

include("types.jl")

export table
include("table.jl")

export colored
include("colored.jl")

end # module


text:
###############
# API methods #
###############

"""
    ForwardDiff.hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x))

Return `H(f)` (i.e. `J(∇(f))`) evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Real)`.
"""
function hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x))
    ∇f = y -> gradient(f, y, cfg.gradient_config)
    return jacobian(∇f, x, cfg.jacobian_config)
end

"""
    ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x))

Compute `H(f)` (i.e. `J(∇(f))`) evaluated at `x` and store the result(s) in `result`,
assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Real)`.
"""
function hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x))
    ∇f = y -> gradient(f, y, cfg.gradient_config)
    jacobian!(result, ∇f, x, cfg.jacobian_config)
    return result
end

"""
    ForwardDiff.hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, result, x))

Exactly like `ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig)`, but
because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, result, x)` instead of
`HessianConfig(f, x)`.
"""
function hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, result, x))
    ∇f! = (y, z) -> begin
        inner_result = DiffResult(zero(eltype(y)), y)
        gradient!(inner_result, f, z, cfg.gradient_config)
        result = DiffResults.value!(result, value(DiffResults.value(inner_result)))
        return y
    end
    jacobian!(DiffResults.hessian(result), ∇f!, DiffResults.gradient(result), x, cfg.jacobian_config)
    return result
end

hessian(f, x::SArray) = jacobian(y -> gradient(f, y), x)

hessian(f, x::SArray, cfg::HessianConfig) = hessian(f, x)

hessian!(result::AbstractArray, f, x::SArray) = jacobian!(result, y -> gradient(f, y), x)

hessian!(result::MutableDiffResult, f, x::SArray) = hessian!(result, f, x, HessianConfig(f, result, x))

hessian!(result::ImmutableDiffResult, f, x::SArray, cfg::HessianConfig) = hessian!(result, f, x)

function hessian!(result::ImmutableDiffResult, f::F, x::SArray{S,V}) where {F,S,V}
    TJ = typeof(Tag((f,gradient),V))
    d1 = dualize(TJ, x)
    TG = typeof(Tag(f, eltype(d1)))
    d2 = dualize(TG, d1)
    fd2 = f(d2)
    val = value(TJ,value(TG,fd2))
    grad = extract_gradient(TJ,value(TG,fd2), x)
    hess = extract_jacobian(TJ,partials(TG,fd2), x)
    result = DiffResults.hessian!(result, hess)
    result = DiffResults.gradient!(result, grad)
    result = DiffResults.value!(result, val)
    return result
end


text:
# modified from https://github.com/SciML/DiffEqFlux.jl/blob/ccc3dcd0ccb7c18a176fc30dc9cb964acf98d50b/src/neural_de.jl#L38
# solvers only know how to work with arrays, so modifies to let solver
# know we're interested in error tolerances on the node feature vectors

struct GraphNeuralODE{M,P,RE,T,A,K} <: DiffEqFlux.NeuralDELayer
    model::M
    p::P
    re::RE
    tspan::T
    args::A
    kwargs::K

    function GraphNeuralODE(model,tspan,args...;p = nothing,kwargs...)
        _p,re = Flux.destructure(model)
        if p === nothing
            p = _p
        end
        new{typeof(model),typeof(p),typeof(re),
            typeof(tspan),typeof(args),typeof(kwargs)}(
            model,p,re,tspan,args,kwargs)
    end
end

function (n::GraphNeuralODE)(g::AbstractGraphTuple,p=n.p,reltol=1e-3,abstol=1e-3,save_everystep=false,save_start=false)
# function (n::GraphNeuralODE)(g::AbstractGraphTuple,p=n.p,save_everystep=false,save_start=false)
    function dudt_(u,p,t)
		g = (@set g.nodes = u)
		(n.re(p)(g)).nodes
	end

    ff = ODEFunction{false}(dudt_,tgrad=x->zero(x))
    prob = ODEProblem{false}(ff,g.nodes,n.tspan,p)
    sense = InterpolatingAdjoint(autojacvec=ZygoteVJP())
    @set g.nodes = solve(prob,n.args...;sense=sense,reltol=reltol,abstol=abstol,save_everystep=save_everystep,save_start=save_start,n.kwargs...).u[end]
    # @set g.nodes = solve(prob,n.args...;sense=sense,save_everystep=save_everystep,save_start=save_start,n.kwargs...).u[end]
end


text:
#
println("Cb) Test of the Hfs module for HFS A,B parameters and hyperfine representation with ASF from an internally generated multiplet.")
#
setDefaults("unit: energy", "Hz")
setDefaults("print summary: open", "zzz-Hfs.sum")
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=JAC.Radial.Grid(true),
                        nuclearModel=Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0), properties=[HFS()],
                        # grid=JAC.Radial.Grid(true),
                        configs=[Configuration("[Ne] 3s"), Configuration("[Ne] 3p"), Configuration("[Ne] 3d")],
                        hfsSettings=Hfs.Settings(true, true, false, false, false, false, false, Int64[] ) )

wb = perform(wa)
setDefaults("print summary: close", "")

#==
wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0); properties=[HFS()],
                        # grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-6, h = 5.0e-2, hp = 0., NoPoints = 600),
                        configs=[Configuration("[Ne] 3s"), Configuration("[Ne] 3p"), Configuration("[Ne] 3d")],
                        hfsSettings=Hfs.Settings(true, true, false, false, false, false, false, Int64[] ) )

wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0); properties=[HFS()],
                        # grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-6, h = 5.0e-2, hp = 0., NoPoints = 600),
                        configs=[Configuration("[Ne] 3s^2"), Configuration("[Ne] 3s 3p"), Configuration("[Ne] 3p^2")],
                        hfsSettings=Hfs.Settings(true, true, false, false, false, false, false, Int64[] ) )

wa = Atomic.Computation("xx",  Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0); properties=[HFS()],
                        # grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-6, h = 5.0e-2, hp = 0., NoPoints = 600),
                        configs=[Configuration("[Ne] 3s^2 3d"), Configuration("[Ne] 3d^3")],
                        hfsSettings=Hfs.Settings(true, true, false, false, false, false, false, Int64[] ) )
==#


text:
# The format of simplegraph files is as follows:
# a one line header: <num_vertices>, <num_edges>, {"d" | "u"}, <name>[, <ver>, <datatype>, <graphcode>]
#   - num_vertices is an integer
#   - num_edges is an integer
#   - "d" for directed graph, "u" for undirected. Note that this
#       option does not perform any additional edge construction; it's
#       merely used to return the correct type of graph.
#   - name is a string
#   - ver is optional and is an int
#   - datatype is mandatory if version is set and is a string ("UInt8", etc.)
#   - graphcode is mandatory if version is set and is a string.
# header followed by a list of (comma-delimited) edges - src,dst.
# Multiple graphs may be present in one file.

struct LGFormat <: AbstractGraphFormat end

struct LGHeader
    nv::Int
    ne::Int
    is_directed::Bool
    name::String
    ver::Int64
    dtype::DataType
    code::String
end
function show(io::IO, h::LGHeader)
    isdir = h.is_directed ? "d" : "u"
    print(io, "$(h.nv),$(h.ne),$isdir,$(h.name),$(h.ver),$(h.dtype),$(h.code)")
end

LGHeader(nv::Int, ne::Int, is_directed::Bool, name::AbstractString) =
    LGHeader(nv, ne, is_directed, name, 1, Int64, "simplegraph")

function _lg_read_one_graph(f::IO, header::LGHeader)
    T = header.dtype
    if header.is_directed
        g = DiGraph{T}(header.nv)
    else
        g = Graph{T}(header.nv)
    end
    for i = 1:header.ne
        line = chomp(readline(f))
        if length(line) > 0
            src_s, dst_s = split(line, r"\s*,\s*")
            src = parse(T, src_s)
            dst = parse(T, dst_s)
            add_edge!(g, src, dst)
        end
    end
    return g
end

function _lg_skip_one_graph(f::IO, n_e::Integer)
    for _ in 1:n_e
        readline(f)
    end
end

function _parse_header(s::AbstractString)
    addl_info = false
    nvstr, nestr, dirundir, graphname = split(s, r"s*,s*", limit=4)
    if occursin(",", graphname) # version number and type
        graphname, _ver, _dtype, graphcode = split(graphname, r"s*,s*")
        ver = parse(Int, _ver)
        dtype = eval(Symbol(_dtype))
        addl_info = true
    end
    n_v = parse(Int, nvstr)
    n_e = parse(Int, nestr)
    dirundir = strip(dirundir)
    directed = !(dirundir == "u")
    graphname = strip(graphname)
    if !addl_info
        header = LGHeader(n_v, n_e, directed, graphname)
    else
        header = LGHeader(n_v, n_e, directed, graphname, ver, dtype, graphcode)
    end
    return header
end

"""
    loadlg_mult(io)

Return a dictionary of (name=>graph) loaded from IO stream `io`.
"""
function loadlg_mult(io::IO)
    graphs = Dict{String,AbstractGraph}()
    while !eof(io)
        line = strip(chomp(readline(io)))
        if !(startswith(line, "#") || line == "")
            header = _parse_header(line)
            g = _lg_read_one_graph(io, header)
            graphs[header.name] = g
        end
    end
    return graphs
end

function loadlg(io::IO, gname::String)
    while !eof(io)
        line = strip(chomp(readline(io)))
        (startswith(line, "#") || line == "") && continue
        header = _parse_header(line)
        if gname == header.name
            return _lg_read_one_graph(io, header)
        else
            _lg_skip_one_graph(io, header.ne)
        end
    end
    throw(ArgumentError("graph $gname not found"))
end

"""
    savelg(io, g, gname)

Write a graph `g` with name `gname` in a proprietary format
to the IO stream designated by `io`. Return 1 (number of graphs written).
"""
function savelg(io::IO, g::AbstractGraph{T}, gname::String) where T
    header = LGHeader(nv(g), ne(g), is_directed(g), gname, 2, T, "simplegraph")
    # write header line
    line = string(header)
    write(io, "$line\n")
    # write edges
    for e in edges(g)
        write(io, "$(src(e)),$(dst(e))\n")
    end
    return 1
end

"""
    savelg_mult(io, graphs)

Write a dictionary of (name=>graph) to an IO stream `io`,
with default `GZip` compression. Return number of graphs written.
"""
function savelg_mult(io::IO, graphs::Dict)
    ng = 0
    for (gname, g) in graphs
        ng += savelg(io, g, gname)
    end
    return ng
end

loadgraph(io::IO, gname::String, ::LGFormat) = loadlg(io, gname)
loadgraphs(io::IO, ::LGFormat) = loadlg_mult(io)
savegraph(io::IO, g::AbstractGraph, gname::String, ::LGFormat) = savelg(io, g, gname)
savegraph(io::IO, g::AbstractGraph, ::LGFormat) = savelg(io, g, "graph")
savegraph(io::IO, d::Dict, ::LGFormat) = savelg_mult(io, d)


text:
module BooleanAD

using ChainRulesCore
using Lazy
using MacroTools
using Random
using Lazy: @forward
import Base: &, |, !, ==
using Flux
using Flux: @adjoint
using Zygote

export DBool, majorityN, orN, andN, kjuntaN, weightedMajorityN

include("rules.jl")
include("voting.jl")
include("order.jl")
include("utility.jl")

# @warn "Experimental package BooleanAD"

end # module

# Why is this even a point here.


text:
# Lazy concatenation of AbstractVector's.
# Similar to Iterators.Flatten and some code has been reused from julia/base/iterators.jl

const Vcat{T,N,I<:Tuple} = ApplyArray{T,N,typeof(vcat),I}

Vcat(A...) = ApplyArray(vcat, A...)
Vcat{T}(A...) where T = ApplyArray{T}(vcat, A...)
Vcat() = Vcat{Any}()

Vcat(A::AbstractVector...) = ApplyVector(vcat, A...)
Vcat{T}(A::AbstractVector...) where T = ApplyVector{T}(vcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(vcat)})
    isempty(A.args) && return A
    m = size(A.args[1],2)
    for k=2:length(A.args)
        size(A.args[k],2) == m || throw(ArgumentError("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(vcat)}) = promote_type(map(eltype,A.args)...)
@inline eltype(A::Applied{<:Any,typeof(vcat),Tuple{}}) = Any
@inline ndims(A::Applied{<:Any,typeof(vcat),I}) where I = max(1,maximum(map(ndims,A.args)))
@inline ndims(A::Applied{<:Any,typeof(vcat),Tuple{}}) = 1
@inline axes(f::Vcat{<:Any,1,Tuple{}}) = (OneTo(0),)
@inline axes(f::Vcat{<:Any,1}) = tuple(OneTo(+(map(length,f.args)...)))
@inline axes(f::Vcat{<:Any,2}) = (OneTo(+(map(a -> size(a,1), f.args)...)), axes(f.args[1],2))
Base.IndexStyle(::Type{<:Vcat{T,1}}) where T = Base.IndexLinear()
Base.IndexStyle(::Type{<:Vcat{T,2}}) where T = Base.IndexCartesian()

@propagate_inbounds @inline function vcat_getindex(f, k::Integer)
    T = eltype(f)
    κ = k
    for A in f.args
        n = length(A)
        κ ≤ n && return convert(T,A[κ])::T
        κ -= n
    end
    throw(BoundsError(f, k))
end

@propagate_inbounds @inline function vcat_getindex(f, k::Integer, j::Integer)
    T = eltype(f)
    κ = k
    for A in f.args
        n = size(A,1)
        κ ≤ n && return convert(T,A[κ,j])::T
        κ -= n
    end
    throw(BoundsError(f, (k,j)))
end

@propagate_inbounds @inline getindex(f::Vcat{<:Any,1}, k::Integer) = vcat_getindex(f, k)
@propagate_inbounds @inline getindex(f::Vcat{<:Any,2}, k::Integer, j::Integer) = vcat_getindex(f, k, j)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer)= vcat_getindex(f, k)
getindex(f::Applied{<:Any,typeof(vcat)}, k::Integer, j::Integer)= vcat_getindex(f, k, j)

@propagate_inbounds @inline function setindex!(f::Vcat{T,1}, v, k::Integer) where T
    κ = k
    for A in f.args
        n = length(A)
        κ ≤ n && return setindex!(A, v, κ)
        κ -= n
    end
    throw(BoundsError(f, k))
end

@propagate_inbounds @inline function setindex!(f::Vcat{T,2}, v, k::Integer, j::Integer) where T
    κ = k
    for A in f.args
        n = size(A,1)
        κ ≤ n && return setindex!(A, v, κ, j)
        κ -= n
    end
    throw(BoundsError(f, (k,j)))
end

reverse(f::Vcat{<:Any,1}) = Vcat((reverse(itr) for itr in reverse(f.args))...)

const Hcat{T,I<:Tuple} = ApplyArray{T,2,typeof(hcat),I}

Hcat(A...) = ApplyArray(hcat, A...)
Hcat{T}(A...) where T = ApplyArray{T}(hcat, A...)

function instantiate(A::Applied{DefaultApplyStyle,typeof(hcat)})
    isempty(A.args) && return A
    m = size(A.args[1],1)
    for k=2:length(A.args)
        size(A.args[k],1) == m || throw(ArgumentError("number of rows of each array must match (got $(map(x->size(x,1), A)))"))
    end
    Applied{DefaultApplyStyle}(A.f,map(instantiate,A.args))
end

@inline eltype(A::Applied{<:Any,typeof(hcat)}) = promote_type(map(eltype,A.args)...)
ndims(::Applied{<:Any,typeof(hcat)}) = 2
size(f::Applied{<:Any,typeof(hcat)}) = (size(f.args[1],1), +(map(a -> size(a,2), f.args)...))
Base.IndexStyle(::Type{<:Hcat}) where T = Base.IndexCartesian()

function hcat_getindex(f, k::Integer, j::Integer)
    T = eltype(f)
    ξ = j
    for A in f.args
        n = size(A,2)
        ξ ≤ n && return T(A[k,ξ])::T
        ξ -= n
    end
    throw(BoundsError(f, (k,j)))
end

getindex(f::Hcat, k::Integer, j::Integer) = hcat_getindex(f, k, j)
getindex(f::Applied{DefaultArrayApplyStyle,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)
getindex(f::Applied{<:Any,typeof(hcat)}, k::Integer, j::Integer)= hcat_getindex(f, k, j)

function setindex!(f::Hcat{T}, v, k::Integer, j::Integer) where T
    ξ = j
    for A in f.args
        n = size(A,2)
        ξ ≤ n && return setindex!(A, v, k, ξ)
        ξ -= n
    end
    throw(BoundsError(f, (k,j)))
end

## copyto!
# based on Base/array.jl, Base/abstractarray.jl
copyto!(dest::AbstractArray, V::Vcat) = vcat_copyto!(dest, arguments(V)...)
function vcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest,1)
    nrows == sum(a->size(a, 1), arrays) || throw(DimensionMismatch("sum of rows each matrix must equal $nrows"))
    ncols = size(dest, 2)
    for a in arrays
        if size(a, 2) != ncols
            throw(DimensionMismatch("number of columns of each array must match (got $(map(x->size(x,2), A)))"))
        end
    end
    pos = 1
    for a in arrays
        p1 = pos+size(a,1)-1
        copyto!(view(dest,pos:p1, :), a)
        pos = p1+1
    end
    return dest
end

function vcat_copyto!(arr::AbstractVector, arrays...)
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    i = 0
    @inbounds for a in arrays
        m = length(a)
        copyto!(view(arr,i+1:i+m), a)
        i += m
    end
    arr
end

function vcat_copyto!(arr::Vector{T}, arrays::Vector{T}...) where T
    n = 0
    for a in arrays
        n += length(a)
    end
    n == length(arr) || throw(DimensionMismatch("destination must have length equal to sums of concatenated vectors"))

    ptr = pointer(arr)
    if isbitstype(T)
        elsz = Core.sizeof(T)
    elseif isbitsunion(T)
        elsz = bitsunionsize(T)
        selptr = convert(Ptr{UInt8}, ptr) + n * elsz
    else
        elsz = Core.sizeof(Ptr{Cvoid})
    end
    t = @_gc_preserve_begin arr
    for a in arrays
        na = length(a)
        nba = na * elsz
        if isbitstype(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
        elseif isbitsunion(T)
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  ptr, a, nba)
            # copy selector bytes
            ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt),
                  selptr, convert(Ptr{UInt8}, pointer(a)) + nba, na)
            selptr += na
        else
            ccall(:jl_array_ptr_copy, Cvoid, (Any, Ptr{Cvoid}, Any, Ptr{Cvoid}, Int),
                  arr, ptr, a, pointer(a), na)
        end
        ptr += nba
    end
    @_gc_preserve_end t
    return arr
end

copyto!(dest::AbstractMatrix, H::Hcat) = hcat_copyto!(dest, arguments(H)...)
function hcat_copyto!(dest::AbstractMatrix, arrays...)
    nargs = length(arrays)
    nrows = size(dest, 1)
    ncols = 0
    dense = true
    for a in arrays
        dense &= isa(a,Array)
        nd = ndims(a)
        ncols += (nd==2 ? size(a,2) : 1)
    end

    nrows == size(first(arrays),1) || throw(DimensionMismatch("Destination rows must match"))
    ncols == size(dest,2) || throw(DimensionMismatch("Destination columns must match"))

    pos = 1
    if dense
        for a in arrays
            n = length(a)
            copyto!(dest, pos, a, 1, n)
            pos += n
        end
    else
        for a in arrays
            p1 = pos+(isa(a,AbstractMatrix) ? size(a, 2) : 1)-1
            copyto!(view(dest,:, pos:p1), a)
            pos = p1+1
        end
    end
    return dest
end

function hcat_copyto!(dest::AbstractMatrix, arrays::AbstractVector...)
    height = size(dest, 1)
    for j = 1:length(arrays)
        if length(arrays[j]) != height
            throw(DimensionMismatch("vectors must have same lengths"))
        end
    end
    for j=1:length(arrays)
        copyto!(view(dest,:,j), arrays[j])
    end

    dest
end

#####
# adjoint/transpose
#####

for adj in (:adjoint, :transpose)
    @eval begin
        $adj(A::Hcat{T}) where T = Vcat{T}(map($adj,A.args)...)
        $adj(A::Vcat{T}) where T = Hcat{T}(map($adj,A.args)...)
    end
end

_vec(a) = a
_vec(a::AbstractArray) = vec(a)
_vec(a::Adjoint{<:Number,<:AbstractVector}) = _vec(parent(a))
vec(A::Hcat) = Vcat(_vec.(A.args)...)

_permutedims(a) = a
_permutedims(a::AbstractArray) = permutedims(a)

permutedims(A::Hcat{T}) where T = Vcat{T}(map(_permutedims,A.args)...)
permutedims(A::Vcat{T}) where T = Hcat{T}(map(_permutedims,A.args)...)

#####
# broadcasting
#
# We want broadcasting for numbers with concaenations to pass through
# to take advantage of special implementations of the sub-components
######

BroadcastStyle(::Type{<:Vcat{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(::Type{<:Hcat{<:Any}}) where N = LazyArrayStyle{2}()

broadcasted(::LazyArrayStyle, op, A::Vcat) =
    Vcat(broadcast(x -> broadcast(op, x), A.args)...)

for Cat in (:Vcat, :Hcat)
    @eval begin
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Number) =
            $Cat(broadcast((x,y) -> broadcast(op, x, y), A.args, c)...)
        broadcasted(::LazyArrayStyle, op, c::Number, A::$Cat) =
            $Cat(broadcast((x,y) -> broadcast(op, x, y), c, A.args)...)
        broadcasted(::LazyArrayStyle, op, A::$Cat, c::Ref) =
            $Cat(broadcast((x,y) -> broadcast(op, x, Ref(y)), A.args, c)...)
        broadcasted(::LazyArrayStyle, op, c::Ref, A::$Cat) =
            $Cat(broadcast((x,y) -> broadcast(op, Ref(x), y), c, A.args)...)
    end
end

# determine indices of components of a vcat
_vcat_axes(::Tuple{}) = (1,)
_vcat_axes(a::Tuple{<:AbstractUnitRange}) = (first(a),)
_vcat_axes(::Tuple{}, b, c...) = tuple(1, broadcast(x -> broadcast(+, 1, x), _vcat_axes(b, c...))...)
_vcat_axes(a::Tuple{<:AbstractUnitRange}, b, c...) = tuple(first(a), broadcast((α,x) -> broadcast(+, α, x), last(first(a)),
                                                            _vcat_axes(b, c...))...)

_vcat_getindex_eval(y) = ()
_vcat_getindex_eval(y, a, b...) = tuple(y[a], _vcat_getindex_eval(y, b...)...)

function broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::AbstractVector)
    kr = _vcat_axes(map(axes,A.args)...)  # determine how to break up B
    B_arrays = _vcat_getindex_eval(B,kr...)    # evaluate B at same chunks as A
    ApplyVector(vcat, broadcast((a,b) -> broadcast(op,a,b), A.args, B_arrays)...)
end

function broadcasted(::LazyArrayStyle, op, A::AbstractVector, B::Vcat{<:Any,1})
    kr = _vcat_axes(axes.(B.args)...)
    A_arrays = _vcat_getindex_eval(A,kr...)
    Vcat(broadcast((a,b) -> broadcast(op,a,b), A_arrays, B.args)...)
end

# Cannot broadcast Vcat's in a lazy way so stick to BroadcastArray
broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::Vcat{<:Any,1}) =
    Broadcasted{LazyArrayStyle}(op, (A, B))

# ambiguities
broadcasted(::LazyArrayStyle, op, A::Vcat{<:Any,1}, B::CachedVector) = cache_broadcast(op, A, B)
broadcasted(::LazyArrayStyle, op, A::CachedVector, B::Vcat{<:Any,1}) = cache_broadcast(op, A, B)

broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Vcat{<:Any,1}, b::Zeros{<:Any,1})=
    broadcast(DefaultArrayStyle{1}(), *, a, b)

function +(A::Vcat, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::Vcat, B::AbstractArray)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end
function +(A::AbstractArray, B::Vcat)
    size(A) == size(B) || throw(DimensionMismatch("dimensions must match."))
    A .+ B
end

######
# Special Vcat broadcasts
#
# We use Vcat for infinite padded vectors, so we need to special case
# two arrays. This may be generalisable in the future
######

function _vcat_broadcasted(::Type{T}, op, (Ahead, Atail)::Tuple{<:AbstractVector,<:AbstractFill},
                               (Bhead, Btail)::Tuple{<:AbstractVector,<:AbstractFill}) where T
    if length(Ahead) ≥ length(Bhead)
        M,m = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(view(Ahead,1:m), Bhead)
        view(Chead,m+1:M) .= op.(view(Ahead,m+1:M),Btail[1:M-m])

        Ctail = op.(Atail, Btail[M-m+1:end])
    else
        m,M = length(Ahead), length(Bhead)
        Chead = Vector{T}(undef,M)
        view(Chead,1:m) .= op.(Ahead, view(Bhead,1:m))
        view(Chead,m+1:M) .= op.(Atail[1:M-m],view(Bhead,m+1:M))

        Ctail = op.(Atail[M-m+1:end], Btail)
    end

    Vcat(Chead, Ctail)
end

_vcat_broadcasted(::Type{T}, op, (Ahead, Atail)::Tuple{<:Number,<:AbstractFill},
                           (Bhead, Btail)::Tuple{<:Number,<:AbstractFill}) where {M,T} =
   Vcat(op.(Ahead,Bhead), op.(Atail,Btail))

_vcat_broadcasted(::Type{T}, op, (Ahead, Atail)::Tuple{<:SVector{M},<:AbstractFill},
                           (Bhead, Btail)::Tuple{<:SVector{M},<:AbstractFill}) where {M,T} =
   Vcat(op.(Ahead,Bhead), op.(Atail,Btail))

# default is BroadcastArray
_vcat_broadcasted(::Type{T}, op, A, B) where T =
    Broadcasted{LazyArrayStyle}(op, (Vcat(A...), Vcat(B...)))

broadcasted(::LazyArrayStyle{1}, op, A::Vcat{T, 1, <:Tuple{<:Any,<:Any}},
                                     B::Vcat{V, 1, <:Tuple{<:Any,<:Any}}) where {T,V} =
  _vcat_broadcasted(promote_type(T,V), op, A.args, B.args)

####
# Cumsum
####

sum(V::Vcat) = mapreduce(sum, +, V.args)
all(V::Vcat) = all(all.(V.args))
any(V::Vcat) = any(any.(V.args))
all(f::Function, V::Vcat) = all(all.(f, V.args))
any(f::Function, V::Vcat) = any(any.(f, V.args))

_dotplus(a,b) = broadcast(+, a, b)

@inline _cumsum(x::Number) = x
@inline _cumsum(x) = cumsum(x)
@generated function _vcat_cumsum(x...)
    N = length(x)
    ret = quote
        @nexprs $N d->(c_d = _cumsum(x[d]))
        d_1 = c_1
        @nexprs $(N-1) k->(d_{k+1} = broadcast(+, last(d_k), c_{k+1}))
        @ntuple $N d
    end
end

@inline cumsum(V::Vcat{<:Any,1}) = ApplyVector(vcat,_vcat_cumsum(V.args...)...)

_vcat_diff(x::Number) = ()
_vcat_diff(x) = (diff(x),)

_vcat_diff(a::Number, b, c...) = (first(b)-a, _vcat_diff(b,c...)...)
_vcat_diff(a, b, c...) = (diff(a), first(b)-last(a), _vcat_diff(b,c...)...)
@inline diff(V::Vcat{T,1}) where T = ApplyVector{T}(vcat,_vcat_diff(V.args...)...)

####
# maximum/minimum
####

for op in (:maximum, :minimum)
    @eval $op(V::Vcat) = $op($op.(V.args))
end

function in(x, V::Vcat)
    for a in V.args
        in(x, a) && return true
    end
    false
end

_fill!(a, x) = fill!(a,x)
function _fill!(a::Number, x)
    a == x || throw(ArgumentError("Cannot set $a to $x"))
    a
end

function fill!(V::Union{Vcat,Hcat}, x)
    for a in V.args
        _fill!(a, x)
    end
    V
end

###
# *
###

function materialize!(M::MatMulVecAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 function materialize!(M::MatMulMatAdd{ApplyLayout{typeof(hcat)},ApplyLayout{typeof(vcat)}})
    α,A,B,β,C =  M.α,M.A,M.B,M.β,M.C
    T = eltype(C)
    _fill_lmul!(β,C) # this is temporary until strong β = false is supported
    for (a,b) in zip(A.args,B.args)
        materialize!(MulAdd(α,a,b,one(T),C))
    end
    C
 end

 ####
 # col/rowsupport
 ####

most(a) = reverse(tail(reverse(a)))
colsupport(M::Vcat, j) = first(colsupport(first(M.args),j)):(size(Vcat(most(M.args)...),1)+last(colsupport(last(M.args),j)))

###
# padded
####

struct PaddedLayout{L} <: MemoryLayout end
applylayout(::Type{typeof(vcat)}, ::A, ::ZerosLayout) where A = PaddedLayout{A}()
cachedlayout(::A, ::ZerosLayout) where A = PaddedLayout{A}()

paddeddata(A::CachedArray) = view(A.data,OneTo.(A.datasize)...)
paddeddata(A::Vcat) = A.args[1]

function ==(A::CachedVector{<:Any,<:Any,<:Zeros}, B::CachedVector{<:Any,<:Any,<:Zeros})
    length(A) == length(B) || return false
    n = max(A.datasize[1], B.datasize[1])
    resizedata!(A,n); resizedata!(B,n)
    view(A.data,OneTo(n)) == view(B.data,OneTo(n))
end

# special copyto! since `similar` of a padded returns a cached
for Typ in (:Number, :AbstractVector)
    @eval function copyto!(dest::CachedVector{T,Vector{T},<:Zeros{T,1}}, src::Vcat{<:Any,1,<:Tuple{<:$Typ,<:Zeros}}) where T
        length(src) ≤ length(dest)  || throw(BoundsError())
        a,_ = src.args
        n = length(a)
        resizedata!(dest, n) # make sure we are padded enough
        copyto!(view(dest.data,OneTo(n)), a)
        dest
    end
end

function copyto!(dest::CachedVector{T,Vector{T},<:Zeros{T,1}}, src::CachedVector{V,Vector{V},<:Zeros{V,1}}) where {T,V}
    length(src) ≤ length(dest)  || throw(BoundsError())
    n = src.datasize[1]
    resizedata!(dest, n)
    copyto!(view(dest.data,OneTo(n)), view(src.data,OneTo(n)))
    dest
end

struct Dot{StyleA,StyleB,ATyp,BTyp}
    A::ATyp
    B::BTyp
end

Dot(A::ATyp,B::BTyp) where {ATyp,BTyp} = Dot{typeof(MemoryLayout(ATyp)), typeof(MemoryLayout(BTyp)), ATyp, BTyp}(A, B)
materialize(d::Dot{<:Any,<:Any,<:AbstractArray,<:AbstractArray}) = Base.invoke(dot, Tuple{AbstractArray,AbstractArray}, d.A, d.B)
function materialize(d::Dot{<:PaddedLayout,<:PaddedLayout,<:AbstractVector{T},<:AbstractVector{V}}) where {T,V}
    a,b = paddeddata(d.A), paddeddata(d.B)
    m = min(length(a), length(b))
    convert(promote_type(T,V), dot(view(a,1:m), view(b,1:m)))
end

dot(a::CachedArray, b::AbstractArray) = materialize(Dot(a,b))
dot(a::LazyArray, b::AbstractArray) = materialize(Dot(a,b))

###
# norm
###

for Cat in (:Vcat, :Hcat)
    for (op,p) in ((:norm1,1), (:norm2,2), (:normInf,Inf))
        @eval $op(a::$Cat) = $op(norm.(a.args,$p))
    end
    @eval normp(a::$Cat, p) = norm(norm.(a.args, p), p)
end

###
# subarrays
###

sublayout(::ApplyLayout{typeof(vcat)}, _) = ApplyLayout{typeof(vcat)}()
sublayout(::ApplyLayout{typeof(hcat)}, _) = ApplyLayout{typeof(hcat)}()
# a row-slice of an Hcat is equivalent to a Vcat
sublayout(::ApplyLayout{typeof(hcat)}, ::Type{<:Tuple{Number,AbstractVector}}) = ApplyLayout{typeof(vcat)}()

arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Slice,<:Any}}) =
    view.(arguments(parent(V)), Ref(:), Ref(parentindices(V)[2]))
arguments(::ApplyLayout{typeof(hcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{<:Any,<:Slice}}) =
    view.(arguments(parent(V)), Ref(parentindices(V)[1]), Ref(:))

copyto!(dest::AbstractArray{T,N}, src::SubArray{T,N,<:Vcat{T,N}}) where {T,N} =
    vcat_copyto!(dest, arguments(ApplyLayout{typeof(vcat)}(), src)...)
copyto!(dest::AbstractMatrix{T}, src::SubArray{T,2,<:Hcat{T}}) where T =
    hcat_copyto!(dest, arguments(ApplyLayout{typeof(hcat)}(), src)...)

_vcat_lastinds(sz) = _vcat_cumsum(sz...)
_vcat_firstinds(sz) = (1, (1 .+ most(_vcat_lastinds(sz)))...)

_argsindices(sz) = broadcast(:, _vcat_firstinds(sz), _vcat_lastinds(sz))

_view_vcat(a::Number, kr) = Fill(a,length(kr))
_view_vcat(a::Number, kr, jr) = Fill(a,length(kr), length(jr))
_view_vcat(a, kr...) = view(a, kr...)

function _vcat_sub_arguments(::ApplyLayout{typeof(vcat)}, A, V)
    kr = parentindices(V)[1]
    sz = size.(arguments(A),1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(arguments(A), skr2)
end
_vcat_sub_arguments(::ApplyLayout{typeof(hcat)}, A, V) = arguments(ApplyLayout{typeof(hcat)}(), V)
_vcat_sub_arguments(::PaddedLayout, A, V) = _vcat_sub_arguments(ApplyLayout{typeof(vcat)}(), A, V)

_vcat_sub_arguments(A, V) = _vcat_sub_arguments(MemoryLayout(typeof(A)), A, V)
arguments(::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,1}) = _vcat_sub_arguments(parent(V), V)

function arguments(L::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2})
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,1)
    skr = intersect.(_argsindices(sz), Ref(kr))
    skr2 = broadcast((a,b) -> a .- b .+ 1, skr, _vcat_firstinds(sz))
    _view_vcat.(args, skr2, Ref(jr))
end

_view_hcat(a::Number, kr, jr) = Fill(a,length(kr),length(jr))
_view_hcat(a, kr, jr) = view(a, kr, jr)

function arguments(L::ApplyLayout{typeof(hcat)}, V::SubArray)
    A = parent(V)
    args = arguments(L, A)
    kr,jr = parentindices(V)
    sz = size.(args,2)
    sjr = intersect.(_argsindices(sz), Ref(jr))
    sjr2 = broadcast((a,b) -> a .- b .+ 1, sjr, _vcat_firstinds(sz))
    _view_hcat.(args, Ref(kr), sjr2)
end

function sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix)
    ret = similar(V)
    n = 0
    _,jr = parentindices(V)
    for a in arguments(V)
        m = size(a,1)
        copyto!(view(ret,n+1:n+m,:), a)
        n += m
    end
    ret
end

function sub_materialize(::ApplyLayout{typeof(hcat)}, V)
    ret = similar(V)
    n = 0
    kr,_ = parentindices(V)
    for a in arguments(V)
        m = size(a,2)
        copyto!(view(ret,:,n+1:n+m), a)
        n += m
    end
    ret
end

# temporarily allocate. In the future, we add a loop over arguments
materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))
materialize!(M::MatMulVecAdd{<:AbstractColumnMajor,<:ApplyLayout{typeof(vcat)}}) =
    materialize!(MulAdd(M.α,M.A,Array(M.B),M.β,M.C))

sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()
sublayout(::PaddedLayout{L}, ::Type{I}) where {L,I<:Tuple{AbstractUnitRange,AbstractUnitRange}} =
    PaddedLayout{typeof(sublayout(L(), I))}()

_lazy_getindex(dat, kr2) = lazy_getindex(dat, kr2)
_lazy_getindex(dat::Number, kr2) = dat

function sub_materialize(::PaddedLayout, v::AbstractVector{T}) where T
    A = parent(v)
    dat = paddeddata(A)
    (kr,) = parentindices(v)
    kr2 = kr ∩ axes(dat,1)
    Vcat(_lazy_getindex(dat, kr2), Zeros{T}(length(kr) - length(kr2)))
end

function sub_materialize(::PaddedLayout, v::AbstractMatrix{T}) where T
    A = parent(v)
    dat = paddeddata(A)
    kr,jr = parentindices(v)
    kr2 = kr ∩ axes(dat,1)
    Vcat(lazy_getindex(dat, kr2, jr), Zeros{T}(length(kr) - length(kr2), length(jr)))
end

## print

_replace_in_print_matrix(A::AbstractArray, k, j, s) = replace_in_print_matrix(A, k, j, s)
_replace_in_print_matrix(_, k, j, s) = s

function replace_in_print_matrix(f::Vcat{<:Any,1}, k::Integer, j::Integer, s::AbstractString)
    @assert j == 1
    κ = k
    for A in f.args
        n = length(A)
        κ ≤ n && return _replace_in_print_matrix(A, κ, 1, s)
        κ -= n
    end
    throw(BoundsError(f, k))
end

function replace_in_print_matrix(f::Vcat{<:Any,2}, k::Integer, j::Integer, s::AbstractString)
    κ = k
    for A in f.args
        n = size(A,1)
        κ ≤ n && return _replace_in_print_matrix(A, κ, j, s)
        κ -= n
    end
    throw(BoundsError(f, (k,j)))
end

# searchsorted

# function searchsorted(f::Vcat{<:Any,1}


text:
#####
##### 2LM "Optimizer"
#####

struct Optimizer2LM <: AbstractOptimizer{Rational{Int64}} end
Utils.getratio(::Optimizer2LM) = 0 // 1
name(::Optimizer2LM) = "2lm"

function factory(
        backend::nGraph.Backend{nGraph.CPU},
        func,
        opt::Optimizer2LM;
        # Useful for the GPU case
        adjust_io = false,
        defrag = true,
        just_profile = false,
        use_scratchpad = false,
        threshold = 20,
        profile_kw...
    )

    # Get the function, arguments, and keyword arguments from the provided function
    Profiler.enable_passes()
    A = func()

    #A callback that profiles the ngraph function
    function cb(f::nGraph.NFunction)
        # Do some minor editing the order of nodes in the graph to hopefully yield slightly
        # better memory characteristics
        priority_pass!(f)

        if use_scratchpad
            data = Profiler.FunctionData(f, nGraph.CPU)
            setup_scratchpad!(data; threshold = threshold)
        end
    end

    callbacks = CallbackChain()
    callback!(callbacks, cb)

    fex = nGraph.compile(
        backend,
        A.f,
        A.args...;
        callback = callbacks,
        emit_timing = true,
        A.kw...
    )

    return fex
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
# # How to Make a report

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ### Set up a report template
# In [report_templates](http://github.com/nrel-siip/PowerGraphics.jl/tree/master/report_templates) there are templates that can be copied and re-named to make a report.

# ```julia
# design_template = "generic_report_template.jmd"
# ```

# ### Generate report

# The report can be generated as a LaTeX PDF (default) or an html file.
# To set the document type as HTML, use the key word argument `doctype = "md2html"`

# ```julia
# out_path = joinpath(pwd(), "Generated_Reports/")
# if !isdir(out_path)
#     mkdir(out_path)
# end
# report(results, out_path, design_template; doctype = "md2html")
# ```

# creates an HTML

# ```julia
# report(results, out_path, design_template)
# ```

# creates a PDF


text:
# Rosenbrock distribution — based on https://arxiv.org/abs/1903.09556
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export RosenbrockDistribution, nl_pdf, nl_pdf_grad, rand

import Random: rand

struct RosenbrockDistribution
    μ::Float32
    a::Float32
    b::Array{Float32,2}
    n1::Int64
    n2::Int64
end

function RosenbrockDistribution(μ::Float32, a::Float32, b::Array{Float32,2})
    return RosenbrockDistribution(μ, a, b, size(b, 1), size(b, 2))
end

function nl_pdf_grad(X::AbstractArray{Float32, 4}, RB::RosenbrockDistribution)

    @assert size(X, 3) == 2 "Not implemented for larger dimension"

    nl_pdf_grad = cat(
        2f0*RB.a*(X[:, :, 1:1, :] .- RB.μ) - 4f0*RB.b[2, 1]*X[:, :, 1:1, :].*(X[:, :, 2:2, :] - X[:, :, 1:1, :].^2f0),
        2f0*RB.b[2, 1]*(X[:, :, 2:2, :] - X[:, :, 1:1, :].^2f0); dims=3)
    return nl_pdf_grad
end

function nl_pdf(X::AbstractArray{Float32, 4}, RB::RosenbrockDistribution)

    @assert size(X, 3) == (RB.n1 - 1) * RB.n2 + 1

    nlpdf = RB.a * (X[1, 1, 1, :] .- RB.μ).^2f0

    for j = 1:RB.n2
        for i = 2:RB.n1
            nlpdf += RB.b[i, j]*(
                X[1, 1, i + (j-1)*(RB.n1-1), :]
                - X[1, 1, i - 1 + (j-1)*(RB.n1-1), :].^2f0
            ).^2f0
        end
    end

    nc = sqrt(RB.a)/((1f0π)^(5f-1*((RB.n1 - 1) * RB.n2 + 1)))
    for j = 1:RB.n2
        for i = 2:RB.n1
            nc *= sqrt(RB.b[i, j])
        end
    end
    return -log(nc) .+ nlpdf
end

function rand(RB::RosenbrockDistribution, n_samples::Int64)

    dim = (RB.n1 - 1) * RB.n2 + 1
    X = zeros(Float32, 1, 1, dim, n_samples)

    X[1, 1, 1, :] = randn(Float32, n_samples)/sqrt(2f0*RB.a) .+ RB.μ
    for j = 1:RB.n2
        for i = 2:RB.n1
            X[1, 1, i+(j-1)*(RB.n1-1), :] = randn(Float32, n_samples)/sqrt(2f0*RB.b[i, j])
            X[1, 1, i+(j-1)*(RB.n1-1), :] += X[1, 1, i-1+(j-1)*(RB.n1 - 1), :].^2f0
        end
    end

    return X
end


text:
using Glob, PyCall, SeisNoise, Serialization
# read files
CORRDIR = "/media/four/LASSO/CORR/"
# do them all
files = glob("*/*",CORRDIR)
bins = collect(0.:0.05:55)
bins1 = collect(0.05:0.05:55.05)
counts = zeros(length(bins))
freqmin = 1.
freqmax = 2.
C = deserialize(files[1])
Cmat = zeros(eltype(Cs[1].corr),size(C.corr,1),length(bins))
for ii = 1:length(files)
    C = deserialize(files[ii])
    ind = findall((C.dist .>= bins) .& (C.dist .<= bins1))[1]
    clean_up!(C,freqmin,freqmax)
    abs_max!(C)
    Cmat[:,ind] .+= C.corr[:]
    counts[ind] += 1
    println("Reading file $ii")
end

maxdist = 25.
ind = findall((bins .> 0.15) .& (bins .<= maxdist))
Cmat = Cmat[:,ind]
counts = counts[ind]
bins = bins[ind]
Cmat ./= counts'
lags = -C.maxlag:1/C.fs:C.maxlag
maxlag = 15.
lagind = findall(abs.(lags) .<= maxlag)
Cmat = Cmat[lagind,:]
abs_max!(Cmat)
# plot envolopes at 5 km/s and 1.3 km/s
heatmap(
    lags[lagind],
    bins,
    Cmat',
    c=:balance,
    xlabel="Lag [s]",
    ylabel="Inter-station Distance [km]",
    legend = :none,
    dpi=500,
)
plot!([0.,5.],[0.1,25],line=(1.5,:dash),color=:black,label="")
plot!([0.,-5.],[0.1,25],line=(1.5,:dash),color=:black,label="")
plot!([0.,15.],[0.03,22.5],line=(1.5,:dash),color=:black,label="")
plot!([0.,-15.],[0.03,22.5],line=(1.5,:dash),color=:black,label="")
annotate!(2.5,17, Plots.text("5 km/s", 14, :dark, rotation = 78 ))
annotate!(8,10, Plots.text("1.5 km/s", 14, :dark, rotation = 48 ))
ylims!((minimum(bins),maximum(bins)))
xlims!((-maxlag,maxlag))
png(expanduser("~/SeisNoise_paper/LASSO-moveout"))


text:
using Test
using LogicCircuits
using DataFrames

@testset "MNIST Loader tests" begin

    ####################################
    train, valid, test = sampled_mnist()

    @test num_features(train) == 784
    @test num_examples(train) == 50000
    @test num_features(valid) == 784
    @test num_examples(valid) == 10000
    @test num_features(test) == 784
    @test num_examples(test) == 10000

    @test !isfpdata(train)
    @test isbinarydata(train)
    @test !isfpdata(valid)
    @test isbinarydata(valid)
    @test !isfpdata(test)
    @test isbinarydata(test)

    @test feature_values(train, 1) isa BitVector
    @test feature_values(train, 784) isa BitVector
    @test feature_values(valid, 1) isa BitVector
    @test feature_values(valid, 784) isa BitVector
    @test feature_values(test, 1) isa BitVector
    @test feature_values(test, 784) isa BitVector
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using SafeTestsets

@safetestset "Placeholder Tests" begin include("placeholder.jl") end


text:
@testset "Error Handling{$ST}" begin
  msg = PETSc.PetscErrorMessage(73)
  @test msg == "Object is in wrong state"
  @test_throws PETSc.PetscError PETSc.chk(76)
end


text:
### MUSE result

@doc doc"""

Stores the result of a MUSE run. Can be constructed by-hand as
`MuseResult()` and passed to any of the inplace `muse!`, `get_J!`, or
`get_H!`.

Fields:

* `θ` — The estimate of the $\theta$ parameters.
* `Σ, Σ⁻¹` — The approximate covariance of $\theta$ and its inverse.
* `H, J` — The $H$ and $J$ matrices which form the covariance (see
  [Millea & Seljak, 2021](https://arxiv.org/abs/2112.09354))
* `gs` — The MAP gradient sims used to compute `J`.
* `Hs` — The jacobian sims used to compute `H`.
* `dist` — A `Normal` or `MvNormal` built from `θ` and `Σ`, for
  convenience.
* `history` — Internal diagnostic info from the run.
* `rng` — RNG used to generate sims for this run (so the same sims can
  be reused if resuming later).
* `time` — Total `Millisecond` wall-time spent computing the result.

"""
Base.@kwdef mutable struct MuseResult
    θ = nothing
    H = nothing
    J = nothing
    Σ⁻¹ = nothing
    Σ = nothing
    dist = nothing
    history = []
    gs = []
    Hs = []
    rng = nothing
    time = Millisecond(0)
end

function Base.show(io::IO, result::MuseResult)
    _print(μ) = @sprintf("%.4g", μ)
    _print(μ, σ) = @sprintf("%.4g±%.3g", μ, σ)
    print(io, "MuseResult(")
    if result.θ != nothing && result.Σ != nothing
        σ² = result.θ isa AbstractVector ? diag(result.Σ) : result.Σ
        str = sprint(show, _print.(result.θ, sqrt.(σ²)))
    elseif result.θ != nothing
        str = sprint(show, _print.(result.θ))
    else
        str = ""
    end
    print(io, replace(str, "\"" => ""))
    print(io,")")
end

### MUSE solver

@doc doc"""

    muse(prob::AbstractMuseProblem, θ₀; kwargs...)
    muse!(result::MuseResult, prob::AbstractMuseProblem, [θ₀=nothing]; kwargs...)

Run the MUSE estimate. The `muse!` form resumes an existing result. If
the `muse` form is used instead, `θ₀` must give a starting guess for
$\theta$.

See [`MuseResult`](@ref) for description of return value.

Keyword arguments:

* `rng` — Random number generator to use. Taken from `result.rng` or
  `Random.default_rng()` if not passed.
* `z₀` — Starting guess for the latent space MAP.
* `maxsteps = 50` — Maximum number of iterations.
* `θ_rtol = 1e-2` — Error tolerance on $\theta$ relative to its
  standard deviation.
* `∇z_logLike_atol = 1e-2` — Absolute tolerance on the $z$-gradient at
  the MAP solution.
* `nsims = 100` — Number of simulations.
* `α = 0.7` — Step size for root-finder.
* `progress = false` — Show progress bar.
* `pmap` — Parallel map function.
* `regularize = identity` — Apply some regularization after each step.
* `H⁻¹_like = nothing` — Initial guess for the inverse Jacobian of
  $s^{\rm MUSE}(\theta)$
* `H⁻¹_update` — How to update `H⁻¹_like`. Should be `:sims`,
  `:broyden`, or `:diagonal_broyden`.
* `broyden_memory = Inf` — How many past steps to keep for Broyden
  updates.
* `checkpoint_filename = nothing` — Save result to a file after each
  iteration.
* `get_covariance = false` — Also call `get_H` and `get_J` to get the
  full covariance.

"""
muse(args...; kwargs...) = muse!(MuseResult(), args...; kwargs...)

@doc doc"""
See [`muse`](@ref).
"""
function muse!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = nothing;
    rng = nothing,
    z₀ = nothing,
    maxsteps = 50,
    θ_rtol = 1e-1,
    ∇z_logLike_atol = 1e-2,
    nsims = 100,
    α = 0.7,
    progress = false,
    pmap = _map,
    batch_size = 1,
    regularize = identity,
    H⁻¹_like′ = nothing,
    H⁻¹_update = :sims,
    broyden_memory = Inf,
    checkpoint_filename = nothing,
    get_covariance = false
)

    result.rng = rng = @something(rng, result.rng, copy(Random.default_rng()))
    θunreg  = θ  = θ₀ = standardizeθ(prob, @something(result.θ, θ₀))
    θunreg′ = θ′ = transform_θ(prob, θ)
    history = result.history

    _rng = copy(rng)
    xs_ẑs_sims = map(1:nsims) do i
        (x, z) = sample_x_z(prob, _rng, θ)
        (x, @something(z₀, z))
    end
    xs = [[prob.x];                                      first.(xs_ẑs_sims)]
    ẑs = [[@something(z₀, sample_x_z(prob, _rng, θ).z)]; last.(xs_ẑs_sims)]

    # set up progress bar
    pbar = progress ? RemoteProgress((maxsteps-length(result.history))*(nsims+1)÷batch_size, 0.1, "MUSE: ") : nothing

    try

        for i = (length(result.history)+1):maxsteps

            t₀ = now()

            if i > 1
                _rng = copy(rng)
                xs = [[prob.x]; [sample_x_z(prob, _rng, θ).x for i=1:nsims]]
            end

            if i > 2
                Δθ′ = history[end].θ′ - history[end-1].θ′
                sqrt(-(Δθ′' * history[end].H⁻¹_post′ * Δθ′)) < θ_rtol && break
            end

            # MUSE gradient
            gẑs = pmap(xs, ẑs, fill(θ,length(xs)); batch_size) do x, ẑ_prev, θ
                local ẑ, history = ẑ_at_θ(prob, x, ẑ_prev, θ; ∇z_logLike_atol)
                g  = ∇θ_logLike(prob, x, ẑ, θ,  UnTransformedθ())
                g′ = ∇θ_logLike(prob, x, ẑ, θ′, Transformedθ())
                progress && ProgressMeter.next!(pbar)
                (;g, g′, ẑ, history)
            end
            ẑs = getindex.(gẑs, :ẑ)
            ẑ_history_dat, ẑ_history_sims = peel(getindex.(gẑs, :history))
            g_like_dat′, g_like_sims′ = peel(getindex.(gẑs, :g′))
            _,           g_like_sims  = peel(getindex.(gẑs, :g))
            g_like′ = g_like_dat′ .- mean(g_like_sims′)
            g_prior′ = AD.gradient(AD.ForwardDiffBackend(), θ′ -> logPriorθ(prob, θ′, Transformedθ()), θ′)[1]
            g_post′ = g_like′ .+ g_prior′

            # Jacobian
            h⁻¹_like_sims′ = -1 ./ var(collect(g_like_sims′))
            H⁻¹_like_sims′ = h⁻¹_like_sims′ isa Number ? h⁻¹_like_sims′ : Diagonal(h⁻¹_like_sims′)
            if (H⁻¹_like′ == nothing) || (H⁻¹_update == :sims)
                H⁻¹_like′ = H⁻¹_like_sims′
            elseif i > 2 && (H⁻¹_update in [:broyden, :diagonal_broyden])
                # on subsequent steps, do a Broyden's update using at
                # most the previous `broyden_memory` steps
                j₀ = Int(max(2, i - broyden_memory))
                H⁻¹_like′ = history[j₀-1].H⁻¹_like_sims′
                for j = j₀:i-1
                    Δθ′      = history[j].θ′      - history[j-1].θ′
                    Δg_like′ = history[j].g_like′ - history[j-1].g_like′
                    H⁻¹_like′ = H⁻¹_like′ + ((Δθ′ - H⁻¹_like′ * Δg_like′) / (Δθ′' * H⁻¹_like′ * Δg_like′)) * Δθ′' * H⁻¹_like′
                    if H⁻¹_update == :diagonal_broyden
                        H⁻¹_like′ = Diagonal(H⁻¹_like′)
                    end
                end
            end

            H_prior′ = AD.hessian(AD.ForwardDiffBackend(), θ′ -> logPriorθ(prob, θ′, Transformedθ()), θ′)[1]
            H⁻¹_post′ = inv(inv(H⁻¹_like′) + H_prior′)

            t = now() - t₀
            push!(
                history,
                (;
                    θ, θunreg, θ′, θunreg′,
                    g_like_sims,
                    g_like_dat′, g_like_sims′, g_like′, g_prior′, g_post′,
                    H⁻¹_post′, H_prior′, H⁻¹_like′, H⁻¹_like_sims′,
                    ẑ_history_dat, ẑ_history_sims, t
                )
            )

            # Newton-Rhapson step
            θunreg′ = θ′ .- α .* (H⁻¹_post′ * g_post′)
            θunreg  = inv_transform_θ(prob, θunreg′)
            θ′ = regularize(θunreg′)
            θ  = inv_transform_θ(prob, θ′)

            (checkpoint_filename != nothing) && save(checkpoint_filename, "result", result)

        end

    finally

        progress && ProgressMeter.finish!(pbar)

    end

    result.time += sum(getindex.(history,:t))
    result.θ = θunreg
    result.gs = collect(history[end].g_like_sims)
    if get_covariance
        get_J!(result, prob; rng, nsims)
        get_H!(result, prob; rng, nsims=max(1,nsims÷10), ∇z_logLike_atol)
    end
    result

end

@doc doc"""

    get_H!(result::MuseResult, prob::AbstractMuseProblem, [θ₀=nothing]; kwargs...)

Compute the $H$ matrix, which is part of the MUSE covariance
computation (see [Millea & Seljak,
2021](https://arxiv.org/abs/2112.09354)).

Positional arguments:

* `result` — `MuseResult` into which to store result
* `prob` — `AbstractMuseProblem` being solved
* `θ₀` — the `θ` at which to evaluate $H$ (default: `result.θ` if it
  exists, otherwise `θ₀` must be given)

Keyword arguments:

* `z₀` — Starting guess for the latent space MAPs. Defaults to random
  sample from prior.
* `∇z_logLike_atol = 1e-2` — Absolute tolerance on the $z$-gradient at
  the MAP solution.
* `rng` — Random number generator to use. Taken from `result.rng` or
  `Random.default_rng()` if not passed.
* `nsims` — How many simulations to average over (default: `10`)
* `pmap` — Parallel map function.
* `progress` — Show progress bar (default: `false`),
* `skip_errors` — Ignore any MAP that errors (default: `false`)
* `fdm` — A `FiniteDifferenceMethod` used to compute the finite
  difference Jacobian of AD gradients involved in computing $H$
  (defaults to: `FiniteDifferences.central_fdm(3,1)`)
* `step` — A guess for the finite difference step-size (defaults to
  0.1σ for each parameter using J to estimate σ; for this reason its
  recommended to run `get_J!` before `get_H!`). Is only a guess
  because different choices of `fdm` may adapt this.

"""
function get_H!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = result.θ;
    fdm :: FiniteDifferenceMethod = central_fdm(3,1),
    ∇z_logLike_atol = 1e-2,
    rng = Random.default_rng(),
    nsims = 10,
    step = nothing,
    pmap = _map,
    batch_size = 1,
    pmap_over = :auto,
    progress = false,
    skip_errors = false,
    z₀ = nothing
)

    θ₀ = standardizeθ(prob, @something(θ₀, result.θ))
    𝟘 = zero(θ₀) * zero(θ₀)' # if θ::ComponentArray, helps keep component labels
    nsims_remaining = nsims - length(result.Hs)
    (nsims_remaining <= 0) && return
    pbar = progress ? RemoteProgress(nsims_remaining*(1+length(θ₀))÷batch_size, 0.1, "get_H: ") : nothing
    t₀ = now()

    # default to finite difference step size of 0.1σ with σ roughly
    # estimated from g sims, if we have them
    if step == nothing && !isempty(result.gs)
        step = 0.1 ./ std(result.gs)
    end

    # generate simulation locally, advancing rng, and saving rng state to be reused remotely
    xs_ẑ₀s_rngs = map(1:nsims_remaining) do i
        _rng = copy(rng)
        (x, z) = sample_x_z(prob, rng, θ₀)
        (x, @something(z₀, z), _rng)
    end

    # initial fit at fiducial, used as starting points for finite difference below
    ẑ₀s_rngs = pmap(xs_ẑ₀s_rngs; batch_size) do (x, ẑ₀, rng)
        ẑ, = ẑ_at_θ(prob, x, ẑ₀, θ₀; ∇z_logLike_atol)
        progress && ProgressMeter.next!(pbar)
        (ẑ, rng)
    end

    # finite difference Jacobian
    pmap_sims, pmap_jac = (pmap_over == :jac || (pmap_over == :auto && length(θ₀) > nsims_remaining)) ? (_map, pmap) : (pmap, _map)
    append!(result.Hs, skipmissing(pmap_sims(ẑ₀s_rngs; batch_size) do (ẑ₀, rng)
        try
            return first(pjacobian(fdm, θ₀, step; pmap=pmap_jac, batch_size, pbar) do θ
                x, = sample_x_z(prob, copy(rng), θ)
                ẑ, = ẑ_at_θ(prob, x, ẑ₀, θ₀; ∇z_logLike_atol)
                ∇θ_logLike(prob, x, ẑ, θ₀, UnTransformedθ())
            end)
        catch err
            if skip_errors && !(err isa InterruptException)
                @warn err
                return missing
            else
                rethrow(err)
            end
        end
    end))

    result.H = (θ₀ isa Number) ? mean(first.(result.Hs)) : (mean(result.Hs) .+ 𝟘)
    result.time += now() - t₀
    finalize_result!(result, prob)

end

@doc doc"""

    get_J!(result::MuseResult, prob::AbstractMuseProblem, [θ₀=nothing]; kwargs...)

Compute the $J$ matrix, which is part of the MUSE covariance
computation (see [Millea & Seljak,
2021](https://arxiv.org/abs/2112.09354)).

Positional arguments:

* `result` — `MuseResult` into which to store result
* `prob` — `AbstractMuseProblem` being solved
* `θ₀` — the `θ` at which to evaluate $J$ (default: `result.θ` if it
  exists, otherwise `θ₀` must be given)

Keyword arguments:

* `z₀` — Starting guess for the latent space MAPs. Defaults to random
  sample from prior.
* `∇z_logLike_atol = 1e-2` — Absolute tolerance on the $z$-gradient at
  the MAP solution.
* `rng` — Random number generator to use. Taken from `result.rng` or
  `Random.default_rng()` if not passed.
* `nsims` — How many simulations to average over (default: `100`)
* `pmap` — Parallel map function.
* `progress` — Show progress bar (default: `false`),
* `skip_errors` — Ignore any MAP that errors (default: `false`)
* `covariance_method` — A `CovarianceEstimator` used to compute $J$
  (default: `SimpleCovariance(corrected=true)`)

"""
function get_J!(
    result :: MuseResult,
    prob :: AbstractMuseProblem,
    θ₀ = nothing;
    z₀ = nothing,
    ∇z_logLike_atol = 1e-2,
    rng = Random.default_rng(),
    nsims = 100,
    pmap = _map,
    batch_size = 1,
    progress = false,
    skip_errors = false,
    covariance_method = SimpleCovariance(corrected=true),
)

    θ₀ = standardizeθ(prob, @something(θ₀, result.θ))
    nsims_remaining = nsims - length(result.gs)

    if nsims_remaining > 0

        pbar = progress ? RemoteProgress(nsims_remaining÷batch_size, 0.1, "get_J: ") : nothing

        xs_z₀s = map(1:nsims_remaining) do i
            (x, z) = sample_x_z(prob, rng, θ₀)
            (x, @something(z₀, z))
        end

        append!(result.gs, skipmissing(pmap(xs_z₀s, fill(θ₀,nsims_remaining); batch_size) do (x, z₀), θ₀
            try
                ẑ, = ẑ_at_θ(prob, x, z₀, θ₀; ∇z_logLike_atol)
                g = ∇θ_logLike(prob, x, ẑ, θ₀, UnTransformedθ())
                progress && ProgressMeter.next!(pbar)
                return g
            catch err
                if skip_errors && !(err isa InterruptException)
                    @warn err
                    return missing
                else
                    rethrow(err)
                end
            end
        end))

    end

    result.J = (θ₀ isa Number) ? var(result.gs) : cov(covariance_method, identity.(result.gs))
    finalize_result!(result, prob)

end

function finalize_result!(result::MuseResult, prob::AbstractMuseProblem)
    @unpack H, J, θ = result
    if H != nothing && J != nothing && θ != nothing
        𝟘 = zero(J) # if θ::ComponentArray, helps keep component labels
        H_prior = -AD.hessian(AD.ForwardDiffBackend(), θ -> logPriorθ(prob, θ, UnTransformedθ()), result.θ)[1]
        result.Σ⁻¹ = H' * inv(J) * H + H_prior + 𝟘
        result.Σ = inv(result.Σ⁻¹) + 𝟘
        if length(result.θ) == 1
            result.dist = Normal(result.θ[1], sqrt(result.Σ[1]))
        else
            result.dist = MvNormal(result.θ, Symmetric(Array(result.Σ)))
        end
    end
    result
end


text:
using ValueOrientedRiskManagementInsurance
using DataFrames
using Base.Test

println("start ECModel2.jl...")

lob_strings = ["Fire", "Liability", "Theft"]
df_claims = Vector(length(lob_strings))
claims = Vector(length(lob_strings))
mack = Vector(length(lob_strings))
β = Vector(length(lob_strings))
for i = 1:length(df_claims)
  df_claims[i] =
    readtable("test/ECModel_Input_" * lob_strings[i] * ".csv",
              header = false)
  mack[i] = Mack(df_claims[i])
  β[i] = mack[i].futureclaims / sum(mack[i].futureclaims)
end

function logpar2statpar(m_s::Vector)
  Real[exp(m_s[1] + m_s[2]^2 / 2),
       √(exp(m_s[2]^2-1)) * exp(m_s[1] + m_s[2]^2 / 2)]
end

function statpar2logpar(μ_σ::Vector)
  vc = μ_σ[2]/μ_σ[1]
  Real[log(μ_σ[1]) - 0.5 * (1+log(vc^2)), √(1+log(vc^2))]
end

logpar2statpar([0,1])
x = logpar2statpar([0,1])

statpar2logpar(x)

statpar2logpar(x[1],x[2])
#
# mean, sde -> logmean, logsigma
# R[i] = LogNormal(log_mean[i], log_sigma[i])
# v⋅mack[i].β * R[i]
#
# v: stochastic discont
# β
# R[i] = stochastic total undiscounted future paiments,
#

println("...end ECModel2.jl")


text:
"""Generic combiner"""
combiner_parsejson(input) = [chomp(line) for line in input]


text:
using ..Spaces:
    AbstractSpace, SpectralElementSpace1D, SpectralElementSpace2D, Quadratures
using ..Topologies: Topology2D, IntervalTopology
using ..Fields: Field
using ..DataLayouts
using SparseArrays, LinearAlgebra

struct LinearRemap{T <: AbstractSpace, S <: AbstractSpace, M <: AbstractMatrix}
    target::T
    source::S
    map::M # linear mapping operator
end

"""
    LinearRemap(target::AbstractSpace, source::AbstractSpace)

A remapping operator from the `source` space to the `target` space.
"""
function LinearRemap(target::AbstractSpace, source::AbstractSpace)
    R = linear_remap_op(target, source)
    LinearRemap(target, source, R)
end

"""
    remap(R::LinearRemap, source_field::Field)

Applies the remapping operator `R` to `source_field`. Outputs a new field on the target space mapped to by `R`.
"""
function remap(R::LinearRemap, source_field::Field)
    target_space = R.target
    target_field = similar(source_field, target_space, eltype(source_field))
    remap!(target_field, R, source_field)
end

"""
    remap!(target_field::Field, R::LinearRemap, source_field::Field)

Applies the remapping operator `R` to `source_field` and stores the solution in `target_field`.
"""
function remap!(target_field::Field, R::LinearRemap, source_field::Field)
    mul!(vec(parent(target_field)), R.map, vec(parent(source_field)))
    return target_field
end

"""
    linear_remap_op(target::AbstractSpace, source::AbstractSpace)

Computes linear remapping operator `R` for remapping from `source` to `target` spaces.

Entry `R_{ij}` gives the contribution weight to the target node `i` from
source node `j`; nodes are indexed by their global position, determined
by both element index and nodal order within the element.
"""
function linear_remap_op(target::AbstractSpace, source::AbstractSpace)
    J = 1.0 ./ local_weights(target) # workaround for julia #26561
    W = overlap_weights(target, source)
    return W .* J
end

"""
    overlap_weights(target, source)

Computes local weights of the overlap mesh for `source` to `target` spaces.
"""
function overlap_weights(
    target::T,
    source::S,
) where {
    T <: SpectralElementSpace1D{<:IntervalTopology, Quadratures.GL{1}},
    S <: SpectralElementSpace1D{<:IntervalTopology, Quadratures.GL{1}},
}
    # Calculate element overlap pattern
    # X_ov[i,j] = overlap length between target elem i and source elem j
    X_ov = fv_x_overlap(target, source)

    return X_ov
end

function overlap_weights(
    target::T,
    source::S,
) where {
    T <: SpectralElementSpace1D{<:IntervalTopology},
    S <: SpectralElementSpace1D{<:IntervalTopology, Quadratures.GL{1}},
}
    FT = Spaces.undertype(source)
    target_topo = Spaces.topology(target)
    source_topo = Spaces.topology(source)
    nelems_t = Topologies.nlocalelems(target)
    nelems_s = Topologies.nlocalelems(source)
    QS_t = Spaces.quadrature_style(target)
    QS_s = Spaces.quadrature_style(source)
    Nq_t = Quadratures.degrees_of_freedom(QS_t)
    Nq_s = Quadratures.degrees_of_freedom(QS_s)
    J_ov = spzeros(nelems_t * Nq_t, nelems_s * Nq_s)

    # Calculate element overlap pattern
    # X_ov[i,j] = overlap length between target elem i and source elem j
    for i in 1:nelems_t
        vertices_i = Topologies.vertex_coordinates(target_topo, i)
        min_i, max_i = Geometry.component(vertices_i[1], 1),
        Geometry.component(vertices_i[2], 1)
        for j in 1:nelems_s # elems coincide w nodes in FV source
            vertices_j = Topologies.vertex_coordinates(source_topo, j)
            # get interval for quadrature
            min_j, max_j = Geometry.component(vertices_j[1], 1),
            Geometry.component(vertices_j[2], 1)
            min_ov, max_ov = max(min_i, min_j), min(max_i, max_j)
            if max_ov <= min_ov
                continue
            end
            ξ, w = Quadratures.quadrature_points(FT, QS_t)
            x_ov =
                FT(0.5) * (min_ov + max_ov) .+ FT(0.5) * (max_ov - min_ov) * ξ
            x_t = FT(0.5) * (min_i + max_i) .+ FT(0.5) * (max_i - min_i) * ξ

            # column k of I_mat gives the k-th target basis function defined on the overlap element
            I_mat = Quadratures.interpolation_matrix(x_ov, x_t)
            for k in 1:Nq_t
                idx = Nq_t * (i - 1) + k # global nodal index
                # (integral of basis on overlap) / (reference elem length * overlap elem length)
                J_ov[idx, j] = w' * I_mat[:, k] ./ 2 * (max_ov - min_ov)
            end
        end
    end
    return J_ov
end

function overlap_weights(
    target::T,
    source::S,
) where {
    T <: SpectralElementSpace1D{<:IntervalTopology, Quadratures.GL{1}},
    S <: SpectralElementSpace1D{<:IntervalTopology},
}
    J_ov = overlap_weights(source, target)
    return J_ov'
end

function overlap_weights(
    target::T,
    source::S,
) where {
    T <: SpectralElementSpace2D{<:Topology2D, Quadratures.GL{1}},
    S <: SpectralElementSpace2D{<:Topology2D, Quadratures.GL{1}},
}
    # Calculate element overlap pattern in x-dimension
    # X_ov[i,j] = overlap length along x-dimension between target elem i and source elem j
    X_ov = fv_x_overlap(target, source)

    # Calculate element overlap pattern in y-dimension
    Y_ov = fv_y_overlap(target, source)

    return kron(Y_ov, X_ov)
end

function fv_x_overlap(
    target::T,
    source::S,
) where {
    T <: Union{SpectralElementSpace1D, SpectralElementSpace2D},
    S <: Union{SpectralElementSpace1D, SpectralElementSpace2D},
}
    target_topo = Spaces.topology(target)
    source_topo = Spaces.topology(source)
    ntarget, nsource = nxelems(target_topo), nxelems(source_topo)
    v1, v2 = 1, 2

    W_ov = spzeros(ntarget, nsource)
    for i in 1:ntarget
        vertices_i = Topologies.vertex_coordinates(target_topo, i)
        min_i, max_i = xcomponent(vertices_i[v1]), xcomponent(vertices_i[v2])
        for j in 1:nsource
            vertices_j = Topologies.vertex_coordinates(source_topo, j)
            min_j, max_j =
                xcomponent(vertices_j[v1]), xcomponent(vertices_j[v2])
            min_ov, max_ov = max(min_i, min_j), min(max_i, max_j)
            overlap_length = max_ov > min_ov ? max_ov - min_ov : continue
            W_ov[i, j] = overlap_length
        end
    end
    return W_ov
end

function fv_y_overlap(
    target::T,
    source::S,
) where {
    T <: Union{SpectralElementSpace1D, SpectralElementSpace2D},
    S <: Union{SpectralElementSpace1D, SpectralElementSpace2D},
}
    target_topo = Spaces.topology(target)
    source_topo = Spaces.topology(source)
    ntarget, nsource = nyelems(target_topo), nyelems(source_topo)
    nx_target, nx_source = nxelems(target_topo), nxelems(source_topo)
    elem_idx = (i, n) -> 1 + (i - 1) * n
    v1, v2 = 1, 4

    W_ov = spzeros(ntarget, nsource)
    for i in 1:ntarget
        vertices_i =
            Topologies.vertex_coordinates(target_topo, elem_idx(i, nx_target))
        min_i, max_i = ycomponent(vertices_i[v1]), ycomponent(vertices_i[v2])
        for j in 1:nsource
            vertices_j = Topologies.vertex_coordinates(
                source_topo,
                elem_idx(j, nx_source),
            )
            min_j, max_j =
                ycomponent(vertices_j[v1]), ycomponent(vertices_j[v2])
            min_ov, max_ov = max(min_i, min_j), min(max_i, max_j)
            overlap_length = max_ov > min_ov ? max_ov - min_ov : continue
            W_ov[i, j] = overlap_length
        end
    end
    return W_ov
end

nxelems(topology::Topologies.IntervalTopology) =
    Topologies.nlocalelems(topology)
nxelems(topology::Topologies.Topology2D) =
    size(Meshes.elements(topology.mesh), 1)
nyelems(topology::Topologies.Topology2D) =
    size(Meshes.elements(topology.mesh), 1)

xcomponent(x::Geometry.XPoint) = Geometry.component(x, 1)
xcomponent(xy::Geometry.XYPoint) = Geometry.component(xy, 1)
ycomponent(y::Geometry.YPoint) = Geometry.component(y, 1)
ycomponent(xy::Geometry.XYPoint) = Geometry.component(xy, 2)

"""
    local_weights(space::AbstractSpace)

Each degree of freedom is associated with a local weight J_i.
For finite volumes the local weight J_i would represent the geometric area
of the associated region. For nodal finite elements, the local weight
represents the value of the global Jacobian, or some global integral of the
associated basis function.
"""
function local_weights(space::AbstractSpace)
    wj = space.local_geometry.WJ
    return vec(parent(wj))
end

slab_value(data::DataLayouts.IJFH, i) = slab(data, i)[1, 1]
slab_value(data::DataLayouts.IFH, i) = slab(data, i)[1]


text:
module mcmc

using LinearAlgebra
using Flux
using Flux: Tracker, params

function mala(logdensity, h, M, iters, init_proposal_mean; δ_trunc=Inf)

    function log_llh_and_proposal(z)
        par_z = param(z)
        logp_u = logdensity(par_z)
        Tracker.back!(logp_u)
        grad = Tracker.grad(par_z)

        #  ----- Truncated MALA (T-MALA) ----
        #  ----- (aka norm clipping!) -------
        if δ_trunc < Inf
            grad_norm = sum(x -> x^2, grad)
            if grad_norm > δ_trunc
                grad *= δ_trunc / grad_norm
            end
        end
        return Tracker.data(logp_u), z + 0.5*h*M*grad, grad   # prev (unnorm) logp, gradient step, gradient
    end

    # set up
    d = length(init_proposal_mean)
    ztrace = Array{Float64, 2}(undef, d, iters)
    cholΣ = Matrix(cholesky(h*M).L)   # `Matrix` keeps Flux happy for the time being. For d < 20, it really makes no difference.

    function gauss_logpdf(z)
        u = cholΣ \ z
        exponent = -0.5*sum(x->x*x, u)
#         lognormconst = -0.5*logdet(2π * h * M)
        out = exponent #+ lognormconst
        return out
    end
    gauss_logpdf(z, μ) = gauss_logpdf(z - μ)

    # initialise
    accepted = 0
    z = init_proposal_mean
    cur_u_logp = -Inf
    proposal_mean = init_proposal_mean
    prev_proposal_mean = init_proposal_mean

    # Markov Chain loop
    for i = 1:iters
        # transition r.v.
        ϵ_prop = randn(d)
        prop_transition_logp = -0.5*sum(x->x*x, ϵ_prop)

        # deterministic proposal
        z_prop = proposal_mean + cholΣ * ϵ_prop

        # calculate logdensity(z_prop) and accompanying new proposal mean
        prop_u_logp, ν_proposal_mean, _ = log_llh_and_proposal(z_prop)

        # metropolis-hastings
        r = prop_u_logp - cur_u_logp +
                gauss_logpdf(z, ν_proposal_mean) -
                prop_transition_logp   |> exp

        # if <accept>
        if rand() < r
            z = z_prop                                                       # update z
            prev_proposal_mean = proposal_mean                               # store prev proposal mean
            cur_u_logp, proposal_mean = prop_u_logp, ν_proposal_mean
            accepted += 1
        end

        ztrace[:,i] = z      # store z regardless
    end

    return ztrace, accepted/iters
end

end


text:
using Zygote, Test
using Zygote: gradient

@testset "Zygote" begin

@info "Testing compiler features"

@testset "Features" begin
  include("features.jl")
end

@info "Testing data structures"

@testset "Data Structures" begin
  include("structures.jl")
end

@info "Running Gradient Checks"

@testset "Gradients" begin
  include("gradcheck.jl")
end

@testset "Complex" begin
  include("complex.jl")
end

@info "Testing Inference & Debug Info"

@testset "Compiler" begin
  include("compiler.jl")
end

end