text:
x, y, z = activate_github("arnavs/InstantiationTest")
@test z == Dict{AbstractString,Any}("Expectations" => "2fe49d83-0758-5602-8f54-1f90ad0d522b", "Distances" => "b4f34e82-e78d-54a5-968a-f98e89d6e8f7", "Suppressor" => "fd094767-a336-5f1f-9728-57cf17d0bbfb", "StatsFuns" => "4c63d2b9-4356-54db-8cca-17b64c39e42c")


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

x = 1.0
y = 0.54

si = billiard_sinai(y/4, x, y)
re = billiard_rectangle(x,  y)

ps = particlebeam(4x/5, y/5, π/9, 20, 0.01)

colors = COLORS[1:2]

bd = si

for (bd, name) in zip((si, re), ("sinai", "rect"))

ps = particlebeam(4x/5, y/5, π/9, 20, 0.01)

# function billiard_video(file::String, bd::Billiard, ps::Vector{<:AbstractParticle};
reslarge = (1200*x, 1200*y)
reszoom = (800, 800)
dt = 0.0001
speed = 10
tail = 1200
frames = 650
framerate = 60
particle_size = 2

# first video without any adjustments
file = videodir("stretching_$(name)_large.mp4")

billiard_video(
    file, bd, deepcopy(ps);
    frames, colors,
    tailwidth = 3.5, particle_size = 0.8, res = reslarge,
    plot_particles = true,
    dt, speed, tail = round(Int, tail*1.5)
)

vr = InteractiveDynamics._estimate_vr(bd)/8

ps = particlebeam(4x/5, y/5, π/9, 20, 0.01)

fig, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
    res = reszoom, particle_size, tail, vr, add_controls = false, colors,
    displayfigure = false,
)
N = length(ps)
ci = N÷2 # index of chosen particle to zoom at and set axis limits
ax = content(fig[1,1]) # axis we need to be zooming in.
zx = 0.08 # zoom limits

file = videodir("stretching_$(name)_zoom.mp4")
Makie.inline!(true) # to not display figure while recording
record(fig, file, 1:frames; framerate) do j
# for j in 1:frames
    for i in 1:N
        parobs = allparobs[i]
        InteractiveDynamics.animstep!(parobs, bd, dt, true)
        balls[][i] = parobs.p.pos
        vels[][i] = vr*particle_size*parobs.p.vel
    end
    balls[] = balls[]
    vels[] = vels[]
    # Set axis limits
    pos = allparobs[ci].p.pos
    xlims!(ax, pos[1] - zx, pos[1] + zx)
    ylims!(ax, pos[2] - zx, pos[2] + zx)
    # sleep(0.000001)
    # yield()
    for _ in 1:speed
        for i in 1:N
            parobs = allparobs[i]
            InteractiveDynamics.animstep!(parobs, bd, dt, false)
        end
    end
end
Makie.inline!(false)

end

ext: jl

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using ExprSearch, SymbolicRegression
using ExprSearch.GP

using RLESUtils, TreeUtils
using DerivTreeVis

function test_rand(target_depth=5)
    problem = Symbolic("gt_easy.jl")
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    for i = 1:5
        ind = rand(grammar, mda, target_depth)
        @show get_expr(ind.derivtree)
    end
end

function test_randnode()
    problem = Symbolic()
    p = GPESParams(10,5,10,10,0.1,0.4,0.2,0.2,0.0)
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    ind = rand(grammar, mda, 5)
    node = rand_node(ind.derivtree.root)
end

function test_crossover()
    problem = Symbolic()
    p = GPESParams(10,5,10,10,0.1,0.4,0.2,0.2,0.0)
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    ind1 = rand(grammar, mda, 5)
    node1 = rand_node(ind1.derivtree.root)
    derivtreevis(ind1.derivtree, "node1")
    ind2 = rand(grammar, mda, 5)
    node2 = rand_node(ind2.derivtree.root)
    derivtreevis(ind2.derivtree, "node2")
    (ind3, ind4) = crossover(ind1, ind2, grammar,5)
    derivtreevis(ind3.derivtree, "node3")
    derivtreevis(ind4.derivtree, "node4")
end

function test_maxdepth()
    problem = Symbolic()
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    for i = 1:5
        ind = rand(grammar, mda, 5)
        @show max_depth(ind)
        derivtreevis(ind.derivtree, "ind$i")
    end
end

function test_mutate()
    problem = Symbolic()
    p = GPESParams(10,5,10,10,0.1,0.4,0.2,0.2,0.0)
    grammar = get_grammar(problem)
    mda = min_depth_actions(grammar)
    ind1 = rand(grammar, mda, 10)
    derivtreevis(ind1.derivtree, "ind1")
    ind2 = mutate(ind1, grammar, mda, 10)
    derivtreevis(ind2.derivtree, "ind2")
end

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
############################################################
## joAbstractSAparallelLinearOperator(s) - overloaded Base functions

