text:
# TODO
# we need to have nquat directly inside the model

################################################################################
# z indices
################################################################################

"""
	Returns the indices of q2 in z and Δq2 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_q2(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nquat = 0
	iq2 = Vector(1:nq - nquat)
	return iq2
end

"""
	Returns the indices of γ1 in z and Δγ1 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_γ1(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	off = (nq - nquat)
	iγ1 = Vector(off .+ (1:nc))
	return iγ1
end

"""
	Returns the indices of b1 in z and Δb1 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_b1(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc
	ib1 = Vector(off .+ (1:nb))
	return ib1
end

"""
	Returns the indices of ψ1 in z and Δψ1 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_ψ1(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb
	iψ1 = Vector(off .+ (1:nc))
	return iψ1
end

"""
	Returns the indices of s1 in z and Δs1 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_s1(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc
	is1 = Vector(off .+ (1:nc))
	return is1
end

"""
	Returns the indices of η1 in z and Δη1 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_η1(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc + nc
	iη1 = Vector(off .+ (1:nb))
	return iη1
end

"""
	Returns the indices of s2 in z and Δs2 in Δz. Important to note that for the indexing
	in z, we want to set quat = false. For indexing in Δz, we set quat = true, so that nquat =
	number of quaternions.
"""
function index_s2(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc + nc + nb
	is2 = Vector(off .+ (1:nc))
	return is2
end

################################################################################
# θ indices
################################################################################

"""
	Returns the indices of q0 in θ.
"""
function index_q0(model::Model)
	nq = model.nq
	iq0 = Vector(1:nq)
	return iq0
end

"""
	Returns the indices of q1 in θ.
"""
function index_q1(model::Model)
	nq = model.nq
	off = nq
	iq1 = Vector(off .+ (1:nq))
	return iq1
end

"""
	Returns the indices of u1 in θ.
"""
function index_u1(model::Model)
	nq = model.nq
	nu = model.nu
	off = nq + nq
	iu1 = Vector(off .+ (1:nu))
	return iu1
end

"""
	Returns the indices of w1 in θ.
"""
function index_w1(model::Model)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	off = nq + nq + nu
	iw1 = Vector(off .+ (1:nw))
	return iw1
end

"""
	Returns the indices of μ in θ.
"""
function index_μ(model::Model)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	off = nq + nq + nu + nw
	iμ = Vector(off .+ (1:1))
	return iμ
end

"""
	Returns the indices of h in θ.
"""
function index_h(model::Model)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	off = nq + nq + nu + nw + 1
	ih = Vector(off .+ (1:1))
	return ih
end

################################################################################
# Residual indices
################################################################################

"""
	Returns the indices of the dynamics equation in the residual r.
"""
function index_dyn(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nquat = 0
	idyn = Vector(1:nq - nquat)
	return idyn
end

"""
	Returns the indices of the impact equation in the residual r.
"""
function index_imp(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	off = (nq - nquat)
	iimp = Vector(off .+ (1:nc))
	return iimp
end

"""
	Returns the indices of the maximum dissipation principle equation in the residual r.
"""
function index_mdp(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc
	imdp = Vector(off .+ (1:nb))
	return imdp
end

"""
	Returns the indices of the friction cone equation in the residual r.
"""
function index_fri(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb
	ifri = Vector(off .+ (1:nc))
	return ifri
end

"""
	Returns the indices of the bilinear impact equation in the residual r.
"""
function index_bimp(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc
	ibimp = Vector(off .+ (1:nc))
	return ibimp
end

"""
	Returns the indices of the bilinear maximum dissipation principle equation in the residual r.
"""
function index_bmdp(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc + nc
	ibmdp = Vector(off .+ (1:nb))
	return ibmdp
end

"""
	Returns the indices of the bilinear friction equation in the residual r.
"""
function index_bfri(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	off = (nq - nquat) + nc + nb + nc + nc + nb
	ibfri = Vector(off .+ (1:nc))
	return ibfri
end

################################################################################
# Aggregated Indices
################################################################################

"""
	Returns the indices of z and Δz, where 3 groups are formed.

	w1 = (q2)
	w2 = (γ1, b1, ψ1)
	w3 = (s1, η1, s2)

	Δw1 = (Δq2)
	Δw2 = (Δγ1, Δb1, Δψ1)
	Δw3 = (Δs1, Δη1, Δs2)

	Note that Δq2 does not have the same size as q2 when optimizing over
	non-euclidean spaces.
"""
function linearization_var_index(model::Model, env::Environment; quat::Bool = false)
	iq2 = index_q2(model, env, quat = quat)
	iγ1 = index_γ1(model, env, quat = quat)
	ib1 = index_b1(model, env, quat = quat)
	iψ1 = index_ψ1(model, env, quat = quat)
	is1 = index_s1(model, env, quat = quat)
	iη1 = index_η1(model, env, quat = quat)
	is2 = index_s2(model, env, quat = quat)
	iw1 = iq2
	iw2 = [iγ1; ib1; iψ1]
	iw3 = [is1; iη1; is2]
	return iw1, iw2, iw3
end

"""
	Returns the indices of the residual r, where 3 groups are formed.
	r = (rdyn, rrst, rbil)
"""
function linearization_term_index(model::Model, env::Environment; quat::Bool = false)
	# dyn = [dyn]
	# rst = [s1  - ..., ≡ ialt
	#        η1  - ...,
	#        s2  - ...,]
	# bil = [γ1 .* s1 .- κ;
	#        b1 .* η1 .- κ;
	#        ψ1 .* s2 .- κ]
	idyn = index_dyn(model, env, quat = quat)
	iimp = index_imp(model, env, quat = quat)
	imdp = index_mdp(model, env, quat = quat)
	ifri = index_fri(model, env, quat = quat)
	ibimp = index_bimp(model, env, quat = quat)
	ibmdp = index_bmdp(model, env, quat = quat)
	ibfri = index_bfri(model, env, quat = quat)

	irst = [iimp; imdp; ifri]
	ibil = [ibimp; ibmdp; ibfri]
	ialt = iimp
	return idyn, irst, ibil, ialt
end

"""
	Returns the positive orthant indices in z or Δz.
"""
function index_ort(model::Model, env::Environment{<:World,LinearizedCone}; quat::Bool = false)
	iγ1 = index_γ1(model, env, quat = quat)
	ib1 = index_b1(model, env, quat = quat)
	iψ1 = index_ψ1(model, env, quat = quat)
	is1 = index_s1(model, env, quat = quat)
	iη1 = index_η1(model, env, quat = quat)
	is2 = index_s2(model, env, quat = quat)
	return [[iγ1; ib1; iψ1], [is1; iη1; is2]]
end

"""
	Returns the positive orthant indices in z or Δz.
"""
function index_ort(model::Model, env::Environment{<:World,NonlinearCone}; quat::Bool = false)
	iγ1 = index_γ1(model, env, quat = quat)
	is1 = index_s1(model, env, quat = quat)
	return [iγ1, is1]
end

"""
	Returns the second order cone indices in z or Δz.
"""
index_soc(model::Model, env::Environment{<:World,LinearizedCone}; quat::Bool = false) = Vector{Vector{Int}}[]

"""
	Returns the second order cone indices in z or Δz.
"""
function index_soc(model::Model, env::Environment{<:World,NonlinearCone}; quat::Bool = false)
	nc = model.nc
	nf = friction_dim(env)

	ib1 = index_b1(model, env, quat = quat) # primal cones: vector part
	iψ1 = index_ψ1(model, env, quat = quat) # primal cones: scalar part
	iη1 = index_η1(model, env, quat = quat) # dual cones: vector part
	is2 = index_s2(model, env, quat = quat) # dual cones: scalar part

	[[[iψ1[i]; iη1[(i - 1) * nf .+ (1:nf)]], [is2[i]; ib1[(i - 1) * nf .+ (1:nf)]]] for i = 1:nc]
end

function num_var(model::Model, env::Environment; quat::Bool = false)
	nq = model.nq
	nc = model.nc
	nquat = 0
	nb = nc * friction_dim(env)
	(nq - nquat) + nc + nb + nc + nc + nb + nc
end

function num_data(model::Model)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	nq + nq + nu + nw + 1 + 1
end

function num_bilinear(model::Model, env::Environment)
	nc = model.nc
	nb = nc * friction_dim(env)
	nc + nb + nc
end

################################################################################
# Packing and Unpacking Variables
################################################################################

function unpack_θ(model::Model, θ)
	iq0 = index_q0(model)
	iq1 = index_q1(model)
	iu1 = index_u1(model)
	iw1 = index_w1(model)
	iμ  = index_μ(model)
	ih  = index_h(model)

	q0 = θ[iq0]
	q1 = θ[iq1]
	u1 = θ[iu1]
	w1 = θ[iw1]
	μ = θ[iμ]
	h = θ[ih]
	return q0, q1, u1, w1, μ, h
end

function pack_θ(model::Model, q0, q1, u1, w1, μ, h)
	return [q0; q1; u1; w1; μ; h]
end

function unpack_z(model::Model, env::Environment, z)
	iq2 = index_q2(model, env, quat = false)
	iγ1 = index_γ1(model, env, quat = false)
	ib1 = index_b1(model, env, quat = false)
	iψ1 = index_ψ1(model, env, quat = false)
	is1 = index_s1(model, env, quat = false)
	iη1 = index_η1(model, env, quat = false)
	is2 = index_s2(model, env, quat = false)

	# system variables
	q2 = z[iq2]
	γ1 = z[iγ1]
	b1 = z[ib1]
	ψ1 = z[iψ1]
	s1 = z[is1]
	η1 = z[iη1]
	s2 = z[is2]
	return q2, γ1, b1, ψ1, s1, η1, s2
end

function pack_z(model::Model, env::Environment{<:World,LinearizedCone}, q2, γ1, b1, ψ1, η1)
	s1 = ϕ_func(model, env, q2)
	s2 = model.μ_world * γ1 .- E_func(model, env) * b1
	return pack_z(model, env, q2, γ1, b1, ψ1, s1, η1, s2)
end

function pack_z(model::Model, env::Environment{<:World,NonlinearCone}, q2, γ1, b1, ψ1, η1)
	s1 = ϕ_func(model, env, q2)
	s2 = model.μ_world .* γ1
	return pack_z(model, env, q2, γ1, b1, ψ1, s1, η1, s2)
end

function pack_z(model::Model, env::Environment, q2, γ1, b1, ψ1, s1, η1, s2)
	return [q2; γ1; b1; ψ1; s1; η1; s2]
end

################################################################################
# Aggregated indices
################################################################################

function index_equr(model::Model, env::Environment; quat::Bool = false)
	idyn = index_dyn(model, env, quat = quat)
	iimp = index_imp(model, env, quat = quat)
	imdp = index_mdp(model, env, quat = quat)
	ifri = index_fri(model, env, quat = quat)
	iequ = [idyn; iimp; imdp; ifri]
	return iequ
end

function index_ortr(model::Model, env::Environment{<:World,LinearizedCone}; quat::Bool = false)
	ibimp = index_bimp(model, env, quat = quat)
	ibmdp = index_bmdp(model, env, quat = quat)
	ibfri = index_bfri(model, env, quat = quat)
	iortr = [ibimp; ibmdp; ibfri]
	return iortr
end

function index_ortr(model::Model, env::Environment{<:World,NonlinearCone}; quat::Bool = false)
	ibimp = index_bimp(model, env, quat = quat)
	iortr = ibimp
	return iortr
end

function index_socr(model::Model, env::Environment{<:World,LinearizedCone}; quat::Bool = false)
	isocr = Vector{Int64}()
	return isocr
end

function index_socr(model::Model, env::Environment{<:World,NonlinearCone}; quat::Bool = false)
	ibmdp = index_bmdp(model, env, quat = quat)
	ibfri = index_bfri(model, env, quat = quat)
	isocr = [ibmdp; ibfri]
	return isocr
end

function index_socri(model::Model, env::Environment{<:World,LinearizedCone}; quat::Bool = false)
	isocri = Vector{Vector{Int}}()
	return isocri
end

function index_socri(model::Model, env::Environment{<:World,NonlinearCone}; quat::Bool = false)
	isocri = Vector{Vector{Int}}()
	return isocri
end

# struct DynamicsStructure <: IndicesStructure
# 	ny::Int
# 	dynz::Vector{Int}
# 	dynΔ::Vector{Int}
# 	dyn::Vector{Int}
# 	rst::Vector{Int}
# 	alt::Vector{Int}
# end

function IndicesOptimization(model::Model, env::Environment)
	# Dimensions
	nz = num_var(model, env, quat = false)
	nΔ = num_var(model, env, quat = true)
	ny = num_bilinear(model, env)

	# Variables
	dynz = index_q2(model, env, quat = false)
	dynΔ = index_q2(model, env, quat = true)
	ortz = index_ort(model, env, quat = true)
	ortΔ = index_ort(model, env, quat = false)
	socz = index_soc(model, env, quat = true)
	socΔ = index_soc(model, env, quat = false)

	# Residual
	# dyn = [dyn]
	# rst = [s1  - ..., ≡ ialt
	#        η1  - ...,
	#        s2  - ...,]
	# bil = [γ1 .* s1 .- κ;
	#        b1 .* η1 .- κ;
	#        ψ1 .* s2 .- κ]

	equr = index_equr(model, env, quat = true)
	ortr = index_ortr(model, env, quat = true)
	socr = index_socr(model, env, quat = true)
	socri = index_socri(model, env, quat = true)

	dyn = index_dyn(model, env, quat = true)
	imp = index_imp(model, env, quat = true)
	mdp = index_mdp(model, env, quat = true)
	fri = index_fri(model, env, quat = true)
	bimp = index_bimp(model, env, quat = true)
	bmdp = index_bmdp(model, env, quat = true)
	bfri = index_bfri(model, env, quat = true)

	rst = [imp; mdp; fri]
	bil = [bimp; bmdp; bfri]

	IndicesOptimization(
		nz, nΔ,
		ortz, ortΔ, socz, socΔ,
		equr, ortr, socr, socri, bil)
end


text:
include("../searchPatches/searchPatches.jl")

#=
A pass that searches.

This understands that a pass is an iteration over the target values (pixels).
The result of a pass is side effects on the array (image) and other data.
The only non side effect is the returned value.

Return the count of better points that the pass found.

Original name in C: synthesize()
=#

function makeAPass(targetImage, corpusImage, targetPoints, synthPatch, synthResult, searchResult, sortedOffsets)::Int64
    #=
    Retry a target point even if a perfect match was found on prior pass
    since its patch might have changed.
    That is, we don't retain patch diffs between passes.
    All we keep is the best matching point.
    =#
    countBetterPixels::Int64 = 0

    for targetPoint in targetPoints

        synthPatchCenterPoint = targetPoint

        #=
        The scatterPatch was passed in, so it was allocated once.)
        Mutate it, instead of allocating a new one.
        =#
        patchSize = prepareScatterPatch(
            synthPatch,
            targetImage,           # frame
            synthPatchCenterPoint, # framed point
            synthResult,
            sortedOffsets)
        #=
        Assert for certain edge-case inputs, the synthPatch is not full size.
        TODO pass patchSize to the search, and iterate over only the valid patch.
        =#

        #=
        Set searchResult to the starting value of (NotBetter, Inf)
        We can't use diff for this targetpoint/patch from prior passes,
        since, as discussed below, the patch changes.
        =#
        setStartSearchResult(searchResult, ProbeResult())

        # !!! side effects on targetImage, and synthResult
        searchResult = searchForPatchMatches(
            synthPatch,
            targetImage, synthPatchCenterPoint, # framed point
            corpusImage,
            synthResult,
            searchResult)

        if searchResult.bestProbeResult.betterment != NotBetter
            # TODO this is too crude, for equal to the old result ??
            # See original algorithm
            countBetterPixels += 1

            #@debug "Corpus and target points"
            #@debug searchResult.bestMatchPointInCorpus
            #@debug targetPoint

            #=
            TODO When parameters.withReplacement.
            Side effect on targetImage.
            Change the value at the target point for which we found a better match.
            =#
            targetImage.image[targetPoint] = corpusImage.image[searchResult.bestMatchPointInCorpus]

            #=
            All we save from the searchResult is bestMatchPointInCorpus.
            !!! Not save the diff, since patch may change.
            That is, the existing diff for this point was for a patch that is now different
            because the neighbors changed.

            !!! This also records the boolean event that we synthesized targetPoint,
            so we now know its value (e.g. color or best matching color)
            even if we did not actually replace the value in the target.
            =#
            setBetterSearchResult(synthResult, targetPoint, searchResult.bestMatchPointInCorpus)
        end
    end
    return countBetterPixels
end

ext: jl

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module Conduit

import ..API

struct ConduitNode
    ptr::Ptr{API.conduit_node}

    function ConduitNode()
	    nptr = API.conduit_node_create()
	    return new(nptr)
    end
end
Base.unsafe_convert(::Type{Ptr{API.conduit_node}}, node::ConduitNode) = node.ptr

function ConduitNode(f::Function)
    node = ConduitNode()
    try
	    f(node)
    finally
	    node_destroy!(node)
    end
    return nothing
end

Base.setindex!(node::ConduitNode, val, path::String) = node_set!(node, path, val)

for numtype in (:UInt8, :Int8, :UInt16, :Int16, :UInt32, :Int32, :UInt64, :Int64, :Float32, :Float64)
    cnumtype = Symbol(lowercase("$(numtype)"))
    node_set_path = Symbol("conduit_node_set_path_$(cnumtype)")
    node_set_path_ptr = Symbol("conduit_node_set_path_$(cnumtype)_ptr")
    @eval begin
        function node_set!(node::ConduitNode, path::String, val::$numtype)
            API.$node_set_path(node, path, val)
            return node
        end
        function node_set!(node::ConduitNode, path::String, val::Array{$numtype})
            API.$node_set_path_ptr(node, path, pointer(val), length(val))
            return node
        end
    end
end

function node_set!(node::ConduitNode, path::String, val::String)
    API.conduit_node_set_path_char8_str(node, path, val)
    return node
end

function node_set!(node::ConduitNode, path::String, val::ConduitNode)
    API.conduit_node_set_path_node(node, path, val)
    return node
end

function node_has_path(node::ConduitNode, path::String)
    return Bool(API.conduit_node_has_path(node, path))
end

function node_path(node::ConduitNode)
    return unsafe_string(API.conduit_node_path(node))
end

function node_remove_path!(node::ConduitNode, path::String)
    API.conduit_node_remove_path(node, path)
    return node
end

function node_update!(dest_node::ConduitNode, src_node::ConduitNode)
    API.conduit_node_update(dest_node, src_node)
    return dest_node
end

function node_update_compatible!(dest_node::ConduitNode, src_node::ConduitNode)
    API.conduit_node_update_compatible(dest_node, src_node)
    return dest_node
end

function node_update_external!(dest_node::ConduitNode, src_node::ConduitNode)
    API.conduit_node_update_external(dest_node, src_node)
    return dest_node
end

function node_info!(dest_node::ConduitNode, src_node::ConduitNode)
    API.conduit_node_info(src_node, dest_node)
    return dest_node
end
function node_info(f::Function, node::ConduitNode)
    info = ConduitNode()
    try
        node_info!(info, node)
        f(info)
    finally
	    node_destroy!(info)
    end
    return
end

function node_destroy!(node::ConduitNode)
    API.conduit_node_destroy(node)
    return nothing
end

function node_print(node::ConduitNode; detailed = false)
    if detailed
	    API.conduit_node_print_detailed(node)
    else
	    API.conduit_node_print(node)
    end
    return
end

function node_print_string(node::ConduitNode; detailed = false)
    output_buffer = IOBuffer()
    default_stdout = stdout

    pipe = Pipe()
    Base.link_pipe!(pipe; reader_supports_async = true, writer_supports_async = true)
    pe_stdout = pipe.in
    redirect_stdout(pe_stdout)

    buffer_redirect_task = @async write(output_buffer, pipe)
    try
        node_print(node; detailed);
    finally
        redirect_stdout(default_stdout)
        close(pe_stdout)
        wait(buffer_redirect_task)
    end
    return String(take!(output_buffer))
end

function about()
    ConduitNode() do node
        API.conduit_about(node)
	    node_print(node)
    end
end

function example_basic_mesh(;mesh = "hexs", nx = 10, ny = 10, nz = 10)
    node = ConduitNode()
    GC.@preserve node begin
        API.conduit_blueprint_mesh_examples_basic(mesh, nx, ny, nz, node)
    end
    return node
end

function example_braid_mesh(;mesh = "hexs", nx = 50, ny = 50, nz = 50)
    node = ConduitNode()
    GC.@preserve node begin
        API.conduit_blueprint_mesh_examples_braid(mesh, nx, ny, nz, node)
    end
    return node
end

for example in (:example_basic_mesh, :example_braid_mesh)
    @eval begin
        function ($example)(f::Function; kw...)
            node = ($example)(; kw...)
            try
                f(node)
            finally
                node_destroy!(node)
            end
        end
    end
end

end


text:
module SummitBackgroundManager

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Summit Background Manager" => Ahorn.EntityPlacement(
        Maple.SummitBackgroundManager
    )
)

function Ahorn.selection(entity::Maple.SummitBackgroundManager)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle[Ahorn.Rectangle(x - 12, y - 12, 24, 24)]
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SummitBackgroundManager, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.summitBackgroundManager, -12, -12)

end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# This file was generated, do not modify it. # hide
ẑ = (W, z) |> RidgeRegressor(lambda=0.1)

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using Clustering
using LinearAlgebra
using Plots
using Plots.PlotMeasures

mutable struct Posterior_similarity_matrix
    psm::Vector
    names::Vector{String}
    Posterior_similarity_matrix(K, n_obs) = new([Matrix(1.0I, n_obs, n_obs) for k in 1:(K + (K > 1))],
                                                Vector{String}(undef, K + (K > 1)))
end

"""
`generate_psm(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)`

Generates a posterior similarity matrix from pmdi output
## Input
- `outputFile::S    tring` a string referring to the location on disk of pmdi output
- `burnin::Int64` an integer of the number of initial iterations to discard as burn-in
- `thin::Int64` an integer for the rate at which to thin output, `thin = 2` discards
every second iteration

## Output
A struct of class `Posterior_similarity_matrix` containing:
    - `psm` K n × n matrices where element (i, j) measures how frequently observation
    i and j are co-clustered in dataset k. If `K > 1` an overall consensus matrix
    is also included which is an element-wise average of the K other matrices
    - `names` a vector of the dataset names
"""
function generate_psm(outputFile::String, burnin::Int64 = 0, thin::Int64 = 1)
    outputNames = split(readline(outputFile), ',')
    output = readdlm(outputFile, ',', header = false, skipstart = burnin + 1)
    K = sum(map(outputNames) do str
        occursin(r"MassParameter", str)
    end)

    output = output[1:thin:end, (K + binomial(K, 2) + (K == 1) + 2):end]

    n_obs = size(output, 2) / K
    @assert mod(n_obs, 1) == 0 "Error: Datasets have different number of observations"
    n_obs = Int64(n_obs)
    n_iter = size(output, 1)

    psm = Posterior_similarity_matrix(K, n_obs)
    psm.names[1:K] = unique(map(x -> split(x, '_')[1], outputNames[(K + (K == 1) + binomial(K, 2) + 2):end]))
    if K > 1
        psm.names[K + 1] = "Overall"
    end
    for k = 1:K
        for j = 1:(n_obs - 1)
            for i = (j + 1):n_obs
                psm.psm[k][i, j] = sum(output[:, i + n_obs * (k - 1)] .== output[:, j + n_obs * (k - 1)]) / n_iter
            end
        end
    end
    if K > 1
        for k = 1:K
            psm.psm[K + 1] += psm.psm[k] / K
        end
            psm.psm[K + 1][diagind(psm.psm[K + 1])] .= 1.0
    end

    return psm
end

function generate_psm2(output, K, n_obs)
    output = output[(K + binomial(K, 2) + (K == 1) + 2):end]

    n_iter = size(output, 1)

    psm = Posterior_similarity_matrix(K, n_obs)
    for k = 1:K
        for j = 1:(n_obs - 1)
            for i = (j + 1):n_obs
                psm.psm[k][i, j] = sum(output[i + n_obs * (k - 1)] .== output[j + n_obs * (k - 1)]) / n_iter
            end
        end
    end
    if K > 1
        for k = 1:K
            psm.psm[K + 1] += psm.psm[k] / K
        end
            psm.psm[K + 1][diagind(psm.psm[K + 1])] .= 1.0
    end
    return psm
end

"""
`get_consensus_allocations(psm::Posterior_similarity_matrix, nclust::Int64, orderby::Int64 = 0)`
"""
function get_consensus_allocations(psm::Posterior_similarity_matrix; k::Union{Int64, Nothing}=nothing, h::Union{Float64, Nothing} = nothing,
    linkage = :ward, orderby::Int64 = 0)
    @assert (k != nothing) | (h != nothing) "You must specify either k (number of clusters) or h (height to cut dendrogram)"
    if orderby == 0
        orderby = size(psm.psm, 1)
    end
    hc = hclust(1 .- Symmetric(psm.psm[orderby], :L), linkage = linkage)
    if k != nothing
        return cutree(hc, k = k)
    elseif h != nothing
        return cutree(hc, h = h)
    end

end

"""
`consensus_map(psm::Posterior_similarity_matrix; k::Union{Int64, Nothing} = nothing,
h::Union{Float64, Nothing} = nothing,
orderby::Int64 = 0)`

Generates a consensus map plot for a given posterior similarity matrix
## Input
- `psm` a posterior similarity matrix struct as output from `generate_psm`
- `k` optional - the number of clusters to highlight in the data
- `h` optional - the distance at which to cut the dendrogram. One of `k` or `h` must be specified.
- `orderby` observations in all plots are ordered in a way to separate clusters as well
as possible. `orderby` specifies which dataset should be used to inform this ordering.
`orderby = 0` will let the overall consensus dictate the ordering.
## Output
- K + 1 (or K if K == 1) consensus maps illustrating the clustering output in each
of K datasets.
"""
function consensus_map(psm::Posterior_similarity_matrix;
                       k::Union{Int64, Nothing} = nothing,
                       h::Union{Float64, Nothing} = nothing,
                       orderby::Int64 = 0,
                       linkage = :ward)
    @assert (k != nothing) | (h != nothing) "You must specify either k (number of clusters) or h (height to cut dendrogram)"
    if (orderby == 0) | (orderby == - 1)
        orderby_main = size(psm.psm, 1)
    else orderby_main = orderby
    end
    hc = hclust(1 .- Symmetric(psm.psm[orderby_main], :L), linkage = linkage)
    if h == nothing
        cuts = cutree(hc, k = k)[hc.order]
    elseif k == nothing
        cuts = cutree(hc, h = h)[hc.order]
    end
    nclust = length(unique(cuts))
    ticks = indexin(1:nclust, cuts) .- 0.5
    sort!(ticks)
    append!(ticks, size(psm.psm[1], 1) + 0.5)
    # Surely a better way of doing this
    ticks1 = [[ticks[i], ticks[i+1], NaN][j] for i = 1:nclust for j = 1:3]
    ticks2 = [[ticks[i], ticks[i], NaN][j] for i = 1:nclust for j = 1:3]
    ticks3 = [[ticks[i], ticks[i], NaN][j] for i = 2:(nclust + 1) for j = 1:3]
    K = length(psm.psm)
    if (K > 1) & (orderby != - 1)
        plots = [Plots.heatmap(Symmetric(psm.psm[k], :L)[hc.order, hc.order],
                        ticks = false,
                        yflip = false,
                        title = psm.names[k],
                        legend = false,
                        aspect_ratio = 1,
                        widen = false,
                        c = :viridis,
                        clim = (0, 1)) for k in [K; 1:max(K - 1, 1)]]
        l = @layout [a{0.85w} grid(K - 1, 1)]
    elseif (K > 1) & (orderby == - 1)
        plot_order = [hclust(1 .- Symmetric(psm.psm[k], :L), linkage = :ward).order for k in 1:K]
        plots = [Plots.heatmap(Symmetric(psm.psm[k], :L)[plot_order[k], plot_order[k]],
                        ticks = false,
                        yflip = false,
                        title = psm.names[k],
                        legend = false,
                        aspect_ratio = 1,
                        c = :viridis,
                        clim = (0, 1),
                        widen = false) for k in [K; 1:max(K - 1, 1)]]
        l = @layout [a{0.85w} grid(max(K - 1, 1), 1)]
    elseif K == 1
        plots = [Plots.heatmap(Symmetric(psm.psm[1], :L)[hc.order, hc.order],
                        ticks = false,
                        yflip = false,
                        title = psm.names[1],
                        legend = false,
                        aspect_ratio = 1,
                        widen = false,
                        # titlefont = Plots.font(family = "serif", pointsize = 12),
                        c = :viridis
                        )]
        l = @layout [a;]
    end

    Plots.plot(plots..., layout = l,
               framestyle = :none)
    plot!(ticks1, ticks2, seriestype = :path, c = "#FFFFFF", linestyle = :dash, linewidth = 2.5)
    plot!(ticks1, ticks3, seriestype = :path, c = "#FFFFFF", linestyle = :dash, linewidth = 2.5)
    plot!(ticks2, ticks1, seriestype = :path, c = "#FFFFFF", linestyle = :dash, linewidth = 2.5)
    plot!(ticks3, ticks1, seriestype = :path, c = "#FFFFFF", linestyle = :dash, linewidth = 2.5, widen = false)

end


text:
# julia> ENV["LLVM_CONFIG"]="/usr/bin/llvm-config"
# julia> Pkg.build("Clang")
# julia> using Clang
# julia> include("gen_libcurand.jl")

using Clang

includes = ["/usr/local/include",
            "/usr/include",
            "/home/shindo/local-lemon/cuda-7.5/include",
            "/usr/lib/gcc/x86_64-linux-gnu/4.8/include",
            "/usr/lib/gcc/x86_64-linux-gnu/4.8/include-fixed"]
headers = ["cuda_runtime_api.h", "driver_types.h", "vector_types.h"]
headers = [joinpath("/home/shindo/cuda/include/",h) for h in headers]

# Customize how functions, constants, and structs are written
const skip_expr = [:(const CUDART_DEVICE = __device__)]
const skip_error_check = [:cudaStreamQuery,:cudaGetLastError,:cudaPeekAtLastError]
function rewriter(ex::Expr)
    if in(ex, skip_expr)
        return :()
    end
    # Empty types get converted to Void
    if ex.head == :type
        a3 = ex.args[3]
        if isempty(a3.args)
            objname = ex.args[2]
            return :(typealias $objname Void)
        end
    end
    ex.head == :function || return ex
    decl, body = ex.args[1], ex.args[2]
    # omit types from function prototypes
    for i = 2:length(decl.args)
        a = decl.args[i]
        if a.head == :(::)
            decl.args[i] = a.args[1]
        end
    end
    # Error-check functions that return a cudaError_t (with some omissions)
    ccallexpr = body.args[1]
    if ccallexpr.head != :ccall
        error("Unexpected body expression: ", body)
    end
    rettype = ccallexpr.args[2]
    if rettype == :cudaError_t
        fname = decl.args[1]
        if !in(fname, skip_error_check)
            body.args[1] = Expr(:call, :checkerror, deepcopy(ccallexpr))
        end
    end
    ex
end

rewriter(A::Array) = [rewriter(a) for a in A]
rewriter(s::Symbol) = string(s)
rewriter(arg) = arg

context=wrap_c.init(output_file = "libcudart.jl",
                    common_file = "libcudart_types.jl",
                    header_library = x -> "libcudart",
                    headers = headers,
                    clang_includes = includes,
                    clang_diagnostics = true,
                    # header_wrapped=(x,y)->contains(x,"cuda"),
                    header_wrapped = (x,y) -> true,
                    rewriter=rewriter)

context.options = wrap_c.InternalOptions(true,true)  # wrap structs, too

run(context)


text:
immutable RectangularBinGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement

    function RectangularBinGeometry(
            default_statistic::Gadfly.StatisticElement=Gadfly.Stat.identity())
        new(default_statistic)
    end
end

const rectbin = RectangularBinGeometry

function histogram2d(; xbincount=nothing, xminbincount=3, xmaxbincount=150,
                       ybincount=nothing, yminbincount=3, ymaxbincount=150)
    RectangularBinGeometry(
        Gadfly.Stat.histogram2d(xbincount=xbincount,
                                xminbincount=xminbincount,
                                xmaxbincount=xmaxbincount,
                                ybincount=ybincount,
                                yminbincount=yminbincount,
                                ymaxbincount=ymaxbincount))
end

function default_statistic(geom::RectangularBinGeometry)
    geom.default_statistic
end

function element_aesthetics(::RectangularBinGeometry)
    [:x, :y, :xmin, :xmax, :ymin, :ymax, :color]
end

# Render rectangular bin (e.g., heatmap) geometry.
#
# Args:
#   geom: rectbin geometry
#   theme: the plot's theme
#   aes: some aesthetics
#
# Returns
#   A compose form.
#
function render(geom::RectangularBinGeometry,
                theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    if aes.x === nothing && (aes.xmin === nothing || aes.xmax === nothing)
        error("Geom.rectbin requires either x or both xmin and xmax be defined.")
    end

    if aes.y === nothing && (aes.ymin === nothing || aes.ymax === nothing)
        error("Geom.rectbin requires either y or both ymin and ymax be defined.")
    end

    if aes.xmin != nothing && length(aes.xmin) != length(aes.xmax)
        error("Geom.rectbin requires that xmin and xmax be of equal length.")
    end

    if aes.ymin != nothing && length(aes.ymin) != length(aes.ymax)
        error("Geom.rectbin requires that ymin and ymax be of equal length.")
    end

    nx = aes.xmin === nothing ? length(aes.x) : length(aes.xmin)
    ny = aes.ymin === nothing ? length(aes.y) : length(aes.ymin)

    if nx != ny
        error("""Geom.rectbin requires an equal number of x (or xmin/xmax) and
                 y (or ymin/ymax) values.""")
    end

    if aes.xmin === nothing
        xmin = [x - 0.5 for x in aes.x]
        xmax = [x + 0.5 for x in aes.x]
    else
        xmin = aes.xmin
        xmax = aes.xmax
    end

    if aes.ymin === nothing
        ymin = [y - 0.5 for y in aes.y]
        ymax = [y + 0.5 for y in aes.y]
    else
        ymin = aes.ymin
        ymax = aes.ymax
    end

    n = nx
    forms = Array(Compose.Form, 0)
    for (i, c) in zip(1:n, cycle(aes.color))
        if !isna(c)
            form = compose(rectangle(xmin[i], ymin[i],
                                    (xmax[i] - xmin[i])*cx - theme.bar_spacing,
                                    (ymax[i] - ymin[i])*cy - theme.bar_spacing),
                           fill(c), svgclass("geometry"))

            push!(forms, form)
        end
    end

    compose(combine(forms...),
            stroke(nothing),
            svgattribute("shape-rendering", "crispEdges"))
end


text:
using CSV
using DataFrames
using Plots

data = CSV.read("memory_usage.csv", DataFrame)
println(data[1:5, :])

timestamp = data[:, 1]
# allocated = 4 * data[:, 2] ./ (1024^3)
allocated = 4 * data[:, 2] ./ (1024^2) # kB for rss, bytes otherwise (Rust and FORTRAN), 4 Julia workers

# plot(timestamp, allocated, label = "jemalloc stats.allocated memory [GB]", xlabel = "elapsed time [s]", ylabel = "memory [GB]", legend = :bottomright, title = "memory consumption")
plot(timestamp, allocated, label = "allocated memory [GB]", xlabel = "elapsed time [s]", ylabel = "memory [GB]", legend = :bottomright, title = "Julia memory consumption")

savefig(homedir() * "/mem_allocated_v5.pdf")


text:
module DiffEqParamEstim
using DiffEqBase, LsqFit, PenaltyFunctions,
      RecursiveArrayTools, ForwardDiff, Calculus, Distributions, LinearAlgebra, DiffEqSensitivity, Dierckx

STANDARD_PROB_GENERATOR(prob,p) = remake(prob;u0=eltype(p).(prob.u0),p=p)
STANDARD_PROB_GENERATOR(prob::EnsembleProblem,p) = EnsembleProblem(
                             remake(prob.prob;u0=eltype(p).(prob.prob.u0),p=p),
                             output_func = prob.output_func,
                             prob_func = prob.prob_func,
                             reduction = prob.reduction,
                             u_init = prob.u_init)
include("cost_functions.jl")
include("lm_fit.jl")
include("build_loss_objective.jl")
include("build_lsoptim_objective.jl")
include("kernels.jl")
include("two_stage_method.jl")
include("multiple_shooting_objective.jl")

end # module


text:
function nlp_mi_001_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - mix of variable start values
    # - non-linear objective without constraints
    # - minimization objective
    # - functions ^, exp, cos
    # - integer variable
    # - mix of discrete and continuous variables

    model = Model(optimizer)

    @variable(model, x, start = 1)
    @variable(model, y >= 0.1, start = 3.12, Int)
    @variable(model, z >= 1)

    @NLobjective(model, Min, x * exp(x) + cos(y) + z^3 - z^2)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -1.35787195018718, tol = objective_tol)
    return check_solution([x, y, z], [-1, 3, 1], tol = primal_tol)
end


text:
using BinaryBuilder, SHA

LBT_ROOT = abspath(dirname(dirname(@__DIR__)))

# This script heavily borroed from the OpenBLAS build recipe in Yggdrasil
# X-ref: https://github.com/JuliaPackaging/Yggdrasil/blob/master/O/OpenBLAS/common.jl
version = v"0.3.13"
for openblas32 in (true, false)
    sources = [
        ArchiveSource("https://github.com/xianyi/OpenBLAS/archive/v$(version).tar.gz",
                    "79197543b17cc314b7e43f7a33148c308b0807cd6381ee77f77e15acf3e6459e")
    ]

    script = """
    # We always want threading
    flags=(USE_THREAD=1 GEMM_MULTITHREADING_THRESHOLD=50 NO_AFFINITY=1)

    # We are cross-compiling
    flags+=(CROSS=1 PREFIX=/ "CROSS_SUFFIX=\${target}-")

    if [[ \${nbits} == 64 && "$(openblas32)" != "true" ]]; then
        # We're building an ILP64 BLAS with 64-bit BlasInt
        flags+=("LIBPREFIX=libopenblas64_" INTERFACE64=1 SYMBOLSUFFIX=64_)
    else
        flags+=("LIBPREFIX=libopenblas" INTERFACE64=1)
    fi
    flags+=(BINARY=\${nbits})

    # "build" the project, but only install the headers
    cd \${WORKSPACE}/srcdir/OpenBLAS*/

    # Fake out the buildsystem, make it think the libraries have been built
    touch lib.grd
    make "\${flags[@]}" "PREFIX=\${prefix}" dummy
    cp config.h config_last.h

    # Call `make install`, fully expecting this to blow up in our faces.
    make "\${flags[@]}" "PREFIX=\${prefix}" install || true
    """

    products = [
        FileProduct("include/cblas.h", :cblas_h),
        FileProduct("include/f77blas.h", :f77blas_h),
        FileProduct("include/lapack.h", :lapack_h),
        FileProduct("include/lapacke.h", :lapack_h),
        FileProduct("include/openblas_config.h", :openblas_config_h),
    ]
    name = "OpenBLASHeaders$(openblas32 ? "32" : "64")"
    build_tarballs(ARGS, name, version, sources, script, supported_platforms(;experimental=false), products, Dependency[], lock_microarchitecture=false)
end

# Unpack each set of headers into its own directory inside of `include/`:
rm(joinpath(LBT_ROOT, "include"), recursive=true, force=true)
for tarball_path in readdir("products")
    interface = startswith(tarball_path, "OpenBLASHeaders32") ? "LP64" : "ILP64"
    triplet = split(tarball_path, ".v$(version).")[end][1:end-7]
    tarball_path = joinpath("products", tarball_path)
    target_path = joinpath(LBT_ROOT, "include", interface, triplet)

    mkpath(target_path)
    run(`tar -C $(target_path) -xzf $(tarball_path) --strip-components=1 include/`)
end

# de-duplicate headers that share the exact same contents
targets = unique(vcat([readdir(joinpath(LBT_ROOT, "include", interface)) for interface in ("LP64", "ILP64")]...))
function collect_hashes(fname, interfaces=("LP64", "ILP64"))
    return vec([open(sha256, joinpath(LBT_ROOT, "include", interface, target, fname)) for interface in interfaces, target in targets])
end
function is_file_identical_across_targets(fname; interfaces=("LP64", "ILP64"))
    hashes = collect_hashes(fname, interfaces)
    return all(Ref(hashes[1]) .== hashes)
end
function symlink_common_file(fname; interfaces = ("LP64", "ILP64"), common_dir = joinpath(LBT_ROOT, "include", "common"))
    mkpath(common_dir)
    # Copy the file out to `common`
    common_path = joinpath(common_dir, fname)
    rm(common_path; force=true)
    mv(joinpath(LBT_ROOT, "include", interfaces[1], targets[1], fname), common_path)

    for interface in interfaces
        for target in targets
            src_path = joinpath(LBT_ROOT, "include", interface, target, fname)
            if islink(src_path)
                continue
            end
            rm(src_path; force=true)
            symlink(relpath(common_path, dirname(src_path)), src_path)
        end
    end
end

function cleanup_common_headers()
    @info("Cleaning up identical headers across targets/interfaces")
    common_filenames = unique(vcat([readdir(joinpath(LBT_ROOT, "include", interface, target)) for target in targets, interface in ("LP64", "ILP64")]...))
    for fname in common_filenames
        if is_file_identical_across_targets(fname; interfaces=("LP64", "ILP64"))
            @info(" -> $(fname)")
            symlink_common_file(fname; interfaces=("LP64", "ILP64"), common_dir=joinpath(LBT_ROOT, "include", "common"))
        end
    end
    @info("Cleaning up identical headers across target but not interface")
    for interface in ("LP64", "ILP64")
        common_filenames = unique(vcat([readdir(joinpath(LBT_ROOT, "include", interface, target)) for target in targets]...))
        for fname in common_filenames
            if is_file_identical_across_targets(fname; interfaces=(interface,))
                @info(" -> $(interface)/$(fname)")
                symlink_common_file(fname; interfaces=(interface,), common_dir=joinpath(LBT_ROOT, "include", interface, "common"))
            end
        end
    end
end

cleanup_common_headers()

max_stars_repo_name: Gnimuc/libblastrampoline

max_issues_repo_name: Gnimuc/libblastrampoline

max_forks_repo_name: Gnimuc/libblastrampoline


text:
### TRAITS CHECKS ###########

function test_inmemory_supervised_table_dataset(d::D;
        n_obs, n_features, n_targets,
        feature_names=nothing, target_names=nothing,
        Tx=Any, Ty=Any) where {D<:SupervisedDataset}

    @test hasfield(D, :metadata)
    @test hasfield(D, :features)
    @test hasfield(D, :targets)
    @test hasfield(D, :dataframe)
    @test d.metadata isa Dict{String, Any}

    as_df = d.features isa DataFrame

    if as_df
        @test d.features isa DataFrame
        @test d.targets isa DataFrame
        @test d.dataframe isa DataFrame
        @test isempty(intersect(names(d.features), names(d.targets)))
        @test size(d.dataframe) == (n_obs, n_features + n_targets)
        @test size(d.features) == (n_obs, n_features)
        @test size(d.targets) == (n_obs, n_targets)

        # check that dataframe shares the same storage of features and targets
        for c in names(d.dataframe)
            if c in names(d.targets)
                @test d.dataframe[!, c] === d.targets[!,c]
            else
                @test d.dataframe[!, c] === d.features[!,c]
            end
        end

        if feature_names !== nothing
            @test names(d.features) == feature_names
        end
        if target_names !== nothing
            @test names(d.targets) == target_names
        end
    else
        @test d.features isa Matrix{<:Tx}
        @test d.targets isa Matrix{<:Ty}
        @test d.dataframe === nothing
        @test size(d.features) == (n_features, n_obs)
        @test size(d.targets) == (n_targets, n_obs)

        if feature_names !== nothing
            @test d.metadata["feature_names"] == feature_names
        end
        if target_names !== nothing
            @test d.metadata["target_names"] == target_names
        end
    end

    @test d[] === (; d.features, d.targets)
    @test length(d) == n_obs
    @test numobs(d) == n_obs
    idx = rand(1:n_obs)
    @test isequal(d[idx], getobs((; d.features, d.targets), idx))
    @test isequal(d[idx], getobs(d, idx))
    idxs = rand(1:n_obs, 2)
    @test isequal(d[idxs], getobs((; d.features, d.targets), idxs))
    @test isequal(d[idxs], getobs(d, idxs))
end

function test_supervised_array_dataset(d::D;
        n_obs, n_features, n_targets,
        Tx=Any, Ty=Any,
        conv2img=false) where {D<:SupervisedDataset}

    if n_features isa Int
        @assert n_features != 0 "use n_features=() if you don't want features dimensions"
        Nx == 2
    else # tuple
        Nx = length(n_features) + 1
    end
    Ny = map(x -> x == 1 ? 1 : 2, n_targets)

    @test d.features isa Array{Tx, Nx}
    @test size(d.features) == (n_features..., n_obs)

    if Ny isa NamedTuple
        for k in keys(Ny)
            ny = Ny[k]
            @test d.targets[k] isa Array{Ty, ny}
            if ny == 1
                @test size(d.targets[k]) == (n_obs,)
            else
                @test size(d.targets[k]) == (ny, n_obs)
            end
        end
    else
        @assert Ny isa Int
        if Ny == 1
            @test size(d.targets) == (n_obs,)
        else
            @test size(d.targets) == (Ny, n_obs)
        end
    end

    @test length(d) == n_obs
    @test numobs(d) == n_obs
    X, y = d[]
    @test X === d.features
    @test y === d.targets

    idx = rand(1:n_obs)
    @test isequal(d[idx], getobs((; d.features, d.targets), idx))
    @test isequal(d[idx], getobs(d, idx))
    idxs = rand(1:n_obs, 2)
    @test isequal(d[idxs], getobs((; d.features, d.targets), idxs))
    @test isequal(d[idxs], getobs(d, idxs))

    if conv2img
        img = convert2image(d, 1)
        @test img isa AbstractArray{<:Color}
        x = d[1].features
        @test convert2image(D, x) == img
        @test convert2image(d, x) == img
    end
end

function test_unsupervised_array_dataset(d::D;
        n_obs, n_features,
        Tx=Any,
        conv2img=false) where {D<:UnsupervisedDataset}

    n_features = n_features === nothing ? () : n_features
    if n_features isa Int
        @assert n_features != 0 "use n_features = () if you don't want features dimensions"
        Nx == 2
    else # tuple
        Nx = length(n_features) + 1
    end

    @test d.features isa Array{Tx, Nx}
    @test size(d.features) == (n_features..., n_obs)

    @test length(d) == n_obs
    @test numobs(d) == n_obs
    X = d[]
    @test X === d.features

    idx = rand(1:n_obs)
    @test isequal(d[idx], getobs(d.features, idx))
    @test isequal(d[idx], getobs(d, idx))
    idxs = rand(1:n_obs, 2)
    @test isequal(d[idxs], getobs(d.features, idxs))
    @test isequal(d[idxs], getobs(d, idxs))

    if conv2img
        img = convert2image(d, 1)
        @test img isa AbstractArray{<:Color}
        x = d[1].features
        @test convert2image(D, x) == img
        @test convert2image(d, x) == img
    end
end


text:
using MDToolbox

include("AfmizeConfig.jl")
include("Probe.jl")
include("Sphere.jl")

# 底を0にする
function moveBottom(atoms::Array{Sphere})
    bottom = atoms[1].z - atoms[1].r
    for atom in atoms
        bottom = min(bottom, atom.z - atom.r)
    end
    for atom in atoms
        atom.z -= bottom
    end
end

# 針を球として見たとき、どこで衝突するかの計算
function calcCollisionAsSphere(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    dr = probe.r + atom.r
    if 0 < dr^2 - distXY^2
        return atom.z + sqrt(dr^2 - distXY^2) - probe.r
    else
        return 0.0
    end
end

# 針を円形推台として見たとき、どこで衝突するかの計算
function calcCollisionAsCircularThrusters(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    collisionDist = probe.r + atom.r * cos(probe.angle)
    if collisionDist < distXY
        return atom.z + atom.r * sin(probe.angle) - (distXY - collisionDist) / tan(probe.angle) - probe.r
    elseif probe.r < distXY
        return atom.z + sqrt(atom.r^2 - (distXY - probe.r)^2) - probe.r
    else
        return 0.0
    end
end

# 衝突計算をする範囲を返す
function calRectangle(min_point::Point2D, max_point::Point2D, config::AfmizeConfig)
    if max_point.x <= config.range_min.x || max_point.y <= config.range_min.y
        return nothing
    end
    if config.range_max.x <= min_point.x || config.range_max.y <= min_point.y
        return nothing
    end

    min_point = Point2D(max(min_point.x, config.range_min.x),
                        max(min_point.y, config.range_min.y))
    max_point = Point2D(min(max_point.x, config.range_max.x),
                        min(max_point.y, config.range_max.y))
    width  = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    return (max(1, floor(Int, (min_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            max(1, floor(Int, (min_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)),
            min(width , floor(Int, (max_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            min(height, floor(Int, (max_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)))
end

function calRectangle_sphere(atom::Sphere, config::AfmizeConfig)
    min_point = Point2D(atom.x - atom.r - config.probeRadius,
                        atom.y - atom.r - config.probeRadius)
    max_point = Point2D(atom.x + atom.r + config.probeRadius,
                        atom.y + atom.r + config.probeRadius)

    return calRectangle(min_point, max_point, config)
end

function calRectangle_circularThrusters(atom::Sphere, config::AfmizeConfig)
    dist_xy = (tan(config.probeAngle) * (atom.z + atom.r * sin(config.probeAngle) - config.probeRadius)
                + config.probeRadius + atom.r * cos(config.probeAngle))
    min_point = Point2D(atom.x - dist_xy,
                        atom.y - dist_xy)
    max_point = Point2D(atom.x + dist_xy,
                        atom.y + dist_xy)

    return calRectangle(min_point, max_point, config)
end

function afmize_beta(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.x[i], tra.y[i], tra.z[i],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)
    probes = [Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                    config.range_min.y + (h-0.5) * config.resolution.y,
                    config.probeRadius, config.probeAngle)
             for h in 1:height, w in 1:width]

    # 各原子事に計算する価値のある矩形を求めて、その中で探索を行う
    for atom in atoms
        rectangle = calRectangle_sphere(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
        end
    end

    for atom in atoms
        rectangle = calRectangle_circularThrusters(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end

"""
高速化前

function afmize_beta(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.x[i], tra.y[i], tra.z[i],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)

    for h in 1:height, w in 1:width
        probe = Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                      config.range_min.y + (h-0.5) * config.resolution.y,
                      config.probeRadius, config.probeAngle)
        for atom in atoms
            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end
"""

max_stars_repo_name: matsunagalab/BioAFM

max_issues_repo_name: matsunagalab/BioAFM

max_forks_repo_name: matsunagalab/BioAFM


text:
using PlaneGeometry
using PlaneGeometry.Theorems.Napoleon
using SymPy
using Test

@testset "PlaneGeometry.jl" begin
    A = Point(0,0)
    B = Point(1, 3)
    C = Point(4,2)
    tri = Triangle(A, B, C)

    @test vertices(tri) == [A, B, C]
    @test edges(tri) == [Edge(A, B), Edge(B, C), Edge(C, A)]

    c = Circle(A, 2)
    @test PlaneGeometry.center(c) == A

    @test squaredist(A, B) == 10
    @test distance(A, B) == SymPy.sqrt(Sym(10))

    @test midpoint(Point(0, 0), Point(2,4)) == Point(1, 2)
    @test midpoint(A, B) == Point(1//2, 3//2)

    @test median(A, B, C) == Edge(Point(1//2, 3//2), C)

    @test concurrent(edges(tri)) == nothing

    @test concurrent(Edge(A, B), Edge(A, C)) == A

    elist = map(i->median(circshift([A, B, C], i)...), 0:2)
    @test concurrent(elist) == Point(5//3, 5//3)

    @test equipoints(A, B) ==
    [
    Point(1//2 - 3*sympy.sqrt(3)//2, sympy.sqrt(3)//2 + 3//2),
    Point(1//2 + 3*sympy.sqrt(3)//2, 3//2 - sympy.sqrt(3)//2)
    ]

    outer_tri = outer_equitriangles(tri)

    @test outer_tri ==
    [
    Triangle(A, B, Point(1//2 - 3*sympy.sqrt(3)//2, sympy.sqrt(3)//2 + 3//2)),
    Triangle(C, A, Point(sympy.sqrt(3) + 2, 1 - 2*sympy.sqrt(3))),
    Triangle(B, C, Point(sympy.sqrt(3)//2 + 5//2, 5//2 + 3*sympy.sqrt(3)//2))
    ]

    @test_throws ArgumentError circumcircle([A, B])

    outer_circles = circumcircle.(vertices.(outer_tri))

    @test outer_circles ==
    [
    Circle(Point(1//2 - sympy.sqrt(3)//2, sympy.sqrt(3)//6 + 3//2), sympy.sqrt(30)//3),
    Circle(Point(sympy.sqrt(3)//3 + 2, 1 - 2*sympy.sqrt(3)//3), 2*sympy.sqrt(15)//3),
    Circle(Point(sympy.sqrt(3)//6 + 5//2, sympy.sqrt(3)//2 + 5//2), sympy.sqrt(30)//3)
    ]

    @test isequilateral(tri) == false

    np =  napoleon_tri(tri)

    @test isequilateral(np)

    @vars a1 a2 b1 b2 c1 c2
    tri1 = Triangle(a1, a2, b1, b2, c1, c2)

    # This needs to fail.
    @test !isequilateral(tri1)

    @test np == Triangle(
        Point(1//2 - sympy.sqrt(3)//2, sympy.sqrt(3)//6 + 3//2),
        Point(sympy.sqrt(3)//3 + 2, 1 - 2*sympy.sqrt(3)//3),
        Point(sympy.sqrt(3)//6 + 5//2, sympy.sqrt(3)//2 + 5//2))

    @test ceva(tri, O) == (
        [Point(0, 0), Point(1, 3), Point(4, 2), Point(5//2, 5//2), Point(4//9, 2//9), Point(1//9, 1//3)],
        [Edge(Point(1, 3), Point(5//2, 5//2)),
            Edge(Point(5//2, 5//2), Point(4, 2)),
            Edge(Point(4, 2), Point(4//9, 2//9)),
            Edge(Point(4//9, 2//9), Point(0, 0)),
            Edge(Point(0, 0), Point(1//9, 1//3)),
            Edge(Point(1//9, 1//3), Point(1, 3))])

    @test napoleon_proof()

    @test centroid_proof()

    @test ceva_proof()
end


text:
using Colors
function handle_cli(channels)
    blue_ansi = "\033[94m"
    end_ansi = "\033[0m"
    bold_ansi = "\033[1m"
    signal_channel = channels[1]
    config_channel = channels[2]
    stopped = false
    valid_modes = r"mode\s*(?<mode>\d{1,2})\s*"i
    config_command = r"config\s*(?<args>.*)"i
    while !stopped
        print("$(bold_ansi)$(blue_ansi)irradiance>$(end_ansi)")
        temp = readline(STDIN)
        if typeof(temp) == Char
            temp = convert(String, [temp])
        end
        if ismatch(valid_modes, temp)
            mat = match(valid_modes, temp)
            if length(signal_channel.data) > 0
                take!(signal_channel)
            end
            put!(signal_channel, mat[:mode])
        elseif ismatch(config_command, temp)
            handle_config(config_channel, match(config_command, temp)[:args])
        elseif (
                ismatch(r"shutdown|quit"six, temp) ||
                ismatch(r"\x03"six, temp) ||
                ismatch(r"\x02"six, temp)
            )
            if length(signal_channel.data) > 0
                take!(signal_channel)
            end
            put!(signal_channel, "shutdown")
            stopped = true
        end
    end
end
function handle_config(config_channel::Channel, ipt)
    valid_subcommands = r"\s*(?<subcommand>primary_color|secondary_color|scaling|speed)\s*(?<arg>.*)\s*"i
    valid_arguments = Dict{String, Regex}(
        "primary_color"=>r"(?P<arg>#{0,1}[0-9A-F]{6})\s*"i,
        "secondary_color"=>r"(?P<arg>#{0,1}[0-9A-F]{6})\s*"i,
        "scaling"=>r"(?P<arg>[\d+, \d*.\d+])\s*"i,
        "speed"=>r"(?P<arg>[\d+, \d*.\d+])\s*"i
    )
    handlers = Dict{String, Function}(
        "primary_color"=>handle_primary_color,
        "secondary_color"=>handle_secondary_color,
        "scaling"=>handle_scaling_config,
        "speed"=>handle_speed_config
    )

    if ismatch(valid_subcommands, ipt)
        mat1 = match(valid_subcommands, ipt)
        subcom = lowercase(mat1[:subcommand])
        if ismatch(valid_arguments[subcom], mat1[:arg])
            handlers[subcom](config_channel, mat1[:arg])
        else
            println("Sorry, that was an invalid argument.")
        end
    else
        println("Sorry, that subcommand was unrecognized.")
    end
end

function swap_configs(config_channel::Channel{EffectConfig}, new_config::EffectConfig)
    put!(config_channel, new_config)
    take!(config_channel)
end

function handle_primary_color(config_channel::Channel{EffectConfig}, argument)
    letters_only_regex = r"$[0-9A-F]{6}\s*"i
    if ismatch(letters_only_regex, argument)
        argument = "#$argument"
    end
    cur_con = fetch(config_channel)
    new_prim_color = cur_con.primary_color
    try
        new_prim_color = convert(HSL, parse(Colorant, argument))
    catch
        println("Sorry, failed to parse that input.")
    end
    new_config = EffectConfig(new_prim_color, cur_con.secondary_color, cur_con.scaling, cur_con.speed, cur_con.special)
    swap_configs(config_channel, new_config)
end

function handle_secondary_color(config_channel::Channel{EffectConfig}, argument)
    cur_con = fetch(config_channel)
    new_sec_color = cur_con.secondary_color
    try
        new_sec_color = convert(HSL, parse(Colorant, argument))
    catch
        println("Sorry, failed to parse that input.")
    end
    new_config = EffectConfig(cur_con.primary_color, new_sec_color, cur_con.scaling, cur_con.speed, cur_con.special)
    swap_configs(config_channel, new_config)
end

function handle_scaling_config(config_channel::Channel{EffectConfig}, argument)
    cur_con = fetch(config_channel)
    new_scaling = cur_con.scaling
    try
        new_scaling = parse(Float64, argument)
    catch
        println("Sorry, failed to parse that input.")
    end
    new_config = EffectConfig(cur_con.primary_color, cur_con.secondary_color, new_scaling, cur_con.speed, cur_con.special)
    swap_configs(config_channel, new_config)
end

function handle_speed_config(config_channel::Channel{EffectConfig}, argument)
    cur_con = fetch(config_channel)
    new_speed = cur_con.speed
    try
        new_speed = parse(Float64, argument)
    catch
        println("Sorry, failed to parse that input.")
    end
    new_config = EffectConfig(cur_con.primary_color, cur_con.secondary_color, cur_con.scaling, new_speed, cur_con.special)
    swap_configs(config_channel, new_config)
end


text:
"""This module builds a LALR(1) transition-table for lalr_parser.jl

For now, shift/reduce conflicts are automatically resolved as shifts.
"""

# Author: Erez Shinan (2017)
# Email : erezshin@gmail.com

# Adapted from Python by James Hester (2019)
# Email: james.r.hester@gmail.com

#==
struct Action
    name::Symbol
end
abstract type Action end
struct Shift <: Action end
struct Reduce <: Action end

#const Shift = Action("Shift")
#const Reduce = Action("Reduce")
==#

# States_convert matches integer state with a Set{RulePtr}, stored
# for debugging purposes
#
struct ParseTable
    states::Array{Dict,1}
    start_states::Dict{String,Int64}
    end_states::Dict{String,Int64}
    states_convert::Dict{Set{RulePtr},Int64}
end

ParseTable(rawstates,start_states,end_states) = begin
    enum = collect(keys(rawstates))
    state_to_idx = Dict((s=>i) for (i,s) in enumerate(enum))
    int_states = Array{Dict{String,Tuple{Symbol,Any}},1}(undef,length(rawstates))
    #==println("\nState table:")
    for s in keys(state_to_idx)
        print("\n$s")
        print("$(state_to_idx[s])")
    end
    ==#

    for (s,la) in rawstates
        la = Dict([(k=>if v[1] == :shift
                       (v[1],state_to_idx[v[2]])
                       else
                        v end) for (k,v) in la])
        int_states[state_to_idx[s] ] = la
    end
    start_states = Dict([start=>state_to_idx[s] for (start,s) in start_states])
    end_states = Dict([start=>state_to_idx[s] for (start,s) in end_states])
    return ParseTable(int_states,start_states,end_states,state_to_idx)
end

# digraph and traverse, see The Theory and Practice of Compiler Writing

# computes F(x) = G(x) union (union { G(y) | x R y })
# X: nodes
# R: relation (function mapping node -> list of nodes that satisfy the relation)
# G: set valued function
digraph(X, R, G) = begin
    F = Dict()
    S = []
    N = Dict{Any,Int64}()
    for x in X
        N[x] = 0
    end
    for x in X
        # this is always true for the first iteration, but N[x] may be updated in traverse below
        if N[x] == 0
            traverse(x, S, N, X, R, G, F)
        end
    end
    return F
end

# x: single node
# S: stack
# N: weights
# X: nodes
# R: relation (see above)
# G: set valued function
# F: set valued function we are computing (map of input -> output)
traverse(x, S, N, X, R, G, F) = begin
    push!(S,x)
    d = length(S)
    N[x] = d
    F[x] = G[x]
    for y in R[x]
        if N[y] == 0
            traverse(y, S, N, X, R, G, F)
        end
        n_x = N[x]
        @assert n_x > 0
        n_y = N[y]
        @assert n_y != 0
        if (n_y > 0) && (n_y < n_x)
            N[x] = n_y
        end
        union!(F[x],F[y])
    end
    if N[x] == d
        f_x = F[x]
        while true
            z = pop!(S)
            N[z] = -1
            F[z] = f_x
            if z == x
                break
            end
        end
    end
end

# Grammar analysis

mutable struct LALR_Analyzer <: GrammarAnalyzer
    debug
    rules_by_origin
    start_states
    end_states
    parse_table
    FIRST
    FOLLOW
    NULLABLE
    nonterminal_transitions
    directly_reads::DefaultDict
    reads::DefaultDict
    includes::DefaultDict
    lookback::DefaultDict
    lr0_states::Set
    lr0_start_states::Dict
    lr0_rules_by_origin
end

LALR_Analyzer(parser_conf;debug=true) = begin
    l = LALR_Analyzer(debug,nothing,nothing,nothing,nothing,nothing,nothing,nothing,[],DefaultDict{Any,Set}(),DefaultDict{Any,Set}(),DefaultDict{Any,Set}(),DefaultDict{Any,Set}(),Set(),Dict(),[])
    init_analyser!(l,parser_conf,debug=debug)
    return l
end

compute_lr0_states(l::LALR_Analyzer) = begin
    l.lr0_states = Set()
    cache = Dict()
    step(state) = Channel() do state_chan
        _, unsat = classify_bool(state.closure,rp->is_satisfied(rp))
        d = classify(unsat,key=rp->next(rp))
        for (sym,rps) in d
            kernel = Set([advance(rp,sym) for rp in rps])
            new_state = get(cache,kernel,nothing)
            if new_state === nothing
                closure = copy(kernel)
                for rp in kernel
                    if !is_satisfied(rp) && !is_terminal(next(rp))
                        union!(closure,expand_rule(l,next(rp),l.lr0_rules_by_origin))
                    end
                end
                new_state = LR0ItemSet(kernel,closure)
                cache[kernel] = new_state
            end
            state.transitions[sym] = new_state
            put!(state_chan,new_state)
        end
        push!(l.lr0_states,state)
    end
    collect(bfs(values(l.lr0_start_states),step))
    #println("lr0_states after: $(l.lr0_states)")
end

compute_reads_relations(l::LALR_Analyzer) = begin
    # handle start state
    for root in values(l.lr0_start_states)
        @assert length(root.kernel) == 1
        for rp in root.kernel
            @assert rp.index == 0
            l.directly_reads[(root, next(rp))] = Set([ Terminal("\$END") ])
        end
    end
    for state in l.lr0_states
        seen = Set()
        for rp in state.closure
            if is_satisfied(rp)
                continue
            end
            s = next(rp)
            # if s is a not a nonterminal
            if !(s in keys(l.lr0_rules_by_origin))
                continue
            end
            if s in seen
                continue
            end
            push!(seen,s)
            nt = (state, s)
            push!(l.nonterminal_transitions,nt)
            dr = l.directly_reads[nt]
            r = l.reads[nt]
            next_state = state.transitions[s]
            for rp2 in next_state.closure
                if is_satisfied(rp2) continue end
                s2 = next(rp2)
                # if s2 is a terminal
                if !(s2 in keys(l.lr0_rules_by_origin))
                    push!(dr,s2)
                end
                if s2 in l.NULLABLE
                    push!(r,(next_state, s2))
                end
            end
        end
    end
    #println("Nonterminal transitions: $(length(l.nonterminal_transitions))")
    #__s = sort(l.nonterminal_transitions)
    #for (i,s1) in enumerate(__s) println("$(i-1): $s1") end
end

compute_includes_lookback(l::LALR_Analyzer) = begin
    for nt in l.nonterminal_transitions
        state, nonterminal = nt
        #println("\n===\nState $state\n")
        includes = []
        lookback = l.lookback[nt]
        for rp in state.closure
            if rp.rule.origin != nonterminal
                continue
            end
            #println("Rule $rp")
            # traverse the states for rp(.rule)
            state2 = state
            for i in rp.index:(length(rp.rule.expansion)-1)
                s = rp.rule.expansion[i+1]
                nt2 = (state2, s)
                state2 = state2.transitions[s]
                if !(nt2 in keys(l.reads))
                    continue
                end
                _loop_success = true #emulate for..else in Python
                for j in (i + 1):(length(rp.rule.expansion)-1)
                    if !(rp.rule.expansion[j+1] in l.NULLABLE)
                        _loop_success = false
                        break
                    end
                end
                # mimic for...else in Python code
                if _loop_success
                    push!(includes,nt2)
                end
            end
            # state2 is at the final state for rp.rule
            if rp.index == 0
                for rp2 in state2.closure
                    if (rp2.rule == rp.rule) && is_satisfied(rp2)
                        push!(lookback,(state2, rp2.rule))
                    end
                end
            end
        end

        for nt2 in includes
            push!(l.includes[nt2],nt)
        end
    end
end

compute_lookaheads(l::LALR_Analyzer) = begin
    #println("$(length(l.nonterminal_transitions)), $(length(l.reads)), $(length(l.directly_reads))")
    read_sets = digraph(l.nonterminal_transitions, l.reads, l.directly_reads)

    follow_sets = digraph(l.nonterminal_transitions, l.includes, read_sets)

    for (nt, lookbacks) in l.lookback
        for (state, rule) in lookbacks
            for s in follow_sets[nt]
                push!(state.lookaheads[s],rule)
            end
        end
    end
end

compute_lalr1_states(l::LALR_Analyzer) = begin
    m = IdDict()
    reduce_reduce = []
    for state in l.lr0_states
        actions = IdDict()
        for (la, next_state) in state.transitions
            actions[la] = (:shift, next_state.closure)
        end

        for (la, rules) in state.lookaheads
            if length(rules) > 1
                # Try to resolve conflict based on priority
                p = [(r.options.priority || 0, r) for r in rules]
                sort!(p,by=r -> r[1], rev=true)
                best, second_best = p[1:2]
                if best[1] > second_best[1]
                    rules = [best[2]]
                else
                    push!(reduce_reduce,(state, la, rules))
                end
            end
            if la in keys(actions)
                if l.debug
                    @warn "Shift/Reduce conflict for terminal $(la.name): (resolving as shift)"
                    @warn " * $(first(rules))"
                end
            else
                actions[la] = (:reduce, first(rules))
            end
        end

        m[state] = IdDict([k.name => v for (k, v) in actions])
    end

    if !isempty(reduce_reduce)
        msgs = []
        for (state, la, rules) in reduce_reduce
            msg = "Reduce/Reduce collision in $la between the following rules: "*join([ "\n\t- $r"  for r in rules ])
            if l.debug
                msg *= "\n    collision occurred in state: {" * join(["\n\t$x" for x in state.closure]) * "\n    }"
            end
            push!(msgs,msg)
        end
        throw(GrammarError(join(msgs,"\n\n")))
    end

    states = IdDict([k.closure => v for (k, v) in m ])

    # compute end states
    end_states = Dict()
    for state in keys(states)
        for rp in state
            for start in keys(l.lr0_start_states)
                if rp.rule.origin.name == ("\$root_" * start) && is_satisfied(rp)
                    @assert !(start in keys(end_states))
                    end_states[start] = state
                end
            end
        end
    end

    l.parse_table = ParseTable(states, Dict([start=> state.closure for (start, state) in l.lr0_start_states]), end_states)

end

compute_lalr!(l::LALR_Analyzer) = begin
    compute_lr0_states(l)
    compute_reads_relations(l)
    compute_includes_lookback(l)
    compute_lookaheads(l)
    compute_lalr1_states(l)
end


text:
"""`duplicated(x)` returns whether `x` contains any duplicate values"""
function duplicated(x::AbstractArray)
    xsort = sort(x)
    for i = 1:(length(xsort)-1)
        if xsort[i] == xsort[i+1]
            return true
        end
    end
    false
end


text:
import Base.Test: @testset, @test, @test_throws
import StaticArrays: SVector

using MultipleScattering

sound_p = Acoustic(1., 4. + 0.0im,2)
particles = [Particle(sound_p,Sphere([x,0.0], .5)) for x =0.:2.:10.]

sound_sim = Acoustic(1., 0.5 + 0.0im,2)
source = plane_source(sound_sim, SVector(0.0,0.0), SVector(1.0,0.0), 1.)
sim = FrequencySimulation(particles, source)

ω = 0.5

bounds = Box([[-0.,-1.], [10.,1.]])
ωs = 0.:0.1:1.
simres = run(sim, bounds, ωs)
timres = run(sim, bounds, ωs; ts = [30.], result_in_time=true)

using Plots ; #pyplot()

plot(sim,ω)
plot!(sim)

plot(simres; x = [[1,1]])
plot(simres,0.3)
plot(simres, 0.3, seriestype=:contour)

# include("plot_field.jl")
# plot_field(simres, linetype=:contour)
# plot_field(timres, linetype=:contour)
# plot!(sim)

amplitude = 1.0
# Create new souce as a linear combination of two other sources
s1 = point_source(sound_sim, SVector(8.0,0.6), amplitude)
s2 = point_source(sound_sim, SVector(2.0,-0.6), amplitude)
source = s1 + s2

sound_sim = Acoustic(1., 0.5 + 0.0im,2)
sim = FrequencySimulation(particles, source)

ω = 1.
# pyplot()
plot(sim,ω;res=60)
plot!(sim)


text:
using Libdl

## Collected pieces of metadata that we will track over time
# This value will be overridden in __init__()
time_start = now()

# Last time we updated the registry
last_registry_update = now()

# Total served requests: incremented in serve_file()
total_hits = Int64(0)
# Total cache hits: incremented in fetch()
cached_hits = Int64(0)
# Total fetch successes: incremented in fetch()
fetch_hits = Int64(0)
# Total misses: incremented in start()
total_misses = Int64(0)
# Total payload bytes received from StorageServers
payload_bytes_received = Int64(0)
# Total payload bytes transmitted to clients
payload_bytes_transmitted = Int64(0)

## Functions to serve metadata about the PkgServer itself

const pkgserver_version = Ref{Union{Nothing,String}}(nothing)
function get_pkgserver_version()
    if pkgserver_version[] === nothing
        # Get PkgServer.jl's version and git sha
        version = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
        try
            repo = LibGit2.GitRepo(dirname(@__DIR__))
            gitsha = string(LibGit2.GitHash(LibGit2.GitCommit(repo, "HEAD")))
            pkgserver_version[] = "$(version)-$(gitsha)"
        catch
            pkgserver_version[] = "$(version)"
        end
    end
    return pkgserver_version[]
end

# Return the fully-qualified domain name of this server, if set.
const pkgserver_url = Ref{Union{Nothing,String}}(nothing)
function get_pkgserver_url()
    if pkgserver_url[] === nothing
        # Default to `https://pkg.julialang.org`
        fqdn = get(ENV, "JULIA_PKG_SERVER_FQDN", "")
        if isempty(fqdn)
            fqdn = "pkg.julialang.org"
        end
        pkgserver_url[] = string("https://", fqdn)
    end
    return pkgserver_url[]
end

# Return all officially-run PkgServer instances we know about
function get_pkgserver_siblings()
    # Note that some of these (`us-east`, `us-west` and `eu-central` at the
    # time of this writing) are loadbalancers with multiple PkgServer instances
    # behind the LB, called things like `us-west1` and `us-east3`, etc...
    # We surface the loadbalancer URL as that is the primary endpoint that
    # users should hit when making requests.
    regions = [
        # North America
        "us-west",
        "us-east",
        # South America
        "sa",
        # Europe
        "eu-central",
        # Asia/Oceania
        "in",
        "kr",
        "sg",
        # Special Chinese servers, serving from within-country
        # Note that to connect directly, you should connect to
        # `cn-X.pkg.juliacn.com`; `cn-X.pkg.julialang.org` will
        # redirect you.  You can always discover the canonical
        # address of a server through the `/meta` endpoint.
        "cn-southeast",
        "cn-east",
        "cn-northeast",
        # Australia
        "au",
    ]
    return [string("https://", region, ".pkg.julialang.org") for region in regions]
end

function get_num_hashnamed_files(dir)
    # If this directory doesn't exist, then we haven't cached anything!
    if !isdir(dir)
        return 0
    end

    num_files = 0
    for f in readdir(dir)
        # only consider things with a name length that is exactly a tree-hash (filters out in-progress downloads, etc...)
        # TODO: if we change hashes (e.g. SHA256), add that hash length here as well
        if length(f) != 40
            continue
        end

        # Only consider actual files
        f_path = joinpath(dir, f)
        if !isfile(f_path)
            continue
        end

        num_files += 1
    end
    return num_files
end

function get_num_pkgs_cached()
    cache_root = config.cache.root
    # If we don't even have a `package` directory, then don't worry about it
    if !isdir(joinpath(cache_root, "package"))
        return 0
    end

    num_pkgs_cached = 0
    for d in readdir(joinpath(cache_root, "package"))
        # Only consider things with a name that is exactly a UUID (filters out in-progress downloads, etc...)
        if length(d) != 36
            continue
        end

        # Only consider directories
        pkg_dir = joinpath(cache_root, "package", d)
        if !isdir(pkg_dir)
            continue
        end

        num_pkgs_cached += get_num_hashnamed_files(pkg_dir)
    end
    return num_pkgs_cached
end

function get_num_artifacts_cached()
    num_artifacts_cached = get_num_hashnamed_files(joinpath(config.cache.root, "artifact"))
end

function serve_data(http::HTTP.Stream, msg, content_type)
    HTTP.setheader(http, "Content-Length" => string(length(msg)))
    HTTP.setheader(http, "Content-Type" => content_type)
    startwrite(http)
    if http.message.method == "GET"
        write(http, msg)
    end
    return
end

function serve_json(http::HTTP.Stream, data)
    return serve_data(http, JSON3.write(data), "application/json")
end

libjulia_internal = C_NULL
function get_num_live_tasks()
    global libjulia_internal
    if libjulia_internal == C_NULL
        libjulia_internals = filter(f -> occursin("libjulia-internal", basename(f)), Libdl.dllist())
        if isempty(libjulia_internals)
            # e.g. we're running on Julia v1.5
            return nothing
        end
        libjulia_internal = dlopen(first(libjulia_internals))
    end
    jl_live_tasks_addr = dlsym(libjulia_internal, "jl_live_tasks"; throw_error=false)
    if jl_live_tasks_addr === nothing
        return nothing
    end
    live_tasks = ccall(jl_live_tasks_addr, Vector, ())
    return length(filter(t -> t.state == :runnable, live_tasks))
end

function serve_meta(http::HTTP.Stream)
    # We serve a JSON representation of some metadata about this PkgServer
    task_state(t) = istaskfailed(t) ? "failed" :
                    istaskdone(t) ? "done" :
                    istaskstarted(t) ? "started" : "not-started"
    metadata = Dict(
        "pkgserver_version" => get_pkgserver_version(),
        "pkgserver_url" => get_pkgserver_url(),
        "julia_version" => string(VERSION),
        "start_time" => string(time_start),
        "last_registry_update" => string(last_registry_update),
        "registry_update_task" => task_state(registry_update_task),
        "registry_watchdog_task" => task_state(registry_watchdog_task),
        "maxrss" => Int(Sys.maxrss()),
    )
    live_tasks = get_num_live_tasks()
    if live_tasks !== nothing
        metadata["live_tasks"] = string(live_tasks)
    end
    return serve_json(http, metadata)
end

# Because we need to JSON-serialize cache entries, we must define the serialization strategy
StructTypes.StructType(::Type{FilesystemDatastructures.CacheEntry}) = StructTypes.Struct()
function serve_meta_stats(http::HTTP.Stream)
    # Serve a JSON representation of the top 100 most commonly requested resources
    stats = Dict(
        "lru" => config.cache.entries,
        "packages_cached" => get_num_pkgs_cached(),
        "artifacts_cached" => get_num_artifacts_cached(),

        ## Global response statistics:
        # Number of requests we've had total.  Should equal `cached_hits + fetch_hits`.
        "total_hits" => total_hits,
        # Number of requests served through the happy path; file was already present on disk.
        "cached_hits" => cached_hits,
        # Number of requests that were served by fetching it from the server.
        # Note; if three simultaneous requests come in for a file that we need to fetch,
        # we'll only fetch it once, but we'll stream the content back for all three, and
        # account all three as a `fetch_hit`.
        "fetch_hits" => fetch_hits,
        # Number of requests that could not be served
        "total_misses" => total_misses,
        "payload_bytes_received" => payload_bytes_received,
        "payload_bytes_transmitted" => payload_bytes_transmitted,
    )
    return serve_json(http, stats)
end

function serve_siblings(http::HTTP.Stream)
    return serve_json(http, get_pkgserver_siblings())
end

function serve_children(http::HTTP.Stream)
    return serve_json(http, String[])
end

function serve_parents(http::HTTP.Stream)
    return serve_json(http, config.storage_servers)
end

function serve_robots_txt(http::HTTP.Stream)
    return serve_data(http, "User-agent: * Disallow: /\n", "text/plain")
end


text:
module QSharp
    using JSON
    using UUIDs

    struct Message
        context :: Union{String, Nothing}
        data :: Any
    end

    mutable struct QSharpServer
        Server :: Base.Process
        MessageBuffer :: Vector{Message}

        QSharpServer(cmd :: Cmd) = begin
            new(
                open(cmd, read = true, write = true),
                Vector{Dict{Any}}()
            )
        end
    end

    # Internal methods

    function start() :: QSharpServer
        QSharpServer(`dotnet run --no-build --project ../../qsharp-server/src/qsharp-server.csproj`)
    end

    function send_command(server :: QSharpServer, command :: String, payload :: Any = Dict()) :: String
        context = string(UUIDs.uuid4())
        JSON.print(server.Server, Dict(
            :command => command,
            :payload => payload,
            :context => context
        ))
        write(server.Server, "\n")
        context
    end

    function stop(server :: QSharpServer)
        send_command(server, "exit")
    end

    # TODO: Make async!
    function readnext(server :: QSharpServer, filter)
        next_msg_json = JSON.parse(server.Server)
        context = haskey(next_msg_json, "context") ? pop!(next_msg_json, "context") :  nothing
        next_msg = Message(context, next_msg_json)
        if filter(next_msg)
            next_msg
        else
            push!(server.MessageBuffer, next_msg)
            nothing
        end
    end

    # TODO: Make async!
    function listen(server :: QSharpServer, context :: Union{String, Nothing})
        filter = if isnothing(context)
            msg -> true
        else
            msg -> msg.context == something(context)
        end
        # Do we already have a message in the buffer?
        idx_msg = findfirst(filter, server.MessageBuffer)
        if isnothing(idx_msg)
            # Need to read the next one.
            while true
                next = readnext(server, filter)
                if !isnothing(next)
                    return next.data
                end
            end
        else
            return popat!(server.MessageBuffer, idx_msg).data
        end

    end

    struct Callable
        full_name :: String
    end

    # TODO: allow passing input
    function invoke(callable :: Callable; simulator :: Union{String, Nothing} = nothing)
        global server
        context = send_command(server, "simulate", Dict(
            "operation" => callable.full_name,
            "input" => Dict(),
            "simulator" => simulator
        ))
        while true
            next_msg = listen(server, context)
            if haskey(next_msg, "data") && haskey(next_msg["data"], "output")
                # Actual output, so we can return.
                return next_msg["data"]["output"]
            elseif haskey(next_msg, "display_output")
                # TODO: Process as Jupyter display!
                println("Displayable output: $(next_msg["display_output"])")
            elseif haskey(next_msg, "console_message")
                stream = get(next_msg, "stream", "StandardOut")
                if stream == "StandardOut"
                    println(next_msg["console_message"])
                elseif stream == "StandardError"
                    println(stderr, next_msg["console_message"])
                else
                    println("Unrecognized stream $stream.")
                end
            else
                println("Unrecognized message $next_msg")
            end
        end
    end

    # Public API

    server = start()

    function compile(code :: String)
        global server
        context = send_command(server, "compile", Dict(
            :source => code
        ))
        response = listen(server, context)
        # TODO: Handle warnings
        [Callable(callable_name) for callable_name in response["data"]["compiled_callables"]]
    end

end

max_stars_repo_name: dfinke/qsharp-server

max_issues_repo_name: dfinke/qsharp-server

max_forks_repo_name: dfinke/qsharp-server


text:
function translate(cursor::CLDeclStmt)
    child_cursors = children(cursor)
    decl_exprs = MetaExpr[]
    for c in child_cursors
        meta = translate(c)
        if meta.expr isa Symbol
            meta.expr = Expr(:local, meta.expr)
            push!(decl_exprs, meta)
        else
            push!(decl_exprs, meta)
        end
    end
    return decl_exprs
end

function translate(cursor::CLCompoundStmt)
    child_cursors = children(cursor)
    meta = MetaExpr[]
    block = Expr(:block)
    for c in child_cursors
        append!(meta, translate(c))
    end
    for m in meta
        push!(block.args, m.expr)
        # TODO: add space between stmt
    end
    return MetaExpr(block)
end

function translate(cursor::CLIfStmt)
    child_cursors = children(cursor)
    if_expr = Expr(:if)
    for c in child_cursors
        if kind(c) == CXCursor_IfStmt
            elseif_expr = translate(c).expr
            elseif_expr.head = :elseif
            push!(if_expr.args, elseif_expr)
        else
            push!(if_expr.args, translate(c).expr)
        end
    end
    return MetaExpr(if_expr)
end

translate(cursor::CLBreakStmt) = MetaExpr(Expr(:break), cursor)
translate(cursor::CLContinueStmt) = MetaExpr(Expr(:continue), cursor)

function translate(cursor::CLGotoStmt)
    meta = translate(first(children(cursor)))
    meta.expr = Expr(:macrocall, Symbol("@goto"), nothing, meta.expr)
    return meta
end

function translate(cursor::CLLabelStmt)
    ex = Expr(:macrocall, Symbol("@label"), nothing, Symbol(spelling(cursor)))
    return MetaExpr(ex, cursor)
end

translate(cursor::CLLabelRef) = MetaExpr(Symbol(spelling(cursor)), cursor)

function translate(cursor::CLWhileStmt)
    child_cursors = children(cursor)
    condition_meta = translate(first(child_cursors))
    body_meta = translate(last(child_cursors))
    return MetaExpr(Expr(:while, condition_meta.expr, body_meta.expr))
end

function translate(cursor::CLDoStmt)
    child_cursors = children(cursor)
    condition_meta = translate(last(child_cursors))
    body_meta = translate(first(child_cursors))
    push!(body_meta.expr.args, Expr(:call, :||, condition_meta.expr, Expr(:break)))
    return MetaExpr(Expr(:while, true, body_meta.expr))
end

# this is a workaround since libclang somehow striped those null objects in `for( ; ; )`
function get_init_cond_inc_body(cursor::CLForStmt)
    child_cursors = children(cursor)
    toks = tokenize(cursor) |> collect
    left_paren_idx = 2
    right_paren_idx = findfirst(x->x.text==")", toks)
    idxs = findall(x->x.text==";", toks[1:right_paren_idx])

    idx_count = 1
    if idxs[1] == left_paren_idx+1
        init = MetaExpr(Expr(:null), cursor)
    else
        init = translate(child_cursors[idx_count])
        idx_count += 1
    end

    if idxs[2] == idxs[1]+1
        condition = MetaExpr(Expr(:null), cursor)
    else
        condition = translate(child_cursors[idx_count])
        idx_count += 1
    end

    if right_paren_idx == idxs[2]+1
        increment = MetaExpr(Expr(:null), cursor)
    else
        increment = translate(child_cursors[idx_count])
        idx_count += 1
    end

    @assert idx_count == length(child_cursors) "unknown for-loop structure: $cursor"
    body = translate(last(child_cursors))

    return init, condition, increment, body
end

function translate(cursor::CLForStmt)
    init, cond, inc, body = get_init_cond_inc_body(cursor)
    if init isa Vector{MetaExpr}
        init = init[]
    end
    for_expr = Expr(:macrocall, Symbol("@cfor"), nothing, init.expr, cond.expr, inc.expr, body.expr)
    return MetaExpr(for_expr)
end

function translate(cursor::CLReturnStmt)
    child_cursors = children(cursor)
    if isempty(child_cursors)
        return MetaExpr(Expr(:return, nothing), cursor)
    else
        meta = translate(first(child_cursors))
        meta.expr = Expr(:return, meta.expr)
        return meta
    end
end

translate(cursor::CLNullStmt) = MetaExpr(Expr(:null), cursor)

function translate(cursor::CLCaseStmt)
    child_cursors = children(cursor)
    case_meta = translate(first(child_cursors))
    case_meta.expr = Expr(:macrocall, Symbol("@case"), nothing, case_meta.expr)
    exprs = MetaExpr[case_meta]
    append!(exprs, translate(last(child_cursors)))
    return exprs
end

function translate(cursor::CLDefaultStmt)
    child_cursors = children(cursor)
    default_meta = MetaExpr(Expr(:macrocall, Symbol("@default"), nothing))
    exprs = MetaExpr[default_meta]
    append!(exprs, translate(last(child_cursors)))
    return exprs
end

function translate(cursor::CLSwitchStmt)
    child_cursors = children(cursor)
    constexpr_meta = translate(first(child_cursors))
    body_meta = translate(last(child_cursors))
    ex = Expr(:macrocall, Symbol("@switch"), nothing, constexpr_meta.expr, body_meta.expr)
    return MetaExpr(ex)
end


text:
"""
DuckDB validity mask
"""
mutable struct ValidityMask
    data::Vector{UInt64}

    function ValidityMask(data::Vector{UInt64})
        result = new(data)
        return result
    end
end

const BITS_PER_VALUE = 64;

function GetEntryIndex(row_idx)
    return ((row_idx - 1) ÷ BITS_PER_VALUE) + 1
end

function GetIndexInEntry(row_idx)
    return (row_idx - 1) % BITS_PER_VALUE
end

function SetInvalid(mask::ValidityMask, index)
    entry_idx = GetEntryIndex(index)
    index_in_entry = GetIndexInEntry(index)
    mask.data[entry_idx] &= ~(1 << index_in_entry)
    return
end

function IsValid(mask::ValidityMask, index)::Bool
    entry_idx = GetEntryIndex(index)
    index_in_entry = GetIndexInEntry(index)
    return (mask.data[entry_idx] & (1 << index_in_entry)) != 0
end


text:
function _do_serial_mem(::Type{M}, ::Type{D}, filename::String, count_mode::CountMode, i) where {M<:AbstractMer,D<:ReadDatastore}
    ds = open(D, filename)
    part = i:nworkers():length(ds)
    counts = serial_mem(M, buffer(ds), count_mode, part)
    close(ds)
    return counts
end

function dist_mem(::Type{M}, ::Type{D}, filename::String, count_mode::CountMode) where {M<:AbstractMer,D<:ReadDatastore}
    nw = nworkers()
    @info "Splitting all reads in $filename, across $nw worker processes and counting kmers"
    if nw === 1
        @warn "There are not any extra worker processes available. This method will not provide any speed benefit over serial_mem."
    end
    local_counts = Vector{Vector{MerCount{M}}}(undef, nw)
    @sync begin
        @inbounds for (i, w) in enumerate(workers())
            @async local_counts[i] = remotecall_fetch(_do_serial_mem, w, M, D, filename, count_mode, i)
        end
        @info "Waiting for workers to finish"
    end
    @info "Merging all counts from $nw worker processes"
    final_count = local_counts[1]
    @inbounds for i in 2:lastindex(local_counts)
        unsafe_merge_into!(final_count, local_counts[i])
    end
    return final_count
end

function dist_mem(::Type{M}, input::ReadDatastore, count_mode::CountMode) where {M<:AbstractMer}
    return dist_mem(M, typeof(input), input.filename, count_mode)
end

function dist_mem(::Type{M}, input::DatastoreBuffer{<:ReadDatastore}, count_mode::CountMode) where {M<:AbstractMer}
    return dist_mem(M, ReadDatastores.datastore(input), count_mode)
end


text:
using Pkg
Pkg.build() # Pkg.build(; verbose = true) for Julia 1.1 and up
Pkg.test(coverage=true)

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
########################################################################
# Observed info -  for 2PL pars and 1-D latent                         #
########################################################################
"""
```julia
_item_observed_information(
    parameters::Parameters1PL,
    latent::Latent1D,
    response_val::Float64,
)
```

# Description

It is equal to the item expected information.

# Arguments

- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`response_val::Float64`** : Required. A scalar response.

# Output
A `Float64` scalar observed (expected) informations.
"""
function _item_observed_information(
    parameters::Parameters1PL,
    latent::Latent1D,
    response_val::Float64
)
    return _item_expected_information(parameters, latent)::Float64
end

########################################################################
# Observed info -  for 2PL pars and 1-D latent                         #
########################################################################
"""
```julia
_item_observed_information(
    parameters::Parameters2PL,
    latent::Latent1D,
    response_val::Float64,
)
```

# Description

Theoretically, it should be equal to its expected counterpart.
However, it is used to compute the standard errors of item parameter estimates which are based on the log-likelihood (which is based on the responses).
Look at the example in the documentation of the exported method.

# Arguments

- **`parameters::Parameters2PL`** : Required. A 2-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`response_val::Float64`** : Required. A scalar response.

# Output
A ``2 \times 2`` matrix of the observed informations.
"""
function _item_observed_information(
    parameters::Parameters2PL,
    latent::Latent1D,
    response_val::Float64
)
    #return _item_expected_information(parameters, latent)::Matrix{Float64}
    p = _probability(latent, parameters)
    p_1_p = p * (1 - p)
    i_aa = p_1_p * (latent.val - parameters.b)^2 # v p.49 Kim, Baker: - L_11
    i_ab =
        (
             ((parameters.a * (latent.val - parameters.b) * (- p) * p_1_p) +
             (response_val - p))
        )  # X p.50 Kim, Baker: - L_12 #TODO
    i_bb = - parameters.a^2 * p_1_p # v p.49 Kim, Baker: - L_22
    return [i_aa i_ab; i_ab i_bb]::Matrix{Float64}
end

########################################################################
# Observed info -  for 3PL pars and 1-D latent                         #
########################################################################
"""
```julia
_item_observed_information(
    parameters::Parameters3PL,
    latent::Latent1D,
    response_val::Float64,
)
```

# Description

It computes the observed information (-second derivative of the likelihood) with respect to the 3 parameters of the 3PL model.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`parameters::Parameters3PL`** : Required. A 3-parameter logistic parameters object.
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`response_val::Float64`** : Required. A scalar response.

# Output
A ``3 \times 3`` matrix of the observed informations.
"""
function _item_observed_information(
    parameters::Parameters3PL,
    latent::Latent1D,
    response_val::Float64,
)
    p = _probability(latent, parameters)
    i = (1 - p) * (p - parameters.c)
    h = (response_val * parameters.c - p^2) * i
    j = response_val * i
    den = ((1 - parameters.c) * p)^2
    i_aa = - h * (latent.val - parameters.b)^2 / den # v p.49 Kim, Baker: - L_11
    i_ab =
        (
            (p - parameters.c) * ((parameters.a * (latent.val - parameters.b) * h) +
            (p * (response_val - p) * (1 - parameters.c)))
        ) / den # X p.50 Kim, Baker: - L_12 #TODO
    i_ac = j * (latent.val - parameters.b) / den # v p.51 Kim, Baker: - L_13
    i_bc = - parameters.a * j / den # v p.51 Kim, Baker: - L_23
    i_bb = - parameters.a^2 * h / den # v p.49 Kim, Baker: - L_22
    i_cc = (response_val - 2 * response_val * p + p^2) / den # X p.50 Kim, Baker: - L_33
    return [i_aa i_ab i_ac; i_ab i_bb i_bc; i_ac i_bc i_cc]::Matrix{Float64}
end