# eltype(jo)
eltype(A::joAbstractSAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = promote_type(DDT,RDT)

# deltype(jo)
deltype(A::joAbstractSAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = DDT

# reltype(jo)
reltype(A::joAbstractSAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = RDT

# show(jo)
function show(A::joAbstractSAparallelLinearOperator)
    println("Type: $(typeof(A).name)")
    println("Name: $(A.name)")
    println("Size: $(size(A))")
    println(" NVC: $(A.nvc)")
    println(" DDT: $(deltype(A))")
    println(" RDT: $(reltype(A))")
    return nothing
end

# display(jo)
display(A::joAbstractSAparallelLinearOperator) = println((typeof(A),((A.m,A.n),A.nvc),A.name))

# size(jo)
size(A::joAbstractSAparallelLinearOperator) = A.m,A.n

# size(jo,1/2)
function size(A::joAbstractSAparallelLinearOperator,ind::Integer)
    if ind==1
		return A.m
	elseif ind==2
		return A.n
	elseif ind==3
		return A.nvc
	else
		throw(joAbstractSAparallelLinearOperatorException("invalid index"))
	end
end

# length(jo)
length(A::joAbstractSAparallelLinearOperator) = A.m*A.n

# jo_full(jo)

# norm(jo)

# real(jo)

# imag(jo)

# conj(jo)
conj(A::joSALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSALinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.PAs_out, A.fclean, A.rclean
        )
conj(A::joSAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributedLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.PAs_out, A.fclean, A.rclean
        )
conj(A::joSAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributingLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_out, A.gclean
        )
conj(A::joSAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAgatheringLinearOperator{DDT,RDT,N}("conj("*A.name*")",A.m,A.n,A.nvc,
        get(A.fop_C), A.fop_A, A.fop_T, A.fop,
        A.iop_C, A.iop_A, A.iop_T, A.iop,
        A.PAs_in, A.gclean
        )

# transpose(jo)
transpose(A::joSALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSALinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
transpose(A::joSAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributedLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
transpose(A::joSAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAgatheringLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_out, A.gclean
        )
transpose(A::joSAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributingLinearOperator{RDT,DDT,N}("transpose("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_T), A.fop, A.fop_C, A.fop_A,
        A.iop_T, A.iop, A.iop_C, A.iop_A,
        A.PAs_in, A.gclean
        )

# adjoint(jo)
adjoint(A::joSALinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSALinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
adjoint(A::joSAdistributedLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributedLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.PAs_in, A.rclean, A.fclean
        )
adjoint(A::joSAdistributingLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAgatheringLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_out, A.gclean
        )
adjoint(A::joSAgatheringLinearOperator{DDT,RDT,N}) where {DDT,RDT,N} =
    joSAdistributingLinearOperator{RDT,DDT,N}("adjoint("*A.name*")",A.n,A.m,A.nvc,
        get(A.fop_A), A.fop_C, A.fop, A.fop_T,
        A.iop_A, A.iop_C, A.iop, A.iop_T,
        A.PAs_in, A.gclean
        )

# isreal(jo)
isreal(A :: joAbstractSAparallelLinearOperator{DDT,RDT}) where {DDT,RDT} = (DDT<:Real && RDT<:Real)

# issymmetric(jo)

# ishermitian(jo)

# getindex(jo,...)

############################################################
## overloaded Base *(...jo...)

# *(jo,jo)
function *(A::joSALinOpsUnion{CDT,ARDT,2},B::joSALinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joSALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joSALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.fclean,A.rclean)
end
function *(A::joSALinOpsUnion{CDT,ARDT,2},B::joSAdistributeLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joSALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joSAdistributingLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               A.PAs_out,B.gclean)
end
function *(A::joSAgatherLinOpsUnion{CDT,ARDT,2},B::joSALinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joSAgatheringLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joSAgatheringLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joSAgatheringLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joSAgatheringLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,B.fclean)
end
function *(A::joSAdistributeLinOpsUnion{CDT,ARDT,2},B::joSAgatherLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joSALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    return joSALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.gclean,A.gclean)
end
function *(A::joSAgatherLinOpsUnion{CDT,ARDT,2},B::joSAdistributeLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.nvc==B.nvc || throw(joLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    return joLinearOperator{BDDT,ARDT}("($(A.name)*$(B.name))",size(A,1),size(B,2),
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF)
end
# extras with warnings
function *(A::joSAdistribute{CDT,ARDT,2},B::joSAgather{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joSALinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.PAs_out.dims[2]==B.PAs_in.dims[2] || throw(joSALinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    @warn "*($(typeof(A)),$(typeof(B))) is a senseless operation. Get rid of it!"
    return joSALinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.PAs_out.dims[2],
               v1->v1, v2->v2, v3->v3, v4->v4,
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,A.PAs_out,B.gclean,A.gclean)
end
function *(A::joSAgather{CDT,ARDT,2},B::joSAdistribute{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    A.PAs_in.dims[2]==B.PAs_out.dims[2] || throw(joLinearOperatorException("*($(A.name),$(B.name)): nvc mismatch"))
    isapprox(A.PAs_in,B.PAs_out) || throw(joLinearOperatorException("*($(A.name),$(B.name)): distributor mismatch"))
    @warn "*($(typeof(A)),$(typeof(B))) is a senseless operation. Get rid of it!"
    return joLinearOperator{BDDT,ARDT}("($(A.name)*$(B.name))",size(A,1),size(B,2),
               v1->v1, v2->v2, v3->v3, v4->v4,
               @joNF, @joNF, @joNF, @joNF)
end
# joSA with joAbstractLinearOperator
function *(A::joSAdistributeLinOpsUnion{CDT,ARDT,2},B::joAbstractLinearOperator{BDDT,CDT}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joAbstractLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    return joSAdistributingLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),A.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               A.PAs_out,A.gclean)
end
function *(A::joAbstractLinearOperator{CDT,ARDT},B::joSAgatherLinOpsUnion{BDDT,CDT,2}) where {ARDT,BDDT,CDT}
    size(A,2) == size(B,1) || throw(joAbstractLinearOperatorException("*($(A.name),$(B.name)): shape mismatch"))
    return joSAgatheringLinearOperator{BDDT,ARDT,2}("($(A.name)*$(B.name))",size(A,1),size(B,2),B.nvc,
               v1->A*(B*v1),
               v2->transpose(B)*(transpose(A)*v2),
               v3->adjoint(B)*(adjoint(A)*v3),
               v4->conj(A)*(conj(B)*v4),
               @joNF, @joNF, @joNF, @joNF,
               B.PAs_in,B.gclean)
end

# *(jo,mvec)
function *(A::joSALinearOperator{ADDT,ARDT,2},mv::SharedArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joSALinearOperatorException("shape mismatch"))
    A.nvc == size(mv,2) || throw(joSALinearOperatorException("shape mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joSALinearOperatorException("*($(A.name),mv::SharedArray): input distributor mismatch"))
    MV=A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joSAdistributedLinearOperator{ADDT,ARDT,2},mv::SharedArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joSAdistributedLinearOperatorException("shape mismatch"))
    A.nvc == size(mv,2) || throw(joSAdistributedLinearOperatorException("shape mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joSAdistributedLinearOperatorException("*($(A.name),mv::SharedArray): input distributor mismatch"))
    MV=salloc(A.PAs_out)
    joSAutils.jo_x_mv!(A.fop,A.PAs_in,A.PAs_out,mv,MV)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joSAdistributingLinearOperator{ADDT,ARDT,2},mv::LocalMatrix{mvDT}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joSAdistributingLinearOperatorException("shape mismatch in A$(size(A))*v$(size(mv))"))
    A.nvc == size(mv,2) || throw(joSAdistributingLinearOperatorException("nvc size mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    MV = A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end
function *(A::joSAgatheringLinearOperator{ADDT,ARDT,2},mv::SharedArray{mvDT,2}) where {ADDT,ARDT,mvDT<:Number}
    A.n == size(mv,1) || throw(joSAgatheringLinearOperatorException("shape mismatch in A$(size(A))*v$(size(mv))"))
    A.nvc == size(mv,2) || throw(joSAgatheringLinearOperatorException("nvc size mismatch"))
    jo_check_type_match(ADDT,mvDT,join(["DDT for *(jo,mvec):",A.name,typeof(A),mvDT]," / "))
    isequiv(A.PAs_in,mv) || throw(joSAgatheringLinearOperatorException("*($(A.name),mv::SharedArray): input distributor mismatch"))
    MV = A.fop(mv)
    jo_check_type_match(ARDT,eltype(MV),join(["RDT from *(jo,mvec):",A.name,typeof(A),eltype(MV)]," / "))
    return MV
end

# *(mvec,jo)

# *(jo,vec)

# *(vec,jo)

# *(num,jo)

# *(jo,num)

############################################################
## overloaded Base \(...jo...)

# \(jo,jo)

# \(jo,mvec)

# \(mvec,jo)

# \(jo,vec)

# \(vec,jo)

# \(num,jo)

# \(jo,num)
#\{ADDT,ARDT}(A::joAbstractSAparallelLinearOperator{ADDT,ARDT},a::Number) = inv(a)*A
#\{ADDT,ARDT}(A::joAbstractSAparallelLinearOperator{ADDT,ARDT},a::joNumber{ADDT,ARDT}) = inv(a)*A

############################################################
## overloaded Base +(...jo...)

# +(jo)
+(A::joAbstractSAparallelLinearOperator) = A

# +(jo,jo)

# +(jo,mvec)

# +(mvec,jo)

# +(jo,vec)

# +(vec,jo)

# +(jo,num)

# +(num,jo)
+(b::Number,A::joAbstractSAparallelLinearOperator{ADDT,ARDT}) where {ADDT,ARDT} = A+b
+(b::joNumber{ADDT,ARDT},A::joAbstractSAparallelLinearOperator{ADDT,ARDT}) where {ADDT,ARDT} = A+b

############################################################
## overloaded Base -(...jo...)

# -(jo)

# -(jo,jo)
-(A::joAbstractSAparallelLinearOperator,B::joAbstractSAparallelLinearOperator) = A+(-B)

# -(jo,mvec)

# -(mvec,jo)

# -(jo,vec)

# -(vec,jo)

# -(jo,num)
-(A::joAbstractSAparallelLinearOperator,b::Number) = A+(-b)
-(A::joAbstractSAparallelLinearOperator,b::joNumber) = A+(-b)

# -(num,jo)
-(b::Number,A::joAbstractSAparallelLinearOperator) = -A+b
-(b::joNumber,A::joAbstractSAparallelLinearOperator) = -A+b

############################################################
## overloaded Base .*(...jo...)

# .*(num,jo)
Base.Broadcast.broadcasted(::typeof(*),a::Number,A::joAbstractSAparallelLinearOperator) = a*A
Base.Broadcast.broadcasted(::typeof(*),a::joNumber,A::joAbstractSAparallelLinearOperator) = a*A

# .*(jo,num)
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractSAparallelLinearOperator,a::Number) = a*A
Base.Broadcast.broadcasted(::typeof(*),A::joAbstractSAparallelLinearOperator,a::joNumber) = a*A

############################################################
## overloaded Base .\(...jo...)

############################################################
## overloaded Base .+(...jo...)

# .+(jo,num)
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractSAparallelLinearOperator,b::Number) = A+b
Base.Broadcast.broadcasted(::typeof(+),A::joAbstractSAparallelLinearOperator,b::joNumber) = A+b

# .+(num,jo)
Base.Broadcast.broadcasted(::typeof(+),b::Number,A::joAbstractSAparallelLinearOperator) = A+b
Base.Broadcast.broadcasted(::typeof(+),b::joNumber,A::joAbstractSAparallelLinearOperator) = A+b

############################################################
## overloaded Base .-(...jo...)

# .-(jo,num)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractSAparallelLinearOperator,b::Number) = A+(-b)
Base.Broadcast.broadcasted(::typeof(-),A::joAbstractSAparallelLinearOperator,b::joNumber) = A+(-b)

# .-(num,jo)
Base.Broadcast.broadcasted(::typeof(-),b::Number,A::joAbstractSAparallelLinearOperator) = -A+b
Base.Broadcast.broadcasted(::typeof(-),b::joNumber,A::joAbstractSAparallelLinearOperator) = -A+b

############################################################
## overloaded Base block methods

# hcat(...jo...)

# vcat(...jo...)

# hvcat(...jo...)

############################################################
## overloaded LinearAlgebra functions

# mul!(...,jo,...)

# ldiv!(...,jo,...)


text:
function normalize(array, target_range)
    target_min, target_max = target_range
    orig_min, orig_max = minimum(array), maximum(array)
    normalized_array = target_min .+ ((target_max-target_min).*(array.-orig_min))./(orig_max-orig_min)
    return normalized_array
end
a = [1,2,3,4]
print(normalize(a, [-1,1]))


text:
struct StorePutKey{T} <: ResourceKey
  priority :: Int
  id :: UInt
  item :: T
  StorePutKey{T}(priority, id, item) where T = new(priority, id, item)
end

struct StoreGetKey <: ResourceKey
  priority :: Int
  id :: UInt
  filter :: Function
end

mutable struct Store{T} <: AbstractResource
  env :: Environment
  capacity :: UInt
  load :: UInt
  items :: Dict{T, UInt}
  seid :: UInt
  put_queue :: DataStructures.PriorityQueue{Put, StorePutKey{T}}
  get_queue :: DataStructures.PriorityQueue{Get, StoreGetKey}
  function Store{T}(env::Environment; capacity::UInt=typemax(UInt)) where {T}
    new(env, capacity, zero(UInt), Dict{T, UInt}(), zero(UInt), DataStructures.PriorityQueue{Put, StorePutKey{T}}(), DataStructures.PriorityQueue{Get, StoreGetKey}())
  end
end

function put(sto::Store{T}, item::T; priority::Int=0) where T
  put_ev = Put(sto.env)
  sto.put_queue[put_ev] = StorePutKey{T}(priority, sto.seid+=one(UInt), item)
  @callback trigger_get(put_ev, sto)
  trigger_put(put_ev, sto)
  put_ev
end

get_any_item(::T) where T = true

function get(sto::Store{T}, filter::Function=get_any_item; priority::Int=0) where T
  get_ev = Get(sto.env)
  sto.get_queue[get_ev] = StoreGetKey(priority, sto.seid+=one(UInt), filter)
  @callback trigger_put(get_ev, sto)
  trigger_get(get_ev, sto)
  get_ev
end

function do_put(sto::Store{T}, put_ev::Put, key::StorePutKey{T}) where {T}
  if sto.load < sto.capacity
    sto.load += one(UInt)
    sto.items[key.item] = get(sto.items, key.item, zero(UInt)) + one(UInt)
    schedule(put_ev)
  end
  false
end

function do_get(sto::Store{T}, get_ev::Get, key::StoreGetKey) where {T}
  for (item, number) in sto.items
    if key.filter(item)
      sto.load -= one(UInt)
      if number === one(UInt)
        delete!(sto.items, item)
      else
        sto.items[item] = number - one(UInt)
      end
      schedule(get_ev; value=item)
      break
    end
  end
  true
end


text:
module MuttsDicts

export MuttsDict
export branch!, double_branch!, mark_immutable!, is_mutable, get_mutable_version

# TODO: make the default constructor mutable, just too weird for newbies
# TODO: explicitly provide overloads for Mutts.branch! etc.?
# TODO: NodeConfig mutable for Serialization.jl sharing? or just ugh.
# TODO: with xor probing, do not need modular, just bound test+subtract.
# TODO: add functional insert/delete?
# TODO: add efficient == operation on dicts (factor out some common code with setdiff.)
# TODO: commented-out asserts should be somehow made optional like @dassert1.
# TODO: hash and equality functions for MuttsDict instances themselves.
# TODO: use split K[]/V[] arrays in Leaf to avoid allocations for isbitstypes?
# TODO: pretty printing
# TODO: unit test for dict[k1, k2, ks...]
# TODO: MuttsDict{Int,Int}(2=>3, 1=>5, ...).
# TODO: size_hint! to preemptively set NodeConfig before populating

include("debug.jl")
include("config.jl")
include("leaf.jl")
include("inode.jl")

# MuttsDict is implemented using a variation on the hashtrie. A MuttsDict object contains
# a root pointer, a population count n, and a NodeConfig that identifies the structure of
# the hashtrie. When a new MuttsDict is created, it initially consists of a MuttsDict
# object pointing to a single Leaf object:
#
# MuttsDict -----> Leaf
#
# The Leaf is a tiny hash table that can accommodate a small number of entries.
# The leaf hash table is keyed by bits 0-15 of the 64-bit key hash.
#
# Once the dictionary reaches size n=16, we change the configuration to incorporate
# a single INode:
#
# MuttsDict ------> INode +-----> Leaf
#                         |-----> Leaf
#                         |-----> Leaf
#                         +-----> Leaf
#
# This INode initially has branch factor 4, with the branch selected by bits 16-17 of the
# key hash.  We use this configuration until we reach n=64, at which point we add another
# INode level with branch factor 4, and the branch selected by bits 24-25 of the key hash:
#
# MuttsDict ------> INode +-----> INode +-----> Leaf
#                         |             |-----> Leaf
#                         |             |-----> Leaf
#                         |             +-----> Leaf
#                         +-----> INode ... (4 Leafs)
#                         |-----> INode ... (4 Leafs)
#                         +-----> INode ... (4 Leafs)
#
# This process continues until there are 6 levels of INodes.  At this point we stop
# adding INode levels, and instead begin increasing the branch factors of the INodes.
# The NodeConfig.inode_capacity[] array contains the branch factors for the INodes:
# initially it is [] (no inodes), then [4,] (a single inode with branch factor 4),
# then [4,4] (two inode levels with branch factor 4), etc. Once we reach [4,4,4,4,4,4],
# we go to [4,4,4,4,4,8], then [4,4,4,4,8,8], etc. Branch factors are always a power of 2.
#
# The increases in bits per INode level are arranged so the branch factor at each INode
# level is Θ(n^(1/7)), with leaf hash tables containing Θ(n^(1/7)) entries on average.
#
# For each Leaf and INode we track whether it is mutable. If a Leaf is mutable, we can
# usually insert a new (k,v) pair by finding a free slot and placing it there. If a Leaf
# is immutable, we copy the Leaf to make a mutable version, then insert the new entry
# in the mutable version. The new leaf pointer is stored in the parent INode, which itself
# may first have to be copied if it is immutable, and so on up to the root.
#
# The 64-bit key hash is partitioned among the INode levels and Leaf like this:
#
# 0000000000000000000000000000000000000000000000000000000000000000
# +-------+-------+-------+-------+-------+-------+--------------+
# inode1  inode2  inode3  inode4  inode5  inode6  leaf
#
# That is, bits 0-15 are for the leaf, bits 16-23 for the last inode level, bits 24-31
# for the second last inode level, etc.
#
# As the NodeConfig changes with increasing n, we use a laziness technique to avoid
# expensive restructuring of the entire hashtrie. To increase the capacity of an INode
# from, say, 4 entries to 8, we mark the subtrees immutable and then have two paths to
# each child:
#
# INode +-----> child0                INode +-----> child0
#       |-----> child1                      |-----> child1
#       |-----> child2    becomes           |-----> child2
#       +-----> child3                      |-----> child3
#                                           |-----> child0
#                                           |-----> child1
#                                           |-----> child2
#                                           +-----> child3
#
# Suppose some key k lives in the subtree child0. Before resizing, the child0 branch
# was taken when the single byte of the hash for this INode matched the pattern
# xxxxxx00. After resizing, there are two paths to child0: xxxxx000 and xxxxx100.
# So whether the next bit of hash(k) is a 0 or 1, we will still reach child0 and
# find key k there. When the next insert occurs to child0, we will do a copy-on-write
# of a leaf, but only retain those keys that match the path used to reach that child.
# We track the path using pathhash/pathmask parameters; an example of a hash path:
#      xxxxxx00xxxxxx11xxxxxx01xxxxxx10xxxxx110xxxxx101xxxxxxxxxxxxxxxx
# which indicates we took child 00 at the first INode, child 11 at the second, etc.
# The hash path is encoded by the pathhash/pathmask parameters:
#      0000000000000011000000010000001000000110000001010000000000000000
#      0000001100000011000000110000001100000111000001110000000000000000
# A 64-bit hash code matches the hash path if (hash & pathmask)==pathhash.
#
# The lazy restructuring means that the hashtrie is really a hash *graph*, with multiple
# paths to reach some leaf nodes. The benefit is that we achieve worst case Θ(n^(1/7)) for
# inserts. This improves on Base.Dict's Θ(n) worst case insert time.

"""
MuttsDict{K,V} provides a dictionary similar to Julia's Dict{K,V}, but with versioning
and improved worst-case asymptotics.  Lookups (getindex) take O(1) time.  Inserts
and updates (`setindex!`) take O(1) amortized time on unversioned dictionaries. You can
create a mutable copy of a dictionary (branch a new version) in O(1) amortized time,
and `setindex!` on a fresh branch is Θ(n^(1/7)).

MuttsDict has been optimized for both space use and cost of immutable inserts. Space use
for small dictionaries is low. For real-time or low-latency applications, MuttsDict
offers a Θ(n^(1/7)) worst case insert/delete time, which is an improvement over
Base.Dict's Θ(n) worst case.

MuttsDict implements the standard dictionary methods `getindex`, `setindex!`, `length`,
`setdiff`, `delete!`, and `iterate`. It is not a subtype of `AbstractDict`, due to
differences between the MuttsDict semantics and those expected by standard library
functions on `AbstractDict`.

MuttsDict implements the "mutable until shared" philosophy for multithreaded programming.
An object local to a single thread can be mutable, and enjoy the more efficient O(1)
amortized insert/delete performance. When objects become shared, they should first be
made immutable (by calling `mark_immutable!` or `branch!`) before sharing. Immutable
MuttsDicts are read-only, and can be safely shared among multiple threads. For lockfree
operation, a thread can branch a shared dictionary, modify it, mark it immutable, then
write the pointer back with an atomic compare-and-swap.

Branching and versioning is handled by the mutts functions `branch!`, `mark_immutable!`,
`is_mutable`, and `get_mutable_version`:
- `dict2=branch!(dict1)` returns a new mutable copy of `dict1`, and marks `dict1`
immutable.
- `mark_immutable!(dict)` marks `dict` as immutable. Once immutable, a MuttsDict can
safely be shared by multiple threads.
- `is_mutable(dict)` returns true when `dict` is mutable.
- `get_mutable_version(dict)` returns `dict` if it is mutable, otherwise returns
`branch!(dict)`.
"""
mutable struct MuttsDict{K,V} #<: AbstractDict{K,V}
    # The root is usually an inode, but for tiny dictionaries it may just be a leaf.
    root::Union{INode{K,V},Leaf{K,V}}

    # Exact population count
    n::Int

    # The height of the tree, and branch factor for each inode.  The NodeConfig is
    # adjusted as n increases.
    config::NodeConfig

    # Note: if you have many dictionaries of the same type which are by default empty,
    # it uses less space to use a single empty, immutable instance as the default, e.g.
    # empty_dict = mark_immutable!(MuttsDict{Int,Int}())
    # dicts = [empty_dict for i=1:10]
    function MuttsDict{K,V}() where {K,V}
        return MuttsDict{K,V}(Leaf{K,V}(1), 0, _config_for_size(1))
    end

    function MuttsDict{K,V}(root,n,config) where {K,V}
        new(root, n, config)
    end
end

"""
    branch!(dict)

Mark dict as immutable, and return a new mutable copy. This operation takes O(1)
time. Changes made to the returned mutable dictionary do not affect `dict`.
"""
function branch!(dict::MuttsDict{K,V}) where {K,V}
    root = dict.root
    mark_immutable!(root)

    # Note: throughout the MuttsDict implementation we use the following
    # pattern for Union{INode,Leaf} types:
    #    if isa(x,Leaf)
    #       ... # leaf stuff I guess
    #    else
    #       ... # inode stuff
    #    end
    # This forces an efficient specialization of the first branch with the compiler
    # knowing x::Leaf, and the second branch with the compiler knowing x::INode.
    # In some cases it is the exact same code in each branch, but it compiles to more
    # efficient code with the above pattern. Julia can often do a similar style of
    # specialization left to its own devices, but it is not guaranteed to happen,
    # and in my experience it is prone to stop happening at deeply inlined loops
    # where performance is critical.
    if isa(root, Leaf{K,V})
        pathmask = UInt64(0)
        pathhash = UInt64(0)
        new_root = _get_mutable_version(root, pathmask, pathhash)
        return MuttsDict{K,V}(new_root, dict.n, dict.config)
    else
        height = length(dict.config.inode_capacity)
        new_root = _get_mutable_version(root, dict.config, height)
        return MuttsDict{K,V}(new_root, dict.n, dict.config)
    end
end

"""
    double_branch!(dict)

Returns two mutable branches of dict. This is useful when you want to create
a mutable branch but keep the original dict mutable:
```
(dict,dict_branch) = double_branch!(dict)
```
"""
function double_branch!(dict::MuttsDict{K,V}) where {K,V}
    return (branch!(dict), branch!(dict))
end

"""
    mark_immutable!(dict)

Marks dict as immutable.
"""
function mark_immutable!(dict::MuttsDict{K,V}) where {K,V}
    root = dict.root
    if isa(root, Leaf{K,V})
        mark_immutable!(root)
    else
        mark_immutable!(root)
    end
end

"""
    is_mutable(dict)

Returns true just when dict is mutable.
"""
function is_mutable(dict::MuttsDict{K,V}) where {K,V}
    root = dict.root
    if isa(root, Leaf{K,V})
        # Branch specialized for isa(root,Leaf{K,V})
        return is_mutable(root)
    else
        # Branch specialized for isa(root,INode{K,V})
        return is_mutable(root)
    end
end

"""
    get_mutable_version(dict)

Returns dict if `is_mutable(dict)`, otherwise returns `branch!(dict)`.
"""
function get_mutable_version(dict::MuttsDict{K,V}) where {K,V}
    if is_mutable(dict)
        return dict
    else
        return branch!(dict)
    end
end

function Base.length(dict::MuttsDict{K,V}) where {K,V}
    return dict.n
end

function Base.copy(dict::MuttsDict{K,V}) where {K,V}
    # Base.copy is incompatible with MuttsDict, because it can result in two competing
    # MuttsDict objects sharing the same underlying tree, with the result that the
    # population count can become inconsistent.
    error("Do not use Base.copy() with MuttsDict{K,V}. Instead use branch!(dict).")
end

function Base.show(io::IO, m::MuttsDict{K,V}) where {K,V}
    print(io, "MuttsDict(")
    entries = [(k,v) for (k,v) in m]
    sep=""
    for (k,v) in entries
        print(io, sep * "$(k)=>$(v)")
        sep=","
    end
    print(io, ")")
end

include("get.jl")
include("insert.jl")
include("delete.jl")
include("setdiff.jl")
include("iterate.jl")
include("merge.jl")
include("tools.jl")

end


text:
struct RuntimeDevice{D}
    device::D
end
default_device() = RuntimeDevice(default_device(RUNTIME[]))
default_device(::typeof(HSA_rt)) = get_default_agent()

struct RuntimeQueue{Q}
    queue::Q
end
default_queue() = default_queue(default_device())
default_queue(device::RuntimeDevice) = RuntimeQueue(default_queue(RUNTIME[], device))
default_queue(::typeof(HSA_rt), device::RuntimeDevice) =
    get_default_queue(device.device)
get_device(queue::RuntimeQueue{HSAQueue}) = RuntimeDevice(queue.queue.agent)

default_isa(device::RuntimeDevice{HSAAgent}) =
    first(isas(device.device))

struct RuntimeEvent{E}
    event::E
end
create_event(exe; kwargs...) = RuntimeEvent(create_event(RUNTIME[], exe; kwargs...))
Base.wait(event::RuntimeEvent; kwargs...) = wait(event.event; kwargs...)

struct RuntimeExecutable{E}
    exe::E
end
create_executable(device, entry, obj; globals=()) =
    RuntimeExecutable(create_executable(RUNTIME[], device, entry, obj; globals=globals))
function create_executable(::typeof(HSA_rt), device, entry, obj; globals=())
    # link with ld.lld
    @assert ld_lld_path != "" "ld.lld was not found; cannot link kernel"
    path_exe = mktemp() do path_o, io_o
        write(io_o, obj)
        flush(io_o)
        path_exe = path_o*".exe"
        run(`$ld_lld_path -shared -o $path_exe $path_o`)
        path_exe
    end
    data = read(path_exe)
    rm(path_exe)

    return HSAExecutable(device.device, data, entry; globals=globals)
end
get_global(exe::RuntimeExecutable, sym::Symbol) =
    get_global(exe.exe, sym)

create_event(::typeof(HSA_rt), exe::RuntimeExecutable{<:HSAExecutable}; queue=default_queue(), kwargs...) =
    HSAStatusSignal(HSASignal(), exe.exe, queue.queue)

struct RuntimeKernel{K}
    kernel::K
end
create_kernel(device, exe, entry, args) =
    RuntimeKernel(create_kernel(RUNTIME[], device, exe, entry, args))
create_kernel(::typeof(HSA_rt), device, exe, entry, args) =
    HSAKernelInstance(device.device, exe.exe, entry, args)
launch_kernel(queue, kern, event; kwargs...) =
    launch_kernel(RUNTIME[], queue, kern, event; kwargs...)
function launch_kernel(::typeof(HSA_rt), queue, kern, event;
                       groupsize=nothing, gridsize=nothing)
    signal = event.event isa HSAStatusSignal ? event.event.signal : event.event
    launch!(queue.queue, kern.kernel, signal;
                       workgroup_size=groupsize, grid_size=gridsize)
end
barrier_and!(queue, events::Vector{<:RuntimeEvent}) = barrier_and!(queue, map(x->x.event,events))
barrier_or!(queue, events::Vector{<:RuntimeEvent}) = barrier_or!(queue, map(x->x.event,events))


text:
export replace_math

"""
    replace_math(path::AbstractString)
    replace_math(file::File)
    replace_math(dir::Dir)

Replace markdown inline math such as `\$x = 10\$` with Markua equivalent,
as well as math blocks with double dollar sign.
"""
replace_math(path::AbstractString) = replace_math(Path(path))

function replace_math(file::File)
    out = IOBuffer()
    open(file.path) do io
        while !eof(io)
            s = readuntil(io, '$')
            print(out, s)

            eof(io) && break

            mark(io)
            s = readuntil(io, '$')
            if isempty(s)
                reset(io)
                s = readuntil(io, raw"$$")
                println(out, "[stem]")
                print(out, "++++")
                print(out, s[2:end])
                print(out, "++++")
            else
                # Finally ready to process
                print(out, "stem:[")
                print(out, s)
                print(out, "]")
            end
        end
    end
    seekstart(out)
    s = read(out, String)
    open(file.path, "w") do io
        write(io, s)
    end
end

replace_math(dir::Dir) = visitdir(replace_math, dir.path)

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
using Test
using LinearAlgebra

include("./KTBC.jl")
using .KTBC: CreateKTBC

function CreateK2D(dimSquared=2)
  IdentityFixedDims = I(dimSquared)
  K, T, B, C = CreateKTBC(dimSquared)
  return kron(IdentityFixedDims, K) + kron(K, IdentityFixedDims)
end

IMatrixTestDims = I(2)

@test CreateK2D() == [
   K + 2 * IMatrixTestDims   -IMatrixTestDims
  -IMatrixTestDims            K + 2 * IMatrixTestDims
]

export CreateK2D


text:
function ODMFile()
    return ODMFile(())
end

function get_mu_created(obj::ODMFile)
    return jcall(obj, "getMuCreated", jdouble, ())
end

function get_mu_parsed(obj::ODMFile)
    return jcall(obj, "getMuParsed", jdouble, ())
end

function get_mu_set(obj::ODMFile)
    return jcall(obj, "getMuSet", jdouble, ())
end

function get_mu_used(obj::ODMFile)
    return jcall(obj, "getMuUsed", jdouble, ())
end

function set_mu_set(obj::ODMFile, arg0::jdouble)
    return jcall(obj, "setMuSet", void, (jdouble,), arg0)
end


text:
read_string(x, length_type::Type=UInt32) = Base.read(x, length_type) |> Base.read $ x |> String
function write_string(io::IO, str::String, length_type::Type=UInt32)
    Base.write(io,length_type(ncodeunits(str)))
    Base.write(io,str)
end

function unsafe_write(io::IO, obj)
    @assert Base.isbitstype(typeof(obj))
    Base.unsafe_write(io,Ref(obj),Base.sizeof(obj))
end

function unsafe_read(io::IO,::Type{T}) where{T}
    @assert Base.isbitstype(T)
    val = Ref{T}()
    Base.unsafe_read(io,val, Base.sizeof(T))
    return val[]
end

write(io::IO, val::T) where {T <: Real} = unsafe_write(io,val)
read(io::IO, type::Type{T}) where {T <: Real} = unsafe_real(io,type)

write(io::IO, obj::Tuple{Vararg{Real,N}}) where {N} = unsafe_write(io, obj)
read(io::IO, ::Type{Tuple{Vararg{Real,N}}}) where {N} = unsafe_read(io, Tuple{Vararg{Real,N}})

# function write(io::IO, obj::T ) where {T}
#     names = fieldnames(T)
#     for name in names
#         write(io, getfield(obj,name))
#     end
# end

# "Base implementation of read. Will simply loop over all the fields and call read for each of them"
# @generated function read(io::IO, ::Type{T}) where {T}
#     args = T.types .|> t ->  Expr(:call,:read, io,t)
# 	return :(T($(args...)))
# end

# Ok so this might not work. How does Julia evaluate union types at runtime? Need to investigate

function write(io::IO, obj::Vector{T} ) where {T<: Real}
    Base.write(io,UInt64(length(obj)))
    Base.write(io,obj)
end

function read(io::IO, ::Type{Vector{T}}) where {T<: Real}
    length = Base.read(io,UInt64)
    data = zeros(T, length)
    Base.read!(io,data)
    return data
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Sockets

export
    accept,
    bind,
    connect,
    getaddrinfo,
    getalladdrinfo,
    getnameinfo,
    getipaddr,
    getpeername,
    getsockname,
    listen,
    listenany,
    recv,
    recvfrom,
    send,
    TCPSocket,
    UDPSocket,
    @ip_str,
    IPAddr,
    IPv4,
    IPv6

import Base: isless, show, print, parse, bind, convert, isreadable, iswritable, alloc_buf_hook, _uv_hook_close

using Base: LibuvStream, LibuvServer, PipeEndpoint, @handle_as, uv_error, associate_julia_struct, uvfinalize,
    notify_error, stream_wait, uv_req_data, uv_req_set_data, preserve_handle, unpreserve_handle, UVError,
    eventloop, StatusUninit, StatusInit, StatusConnecting, StatusOpen, StatusClosing, StatusClosed, StatusActive,
    uv_status_string, check_open, wait_connected,
    UV_EINVAL, UV_ENOMEM, UV_ENOBUFS, UV_EAGAIN, UV_ECONNABORTED, UV_EADDRINUSE, UV_EACCES, UV_EADDRNOTAVAIL,
    UV_EAI_ADDRFAMILY, UV_EAI_AGAIN, UV_EAI_BADFLAGS,
    UV_EAI_BADHINTS, UV_EAI_CANCELED, UV_EAI_FAIL,
    UV_EAI_FAMILY, UV_EAI_NODATA, UV_EAI_NONAME,
    UV_EAI_OVERFLOW, UV_EAI_PROTOCOL, UV_EAI_SERVICE,
    UV_EAI_SOCKTYPE, UV_EAI_MEMORY

include("IPAddr.jl")
include("addrinfo.jl")

"""
    TCPSocket(; delay=true)

Open a TCP socket using libuv. If `delay` is true, libuv delays creation of the
socket's file descriptor till the first [`bind`](@ref) call. `TCPSocket` has various
fields to denote the state of the socket as well as its send/receive buffers.
"""
mutable struct TCPSocket <: LibuvStream
    handle::Ptr{Cvoid}
    status::Int
    buffer::IOBuffer
    readnotify::Condition
    connectnotify::Condition
    closenotify::Condition
    sendbuf::Union{IOBuffer, Nothing}
    lock::ReentrantLock
    throttle::Int

    function TCPSocket(handle::Ptr{Cvoid}, status)
        tcp = new(
                handle,
                status,
                PipeBuffer(),
                Condition(),
                Condition(),
                Condition(),
                nothing,
                ReentrantLock(),
                Base.DEFAULT_READ_BUFFER_SZ)
        associate_julia_struct(tcp.handle, tcp)
        finalizer(uvfinalize, tcp)
        return tcp
    end
end

# kw arg "delay": if true, libuv delays creation of the socket fd till the first bind call
function TCPSocket(; delay=true)
    tcp = TCPSocket(Libc.malloc(Base._sizeof_uv_tcp), StatusUninit)
    af_spec = delay ? 0 : 2   # AF_UNSPEC is 0, AF_INET is 2
    err = ccall(:uv_tcp_init_ex, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint),
                eventloop(), tcp.handle, af_spec)
    uv_error("failed to create tcp socket", err)
    tcp.status = StatusInit
    return tcp
end

mutable struct TCPServer <: LibuvServer
    handle::Ptr{Cvoid}
    status::Int
    connectnotify::Condition
    closenotify::Condition

    function TCPServer(handle::Ptr{Cvoid}, status)
        tcp = new(
            handle,
            status,
            Condition(),
            Condition())
        associate_julia_struct(tcp.handle, tcp)
        finalizer(uvfinalize, tcp)
        return tcp
    end
end

# Keyword arg "delay": if true, libuv delays creation of socket fd till bind.
# It can be set to false if there is a need to set socket options before
# further calls to `bind` and `listen`, e.g. `SO_REUSEPORT`.
function TCPServer(; delay=true)
    tcp = TCPServer(Libc.malloc(Base._sizeof_uv_tcp), StatusUninit)
    af_spec = delay ? 0 : 2   # AF_UNSPEC is 0, AF_INET is 2
    err = ccall(:uv_tcp_init_ex, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cuint),
                eventloop(), tcp.handle, af_spec)
    uv_error("failed to create tcp server", err)
    tcp.status = StatusInit
    return tcp
end

isreadable(io::TCPSocket) = isopen(io) || bytesavailable(io) > 0
iswritable(io::TCPSocket) = isopen(io) && io.status != StatusClosing

"""
    accept(server[,client])

Accepts a connection on the given server and returns a connection to the client. An
uninitialized client stream may be provided, in which case it will be used instead of
creating a new stream.
"""
accept(server::TCPServer) = accept(server, TCPSocket())

# UDP
"""
    UDPSocket()

Open a UDP socket using libuv. `UDPSocket` has various
fields to denote the state of the socket.
"""
mutable struct UDPSocket <: LibuvStream
    handle::Ptr{Cvoid}
    status::Int
    recvnotify::Condition
    sendnotify::Condition
    closenotify::Condition

    function UDPSocket(handle::Ptr{Cvoid}, status)
        udp = new(
            handle,
            status,
            Condition(),
            Condition(),
            Condition())
        associate_julia_struct(udp.handle, udp)
        finalizer(uvfinalize, udp)
        return udp
    end
end
function UDPSocket()
    this = UDPSocket(Libc.malloc(Base._sizeof_uv_udp), StatusUninit)
    err = ccall(:uv_udp_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                eventloop(), this.handle)
    uv_error("failed to create udp socket", err)
    this.status = StatusInit
    return this
end

show(io::IO, stream::UDPSocket) = print(io, typeof(stream), "(", uv_status_string(stream), ")")

function _uv_hook_close(sock::UDPSocket)
    sock.handle = C_NULL
    sock.status = StatusClosed
    notify(sock.closenotify)
    notify(sock.sendnotify)
    notify_error(sock.recvnotify,EOFError())
end

# Disables dual stack mode.
const UV_TCP_IPV6ONLY = 1

# Disables dual stack mode. Only available when using ipv6 binf
const UV_UDP_IPV6ONLY = 1

# Indicates message was truncated because read buffer was too small. The
# remainder was discarded by the OS.
const UV_UDP_PARTIAL = 2

# Indicates if SO_REUSEADDR will be set when binding the handle in uv_udp_bind. This sets
# the SO_REUSEPORT socket flag on the BSDs and OS X. On other Unix platforms, it sets the
# SO_REUSEADDR flag. What that means is that multiple threads or processes can bind to the
# same address without error (provided they all set the flag) but only the last one to bind
# will receive any traffic, in effect "stealing" the port from the previous listener.
const UV_UDP_REUSEADDR = 4

##

_bind(sock::TCPServer, host::IPv4, port::UInt16, flags::UInt32 = UInt32(0)) = ccall(:jl_tcp_bind, Int32, (Ptr{Cvoid}, UInt16, UInt32, Cuint),
            sock.handle, hton(port), hton(host.host), flags)

_bind(sock::TCPServer, host::IPv6, port::UInt16, flags::UInt32 = UInt32(0)) = ccall(:jl_tcp_bind6, Int32, (Ptr{Cvoid}, UInt16, Ptr{UInt128}, Cuint),
            sock.handle, hton(port), Ref(hton(host.host)), flags)

_bind(sock::UDPSocket, host::IPv4, port::UInt16, flags::UInt32 = UInt32(0)) = ccall(:jl_udp_bind, Int32, (Ptr{Cvoid}, UInt16, UInt32, UInt32),
            sock.handle, hton(port), hton(host.host), flags)

_bind(sock::UDPSocket, host::IPv6, port::UInt16, flags::UInt32 = UInt32(0)) = ccall(:jl_udp_bind6, Int32, (Ptr{Cvoid}, UInt16, Ptr{UInt128}, UInt32),
            sock.handle, hton(port), Ref(hton(host.host)), flags)

"""
    bind(socket::Union{UDPSocket, TCPSocket}, host::IPAddr, port::Integer; ipv6only=false, reuseaddr=false, kws...)

Bind `socket` to the given `host:port`. Note that `0.0.0.0` will listen on all devices.

* The `ipv6only` parameter disables dual stack mode. If `ipv6only=true`, only an IPv6 stack is created.
* If `reuseaddr=true`, multiple threads or processes can bind to the same address without error
  if they all set `reuseaddr=true`, but only the last to bind will receive any traffic.
"""
function bind(sock::Union{TCPServer, UDPSocket}, host::IPAddr, port::Integer; ipv6only = false, reuseaddr = false, kws...)
    if sock.status != StatusInit
        error("$(typeof(sock)) is not in initialization state")
    end
    flags = 0
    if isa(host,IPv6) && ipv6only
        flags |= isa(sock, UDPSocket) ? UV_UDP_IPV6ONLY : UV_TCP_IPV6ONLY
    end
    if isa(sock, UDPSocket) && reuseaddr
        flags |= UV_UDP_REUSEADDR
    end
    err = _bind(sock, host, UInt16(port), UInt32(flags))
    if err < 0
        if err != UV_EADDRINUSE && err != UV_EACCES && err != UV_EADDRNOTAVAIL
            #TODO: this codepath is not currently tested
            throw(UVError("bind", err))
        else
            return false
        end
    end
    sock.status = StatusOpen
    isa(sock, UDPSocket) && setopt(sock; kws...)
    return true
end

bind(sock::TCPServer, addr::InetAddr) = bind(sock, addr.host, addr.port)

"""
    setopt(sock::UDPSocket; multicast_loop = nothing, multicast_ttl=nothing, enable_broadcast=nothing, ttl=nothing)

Set UDP socket options.

* `multicast_loop`: loopback for multicast packets (default: `true`).
* `multicast_ttl`: TTL for multicast packets (default: `nothing`).
* `enable_broadcast`: flag must be set to `true` if socket will be used for broadcast
  messages, or else the UDP system will return an access error (default: `false`).
* `ttl`: Time-to-live of packets sent on the socket (default: `nothing`).
"""
function setopt(sock::UDPSocket; multicast_loop = nothing, multicast_ttl=nothing, enable_broadcast=nothing, ttl=nothing)
    if sock.status == StatusUninit
        error("Cannot set options on uninitialized socket")
    end
    if multicast_loop !== nothing
        uv_error("multicast_loop",ccall(:uv_udp_set_multicast_loop,Cint,(Ptr{Cvoid},Cint),sock.handle,multicast_loop) < 0)
    end
    if multicast_ttl !== nothing
        uv_error("multicast_ttl",ccall(:uv_udp_set_multicast_ttl,Cint,(Ptr{Cvoid},Cint),sock.handle,multicast_ttl))
    end
    if enable_broadcast !== nothing
        uv_error("enable_broadcast",ccall(:uv_udp_set_broadcast,Cint,(Ptr{Cvoid},Cint),sock.handle,enable_broadcast))
    end
    if ttl !== nothing
        uv_error("ttl",ccall(:uv_udp_set_ttl,Cint,(Ptr{Cvoid},Cint),sock.handle,ttl))
    end
end

"""
    recv(socket::UDPSocket)

Read a UDP packet from the specified socket, and return the bytes received. This call blocks.
"""
function recv(sock::UDPSocket)
    addr, data = recvfrom(sock)
    return data
end

"""
    recvfrom(socket::UDPSocket) -> (address, data)

Read a UDP packet from the specified socket, returning a tuple of `(address, data)`, where
`address` will be either IPv4 or IPv6 as appropriate.
"""
function recvfrom(sock::UDPSocket)
    # If the socket has not been bound, it will be bound implicitly to ::0 and a random port
    if sock.status != StatusInit && sock.status != StatusOpen && sock.status != StatusActive
        error("UDPSocket is not initialized and open")
    end
    if ccall(:uv_is_active, Cint, (Ptr{Cvoid},), sock.handle) == 0
        uv_error("recv_start", ccall(:uv_udp_recv_start, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
                                    sock.handle, Base.uv_jl_alloc_buf::Ptr{Cvoid}, uv_jl_recvcb::Ptr{Cvoid}))
    end
    sock.status = StatusActive
    return stream_wait(sock, sock.recvnotify)::Tuple{Union{IPv4, IPv6}, Vector{UInt8}}
end

alloc_buf_hook(sock::UDPSocket, size::UInt) = (Libc.malloc(size), size)

function uv_recvcb(handle::Ptr{Cvoid}, nread::Cssize_t, buf::Ptr{Cvoid}, addr::Ptr{Cvoid}, flags::Cuint)
    # C signature documented as (*uv_udp_recv_cb)(...)
    sock = @handle_as handle UDPSocket
    if nread < 0
        Libc.free(buf_addr)
        notify_error(sock.recvnotify, UVError("recv", nread))
    elseif flags & UV_UDP_PARTIAL > 0
        Libc.free(buf_addr)
        notify_error(sock.recvnotify, "Partial message received")
    else
        buf_addr = ccall(:jl_uv_buf_base, Ptr{Cvoid}, (Ptr{Cvoid},), buf)
        buf_size = ccall(:jl_uv_buf_len, Csize_t, (Ptr{Cvoid},), buf)
        # need to check the address type in order to convert to a Julia IPAddr
        addrout = if addr == C_NULL
                      IPv4(0)
                  elseif ccall(:jl_sockaddr_in_is_ip4, Cint, (Ptr{Cvoid},), addr) == 1
                      IPv4(ntoh(ccall(:jl_sockaddr_host4, UInt32, (Ptr{Cvoid},), addr)))
                  else
                      tmp = [UInt128(0)]
                      ccall(:jl_sockaddr_host6, UInt32, (Ptr{Cvoid}, Ptr{UInt8}), addr, pointer(tmp))
                      IPv6(ntoh(tmp[1]))
                  end
        buf = unsafe_wrap(Array, convert(Ptr{UInt8}, buf_addr), Int(nread), own = true)
        notify(sock.recvnotify, (addrout, buf))
    end
    ccall(:uv_udp_recv_stop, Cint, (Ptr{Cvoid},), sock.handle)
    sock.status = StatusOpen
    nothing
end

function _send(sock::UDPSocket, ipaddr::IPv4, port::UInt16, buf)
    ccall(:jl_udp_send, Cint, (Ptr{Cvoid}, UInt16, UInt32, Ptr{UInt8}, Csize_t, Ptr{Cvoid}),
          sock.handle, hton(port), hton(ipaddr.host), buf, sizeof(buf), uv_jl_sendcb::Ptr{Cvoid})
end

function _send(sock::UDPSocket, ipaddr::IPv6, port::UInt16, buf)
    ccall(:jl_udp_send6, Cint, (Ptr{Cvoid}, UInt16, Ref{UInt128}, Ptr{UInt8}, Csize_t, Ptr{Cvoid}),
          sock.handle, hton(port), hton(ipaddr.host), buf, sizeof(buf), uv_jl_sendcb::Ptr{Cvoid})
end

"""
    send(socket::UDPSocket, host, port::Integer, msg)

Send `msg` over `socket` to `host:port`.
"""
function send(sock::UDPSocket,ipaddr,port,msg)
    # If the socket has not been bound, it will be bound implicitly to ::0 and a random port
    if sock.status != StatusInit && sock.status != StatusOpen && sock.status != StatusActive
        error("UDPSocket is not initialized and open")
    end
    uv_error("send", _send(sock, ipaddr, UInt16(port), msg))
    stream_wait(sock, sock.sendnotify)
    nothing
end

function uv_sendcb(handle::Ptr{Cvoid}, status::Cint)
    sock = @handle_as handle UDPSocket
    if status < 0
        notify_error(sock.sendnotify, UVError("UDP send failed", status))
    end
    notify(sock.sendnotify)
    Libc.free(handle)
    nothing
end

#from `connect`
function uv_connectcb(conn::Ptr{Cvoid}, status::Cint)
    hand = ccall(:jl_uv_connect_handle, Ptr{Cvoid}, (Ptr{Cvoid},), conn)
    sock = @handle_as hand LibuvStream
    if status >= 0
        if !(sock.status == StatusClosed || sock.status == StatusClosing)
            sock.status = StatusOpen
        end
        notify(sock.connectnotify)
    else
        ccall(:jl_forceclose_uv, Cvoid, (Ptr{Cvoid},), hand)
        err = UVError("connect", status)
        notify_error(sock.connectnotify, err)
    end
    Libc.free(conn)
    nothing
end

function connect!(sock::TCPSocket, host::IPv4, port::Integer)
    if sock.status != StatusInit
        error("TCPSocket is not in initialization state")
    end
    if !(0 <= port <= typemax(UInt16))
        throw(ArgumentError("port out of range, must be 0 ≤ port ≤ 65535, got $port"))
    end
    uv_error("connect", ccall(:jl_tcp4_connect, Int32, (Ptr{Cvoid}, UInt32, UInt16, Ptr{Cvoid}),
                             sock.handle, hton(host.host), hton(UInt16(port)), uv_jl_connectcb::Ptr{Cvoid}))
    sock.status = StatusConnecting
    nothing
end

function connect!(sock::TCPSocket, host::IPv6, port::Integer)
    if sock.status != StatusInit
        error("TCPSocket is not in initialization state")
    end
    if !(0 <= port <= typemax(UInt16))
        throw(ArgumentError("port out of range, must be 0 ≤ port ≤ 65535, got $port"))
    end
    uv_error("connect", ccall(:jl_tcp6_connect, Int32, (Ptr{Cvoid}, Ref{UInt128}, UInt16, Ptr{Cvoid}),
                              sock.handle, hton(host.host), hton(UInt16(port)), uv_jl_connectcb::Ptr{Cvoid}))
    sock.status = StatusConnecting
    nothing
end

connect!(sock::TCPSocket, addr::InetAddr) = connect!(sock, addr.host, addr.port)

# Default Host to localhost

"""
    connect([host], port::Integer) -> TCPSocket

Connect to the host `host` on port `port`.
"""
connect(sock::TCPSocket, port::Integer) = connect(sock, localhost, port)
connect(port::Integer) = connect(localhost, port)

# Valid connect signatures for TCP
connect(host::AbstractString, port::Integer) = connect(TCPSocket(), host, port)
connect(addr::IPAddr, port::Integer) = connect(TCPSocket(), addr, port)
connect(addr::InetAddr) = connect(TCPSocket(), addr)

function connect!(sock::TCPSocket, host::AbstractString, port::Integer)
    if sock.status != StatusInit
        error("TCPSocket is not in initialization state")
    end
    ipaddr = getaddrinfo(host)
    sock.status = StatusInit
    connect!(sock,ipaddr,port)
    sock.status = StatusConnecting
    return sock
end

function connect(sock::LibuvStream, args...)
    connect!(sock, args...)
    wait_connected(sock)
    return sock
end

##

const BACKLOG_DEFAULT = 511

"""
    listen([addr, ]port::Integer; backlog::Integer=BACKLOG_DEFAULT) -> TCPServer

Listen on port on the address specified by `addr`.
By default this listens on `localhost` only.
To listen on all interfaces pass `IPv4(0)` or `IPv6(0)` as appropriate.
`backlog` determines how many connections can be pending (not having
called [`accept`](@ref)) before the server will begin to
reject them. The default value of `backlog` is 511.
"""
function listen(addr; backlog::Integer=BACKLOG_DEFAULT)
    sock = TCPServer()
    !bind(sock, addr) && error("cannot bind to port; may already be in use or access denied")
    listen(sock; backlog=backlog)
    return sock
end
listen(port::Integer; backlog::Integer=BACKLOG_DEFAULT) = listen(localhost, port; backlog=backlog)
listen(host::IPAddr, port::Integer; backlog::Integer=BACKLOG_DEFAULT) = listen(InetAddr(host, port); backlog=backlog)

function listen(callback, server::Union{TCPSocket, UDPSocket})
    @async begin
        local client = TCPSocket()
        while isopen(server)
            err = accept_nonblock(server, client)
            if err == 0
                callback(client)
                client = TCPSocket()
            elseif err != UV_EAGAIN
                uv_error("accept", err)
            else
                stream_wait(server, server.connectnotify)
            end
        end
    end
    return sock
end

function listen(sock::LibuvServer; backlog::Integer=BACKLOG_DEFAULT)
    uv_error("listen", trylisten(sock))
    return sock
end

# from `listen`
function uv_connectioncb(stream::Ptr{Cvoid}, status::Cint)
    sock = @handle_as stream LibuvServer
    if status >= 0
        notify(sock.connectnotify)
    else
        err = UVError("connection", status)
        notify_error(sock.connectnotify, err)
    end
    nothing
end

function trylisten(sock::LibuvServer; backlog::Integer=BACKLOG_DEFAULT)
    check_open(sock)
    err = ccall(:uv_listen, Cint, (Ptr{Cvoid}, Cint, Ptr{Cvoid}),
                sock, backlog, uv_jl_connectioncb::Ptr{Cvoid})
    sock.status = StatusActive
    return err
end

##

function accept_nonblock(server::TCPServer, client::TCPSocket)
    if client.status != StatusInit
        error("client TCPSocket is not in initialization state")
    end
    err = ccall(:uv_accept, Int32, (Ptr{Cvoid}, Ptr{Cvoid}), server.handle, client.handle)
    if err == 0
        client.status = StatusOpen
    end
    return err
end

function accept_nonblock(server::TCPServer)
    client = TCPSocket()
    uv_error("accept", accept_nonblock(server, client))
    return client
end

function accept(server::LibuvServer, client::LibuvStream)
    if server.status != StatusActive
        throw(ArgumentError("server not connected, make sure \"listen\" has been called"))
    end
    while isopen(server)
        err = accept_nonblock(server, client)
        if err == 0
            return client
        elseif err != UV_EAGAIN
            uv_error("accept", err)
        end
        stream_wait(server, server.connectnotify)
    end
    uv_error("accept", UV_ECONNABORTED)
end

## Utility functions

const localhost = ip"127.0.0.1"

"""
    listenany([host::IPAddr,] port_hint) -> (UInt16, TCPServer)

Create a `TCPServer` on any port, using hint as a starting point. Returns a tuple of the
actual port that the server was created on and the server itself.
"""
function listenany(host::IPAddr, default_port)
    addr = InetAddr(host, default_port)
    while true
        sock = TCPServer()
        if bind(sock, addr) && trylisten(sock) == 0
            if default_port == 0
                _addr, port = getsockname(sock)
                return (port, sock)
            end
            return (addr.port, sock)
        end
        close(sock)
        addr = InetAddr(addr.host, addr.port + 1)
        if addr.port == default_port
            error("no ports available")
        end
    end
end

listenany(default_port) = listenany(localhost, default_port)

"""
    getsockname(sock::Union{TCPServer, TCPSocket}) -> (IPAddr, UInt16)

Get the IP address and port that the given socket is bound to.
"""
getsockname(sock::Union{TCPSocket, TCPServer}) = _sockname(sock, true)

"""
    getpeername(sock::TCPSocket) -> (IPAddr, UInt16)

Get the IP address and port of the remote endpoint that the given
socket is connected to. Valid only for connected TCP sockets.
"""
getpeername(sock::TCPSocket) = _sockname(sock, false)

function _sockname(sock, self=true)
    rport = Ref{Cushort}(0)
    raddress = zeros(UInt8, 16)
    rfamily = Ref{Cuint}(0)

    if self
        r = ccall(:jl_tcp_getsockname, Int32,
                (Ptr{Cvoid}, Ref{Cushort}, Ptr{Cvoid}, Ref{Cuint}),
                sock.handle, rport, raddress, rfamily)
    else
        r = ccall(:jl_tcp_getpeername, Int32,
                (Ptr{Cvoid}, Ref{Cushort}, Ptr{Cvoid}, Ref{Cuint}),
                sock.handle, rport, raddress, rfamily)
    end
    uv_error("cannot obtain socket name", r)
    if r == 0
        port = ntoh(rport[])
        af_inet6 = @static if Sys.iswindows() # AF_INET6 in <sys/socket.h>
            23
        elseif Sys.isapple()
            30
        elseif Sys.KERNEL ∈ (:FreeBSD, :DragonFly)
            28
        elseif Sys.KERNEL ∈ (:NetBSD, :OpenBSD)
            24
        else
            10
        end

        if rfamily[] == 2 # AF_INET
            addrv4 = raddress[1:4]
            naddr = ntoh(unsafe_load(Ptr{Cuint}(pointer(addrv4)), 1))
            addr = IPv4(naddr)
        elseif rfamily[] == af_inet6
            naddr = ntoh(unsafe_load(Ptr{UInt128}(pointer(raddress)), 1))
            addr = IPv6(naddr)
        else
            error(string("unsupported address family: ", getindex(rfamily)))
        end
    else
        error("cannot obtain socket name")
    end
    return addr, port
end

# domain sockets

include("PipeServer.jl")

# libuv callback handles

function __init__()
    global uv_jl_getaddrinfocb = @cfunction(uv_getaddrinfocb, Cvoid, (Ptr{Cvoid}, Cint, Ptr{Cvoid}))
    global uv_jl_getnameinfocb = @cfunction(uv_getnameinfocb, Cvoid, (Ptr{Cvoid}, Cint, Cstring, Cstring))
    global uv_jl_recvcb        = @cfunction(uv_recvcb, Cvoid, (Ptr{Cvoid}, Cssize_t, Ptr{Cvoid}, Ptr{Cvoid}, Cuint))
    global uv_jl_sendcb        = @cfunction(uv_sendcb, Cvoid, (Ptr{Cvoid}, Cint))
    global uv_jl_connectioncb  = @cfunction(uv_connectioncb, Cvoid, (Ptr{Cvoid}, Cint))
    global uv_jl_connectcb     = @cfunction(uv_connectcb, Cvoid, (Ptr{Cvoid}, Cint))
end

# deprecations

@deprecate convert(dt::Type{<:Integer}, ip::IPAddr)  dt(ip)

@noinline function getaddrinfo(callback::Function, host::AbstractString)
    Base.depwarn("`getaddrinfo` with a callback function is deprecated, wrap code in `@async` instead for deferred execution.", :getaddrinfo)
    @async begin
        r = getaddrinfo(host)
        callback(r)
    end
    nothing
end

end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
"""
    converge_ensemble!(e; max_iterates, backup, clean, converge_criterion, converge_factor, mc_noise, progargs)

Given 'e <: GMC_NS_Ensemble', recursively execute nested sampling steps ('nested_step!()') until ensemble convergence criterion is obtained, or 'max_iterates' is reached.

Convergence criteria specied by 'converge_criterion <: String': "standard" for evidence convergence, "compression" for ensemble likelihood range compression (can be used to improve MAP parameter estimation after evidence has converged). Factor to use specified by 'converge_factor <: AbstractFloat'. "Standard" defined in Skilling 2006.

For likelihood estimates generated by Monte Carlo model evaluation, 'mc_noise <: AbstractFloat' defines the expected standard deviation of likelihood in the MAP region of parameter space; this will be used as an alternative convergence criterion if the ensemble is converged to within this value.

Keyword arguments to GMC_NS_Progress are passed from progargs.

Skilling, John. “Nested Sampling for Bayesian Computations.” In Proc. Valencia. Benidorm (Alicante, Spain): inference.org.uk, 2006. https://www.inference.org.uk/bayesys/Valencia.pdf.
"""
function converge_ensemble!(e::GMC_NS_Ensemble; max_iterates=typemax(Int64), backup::Tuple{Bool,Integer}=(false,0), clean::Tuple{Bool,Integer,Integer}=(false,0,0),  converge_criterion::String="standard", converge_factor::AbstractFloat=1e-3, mc_noise::AbstractFloat=0., progargs...)
    N = length(e.models); curr_it=length(e.log_Li)

    if curr_it==1 || !isfile(e.path*"/tuner")
        serialize(e.path*"/ens",e)
        tuner_dict=get_tuner_dict(e)
    else
        tuner_dict=deserialize(e.path*"/tuner") #restore tuner from saved if any
    end

    meter = GMC_NS_Progress(e, .1; start_it=curr_it, progargs...)

    converge_check = get_convfunc(converge_criterion)

    reflect_cache = nothing #init cache for models precalculated in galilean reflection

    cln_switch = clean[1] && !e.sample_posterior #ignore clean arguments if posterior samples are to be collected

    while !converge_check(e, converge_factor, mc_noise) && (curr_it <= max_iterates)
        warn, reflect_cache = nested_step!(e, tuner_dict, reflect_cache)
        warn == 1 && (@error "Failed to find new models, aborting at current iterate."; return e)
        curr_it += 1

        backup[1] && curr_it%backup[2] == 0 && e_backup(e,tuner_dict) #every backup interval, serialise the ensemble and tuner
        cln_switch && curr_it%clean[2] == 0 && clean_ensemble_dir(e,clean[3]) #every clean interval, remove old discarded models

        update!(meter, converge_check(e,converge_factor, vals=true)...)
    end

    if converge_check(e,converge_factor, mc_noise)
        final_logZ=measure_evidence(e)
        @info "Job done, sampled to convergence. Final logZ $(final_logZ.val) ± $(final_logZ.err)"

        e_backup(e,tuner_dict) #final backup
        cln_switch && clean_ensemble_dir(e,0) #final clean
        return final_logZ
    elseif curr_it==max_iterates
        @info "Job done, sampled to maximum iterate $max_iterates. Convergence criterion not obtained."

        e_backup(e,tuner_dict) #final backup
        cln_switch && clean_ensemble_dir(e,0) #final clean
        return e.log_Zi[end]
    end
end
                function evidence_converge(e, evidence_fraction, mc_noise=0.; vals=false)
                    mc_noise > 0. && noise_check(e,mc_noise) && return true

                    val=lps(findmax([model.log_Li for model in e.models])[1],  e.log_Xi[end])
                    thresh=lps(log(evidence_fraction),e.log_Zi[end])

                    vals ? (return val, thresh) : (return val<thresh)
                end

                function compress_converge(e, compression_ratio, mc_noise=0.; vals=false)
                    mc_noise > 0. && noise_check(e,mc_noise) && return true

                    val=findmax([model.log_Li for model in e.models])[1]-e.contour
                    thresh=compression_ratio

                    val===-Inf && (val=nextfloat(val))
                    vals ? (return val, thresh) : (return val<thresh)
                end

                function get_convfunc(criterion)
                    if criterion == "standard"
                        return evidence_converge
                    elseif criterion == "compression"
                        return compress_converge
                    else
                        throw(ArgumentError("Convergence criterion $criterion not supported! Try \"standard\" or \"compression\"."))
                    end
                end

                #if the entire ensemble is within 99.7% of the noise distribution of its mean likelihood, it is converged to the extent possible given Monte Carlo noise
                function noise_check(e, mc_noise)
                    mod_lhs=[m.log_Li for m in e.models]
                    std(mod_lhs) <= mc_noise ? (return true) : (return false)
                end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Scenario-decomposition generation #
# ================================= #
function generate!(stochasticprogram::StochasticProgram{N}, structure::ScenarioDecompositionStructure{N}) where N
    # Generate all stages
    for stage in 2:N
        generate!(stochasticprogram, structure, stage)
    end
    return nothing
end

function generate!(stochasticprogram::TwoStageStochasticProgram, structure::ScenarioDecompositionStructure{2}, stage::Integer)
    stage == 2 || error("Stage $stage not available in two-stage model.")
    # Check generators
    has_generator(stochasticprogram, :stage_1) || error("First-stage problem not defined in stochastic program. Consider @stage 1.")
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @stage 2.")
    # Sanity check on scenario probabilities
    if num_scenarios(structure, stage) > 0
        p = stage_probability(structure, stage)
        abs(p - 1.0) <= 1e-6 || @warn "Scenario probabilities do not add up to one. The probability sum is given by $p"
    end
    # Generate
    generate_horizontal!(scenarioproblems(structure, stage),
                         generator(stochasticprogram, :stage_1),
                         generator(stochasticprogram, :stage_2),
                         stage_parameters(stochasticprogram, 1),
                         stage_parameters(stochasticprogram, 2),
                         structure.decisions[stage-1],
                         subproblem_optimizer(stochasticprogram))
    return nothing
end

function generate_horizontal!(scenarioproblems::ScenarioProblems,
                              stage_one_generator::Function,
                              stage_two_generator::Function,
                              stage_one_params::Any,
                              stage_two_params::Any,
                              decision_map::DecisionMap,
                              optimizer)
    for i in num_subproblems(scenarioproblems)+1:num_scenarios(scenarioproblems)
        push!(scenarioproblems.problems, _WS(stage_one_generator,
                                             stage_two_generator,
                                             stage_one_params,
                                             stage_two_params,
                                             scenario(scenarioproblems,i),
                                             decision_map,
                                             DecisionMap(),
                                             optimizer))
    end
    return nothing
end
function generate_horizontal!(scenarioproblems::DistributedScenarioProblems,
                              stage_one_generator::Function,
                              stage_two_generator::Function,
                              stage_one_params::Any,
                              stage_two_params::Any,
                              ::DecisionMap,
                              optimizer)
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                stage_one_generator,
                stage_two_generator,
                stage_one_params,
                stage_two_params,
                scenarioproblems.decisions[w-1],
                optimizer) do sp, gen_one, gen_two, one_params, two_params, decisions, opt
                    generate_horizontal!(fetch(sp),
                                         gen_one,
                                         gen_two,
                                         one_params,
                                         two_params,
                                         fetch(decisions),
                                         opt)
                end
        end
    end
    return nothing
end

function clear!(structure::ScenarioDecompositionStructure{N}) where N
    # Clear decisions
    clear!(structure.decisions)
    # Clear all stages
    for stage in 2:N
        clear_stage!(structure, stage)
    end
    return nothing
end

function clear_stage!(structure::ScenarioDecompositionStructure{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    clear!(scenarioproblems(structure, s))
    return nothing
end