text:
function CreateRandomVacancy(universe::Universe, maxSize::Int64)
    position = [rand(universe.box.boundary[d,1]:universe.box.boundary[d,2]) for d in 1:universe.box.nDimension]
    size = rand(1:maxSize)
    Vacancy(position, size)
end

function CreateRandomInterstitial(universe::Universe, maxSize::Int64)
    position = [rand(universe.box.boundary[d,1]:universe.box.boundary[d,2]) for d in 1:universe.box.nDimension]
    size = rand(1:maxSize)
    direction = sample(universe.box.map.interstitialDirections)
    Interstitial(position, size, direction)
end

function GenerateIntersitialDirections(nDimension::Int64)
    nDimension -= 1
    directions = Vector{Vector{Int64}}()
    direction = Vector{Int64}(undef, nDimension)
    IterIntersitialDirections!(directions, direction, 1, nDimension)
    vcat.([1],directions)
end

function IterIntersitialDirections!(directions::Vector{Vector{Int64}},
    direction::Vector{Int64}, d::Int64, nDimension::Int64)
    for i in [1,-1]
        direction[d] = i
        if d < nDimension
            IterIntersitialDirections!(directions, direction, d+1, nDimension)
        end
        if d == nDimension
            push!(directions, deepcopy(direction))
        end
    end
end

function DeleteCell!(universe::Universe, obj::DefectObj)
    cellIndex = obj.cellIndex
    cell = universe.cellContainer.cells[CartesianIndex(Tuple(cellIndex))]
    filter!(o->o!=obj,cell.objs)
end

function ReCell!(universe::Universe, obj::DefectObj)
    box = universe.box
    cellContainer = universe.cellContainer
    cellIndex = Vector{Int64}(undef, box.nDimension)
    #if universe.thermo.iterTime == 503
    #    println(obj.id,"     ########")
    #end
    #if universe.thermo.iterTime > 0# == 507
    #    println(debug[1] in universe.cellContainer.cells[1,9,7].objs,"  @      ", universe.thermo.iterTime)
    #end
    for d in 1:box.nDimension
        cellIndex[d] = fld(obj.position[d] - box.boundary[d,1], cellContainer.cellCutoff)
        if cellIndex[d] >= size(cellContainer.cells)[d]
            cellIndex[d] -= 2
        end
        if cellIndex[d] < 0
            cellIndex[d] = 0
        end
        cellIndex[d] += 1
    end
    if cellIndex != obj.cellIndex
        if obj.cellIndex != []
            DeleteCell!(universe, obj)
        end
        push!(cellContainer.cells[CartesianIndex(Tuple(cellIndex))].objs, obj)
        obj.cellIndex = cellIndex
    end
    #if universe.thermo.iterTime > 0 # == 507
    #    println(debug[1] in universe.cellContainer.cells[1,9,7].objs,"    @    ", universe.thermo.iterTime)
    #end
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mcke

max_forks_repo_name: permissionx/mcke


text:
module RayTracer

using StaticArrays
using LinearAlgebra

global ray_count = 0

const Vec3 = SArray{Tuple{3}, Float64, 1, 3}
const Point3 = SArray{Tuple{3}, Float64, 1, 3}
const Color = SArray{Tuple{3}, Float64, 1, 3}
import Base.zero
zero(::Type{Point3}) = Point3(0,0,0)
zero(::Type{Color}) = Color(0,0,0)

magnitude(x,y) = sqrt(x^2 + y^2)
magnitude(x,y,z) = sqrt(x^2 + y^2 + z^2)
magnitude(v::Vec3) = magnitude(v...)

unit_vector(v::Vec3) = v / magnitude(v)

near_zero(v) = v.x < 1e-8 && v.y < 1e-8 && v.z < 1e-8

function random_in_unit_disk()
	x,y = randf(-1, 1), randf(-1, 1)
	while magnitude(x,y)^2 >= 1
     	x,y = randf(-1, 1), randf(-1, 1)
	end
	x,y
end

function random_in_unit_sphere()
	x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	while magnitude(x,y,z)^2 >= 1
		x,y,z = randf(-1,1), randf(-1,1), randf(-1,1)
	end
	Point3(x,y,z)
end

random_unit_vector() = unit_vector(random_in_unit_sphere())

function random_in_hemisphere(normal)
    in_unit_sphere = random_in_unit_sphere()
    dot(in_unit_sphere, normal) > 0.0 ? in_unit_sphere : -in_unit_sphere
end

reflect(v, n) =  v - 2dot(v,n)*n

function refract(uv, n, etai_over_etat)
    cos_theta = min(dot(-uv, n), 1.0)
    r_out_perp = etai_over_etat * (uv + cos_theta*n)
    r_out_parallel = -sqrt(abs(1.0 - magnitude(r_out_perp)^2)) * n
    r_out_perp + r_out_parallel
end

struct Ray
	origin::Point3
	direction::Vec3
	udirection::Vec3
	tm::Float64
	Ray(o, d, m) = new(o, d, unit_vector(d), m)
	Ray(o, d) = Ray(o, d, 0)
	Ray() = new(Vec3(Inf, Inf, Inf), Vec3(Inf, Inf, Inf))
end

at(r::Ray, t) = r.origin + t * r.direction

struct Camera
	origin::Point3
	lower_left_corner::Point3
	horizontal::Vec3
	vertical::Vec3
	u::Vec3
	v::Vec3
	w::Vec3
	lens_radius::Float64
	time0::Float64
	time1::Float64
	function Camera(lookfrom, lookat, vup, vfov, aspect_ratio, aperture, focus_dist, _time0=0, _time1=0)
		viewport_height = 2.0 * tan(deg2rad(vfov)/2)
		viewport_width = aspect_ratio * viewport_height

		w = unit_vector(lookfrom - lookat)
		u = unit_vector(cross(vup, w))
		v = cross(w, u)

		origin = lookfrom
		horizontal = focus_dist * viewport_width * u
		vertical = focus_dist * viewport_height * v
		lower_left_corner = origin - horizontal/2 - vertical/2 - focus_dist*w

		lens_radius = aperture / 2
		time0 = _time0
		time1 = _time1
		new(origin, lower_left_corner, horizontal, vertical, u, v, w, lens_radius, time0, time1)
	end
	Camera() = Camera(Point3(0,0,-1), Point3(0,0,0), Vec3(0,1,0), 40, 1, 0, 10)
end

abstract type Material end
abstract type Hitable end

struct Scene
	camera::Camera
	hitables::Vector{Hitable}
	Scene(cam) = new(cam, Vector{Hitable}())
end

randf(fmin, fmax) = fmin + (fmax-fmin)*rand()

get_ray(scene::Scene, s, t) = get_ray(scene.camera, s, t)

function get_ray(cam::Camera, s, t)
	x,y = cam.lens_radius .* random_in_unit_disk()
	offset = cam.u * x + cam.v *y
	Ray(cam.origin + offset, cam.lower_left_corner + s*cam.horizontal + t*cam.vertical - cam.origin - offset, randf(cam.time0, cam.time1))
end

#==
import Base.+
+(c1::Color, c2::Color) = Color(c1.rgb[1]+c2.rgb[1], c1.rgb[2]+c2.rgb[2], c1.rgb[3]+c2.rgb[3])

import Base.*
*(f::Float64, c::Color) = Color(f*c.rgb[1], f*c.rgb[2], f*c.rgb[3])
*(c1::Color, c2::Color) = Color(c1.rgb[1]*c2.rgb[1], c1.rgb[2]*c2.rgb[2], c1.rgb[3]*c2.rgb[3])

import Base./
/(c::Color, f) = Color(c.rgb[1]/f, c.rgb[2]/f, c.rgb[3]/f)
==#

clamp(f, b, t) = f < b ? b : f > t ? t : f

function write(io::IO, c::Color)
	val(rgb) = round(Int, 255clamp(sqrt(rgb), 0, 1))
	println(io, "$(val(c[1])) $(val(c[2])) $(val(c[3]))")
end

struct Hit
	p::Point3
	normal::Vec3
	material::Material
	t::Float64
	front_face::Bool
	function Hit(s::Sphere, t, ray)
		p = at(ray, t)
		outward_normal = (p - s.center) / s.radius
		ff = dot(ray.direction, outward_normal)
		norm = ff < 0 ? outward_normal : -outward_normal
		new(p, norm, s.material, t, ff < 0)
	end
	Hit() = new()
end

function trace(scene::Scene, ray::Ray, t_min, t_max)
	closest_t = t_max
	local struck::Hitable

	for hitable in scene.hitables
		t = trace(hitable, ray, t_min, closest_t)
		if t >= 0
			closest_t = t
			struck = hitable
		end
	end

	if @isdefined struck
		Hit(struck, closest_t, ray)
	else
		Hit()
	end
end

struct Lambertian <: Material
	albedo::Color
	Lambertian(c::Color) = new(c)
	Lambertian(r,g,b) = new(Color(r,g,b))
	Lambertian() = Lambertian(Color(rand()*rand(), rand()*rand(), rand()*rand()))
end

function scatter(l::Lambertian, ray::Ray, hit::Hit)
	scatter_direction = hit.normal + random_unit_vector()
	if near_zero(scatter_direction)
		scatter_direction = hit.normal
	end
	Ray(hit.p, scatter_direction), l.albedo
end

struct Metal <: Material
	albedo::Color
	fuzz::Float64
	Metal(a, f) = new(a, f)
	Metal(r,g,b,f) = Metal(Color(r,g,b), f)
end

function scatter(m::Metal, ray::Ray, hit::Hit)
	reflected = reflect(ray.udirection, hit.normal)
        scattered = Ray(hit.p, reflected + m.fuzz * random_in_unit_sphere())
	dot(scattered.direction, hit.normal) > 0 ? (scattered, m.albedo) : (Ray(), zero(Color))
end

struct Dielectric <: Material
	ir::Float64
end

function reflectance(cosine, ref_idx)
	r0 = ((1-ref_idx) / (1+ref_idx))^2
    r0 + (1-r0)*(1 - cosine)^5
end

function scatter(d::Dielectric, ray::Ray, hit::Hit)
	refraction_ratio = hit.front_face ? (1.0/d.ir) : d.ir

	cos_theta = min(dot(-ray.udirection, hit.normal), 1.0)
	sin_theta = sqrt(1.0 - cos_theta^2)
	cannot_refract = refraction_ratio * sin_theta > 1.0
	direction = if cannot_refract || reflectance(cos_theta, refraction_ratio) > rand()
			reflect(ray.udirection, hit.normal)
		else
			refract(ray.udirection, hit.normal, refraction_ratio)
		end

	Ray(hit.p, direction), Color(1,1,1)
end

abstract type Hitable end

struct Sphere <: Hitable
	center::Point3
	radius::Float64
	material::Material
end

function trace(sphere::Sphere, ray::Ray, t_min, t_max)
	oc = ray.origin - sphere.center
	a = magnitude(ray.direction)^2
	half_b = dot(oc, ray.direction)
	c = magnitude(oc)^2 - sphere.radius^2
	discriminant = half_b^2 - a*c
	if discriminant < 0
		return -1
	end

	sqrtd = sqrt(discriminant)
	root = (-half_b - sqrtd) / a
	if root < t_min || t_max < root
		root = (-half_b + sqrtd) / a
		if root < t_min || t_max < root
			return -1
		end
	end

	root
end

function ray_color(scene::Scene, ray::Ray, depth)::Tuple{Float64, Float64, Float64}
	if depth <= 0
        	return 0,0,0
	end

	hit = trace(scene, ray, 0.001, Inf)
	if isdefined(hit, :material)
		s, a = scatter(hit.material, ray, hit)
		if s.origin.x == Inf
			return 0,0,0
		end
		r,g,b = ray_color(scene, s, depth-1)
		return a[1] * r, a[2] * g, a[3] * b
	end

	t::Float64 = 0.5*(ray.udirection.y + 1.0)
	t1m = 1.0 - t
	t1m + 0.5t, t1m + 0.7t, t1m + t
end

add!(s::Scene, h::Hitable) = push!(s.hitables, h)

function add_random_scene!(scene::Scene)

	add!(scene, Sphere(Point3(0,-1000,0), 1000, Lambertian(Color(0.5, 0.5, 0.5))))

	rand_material(p) = if p < 0.8
				Lambertian()
			elseif p < 0.95
				rf = randf(0.5, 1)
				Metal(Color(rf, rf, rf), 0.5rand())
			else
				Dielectric(1.5)
			end

	for a in -11:10, b in -11:10
		center = Point3(a + 0.9rand(), 0.2, b + 0.9rand())
		if magnitude(center - Point3(4, 0.2, 0)) > 0.9
			add!(scene, Sphere(center, 0.2, rand_material(rand())))
		end
	end

	add!(scene, Sphere(Point3(0, 1, 0), 1.0, Dielectric(1.5)))
	add!(scene, Sphere(Point3(-4, 1, 0), 1.0, Lambertian(0.4,0.2,0.1)))
	add!(scene, Sphere(Point3(4, 1, 0), 1.0, Metal(Color(0.7,0.6,0.5), 0.0)))
end

const Scanline = Vector{Color}

function trace_scanline(world, y, samples, width, height, max_depth)
	scanline = Scanline(undef, width)
	@simd for x in 1:width
		r=g=b=0.0
		for i in 1:samples
			(r,g,b) = (r,g,b) .+ ray_color(world, get_ray(world, (x + rand()) / width, (y + rand()) / height), max_depth)
		end
		@inbounds scanline[x] = Color(r/samples, g/samples, b/samples)
	end
	scanline
end

function write_ppm(io, scanlines, w, h)
	println(io, "P3\n$w $h\n255")
	foreach(scanline->foreach(p->write(io, p), scanline), scanlines[end:-1:1])
end

function main(io; image_width=1200, aspect_ratio=16/9, samples_per_pixel=10, max_depth=50)

	image_height = round(Int, image_width / aspect_ratio)

	world = Scene(Camera(Point3(13.,2.,3.), zero(Point3), Vec3(0,1,0), 20, aspect_ratio, 0.1, 10.0))
	add_random_scene!(world)

	scanlines = Vector{Scanline}(undef, image_height)

	#Threads.@threads
	for y in 1:image_height
		@inbounds scanlines[y] = trace_scanline(world, y, samples_per_pixel, image_width, image_height, max_depth)
    end

	write_ppm(io, scanlines, image_width, image_height)
end

function go()
	open("r.ppm", "w") do io
		@time main(io, image_width=1200, samples_per_pixel=10)
	end
end

function bmark()
	# using BenchmarkTools
	world = Scene(Camera(Point3(13.,2.,3.), zero(Point3), Vec3(0,1,0), 20, 16/9, 0.1, 10.0))
	add_random_scene!(world)
	#@benchmark trace_scanline(world, 400, 10, 1200, 800, 50) setup=(evals=3)
end
#==
includet("main.jl")
world = Scene(Camera(Point3(13.,2.,3.), zero(Point3), Vec3(0,1,0), 20, 16/9, 0.1, 10.0))
@code_warntype ray_color(world, get_ray(world, 0.5, 0.5), 5)
==#

ext: jl


text:
"""
$(TYPEDEF)

Generalized power cone parametrized by powers `α` in the unit simplex and
dimension `d` of the normed variables.

    $(FUNCTIONNAME){T}(α::Vector{T}, d::Int, use_dual::Bool = false)
"""
mutable struct GeneralizedPower{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    α::Vector{T}
    n::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    u_idxs::UnitRange{Int}
    w_idxs::UnitRange{Int}
    z::T
    w2::T
    zw::T
    zwzwi::T
    tempu1::Vector{T}
    tempu2::Vector{T}

    function GeneralizedPower{T}(
        α::Vector{T},
        n::Int;
        use_dual::Bool = false,
        ) where {T <: Real}
        @assert n >= 1
        dim = length(α) + n
        @assert dim >= 3
        @assert all(α_i > 0 for α_i in α)
        @assert sum(α) ≈ 1
        cone = new{T}()
        cone.n = n
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.α = α
        return cone
    end
end

function setup_extra_data!(cone::GeneralizedPower{T}) where {T <: Real}
    m = length(cone.α)
    cone.u_idxs = 1:m
    cone.w_idxs = (m + 1):cone.dim
    cone.tempu1 = zeros(T, m)
    cone.tempu2 = zeros(T, m)
    return cone
end

get_nu(cone::GeneralizedPower) = length(cone.α) + 1

function set_initial_point!(arr::AbstractVector, cone::GeneralizedPower)
    @. arr[cone.u_idxs] = sqrt(1 + cone.α)
    arr[cone.w_idxs] .= 0
    return arr
end

function update_feas(cone::GeneralizedPower{T}) where {T <: Real}
    @assert !cone.feas_updated
    @views u = cone.point[cone.u_idxs]

    if all(>(eps(T)), u)
        cone.z = exp(2 * sum(α_i * log(u_i) for (α_i, u_i) in zip(cone.α, u)))
        @views cone.w2 = sum(abs2, cone.point[cone.w_idxs])
        cone.zw = cone.z - cone.w2
        cone.is_feas = (cone.zw > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::GeneralizedPower{T}) where {T <: Real}
    @views u = cone.dual_point[cone.u_idxs]

    if all(>(eps(T)), u)
        p = exp(2 * sum(α_i * log(u_i / α_i) for (α_i, u_i) in zip(cone.α, u)))
        @views w = cone.dual_point[cone.w_idxs]
        return (p - sum(abs2, w) > eps(T))
    end

    return false
end

function update_grad(cone::GeneralizedPower)
    @assert cone.is_feas
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    @views u = cone.point[u_idxs]
    @views w = cone.point[w_idxs]
    g = cone.grad
    cone.zwzwi = (cone.z + cone.w2) / cone.zw

    @. g[u_idxs] = -(cone.zwzwi * cone.α + 1) / u
    @. g[w_idxs] = 2 * w / cone.zw

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::GeneralizedPower)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    @views u = cone.point[u_idxs]
    H = cone.hess.data
    g = cone.grad
    zw = cone.zw
    aui = cone.tempu1
    auizzwi = cone.tempu2
    @. aui = 2 * cone.α / u
    @. auizzwi = -cone.z * aui / zw
    zzwim1 = -cone.w2 / zw
    zwi = 2 / zw

    @inbounds for j in u_idxs
        auizzwim1 = auizzwi[j] * zzwim1
        for i in 1:j
            H[i, j] = aui[i] * auizzwim1
        end
        H[j, j] -= g[j] / u[j]
    end

    @inbounds for j in w_idxs
        gj = g[j]
        for i in u_idxs
            H[i, j] = auizzwi[i] * gj
        end
        for i in first(w_idxs):j
            H[i, j] = g[i] * gj
        end
        H[j, j] += zwi
    end

    cone.hess_updated = true
    return cone.hess
end

function update_inv_hess(cone::GeneralizedPower{T}) where {T <: Real}
    @assert cone.grad_updated
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    H = cone.inv_hess.data
    α = cone.α
    @views begin
        w = cone.point[w_idxs]
        u = cone.point[u_idxs]
        gw = cone.grad[w_idxs]
        gu = cone.grad[u_idxs]
    end
    z = cone.z
    zw = cone.zw
    w2 = cone.w2
    gww = dot(gw, w)
    guiα = cone.tempu1
    uiα = cone.tempu2
    @. guiα = α / gu
    @. uiα = α / u
    k2 = -dot(uiα, guiα)
    k1 = (1 + w2 / z) / 2 - k2 * gww
    k3 = zw * (-k2 * (gww + 2) / k1 + 1) / (z + w2)
    k4 = -gww / k1
    zw2 = zw / 2

    @inbounds for j in u_idxs
        guiαj = -guiα[j] * k4
        for i in 1:j
            H[i, j] = guiα[i] * guiαj
        end
        H[j, j] -= u[j] / gu[j]
    end

    @inbounds for j in w_idxs
        wj = cone.point[j]
        for i in u_idxs
            H[i, j] = -guiα[i] * wj / k1
        end
        for i in first(w_idxs):j
            H[i, j] = -cone.point[i] * wj * k3
        end
        H[j, j] += zw2
    end

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::GeneralizedPower,
    )
    @assert cone.grad_updated
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    @views u = cone.point[u_idxs]
    @views w = cone.point[w_idxs]
    α = cone.α
    z = cone.z
    zw = cone.zw
    tempu1 = cone.tempu1
    @. tempu1 = 1 + cone.zwzwi * α

    @inbounds for j in 1:size(arr, 2)
        @views begin
            arr_u = arr[u_idxs, j]
            arr_w = arr[w_idxs, j]
            prod_u = prod[u_idxs, j]
            prod_w = prod[w_idxs, j]
        end
        @. prod_u = arr_u / u
        @. prod_w = 2 * arr_w / zw

        dot1 = -4 * dot(α, prod_u) * z / zw
        dot2 = (dot1 + 2 * dot(w, prod_w)) / zw
        dot3 = dot1 - dot2 * z
        prod_u .*= tempu1
        @. prod_u += dot3 * α
        prod_u ./= u
        @. prod_w += dot2 * w
    end

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::GeneralizedPower,
    )
    @assert cone.grad_updated
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    α = cone.α
    @views begin
        w = cone.point[w_idxs]
        u = cone.point[u_idxs]
        gw = cone.grad[w_idxs]
        gu = cone.grad[u_idxs]
    end
    z = cone.z
    zw = cone.zw
    w2 = cone.w2
    gww = dot(gw, w)
    guiα = cone.tempu1
    uiα = cone.tempu2
    @. guiα = α / gu
    @. uiα = α / u
    k2 = -dot(uiα, guiα)
    k1 = (1 + w2 / z) / 2 - k2 * gww
    k3 = zw * (-k2 * (gww + 2) + k1) / (z + w2)

    @inbounds for j in 1:size(arr, 2)
        @views begin
            p = arr[u_idxs, j]
            r = arr[w_idxs, j]
            prod_u = prod[u_idxs, j]
            prod_w = prod[w_idxs, j]
        end
        dot_wr = dot(r, w)
        dot_pu = dot(guiα, p)
        @. prod_w = zw * r / 2 - w * (k3 * dot_wr + dot_pu) / k1
        @. prod_u = p * u / -gu - guiα * (-gww * dot_pu + dot_wr) / k1
    end

    return prod
end

function dder3(cone::GeneralizedPower, dir::AbstractVector)
    @assert cone.grad_updated
    u_idxs = cone.u_idxs
    w_idxs = cone.w_idxs
    @views u = cone.point[u_idxs]
    @views w = cone.point[cone.w_idxs]
    dder3 = cone.dder3
    @views u_dder3 = dder3[u_idxs]
    @views w_dder3 = dder3[cone.w_idxs]
    @views u_dir = dir[u_idxs]
    @views w_dir = dir[cone.w_idxs]
    α = cone.α
    z = cone.z
    zw = cone.zw
    w2 = cone.w2
    zwzwi = cone.zwzwi
    zzwi = 2 * z / zw
    zwi = 2 / zw
    udu = cone.tempu1
    @. udu = u_dir / u

    wwd = 2 * dot(w, w_dir)
    c15 = wwd / zw
    audu = dot(α, udu)
    sumaudu2 = sum(α_i * abs2(udu_i) for (α_i, udu_i) in zip(α, udu))
    c1 = 2 * zwzwi * abs2(audu) + sumaudu2
    c10 = sum(abs2, w_dir) + wwd * c15

    c13 = zzwi * (w2 * c1 - 2 * wwd * zwzwi * audu + c10) / zw
    c14 = zzwi * (2 * audu * w2 - wwd) / zw
    @. u_dder3 = (c13 * α + ((c14 + zwzwi * udu) * α + udu) * udu) / u

    c6 = zwi * (z * (4 * audu * c15 - c1) - c10) / zw
    c7 = zwi * (2 * z * audu - wwd) / zw
    @. w_dder3 = c7 * w_dir + c6 * w

    return cone.dder3
end


text:
## ================================================================================================================================
# Serialization
## ================================================================================================================================

# abstract type JSONManifoldKernelDensity end

import DistributedFactorGraphs: getVariableType

export PackedManifoldKernelDensity

# NOTE, user variables and manifolds will require the same definitions, TODO perhaps add into `@defVariable`
# unusual definitions, but all they need to do is pack and unpack as one-to-one
# this step actually occurs separate from the actual variables or factors (with their own manifolds)
# relies on later use of getManifold to give back the same <:AbstractManifold
# NOTE added to DFG.@defVariable
getVariableType(M::Euclidean{Tuple{N}}) where N = ContinuousEuclid(N)
getVariableType(M::TranslationGroup{Tuple{N}}) where N = ContinuousEuclid(N)

# getVariableType(M::RealCircleGroup) = Circular()
# getVariableType(M::Circle) = error("Circle manifold is deprecated use RealCircleGroup, will come back when we generalize to non-group Riemannian")

##

struct PackedManifoldKernelDensity <: PackedSamplableBelief
  _type::String
  varType::String
  pts::Vector{Vector{Float64}}
  bw::Vector{Float64}
  partial::Vector{Int}
  infoPerCoord::Vector{Float64}
end

# Type converters for MKD
Base.convert(::Type{<:SamplableBelief}, ::Type{<:PackedManifoldKernelDensity}) = ManifoldKernelDensity
Base.convert(::Type{<:PackedSamplableBelief}, ::Type{<:ManifoldKernelDensity}) = PackedManifoldKernelDensity

# Data converters for MKD
function Base.convert(::Type{<:AbstractString},
                      mkd::ManifoldKernelDensity )
  #
  pts = getPoints(mkd)

  packedMKD = PackedManifoldKernelDensity(
    "PackedManifoldKernelDensity",
    # piggy back on InferenceVariable serialization rather than try serialize anything Manifolds.jl
    DFG.typeModuleName(getVariableType(mkd.manifold)),
    [AMP.makeCoordsFromPoint(mkd.manifold, pt) for pt in pts],
    getBW(mkd.belief)[:,1],
    mkd._partial isa Nothing ? collect(1:manifold_dimension(mkd.manifold)) : mkd._partial ,
    mkd.infoPerCoord
  )

  JSON2.write(packedMKD)
end

Base.convert(::Type{<:PackedSamplableBelief}, mkd::ManifoldKernelDensity) = convert(String, mkd)

# make module specific
# good references:
#  https://discourse.julialang.org/t/converting-string-to-datatype-with-meta-parse/33024/2
#  https://discourse.julialang.org/t/is-there-a-way-to-import-modules-with-a-string/15723/6
function Base.convert(::Type{<:ManifoldKernelDensity}, str::AbstractString)
  data = JSON2.read(str, PackedManifoldKernelDensity)

  # piggy back on serialization of InferenceVariable rather than try serialize anything Manifolds.jl
  M = DFG.getTypeFromSerializationModule(data.varType) |> getManifold
  vecP = [AMP.makePointFromCoords(M, pt) for pt in data.pts]
  partial = length(data.partial) == manifold_dimension(M) ? nothing : data.partial

  return manikde!( M, vecP, bw=data.bw, partial=partial, infoPerCoord=data.infoPerCoord )
end

#


text:
module TestQuicktypes
using Test

import Base: ==
import MacroTools

using QuickTypes
using Accessors
import ConstructionBase

# this is a limitation in `MacroTools.splitarg`. If it is fixed
# this test can be removed and our custom splitarg removed.
try
    MacroTools.splitarg(:(x=$nothing))
    println("MacroTools.splitarg can now handle literal nothing in AST.
            Revisit splitarg_no_default workaround")
catch e
    @assert e isa AssertionError
end

# Examples are taken from https://github.com/cstjean/QuickTypes.jl
# Strings are replaced with symbols so that `===` and `==` works
# nicely.

@qstruct Wall(width, height)

@testset "Wall" begin
    x0 = Wall(400, 600)
    x1 = @set x0.width = 300
    @test x1 === Wall(300, 600)
end

abstract type Vehicle end

@qstruct Car{T<:Number, U}(size::T, nwheels::Int=4; manufacturer::U=nothing,
                           brand::String="off-brand") <: Vehicle

@testset "Car" begin
    c = Car(10; manufacturer=("Danone", "Hershey"))
    @test c isa Car
    @test Car <: Vehicle
    c2 = @set c.size = 10
    @test c2.size === 10
    c3 = @set c.manufacturer = 100
    @test c3 === Car(10;manufacturer =100)
end

@qstruct Empty()
@qstruct Cat(name, age::Int, nlegs=4; species=:Siamese)

@testset "Cat" begin
    x0 = Cat(:Tama, 1)

    x1 = @set x0.nlegs = 8
    @test x1 === Cat(:Tama, 1, 8)

    x2 = @set x0.species = :Singapura
    @test x2 === Cat(:Tama, 1, species=:Singapura)
end

@qstruct Pack{T, N}(animals::NTuple{N, T})

@testset "Pack" begin
    x = Pack((Cat(:Tama, 1), Cat(:Pochi, 2)))

    x = @set x.animals[2].nlegs = 5
    @test x.animals == (Cat(:Tama, 1), Cat(:Pochi, 2, 5))
end

abstract type Tree end
@qstruct Maple(qty_syrup::Float64) <: Tree

@testset "Maple" begin
    x0 = Maple(1)
    x1 = @set x0.qty_syrup = 2
    @test x1 === Maple(2)
end

@qmutable Window(height::Float64, width::Float64)
==(x::Window, y::Window) = x.height == y.height && x.width == y.width

@testset "Window" begin
    x0 = Window(1, 2)

    x1 = @set x0.width = 3.0
    @test isequal(x1, Window(1, 3))
    @test x1 == Window(1, 3)
    x2 = @set x0.width = 3
    @test x1 == x2
    @test !(x1 === x2)
end

@qstruct Human(; name=:Alice, height::Float64=170) do
    @assert height > 0    # arbitrary code, executed in the constructor
end

@testset "Human" begin
    x0 = Human()

    x1 = @set x0.name = :Bob
    @test x1 === Human(name=:Bob)

    @test_throws AssertionError @set x0.height = -10
end

@qstruct Group{x}(members::x; _concise_show=true)

@testset "Group" begin
    x = Group((0, 1))
    x = @set x.members[2] = 111
    @test x.members == (0, 111)
end

@qstruct_fp Plane1(nwheels, weight::Number; brand=:zoomba)

@testset "Plane1" begin
    x0 = Plane1(3, 100)
    x1 = @set x0.nwheels = 5
    @test x1 == Plane1(5, 100)
    @test (@set x0.brand = 31).brand === 31
end

# Another way to "support" QuickTypes with type parameters is to use
# QuickTypes.construct.
@qstruct_fp Plane2(nwheels, weight::Number; brand=:zoomba)
ConstructionBase.constructorof(::Type{<: Plane2}) =
    (args...) -> QuickTypes.construct(Plane2, args...)

@testset "Plane2" begin
    x0 = Plane2(3, 100)

    x1 = @set x0.brand = 31
    @test typeof(x1) != typeof(x0)
    @test x1 == Plane2(3, 100, brand=31)
end

end  # module


text:
CachedOperator(::Type{RaggedMatrix},op::Operator;padding::Bool=false) =
    CachedOperator(op,RaggedMatrix{eltype(op)}(undef, 0, Int[]),padding)

## Grow cached operator

function resizedata!(B::CachedOperator{T,RaggedMatrix{T}},::Colon,n::Integer) where T<:Number
    if n > size(B,2)
        throw(ArgumentError("Cannot resize beyond size of operator"))
    end

    if n > B.datasize[2]
        resize!(B.data.cols,n+1)

        if B.padding
            # K is largest colstop.  We get previous largest by looking at precalulated
            # cols
            K = B.datasize[2]==0 ? 0 : B.data.cols[B.datasize[2]+1]-B.data.cols[B.datasize[2]]

            for j = B.datasize[2]+1:n
                K = max(K,colstop(B.op,j))
                B.data.cols[j+1] = B.data.cols[j] + K
            end
        else
            K = B.datasize[2]==0 ? 0 : B.data.m# more robust but slower: maximum(diff(B.data.cols))

            for j = B.datasize[2]+1:n
                cs = colstop(B.op,j)
                K = max(K,cs)
                B.data.cols[j+1] = B.data.cols[j] + cs
            end
        end

        # avoid padding with negative length
        if B.data.cols[n+1] ≤ 0
            return B
        end

        pad!(B.data.data,B.data.cols[n+1]-1)
        B.data.m = K

        jr=B.datasize[2]+1:n
        kr=1:K
        BLAS.axpy!(1.0,view(B.op,kr,jr),view(B.data,kr,jr))

        B.datasize = (K,n)

    end

    B
end

function resizedata!(B::CachedOperator{T,RaggedMatrix{T}},n::Integer,m::Integer) where T<:Number
    resizedata!(B,:,m)
    B.data.m = max(B.data.m,n)   # make sure we have at least n rows

    B
end

## Grow QR

QROperator(R::CachedOperator{T,RaggedMatrix{T}}) where {T} =
    QROperator(R,RaggedMatrix{T}(undef,0,Int[]),0)

function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T},
                                                  MM,DS,RS,BI}},
         ::Colon,col) where {T,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R_cache
    W=QR.H

    if col > MO.datasize[2]
        m = MO.datasize[2]
        resizedata!(MO,:,col+100)  # last rows plus a bunch more

        # apply previous Householders to new columns of R
        for J=1:QR.ncols
            wp=view(W,1:colstop(W,J),J)
            for j=m+1:MO.datasize[2]
                kr=J:J+length(wp)-1
                v=view(MO.data,kr,j)
                dt=BandedMatrices.dot(wp,v)
                LinearAlgebra.axpy!(-2*dt,wp,v)
            end
        end
    end

    if col > size(W,2)
        m=size(W,2)
        resize!(W.cols,col+101)

        for j=m+1:col+100
            cs=colstop(MO.data,j)
            W.cols[j+1]=W.cols[j] + cs-j+1
            W.m=max(W.m,cs-j+1)
        end

        resize!(W.data,W.cols[end]-1)
    end

    for k=QR.ncols+1:col
        cs = colstop(MO.data,k)
        W[1:cs-k+1,k] = view(MO.data,k:cs,k) # diagonal and below
        wp=view(W,1:cs-k+1,k)
        W[1,k]+= flipsign(norm(wp),W[1,k])
        normalize!(wp)

        # scale rows entries
        kr=k:k+length(wp)-1
        for j=k:MO.datasize[2]
            v=view(MO.data,kr,j)
            dt=BandedMatrices.dot(wp,v)
            LinearAlgebra.axpy!(-2*dt,wp,v)
        end
    end
    QR.ncols=col
    QR
end

# BLAS versions, requires BlasFloat

function resizedata!(QR::QROperator{CachedOperator{T,RaggedMatrix{T},
                                       MM,DS,RS,BI}},
::Colon,col) where {T<:BlasFloat,MM,DS,RS,BI}
    if col ≤ QR.ncols
        return QR
    end

    MO=QR.R_cache
    W=QR.H

    sz=sizeof(T)

    w=pointer(W.data)
    R=MO.data
    r=pointer(R.data)

    if col > MO.datasize[2]
        m = MO.datasize[2]
        resizedata!(MO,:,col+100)  # last rows plus a bunch more

        R=MO.data
        r=pointer(R.data)

        # apply previous Householders to new columns of R
        for k=1:QR.ncols
            M=colstop(W,k)  # length of wp
            wp=w+(W.cols[k]-1)*sz  # shift by first index of col J

            for j=m+1:MO.datasize[2]
                v=r+(R.cols[j]+k-2)*sz
                dt=BandedMatrices.dot(M,wp,1,v,1)
                BLAS.axpy!(M,-2*dt,wp,1,v,1)
            end
        end
    end

    if col > size(W,2)
        m=size(W,2)
        resize!(W.cols,col+101)

        for j=m+1:col+100
            cs=colstop(R,j)
            q_len=cs-j+1  # number of entries in j-th column manipulated
            @assert q_len > 0   # Otherwise, diagonal is not included
            W.cols[j+1]=W.cols[j] + q_len
            W.m=max(W.m,q_len)
        end

        resize!(W.data,W.cols[end]-1)
        w=pointer(W.data)
    end

    for k=QR.ncols+1:col
        cs= colstop(R,k)
        M=cs-k+1

        v=r+sz*(R.cols[k]+k-2)    # diagonal entry of R
        wp=w+sz*(W.cols[k]-1)          # k-th column of W
        BLAS.blascopy!(M,v,1,wp,1)
        W.data[W.cols[k]] += flipsign(BLAS.nrm2(M,wp,1),W.data[W.cols[k]])
        normalize!(M,wp)

        # scale rows entries
        for j=k:MO.datasize[2]
            v=r+(R.cols[j]+k-2)*sz
            dt=BandedMatrices.dot(M,wp,1,v,1)
            BLAS.axpy!(M,-2*dt,wp,1,v,1)
        end
    end
    QR.ncols=col
    QR
end

## back substitution
for ArrTyp in (:AbstractVector, :AbstractMatrix)
    @eval function ldiv!(U::UpperTriangular{T, SubArray{T, 2, RaggedMatrix{T}, Tuple{UnitRange{Int}, UnitRange{Int}}, false}},
                             u::$ArrTyp{T}) where T
        n = size(u,1)
        n == size(U,1) || throw(DimensionMismatch())

        V = parent(U)
        @assert parentindices(V)[1][1] == 1
        @assert parentindices(V)[2][1] == 1

        A = parent(V)

        for c=1:size(u,2)
            for k=n:-1:1
                @inbounds ck = A.cols[k]
                @inbounds u[k,c] /= A.data[ck+k-1]
                BLAS.axpy!(-u[k,c], view(A.data,ck:ck+k-2), view(u,1:k-1,c))
            end
        end
        u
    end
end

## Apply Q

function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,RaggedMatrix{T},T},T}},
                      B::AbstractVector{T},tolerance,maxlength) where {RR,T}
    A = parent(Ac)
    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    H=A.QR.H
    M=size(H,1)
    m=length(B)
    Y=pad(B,m+M+10)

    k=1
    yp=view(Y,1:length(B))
    while (k ≤ m || norm(yp) > tolerance )
        if k > maxlength
            @warn "Maximum length $maxlength reached."
            break
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,k+M+50)
            H=A.QR.H
            M=size(H,1)
        end

        cr=colrange(H,k)

        if k+length(cr)-1>length(Y)
            pad!(Y,2*(k+M))
        end

        wp=view(H,cr,k)
        yp=view(Y,k-1+(cr))

        dt=BandedMatrices.dot(wp,yp)
        LinearAlgebra.axpy!(-2*dt,wp,yp)
        k+=1
    end
    resize!(Y,k)  # chop off zeros
end

# BLAS apply Q

function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,RaggedMatrix{T},T},T}},
           B::AbstractVector{T},tolerance,maxlength) where {RR,T<:BlasFloat}
    A = parent(Ac)
    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    H=A.QR.H
    h=pointer(H.data)

    M=size(H,1)
    m=length(B)
    Y=pad(B,m+M+10)

    sz=sizeof(T)

    k=1
    y=pointer(Y)

    yp=y
    while (k ≤ m || BLAS.nrm2(M,yp,1) > tolerance )
        if k > maxlength
            @warn "Maximum length $maxlength reached."
            break
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,k+M+50)
            H=A.QR.H
            h=pointer(H.data)
        end

        M=colstop(H,k)

        if k+M-1>length(Y)
            pad!(Y,2*(k+M))
            y=pointer(Y)
        end

        wp=h + sz*(H.cols[k]-1)
        yp=y+sz*(k-1)

        dt = BandedMatrices.dot(M,wp,1,yp,1)
        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
        k+=1
    end
    resize!(Y,k)  # chop off zeros
end


text:
################################################################################
# Hexagonal
################################################################################

struct Hexagonal <: Style
    links::Int
    io_links::Int
end

#=
          ____          ____
         /    \        /    \
        /      \      /      \
   ____/  (0,1) \____/ (0,3)  \
  /    \        /    \        /
 /      \      /      \      /
/ (0,0)  \____/ (0,2)  \____/
\        /    \        /    \
 \      /      \      /      \
  \____/ (1,1)  \____/ (1,3)  \
  /    \        /    \        /
 /      \      /      \      /
/ (1,0)  \____/ (1,2)  \____/
\        /    \        /
 \      /      \      /
  \____/        \____/
=#

# Encode hexagonal directions using angles.
directions(::Hexagonal) = ("30", "90", "150", "210", "270", "330")
function procrules(style::Hexagonal)

    # Rule to apply if the column address is even.
    even_skeleton = (
        ((0, 1),  "30_out",  "210_in"),
        ((-1, 0), "90_out",  "270_in"),
        ((0, -1), "150_out", "330_in"),
        ((1, -1), "210_out",  "30_in"),
        ((1, 0),  "270_out",  "90_in"),
        ((1, 1),  "330_out", "150_in"),
    )
    even_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in even_skeleton, i in 0:style.links-1
    ])
    even_rule = ConnectionRule(
        even_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> iseven(x[2])
    )

    # Rule to apply if the column address is odd.
    odd_skeleton = (
        ((-1, 1),  "30_out",  "210_in"),
        ((-1, 0),  "90_out",  "270_in"),
        ((-1, -1), "150_out", "330_in"),
        (( 0, -1), "210_out",  "30_in"),
        (( 1, 0),  "270_out",  "90_in"),
        (( 0, 1),  "330_out", "150_in"),
    )
    odd_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in odd_skeleton, i in 0:style.links-1
    ])
    odd_rule = ConnectionRule(
        odd_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> isodd(x[2])
    )

    return (even_rule, odd_rule)
end

function iorules(style::Hexagonal)

    # Connection rule for even columns
    even_skeleton = (
        ( ( 0, 1), "30_out", "210_in"),
        ( ( 1, 1), "330_out", "150_in"),
        ( ( 0,-1), "150_out", "210_in"),
        ( ( 1,-1), "210_out", "30_in"),
    )

    # Processor -> Output_Handler
    even_offsets_1 = squash([
        Offset(a, "$b[$i]", "in[$i]")
        for (a,b,c) in even_skeleton, i in 0:style.io_links-1
    ])
    # Input Handler -> Processor
    even_offsets_2 = squash([
        Offset(a, "out[$i]", "$c[$i]")
        for (a,b,c) in even_skeleton, i in 0:style.io_links-1
    ])
    even_rule = ConnectionRule(
        vcat(even_offsets_1, even_offsets_2),
        source_filter = filter_io,
        dest_filter = filter_io,
        address_filter = x -> iseven(x[2])
    )

    # Connection rule for odd columns
    odd_skeleton = (
        ( (-1, 1), "30_out", "210_in"),
        ( ( 0, 1), "330_out", "150_in"),
        ( (-1,-1), "150_out", "210_in"),
        ( ( 0,-1), "210_out", "30_in"),
    )
    # Processor -> Output_Handler
    odd_offsets_1 = squash([
        Offset(a, "$b[$i]", "in[$i]")
        for (a,b,c) in odd_skeleton, i in 0:style.io_links-1
    ])
    # Input Handler -> Processor
    odd_offsets_2 = squash([
        Offset(a, "out[$i]", "$c[$i]")
        for (a,b,c) in odd_skeleton, i in 0:style.io_links-1
    ])
    odd_rule = ConnectionRule(
        vcat(odd_offsets_1, odd_offsets_2),
        source_filter = filter_io,
        dest_filter = filter_io,
        address_filter = x -> isodd(x[2])
    )

    return (even_rule, odd_rule)
end

function memory_return_rules(::Hexagonal)
    mem2_to_proc = ConnectionRule(
        [
            Offset((-1, 0), "out[0]", "memory_in"),
            Offset((-1, 1), "out[1]", "memory_in")
        ],
        source_filter = filter_memory(2),
        dest_filter = filter_memproc
    )

    mem1_to_proc = ConnectionRule(
        [Offset((-1,0), "out[0]", "memory_in")],
        source_filter = filter_memory(1),
        dest_filter = filter_memory(1),
    )

    return (mem2_to_proc, mem1_to_proc)
end

function memory_request_rules(::Hexagonal)
    proc_to_mem2 = ConnectionRule(
        [
            Offset(( 1, 0), "memory_out", "in[0]"),
            Offset(( 1,-1), "memory_out", "in[1]")
        ],
        source_filter = filter_memproc,
        dest_filter = filter_memory(2)
    )

    proc_to_mem1 = ConnectionRule(
        [Offset(( 1, 0), "memory_out", "in[0]"), ],
        source_filter = filter_memproc,
        dest_filter = filter_memory(1)
    )

    return (proc_to_mem2, proc_to_mem1)
end

# Port orientations
#directions(::Hexagonal) = ("30", "90", "150", "210", "270", "330")
sq() = sqrt(3) / 2
function port_boundaries(::Hexagonal, orientation)
    coords = Dict(
        "30"  => ((   1,      0), ( 1/2, -sq())),
        "90"  => (( 1/2,  -sq()), (-1/2, -sq())),
        "150" => ((-1/2,  -sq()), (  -1,     0)),

        "210" => ((-1/2,   sq()), (   -1,    0)),
        "270" => (( 1/2,   sq()), (-1/2,  sq())),
        "330" => ((   1,      0), ( 1/2,  sq())),
    )

    return coords[orientation]
end

function initial_offset(::Hexagonal, orientation, direction)
    # Select where to start between the start and stop coordinates depending
    # on the direction selected
    starts = Dict(
        ( "30",  Output ) => 0.55,
        ( "90" , Output ) => 0.55,
        ( "150", Output ) => 0.55,
        ( "210",  Input ) => 0.55,
        ( "270" , Input ) => 0.55,
        ( "330" , Input ) => 0.55,

        ( "30",  Input ) => 0.05,
        ( "90" , Input ) => 0.05,
        ( "150", Input ) => 0.05,
        ( "210",  Output ) => 0.05,
        ( "270" , Output ) => 0.05,
        ( "330" , Output ) => 0.05,
    )

    return starts[(orientation, direction)]
end

function cartesian(::Hexagonal, x, y)
    # Determine if this is an odd or even column
    x_out = 2.5 * x
    y_out = 2.5 * y

    if isodd(y)
        x_out = x_out - 1.25
    end

    return (x_out, y_out)
end

hexx(x) = [x+1, x + 1/2, x - 1/2, x - 1, x - 1/2, x + 1/2, x+1]
hexy(y) = [y, y + sq(), y + sq(), y, y - sq(), y - sq(), y]

polygon(::Hexagonal, x, y) = (hexx(x), hexy(y))


text:
import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback, rtmidi_in_cancel_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")
if !(device[].ok)
    @error :rtmidi_open_port "could not found ableton push 2"
    exit(1)
end

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

module OctoStuff

using Octo.Adapters.PostgreSQL
Repo.debug_sql()
Repo.connect(
           adapter = Octo.Adapters.PostgreSQL,
           dbname = "postgresqltest",
           user = "postgres",
       )
Repo.execute([DROP TABLE IF EXISTS :Employee])
Repo.execute(Raw("""
           CREATE TABLE Employee (
               ID SERIAL,
               Name VARCHAR(255),
               Salary FLOAT(8),
               PRIMARY KEY (ID)
           )"""))

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

function midi_callback(msg)
    @info :msg repr(msg)
    if 0x90 == first(msg)
        name = rand(["Jeremy", "Cloris", "John", "Hyunden", "Justin", "Tom"])
        Repo.insert!(Employee, (Name=name,  Salary=10000.50))
    end
end

Repo.query(Employee)

end # module OctoStuff

import .OctoStuff: midi_callback

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        try
            wait(cond)
            msg = codeunits(unsafe_string(r_ecb[].message))
            midi_callback(msg)
        end
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)

#=
midi_callback(msg) = @info :EDIT repr(msg)

rtmidi_in_cancel_callback(device)
close(cond)
istaskdone(cbloop) # true
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)
=#


text:
import Base: step

abstract type Optimizer end

"""
    step(optim)

Run a single optimization step.
"""
step(optim::Optimizer) = throw("unimplemented")

to!(optim::Optimizer, device::Device) = throw("unimplemented")

"""
    zerograd!(optim)

Run zerograd! on every parameters optimized by optim.
"""
function zerograd!(optim::Optimizer)
    for p ∈ optim.params
        zerograd!(p)
    end
end

struct SGD <: Optimizer
    params::Vector{Tensor}
    lr::Real
    momentum::Real
    velocities::Vector{AbstractArray}

    SGD(params, lr) = new(params, lr, 0.0, Vector{Matrix{Real}}(undef, size(params, 1)))
    SGD(params, lr, momentum) = new(params, lr, momentum, [to(zeros(size(p.values)), p.device) for p ∈ params])
end

function step(optim::SGD)
    for (param, prev_vel) ∈ zip(optim.params, optim.velocities)
        vel =  -param.grad
        if optim.momentum > 0
            vel .+= optim.momentum * prev_vel
            prev_vel .= deepcopy(vel)
        end
        param.values += optim.lr * vel
    end
end

function to!(optim::SGD, device::Device)
    for (i, vel) ∈ enumerate(optim.velocities)
        optim.velocities[i] = to(vel, device)
    end
end

max_stars_repo_name: taraspiotr/DeepJulia

max_issues_repo_name: taraspiotr/DeepJulia

max_forks_repo_name: taraspiotr/DeepJulia


text:
#!/usr/bin/env julia

# updates ground truth data for testing

using GraphicalModelLearning

## remove old files
#`rm -rf *.csv`

using Random

include("../common.jl")

for (gm_name, model) in gms
    Random.seed!(0) # fix random number generator
    samples = sample(model, gibbs_test_samples)
    writecsv("$(gm_name)_samples.csv", samples)
end

for (gm_name, model) in gms
    samples = readdlm("$(gm_name)_samples.csv", ",")
    for (form_name, formulation) in formulations
        Random.seed!(0) # fix random number generator
        learned_gm = learn(samples, formulation)
        writecsv("$(gm_name)_$(form_name)_learned.csv", learned_gm)
    end
end


text:
using TrieLab
using Test

@testset "TrieLab.jl" begin
    # Write your tests here.
end


text:
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies: BrickTopology
using Test
using MPI

MPI.Initialized() || MPI.Init()

@testset "2-D Mass matrix" begin
    topology = BrickTopology(MPI.COMM_WORLD, ntuple(_ -> (-1, 1), 2);)

    @testset for N in ((2, 2), (2, 3), (3, 2))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        grid = DiscontinuousSpectralElementGrid(
            topology,
            FloatType = Float64,
            DeviceArray = Array,
            polynomialorder = N,
        )

        @views begin
            ω = grid.ω
            M = grid.vgeo[:, Grids._M, 1]
            MH = grid.vgeo[:, Grids._MH, 1]
            sM = grid.sgeo[Grids._sM, :, :, 1]

            M = reshape(M, Nq[1], Nq[2])
            @test M ≈ [ω[1][i] * ω[2][j] for i in 1:Nq[1], j in 1:Nq[2]]

            MH = reshape(MH, Nq[1], Nq[2])
            @test MH ≈ [ω[1][i] for i in 1:Nq[1], j in 1:Nq[2]]

            sM12 = sM[1:Nfp[1], 1:2]
            @test sM12[:, 1] ≈ [ω[2][j] for j in 1:Nq[2]]
            @test sM12[:, 2] ≈ [ω[2][j] for j in 1:Nq[2]]

            sM34 = sM[1:Nfp[2], 3:4]
            @test sM34[:, 1] ≈ [ω[1][j] for j in 1:Nq[1]]
            @test sM34[:, 2] ≈ [ω[1][j] for j in 1:Nq[1]]
        end
    end
end

@testset "3-D Mass matrix" begin
    topology = BrickTopology(MPI.COMM_WORLD, ntuple(_ -> (-1, 1), 3);)

    @testset for N in ((2, 2, 2), (2, 3, 4), (4, 3, 2), (2, 4, 3))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        grid = DiscontinuousSpectralElementGrid(
            topology,
            FloatType = Float64,
            DeviceArray = Array,
            polynomialorder = N,
        )

        @views begin
            ω = grid.ω
            M = grid.vgeo[:, Grids._M, 1]
            MH = grid.vgeo[:, Grids._MH, 1]
            sM = grid.sgeo[Grids._sM, :, :, 1]

            M = reshape(M, Nq[1], Nq[2], Nq[3])
            @test M ≈ [
                ω[1][i] * ω[2][j] * ω[3][k]
                for i in 1:Nq[1], j in 1:Nq[2], k in 1:Nq[3]
            ]

            MH = reshape(MH, Nq[1], Nq[2], Nq[3])
            @test MH ≈ [
                ω[1][i] * ω[2][j] for i in 1:Nq[1], j in 1:Nq[2], k in 1:Nq[3]
            ]

            sM12 = reshape(sM[1:Nfp[1], 1:2], Nq[2], Nq[3], 2)
            @test sM12[:, :, 1] ≈
                  [ω[2][j] * ω[3][k] for j in 1:Nq[2], k in 1:Nq[3]]
            @test sM12[:, :, 2] ≈
                  [ω[2][j] * ω[3][k] for j in 1:Nq[2], k in 1:Nq[3]]

            sM34 = reshape(sM[1:Nfp[2], 3:4], Nq[1], Nq[3], 2)
            @test sM34[:, :, 1] ≈
                  [ω[1][i] * ω[3][k] for i in 1:Nq[1], k in 1:Nq[3]]
            @test sM34[:, :, 2] ≈
                  [ω[1][i] * ω[3][k] for i in 1:Nq[1], k in 1:Nq[3]]

            sM56 = reshape(sM[1:Nfp[3], 5:6], Nq[1], Nq[2], 2)
            @test sM56[:, :, 1] ≈
                  [ω[1][i] * ω[2][j] for i in 1:Nq[1], j in 1:Nq[2]]
            @test sM56[:, :, 2] ≈
                  [ω[1][i] * ω[2][j] for i in 1:Nq[1], j in 1:Nq[2]]
        end
    end
end


text:
import BinDeps

rm′(f) = (isdir(f) || isfile(f)) && rm(f, recursive = true)

const version = "2.0.5"

folder() = normpath(joinpath(@__FILE__, "../../../deps"))

@static if Sys.isapple()
    uninstall() = map(rm′, filter(x -> !endswith(x, "build.jl"), readdir(folder())))
else
    uninstall() = rm′(joinpath(folder(), "atom"))
end

isinstalled() = Sys.isapple() ?
    isfile(joinpath(folder(), "version")) :
    isdir(joinpath(folder(), "atom"))

function install()
  dir = folder()
  if Sys.isapple()
    _icons = normpath(joinpath(@__FILE__, "../../../res/julia-icns.icns"))
  end
  !isdir(dir) && mkpath(dir)
  uninstall()
  cd(dir) do
    download(x) = run(BinDeps.download_cmd(x, basename(x)))

    download("http://junolab.s3.amazonaws.com/blink/julia.png")

    if Sys.isapple()
      file = "electron-v$version-darwin-x64.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file`)
      rm(file)
      run(`mv Electron.app Julia.app`)
      run(`mv Julia.app/Contents/MacOS/Electron Julia.app/Contents/MacOS/Julia`)
      run(`sed -i.bak 's/Electron/Julia/' Julia.app/Contents/Info.plist`)
      run(`cp $_icons Julia.app/Contents/Resources/electron.icns`)
      run(`touch Julia.app`)  # Apparently this is necessary to tell the OS to double-check for the new icons.
    end

    if Sys.iswindows()
      arch = Int == Int64 ? "x64" : "ia32"
      file = "electron-v$version-win32-$arch.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`7z x $file -oatom`)
      rm(file)
    end

    if Sys.islinux()
      arch = Int == Int64 ? "x64" : "ia32"
      file = "electron-v$version-linux-$arch.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file -d atom`)
      rm(file)
    end
  end
end


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

module LRUExample
# An LRU (Least Recently Used) cache is an associative data structure which
# maintains its contents in an order such that the most recently used item
# is at the beginning of the structure, and the least recently used at the end.
#
# This file specifies two types of LRU caches, both with and without a size
# limit. BoundedLRU has a limit and evicts the LRU item if a new item is added
# after that bound is reached. UnboundedLRU does not have a maximum size, but
# can be used as a basis for more complex LRUs.
#
# LRUs should follow the interfaces for general collections, indexable
# collections, and associative collections.

# The standard implementation of an LRU backs a hash table with a doubly-linked
# list for O(1) operations when reordering on access and eviction. The Julia
# implementation instead backs the table with a Vector. For moderately-sized
# collections, the difference in performance is small, and this implmentation
# is simpler and easier to understand.

import Base.isempty, Base.length, Base.sizeof
import Base.start, Base.next, Base.done
import Base.haskey, Base.get
import Base.setindex!, Base.getindex, Base.delete!, Base.empty!
import Base.show

abstract type LRU{K,V} <: Associative{K,V} end

# Default cache size
const __MAXCACHE = 1024

mutable struct CacheItem{K,V}
    k::K
    v::V
end

mutable struct UnboundedLRU{K,V} <: LRU{K,V}
    ht::Dict
    q::Vector{CacheItem}

    UnboundedLRU{K,V}() where {K,V} = new(Dict(), similar(Array{CacheItem}(1), 0))
end
UnboundedLRU() = UnboundedLRU{Any, Any}()

mutable struct BoundedLRU{K,V} <: LRU{K,V}
    ht::Dict
    q::Vector{CacheItem}
    maxsize::Int

    BoundedLRU{K,V}(m) where {K,V} = new(Dict(), similar(Array{CacheItem}(1), 0), m)
    BoundedLRU{K,V}() where {K,V} = BoundedLRU(__MAXCACHE)
end
BoundedLRU(m) = BoundedLRU{Any, Any}(m)
BoundedLRU() = BoundedLRU{Any, Any}()

## collections ##

isempty(lru::LRU) = isempty(lru.q)
length(lru::LRU) = length(lru.q)

## associative ##

# Should this check count as an access?
haskey(lru::LRU, key) = haskey(lru.ht, key)

get(lru::LRU, key, default) = haskey(lru, key) ? lru[key] : default

function empty!(lru::LRU)
    empty!(lru.ht)
    empty!(lru.q)
end

show(io::IO, lru::UnboundedLRU) = print(io,"UnboundedLRU()")
show(io::IO, lru::BoundedLRU) = print(io,"BoundedLRU($(lru.maxsize))")

## indexable ##

# Method to do the second, slow lookup in the list with early return.
function locate(q, x)
    for i = 1:length(q)
        if q[i] == x
            return i
        end
    end
    error("Item not found.")
end

function getindex(lru::LRU, key)
    item = lru.ht[key]
    idx = locate(lru.q, item)
    splice!(lru.q, idx)
    unshift!(lru.q, item)
    item.v
end

function setindex!(lru::LRU, v, key)
    if haskey(lru, key)
        item = lru.ht[key]
        idx = locate(lru.q, item)
        item.v = v
        splice!(lru.q, idx)
    else
        item = CacheItem(key, v)
        lru.ht[key] = item
    end
    unshift!(lru.q, item)
end

# Eviction
function setindex!{V,K}(lru::BoundedLRU, v::V, key::K)
    invoke(setindex!, Tuple{LRU,V,K}, lru, v, key)
    nrm = length(lru) - lru.maxsize
    for i in 1:nrm
        rm = pop!(lru.q)
        delete!(lru.ht, rm.k)
    end
end

## associative ##

function delete!(lru::LRU, key)
    item = lru.ht[key]
    idx = locate(lru.q, item)
    delete!(lru.ht, key)
    delete!(lru.q, idx)
end

end # module


text:
@doc_str "Calculates x for `0 = a*x^2+b*x+c`, arguments types can be defined
in function definitions"
function quadratic2(a::Float64, b::Float64, c::Float64)
    quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a
    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
    r1, r2
end


text:
#################################################################
# The simulate function runs a simulation of a pomdp and returns
# the accumulated reward. Different behavior is defined by
# creating Simulator types
#################################################################

# Base type for an object defining how a simulation should be carried out
abstract Simulator

# runs a simulation using the specified policy and returns the accumulated reward
@pomdp_func simulate(simulator::Simulator, pomdp::POMDP, policy::Policy, updater::BeliefUpdater, initial_belief::Belief)

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module Utils
include("./problem_set.jl")

import Distributions: Normal, Uniform
import Statistics: quantile
using Random
using Statistics
using DifferentialEquations
using .ProblemSet

Random.seed!(1234)

function transpose_vector(v::Vector{Vector{T}})::Vector{Vector{T}} where T
    [[v[i][j] for i in 1:length(v)] for j in 1:length(v[1])]
end

"""
From bounds, generates a random startign point for
the optimization procedure.
"""
function rand_guess(bounds::Vector)::Vector
    [rand(Uniform(bounds[1][i], bounds[2][i]))
            for i in 1:length(bounds[1])]
end

"""
Adds noise to a Vector of Vector
Used for adding simulated noise to
the states of a system.
Assumes that the values can be negative.
"""
function add_noise(data::Vector{Vector{T}},
                    percent::T,positive::Bool=false)::Vector{Vector{T}} where T
    if percent == 0.0
        return data
    end

    noise_data = deepcopy(data)
    epsilon_arr = [0.01*mean(getindex.(data,i)) for i in 1:length(data[1])]
    sigma = zero(T)
    @inbounds for i in 1:length(noise_data)
        for j in 1:length(noise_data[1])
            sigma = abs(percent*noise_data[i][j] + epsilon_arr[j])
            d = Normal(0,sigma)
            noise_data[i][j] += rand(d)
            noise_data[i][j] = positive ? abs(noise_data[i][j]) : noise_data[i][j]
        end
    end
    return noise_data
end

"""
Adds noise to a Vector of T
Used for generating new reasonable
initial values to an IVP. Assumes
that the values should be positive.
"""
function add_noise(data::Vector{T},
                    percent::T,positive::Bool=true)::Vector{T} where T
    if percent == 0.0
        return data
    end

    noise_data = deepcopy(data)
    epsilon = 0.01 * mean(data)
    sigma = zero(T)
    @inbounds for i in 1:length(noise_data)
        sigma = abs(percent*noise_data[i] + epsilon)
        d = Normal(0,sigma)
        noise_data[i] += rand(d)
        noise_data[i] = positive ? abs(noise_data[i]) : noise_data[i]
    end
    return noise_data
end

function fae(actual::Array, forecast::Array)::Number
    n = length(actual)
    frac = 1/n
    e = abs.(actual-forecast)
    normalize = abs.(actual)+abs.(forecast)
    stat = frac*sum(2*e/normalize)
end

function smape(actual::Array, forecast::Array)::Float16
    100*fae(actual,forecast)
end

function mape(actual::Array, forecast::Array)::Float32
    n = length(actual)
    frac = 1/n
    e = abs.(actual-forecast)
    normalize = abs.(actual)
    stat = 100*frac*sum(e/normalize)
end

function mrae(actual::Array, forecast::Array)::Float32
    n = length(actual)
    e = abs.(actual-forecast)
    stat = (1/n)*sum(e./abs.(actual))
end

"""
Normalized Mean Square Error
"""
function nmse(data::Vector{T}, data_est::Vector{T})::T where T
    normalizer = abs2(maximum(data) - minimum(data))
    res = mean(abs2.(data_est-data))/normalizer
    return res
end

function filter_outlier(arr::Array; p::Float64=2)::Array
    q = quantile(arr)
    iqr = q[3] - q[1]
    maxValue = q[3] + iqr * p
    minValue = q[1] - iqr * p

    return filter(x -> (x >= minValue) && (x <= maxValue), arr)
end

"""
A Trace object is intended to hold
multiple pairs of time and eval.
Therefore, the attributes time and
eval shall be array of arrays, that,
when zipped, recreate the pairs
mentioned above.
"""
struct Trace
    time::Array{<:Array{<:AbstractFloat},1}
    eval::Array{<:Array{<:AbstractFloat},1}
end

function make_trace(trace)::Trace
    time = []
    eval = []
    @inbounds for i in 1:length(trace)
        append!(time, parse(Float64, split(string(trace[i]))[end]))
        append!(eval, parse(Float64, split(string(trace[i]))[2]))
    end
    return Trace(time, eval)
end

"""
Given a trace, make sure that
both time and eval arrays have the same size.
This is done by repeating the last registered
value in each array.
"""
function fill_trace(t::Trace)::Trace
    max_t_array_len = maximum(map(x -> length(x), t.time))
    resized_t = map( x -> append!(x, fill(x[end], max_t_array_len - length(x))), t.time)

    max_eval_array_len = maximum(map(x -> length(x), t.eval))
    resized_eval = map( x -> append!(x, fill(x[end], max_eval_array_len - length(x))), t.eval)

    return Trace(resized_t, resized_eval)
end

function scale_eval(t::Trace)::Trace
    scaled_eval = [e/e[end] for e in t.eval]
    return Trace(t.time, scaled_eval)
end

"""
Given a Trace
select the best trace (the one with the least amount of
time spent to get to the optimal solution) and the worst
trace (the one with the most)
"""
function get_range_traces(trace::Trace)::Trace
    trace_times = [t[end] for t in trace.time]
    (_,min_idx) = findmin(trace_times)
    best = Trace(trace.time[min_idx], trace.eval[min_idx])

    (_,max_idx) = findmax(trace_times)
    worst = Trace(trace.time[max_idx], trace.eval[max_idx])

    if length(trace_times)%2 == 0
        append!(trace_times, 0.0)
        median_idx = findfirst(x -> x == median(trace_times), trace_times)
        median_idx -= 1
    else
        median_idx = findfirst(x -> x == median(trace_times), trace_times)
    end
    med = Trace(trace.time[median_idx], trace.eval[median_idx])

    return best, med, worst
end

function success_rate(x::T)::T where T
    return 2/(1+exp100(x))
end

function step_success_rate(x::T)::Int64 where T
    if x < 0.125
       return 1
    else
        return 0
    end
end

"""
Absolute Log Difference
Since the data is required to be positive,
find the minimum value from both actual and
forecast data and lift up all values to avoid
zeros and negative numbers.
"""
function abs_log_diff(data::Vector{T}, data_est::Vector{T})::T where T
    min_val = minimum(vcat(data,data_est))
    if min_val < 0.0
        data .= data .- (min_val-0.1)
        data_est .= data_est .- (min_val-0.1)
    end
    return abs.(log.(data)-log.(data_est))
end

function diff_calc(problem::DEProblem,
                    estimated::Array{<:AbstractFloat,1},
                    u0::Array{<:AbstractFloat},
                    f::Function)::Float64
    f = problem.fun
    t = problem.t
    known = problem.phi

    de_prob = ODEProblem(f,u0,(t[1],t[end]),known)
    de_sol = DifferentialEquations.solve(de_prob, Tsit5(), saveat=t)
    data = reduce(hcat, de_sol.u)

    de_prob_est = ODEProblem(f,u0,(t[1],t[end]),estimated)
    de_sol_est = DifferentialEquations.solve(de_prob_est, Tsit5(), saveat=t)
    data_est = reduce(hcat, de_sol_est.u)

    fdiff = mean(f(data, data_est))

    return fdiff
end

function max_diff_states(problem::DEProblem,
                        estimated::Array{<:AbstractFloat,1},
                        variance::Float64)::Float64

    len = round(Int64, 1e3*(problem.t[end]-problem.t[1]))
    _t = range(problem.t[1], stop=problem.t[end], length=len)
    _problem = DEProblem(problem.fun, problem.phi,
                        problem.bounds, problem.data, _t)

    u0 = problem.data[1]
    reps = 5
    #u0_arr = collect(Map(x -> abs.(add_noise(x,variance))),eduction(u0 for _ in 1:reps))
    u0_arr = [abs.(add_noise(u0,variance)) for _ in 1:reps]

    median_max_error = mean([diff_calc(_problem,estimated,x) for x in u0_arr])

    return median_max_error
end

function log_scale_time(time::Array{T})::Array{T} where T<:AbstractFloat
	_arr = deepcopy(arr)
	while any(_arr .< 1.0)
		_arr .= _arr .* 10.0
    end
	return _arr
end

end

max_stars_repo_name: atgmello/DEParamEst

max_issues_repo_name: atgmello/DEParamEst

max_forks_repo_name: atgmello/DEParamEst


text:
export regex_escape
## https://github.com/JuliaLang/julia/pull/29643/commits/dfb865385edf19b681bc0936028af23b1f282b1d
"""
        regex_escape(s::AbstractString)

regular expression metacharacters are escaped along with whitespace.
"""
function regex_escape(s)
    res = replace(escape_string(string(s)), r"([()[\]{}?*+\-|^\$.&~#\s=!<>|:])" => s"\\\1")
    replace(res, "\0" => "\\0")
end
include("unicode.jl")
"""
`ValueMatcher` match value at point `c` iif [`ismatch`](@ref)`(c, parser)`.
A `ValueMatcher{T}=NIndexParser{1,T}` and has `state_type` `MatchState`.

See [`AnyValue`](@ref), [`ValueIn`](@ref), and [`ValueNotIn`](@ref).
"""
abstract type ValueMatcher{T} <: NIndexParser{1,T} end
regex_prefix(x::ValueMatcher) = "["
regex_suffix(x::ValueMatcher) = "]"
result_type(::Type{ValueMatcher{T}}) where T = T
_regex_string(x::ValueMatcher) = regex_inner(x)

export regex_string
regex_string(x) = _regex_string(x)

_regex_string(x::AbstractString) = x # regex_escape(x)
_regex_string(x::Char) = ("$x") ##x == '\\' ? "\\\\" : "$x" ## for [] char ranges

_regex_string(x) = "$x"
_regex_string(x::Nothing) = ""
_regex_string(x::ConstantParser) = _regex_string(x.parser)
_regex_string(x::Union{Vector,Set}) = join(_regex_string.(x),",")
_regex_string(x::Union{Vector{<:AbstractChar},Set{<:AbstractChar}}) = join(_regex_string.(x))
_regex_string(x::StepRange) =
    if x.start == x.stop
        x.start
    else
        x.start*"-"*x.stop
    end
_regex_string(x::Tuple) = join([_regex_string(s) for s in x])
_regex_string(x::Function) = "$x(...)"
# _regex_string(x::ValueIn) = ( x.pcre =="" ? _regex_string(x.sets) : x.pcre )

export AnyValue, AnyChar
"""
    AnyValue(T=Char)

Parser matching exactly one `x::T`, returning the value.
```jldoctest
julia> AnyChar()
. AnyValue
::Char
```

"""
struct AnyValue{T} <: ValueMatcher{T} end

AnyValue(T::Type) = AnyValue{T}()
"""
    AnyChar() = AnyValue(Char)
"""
AnyChar() = AnyValue(Char)
@deprecate AnyValue() AnyChar()
regex_inner(x::AnyValue{Char}) = "."
regex_inner(x::AnyValue{T}) where T = "(.::$T)"
regex_prefix(x::AnyValue) = ""
regex_suffix(x::AnyValue) = ""

"""
    iterate_state(parser::ValueMatcher, sequence, till, posi, next_i, state::Nothing)

When implementing a `Custom<:ValueMatcher` it suffices to provide a method [`CombinedParsers._ismatch`](@ref)`(c, parser::Custom)`.
"""
@inline function iterate_state(parser::ValueMatcher, sequence, till, posi, next_i, state::Nothing)
    next_i>till && return nothing
    @inbounds c,ni = sequence[next_i], _nextind(sequence, next_i)
    !ismatch(c,parser) && return nothing
    return ni, MatchState()
end

"""
    _ismatch(x::Char, set::Union{Tuple,Vector})::Bool

Return `_ismatch(x,set...)`.
"""
function _ismatch(x, set::Union{Tuple,Vector})::Bool
    return _ismatch(x,set...)
end

"""
    _ismatch(x, f, r1, r...)

Check if `x` matches any of the options `f, r1,r...`:
If `ismatch(x,f)` return `true`,
otherwise return `_ismatch(x, r1, r...)`.
"""
function _ismatch(x, f, r1, r...)::Bool
    ismatch(x,f) && return true
    return _ismatch(x, r1, r...)
end

"""
    _ismatch(x)

returns `false` (out of options)
"""
function _ismatch(x)::Bool
    return false
end

"""
    _ismatch(x, p)

returns `x==p`
"""
function _ismatch(c,p)::Bool
    c==p
end

"""
    ismatch(c,p)

returns [`_ismatch`](@ref)`(c, p)`
"""
function ismatch(c,p)::Bool
    _ismatch(c, p)
end

"""
    _ismatch(c,p::Function)

returns `p(c)`
"""
_ismatch(c,p::Function)::Bool = p(c)::Bool
"""
    _ismatch(c,p::AnyValue)

`true`
"""
_ismatch(c,p::AnyValue)::Bool = true
"""
    _ismatch(c,p::Union{StepRange,Set})

returns `c in p`
"""
_ismatch(c,p::Union{StepRange,Set})::Bool = c in p

export ValueIn, CharIn
"""
    ValueIn(x)

Parser matching exactly one element `c` (character) in a sequence, iif [`_ismatch`](@ref)`(c,x)`.

```jldoctest
julia> a_z = ValueIn('a':'z')
[a-z] ValueIn
::Char

julia> parse(a_z, "a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> ac = CharIn("ac")
[ac] ValueIn
::Char

julia> parse(ac, "c")
'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)

julia> l = CharIn(islowercase)
[islowercase(...)] ValueIn
::Char

julia> parse(l, "c")
'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)

```
"""
@auto_hash_equals struct ValueIn{T,S} <: ValueMatcher{T}
    pcre::String
    sets::S
    function ValueIn{T}(pcre::AbstractString, x_...) where T
        label, x = flatten_valuepatterns(x_...)
        new{T,typeof(x)}(pcre == "" ? label : pcre,x)
    end
    function ValueIn(pcre::AbstractString, x_...)
        ## @show pcre
        label, x = flatten_valuepatterns(x_...)
        T = valuepattern_type(x)
        ## @show pcre == "" ? label : pcre
        new{T,typeof(x)}(pcre == "" ? label : pcre,x)
    end
end
@inline _ismatch(c,p::ValueIn)::Bool = _ismatch(c,p.sets)
regex_inner(x::ValueIn) = ( x.pcre == "" ? _regex_string(x.sets) : x.pcre )

ValueIn{T}(x_...) where T = ValueIn{T}("",x_...)
ValueIn(x_...) = ValueIn("",x_...)
ValueIn{Char}(chars::AbstractString) = isempty(chars) ? Never() : ValueIn{Char}(chars,chars...)
ValueIn{T}(label::AbstractString=constructor_name(T)) where T =
    ValueIn{T}(label, x-> x isa T)

parser(x::UnicodeClass) = ValueIn{Char}(x)

@deprecate ValueIn(unicode_classes::Symbol...) ValueIn{Char}(UnicodeClass(unicode_classes...))

export ValueNotIn, CharNotIn
"""
    ValueNotIn{T}(label::AbstractString, x)

Parser matching exactly one element (character) in a sequence, iif not in `x`.

    ValueNotIn([label::AbstractString="", ]x...)
    ValueNotIn{T}([label::AbstractString="", ]x...)

Flattens `x` with [`CombinedParsers.flatten_valuepatterns`](@ref), and tries to infer `T` if not provided.

```jldoctest
julia> a_z = CharNotIn('a':'z')
[^a-z] ValueNotIn
::Char

julia> ac = CharNotIn("ca")
[^ca] ValueNotIn
::Char
```

Respects boolean logic:
```jldoctest
julia> CharNotIn(CharNotIn("ab"))("a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```

Respects boolean logic:
```jldoctest
julia> CharIn(CharIn("ab"))("a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> CharIn(CharNotIn("bc"))("a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> parse(CharNotIn(CharIn("bc")), "a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```

"""
@auto_hash_equals struct ValueNotIn{T,S} <: ValueMatcher{T}
    pcre::String
    sets::S
    function ValueNotIn{T}(pcre::String, x_...) where T
        label, x = flatten_valuepatterns(x_...)
        new{T,typeof(x)}(pcre == "" ? label : pcre,x)
    end
    function ValueNotIn(pcre::String, x_...)
        label, x = flatten_valuepatterns(x_...)
        T = valuepattern_type(x)
        new{T,typeof(x)}(pcre == "" ? label : pcre,x)
    end
end
# result_type(::Type{T}) where T = T
@inline _ismatch(c,p::ValueNotIn)::Bool = !_ismatch(c,p.sets)
regex_inner(x::ValueNotIn) = "^"*( x.pcre =="" ? _regex_string(x.sets) : x.pcre )

ValueNotIn{T}(x_...) where T = ValueNotIn{T}("",x_...)
ValueNotIn{T}(label::AbstractString=constructor_name(T)) where T =
    ValueNotIn{T}(label, x-> x isa T)
ValueNotIn(x_...) = ValueNotIn("",x_...)
ValueNotIn{Char}(chars::String) = ValueNotIn{Char}(chars,chars)
# ValueNotIn(chars::StepRange) =
#     ValueNotIn{eltype(chars)}("$(chars.start)-$(chars.stop)",chars)
# ValueNotIn(pcre::String,x::ConstantParser{Char}) =
#     ValueNotIn{Char}(pcre,x.parser)
@deprecate ValueNotIn(unicode_classes::Symbol...) ValueNotIn(UnicodeClass(unicode_classes...))
# ValueIn(x::Tuple{<:ValueNotIn}) = x[1]

_push!(charset::Nothing, ::Nothing) = nothing
_push!(charset::Nothing, x) = x
_push!(charset, x) = Set{Union{typeof(charset),typeof(x)}}(tuple(charset,x))
_push!(charset::Set, x) = push!(charset,x)
_push!(charset::Vector, x) = push!(charset,x)
_push!(charset::Nothing, x::Union{<:Function,<:UnicodeClass,<:ValueNotIn}) = Any[x]

# This case might occur during construction from Regexp brackets.
flatten_valuepatterns!(label,charset,otherstuff,x::Nothing) =
    label, charset, otherstuff

flatten_valuepatterns!(label,charset,otherstuff,x) =
    label*_regex_string(x), _push!(charset,x), otherstuff

flatten_valuepatterns!(label,charset,otherstuff,x::Union{<:Function,<:UnicodeClass,<:ValueNotIn}) =
    label*_regex_string(x), charset, _push!(otherstuff,x)

flatten_valuepatterns!(label,charset,otherstuff,x::ConstantParser) =
    flatten_valuepatterns!(label,charset,otherstuff, x.parser)

# e.g. ValueIn{UnicodeClass}
flatten_valuepatterns!(label,charset,otherstuff,c::ValueIn) =
    label * regex_inner(c), flatten_valuepatterns!(label,charset,otherstuff, c.sets)[2:3]...

ElementIterators = Union{<:Vector,<:Tuple,<:StepRange,<:Set,<:AbstractString,<:AbstractSet}
function flatten_valuepatterns!(label,charset,otherstuff,x::ElementIterators)
    for e in x
        label2, charset, otherstuff = flatten_valuepatterns!(label, charset, otherstuff,e)
    end
    label * _regex_string(x), charset, otherstuff
end

"""
    flatten_valuepatterns(x...)

Used in `ValueMatcher` constructors.

Heuristic is roughly:
- collect `ElementIterators` in a `Set`
- collect everything else in a `Tuple` (`Function`s etc.)
- in the process the `label` is concatenated
- return all that was collected as `Tuple{String, <:Set, <:Tuple}` or `Tuple{String, <:Set}` or `Tuple{String, <:Tuple}`.
"""
function flatten_valuepatterns(x...)
    label, charset,otherstuff = flatten_valuepatterns!("",nothing,nothing,x)
    label, if otherstuff===nothing
        charset === nothing ? nothing : charset
    elseif charset===nothing
        tuple(otherstuff...)
    elseif isempty(otherstuff)
        charset
    else
        tuple(charset,otherstuff...)
    end
end

valuepattern_type(x) =
    if x isa Tuple
        valuepattern_type(x[1])
    elseif x isa AbstractSet
        eltype(x)
    else
        typeof(x) # error()
    end

"""
    CharIn(a...; kw...) = ValueIn{Char}(a...; kw...)
"""
CharIn(a...; kw...) = ValueIn{Char}(a...; kw...)

"""
    CharNotIn(a...; kw...) = ValueNotIn{Char}(a...; kw...)
"""
CharNotIn(a...; kw...) = ValueNotIn{Char}(a...; kw...)


text:
using Base.Test

try
  using ZAdjoin
catch
  info("Testing as a module. Not installed as a package.")
  include("../src/ZAdjoin.jl")
  using ZAdjoin
end

# ZAdjoin.cyclotomic(n) is recursive, requiring cyclotomic
# polynomials for all proper divisors of n (including 1.)
@test ZAdjoin.cyclotomic(24)==[1,0,0,0,-1,0,0,0,1] # 1-x^4+x^8
@test ZAdjoin.cyclotomic(21)==[1,-1,0,1,-1,0,1,0,-1,1,0,-1,1] # 1-x+x^3-x^4+x^6-x^8+x^9-x^11+x^12

nothing

max_stars_repo_name: WilCrofter/ZAdjoin

max_issues_repo_name: WilCrofter/ZAdjoin

max_forks_repo_name: WilCrofter/ZAdjoin


text:
module EvolutionaryPlus

include("DE.jl")
export DE

include("CMA_ES.jl")
export CMA_ES

include("GPs.jl")
export GPs

end # module


text:
macro import_fields(obj)
  quote
    eval(assign_fields($obj))
  end
end

function assign_fields(obj)
  expr = Expr(:block, [], Any)
  for field in fieldnames(obj)
    push!(expr.args, :($field = ($obj).$field))
  end
  return expr
end

export @import_fields

max_stars_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_issues_repo_name: hunghtd/PowerFlowAnalysis-tutorial

max_forks_repo_name: hunghtd/PowerFlowAnalysis-tutorial


text:
using Base: PkgId, loaded_modules, package_callbacks
using Base.Meta: isexpr

export @require

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

loaded(pkg) = haskey(Base.loaded_modules, pkg)

const notified_pkgs = [Base.PkgId(UUID(0x295af30fe4ad537b898300126c2a3abe), "Revise")]

const _callbacks = Dict{PkgId, Vector{Function}}()
callbacks(pkg) = get!(Vector{Function}, _callbacks, pkg)

listenpkg(@nospecialize(f), pkg) =
  loaded(pkg) ? f() : push!(callbacks(pkg), f)

function loadpkg(pkg::Base.PkgId)
  if haskey(_callbacks, pkg)
    fs = _callbacks[pkg]
    delete!(_callbacks, pkg)
    foreach(Base.invokelatest, fs)
  end
end

function withpath(@nospecialize(f), path)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end

function err(@nospecialize(f), listener, mod)
  try
    f()
  catch e
    @warn """
      Error requiring $mod from $listener:
      $(sprint(showerror, e, catch_backtrace()))
      """
  end
end

function parsepkg(ex)
  isexpr(ex, :(=)) || @goto fail
  mod, id = ex.args
  (mod isa Symbol && id isa String) || @goto fail
  return id, String(mod)
  @label fail
  error("Requires syntax is: `@require Pkg=\"uuid\"`")
end

function withnotifications(@nospecialize(args...))
  for id in notified_pkgs
    if loaded(id)
      mod = Base.root_module(id)
      if isdefined(mod, :add_require)
        add_require = getfield(mod, :add_require)::Function
        add_require(args...)
      end
    end
  end
  return nothing
end

macro require(pkg, expr)
  pkg isa Symbol &&
    return Expr(:macrocall, Symbol("@warn"), __source__,
                "Requires now needs a UUID; please see the readme for changes in 0.7.")
  id, modname = parsepkg(pkg)
  pkg = :(Base.PkgId(Base.UUID($id), $modname))
  quote
    if !isprecompiling()
      listenpkg($pkg) do
        $withnotifications($(string(__source__.file)), $__module__, $id, $modname, $(esc(Expr(:quote, expr))))
        withpath($(string(__source__.file))) do
          err($__module__, $modname) do
            $(esc(:(eval($(Expr(:quote, Expr(:block,
                                            :(const $(Symbol(modname)) = Base.require($pkg)),
                                            expr)))))))
          end
        end
      end
    end
  end
end


text:
using Test
using Documenter

using Indexes

@testset "Indexes" begin
    # TODO
    # @test GenomicFeatures.Indexes.Tabix === GenomicFeatures.Indexes.Tabix

# Include doctests.
DocMeta.setdocmeta!(Indexes, :DocTestSetup, :(using Indexes); recursive=true)
doctest(Indexes; manual = false)

end


text:
using MappingRobots
using Base.Dates

timeout = length(ARGS) > 0 ? Second(parse(Int, ARGS[1])) : Second(30)

robot = Robot(Brick("/home/pi/Ev3"))
local_map = run_mapping(robot, timeout=timeout)

open("scan_points.csv", "w") do f
    writedlm(f, local_map.points, ',')
end

open("path.csv", "w") do f
    writedlm(f, [[t.offset[1] for t in local_map.path] [t.offset[2] for t in local_map.path]], ',')
end


text:
export flat_recon

function flat_recon(img, filename, step = 9, Save = True)
    img_cut = img

    (x, y, z) = size(img_cut)

    img_re = zeros(z * step, y, Int(floor(x / step)))

    @showprogress "Image Reconstruction for Record %s..." for i = 1:z
        for j = 1:Int(floor(x / step))
            #for h in 0:step-1
            img_re[1+(i-1)*step:i*step, :, j] = img_cut[(j-1)*step+1:j*step, :, i]
        end
    end

    if Save
        img_save(img_re, "/home/jchang/image/result/", (@sprintf("%s-.tif", filename)))
    end
    return img_re
end


text:
# PROJECT_POS    Project a point onto the nonnegative orthant.
#
#   project_pos(v) is the positive part of v.
function project_pos(v)
    x = max(v, 0)
    return x
end


text:
################################################################################
#                                   SiO₂                                       #
################################################################################
export SiO₂

"""
This code creates a symbolic representation of the Sellmeier Equation model
for the index of refraction of amorphous SiO₂.Equation form is based on:
     Kitamura, et al.
     "Optical constants of silica glass from extreme ultraviolet to far
     infrared at near room temperature." Applied optics 46.33 (2007): 8118-8133.
which references
    Malitson, “Interspecimen comparison of the refractive index of fused
    silica,” J. Opt. Soc. Am.55,1205–1209 (1965)
and has been validated from 0.21-6.7 μm (free space wavelength).
The thermo-optic coefficient (for 300K) is from the literature, but I forgot
the source.
The symbolic index model and its derivatives in turn are used to generate
numerical functions for the SiO₂ index, group index and GVD as a function
of temperature and wavelength.
Variable units are lm in [um] and T in [deg C]
"""

p_n²_SiO₂ = (
    A₀ = 1,
    B₁ = 0.6961663,
    C₁ = (0.0684043)^2,     #                           [μm²]
    B₂ = 0.4079426,
    C₂ = (0.1162414)^2,     #                           [μm²]
    B₃ = 0.8974794,
    C₃ = (9.896161)^2,      #                           [μm²]
    dn_dT = 6.1e-6,         # thermo-optic coefficient  [K⁻¹]
    T₀ = 20,                # reference temperature     [°C]
)

function make_SiO₂(;p_n²=p_n²_SiO₂)
	@variables λ, T
	n = sqrt(n²_sym_fmt1( λ ; p_n²...)) + p_n².dn_dT  *  ( T - p_n².T₀  )
	n² = n^2
	ε 	= diagm([n², n², n²])
	ng = ng_model(n,λ)
	gvd = gvd_model(n,λ)
	models = Dict([
		:n		=>	n,
		:ng		=>	ng,
		:gvd	=>	gvd,
		:ε 		=> 	diagm([n², n², n²]),
	])
	defaults =	Dict([
		:λ		=>		0.8,		# μm
		:T		=>		p_n².T₀,	# °C

	])
	Material(models, defaults, :SiO₂, colorant"aqua")
end

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

SiO₂ = make_SiO₂()
# SiO₂_num =

# struct NumMaterial{T,F1,F2,F3,TC} where {T, F1<:Function, F2<:Function, F3<:Function, TC}
# 	ε::T
# 	fε::F1
# 	fng::F2
# 	fgvd::F3
# 	fχ⁽²⁾::F4
# 	name::Symbol
# 	color::TC
# end


text:
using Test

using PotentialFlow
using LinearAlgebra: norm, ⋅

@test isempty(detect_ambiguities(Vortex))

# Copied verbatim from ssfrr/TestSetExtensions.jl
# The package is not 0.7/1.0 compatible yet, but
# this function is too useful
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(f->endswith(f, ".jl") && f!= basename(rootfile), tests)
        else
            tests = map(f->string(f, ".jl"), tests)
        end
        println();
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

@testset "All tests" begin
    @includetests ARGS
end


text:
abstract type AbstractHashNode{T, N, p} <: AbstractNode{T, N, p} end

Base.IndexStyle(::Type{<: AbstractHashNode}) = IndexLinear()

@inline function Base.getindex(x::AbstractHashNode, i::Int)
    @boundscheck checkbounds(x, i)
    isnull(x) && return null(childtype(x))
    # data can haskey even if the mask is false
    # so trying to return stored data
    # in this case, need to handle stored data very carefully because
    # we can change its data unexpectedly even if the mask is false
    get(x.data, i, null(childtype(x)))
end

@inline function Base.setindex!(x::AbstractHashNode, v, i::Int)
    @boundscheck checkbounds(x, i)
    isnull(v) && throw(ArgumentError("trying to set null node, call `delete!(node, i)` instead"))
    @inbounds x.data[i] = v
    x
end

@inline function Base.setindex!(x::AbstractHashNode, ::Nothing, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        if isactive(x, i)
            x.data[i] = nothing # make `mask` `false` at `i`
            deactivate!(unsafe_getindex(x, i))
        end
    end
    x
end

@inline function deactivate!(x::AbstractHashNode)
    isnull(x) && return x
    fill!(getmask(x), false)
    foreach(deactivate!, values(x.data))
    x
end

function Base.delete!(x::AbstractHashNode, i)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        # only this case allows to set null node in `x`
        x.data[i] = nothing # make mask false
        delete!(x.data, i) # not set null, just delete it
    end
    x
end

function allocate!(x::AbstractHashNode{T}, i::Int) where {T}
    @boundscheck checkbounds(x, i)
    @inbounds begin
        if haskey(x.data, i)
            childnode = unsafe_getindex(x, i)
        else
            childnode = T()
        end
        x[i] = childnode # activated in setindex!
    end
    Activated(x, i)
end

isallocated(x::AbstractHashNode, i::Int) = haskey(x.data, i)

struct HashNode{T <: AbstractNode, N, p} <: AbstractHashNode{T, N, p}
    data::MaskedHashArray{T, N}
    HashNode{T, N, p}(data) where {T, N, p} = new(data)
    HashNode{T, N, p}(::UndefInitializer) where {T, N, p} = new()
end

function HashNode{T, N, p}() where {T, N, p}
    dims = size(HashNode{T, N, p})
    data = MaskedHashArray{T}(undef, dims)
    HashNode{T, N, p}(data)
end

Base.size(x::HashNode) = size(typeof(x))

struct DynamicHashNode{T <: AbstractNode, N} <: AbstractHashNode{T, N, Dynamic()}
    data::MaskedHashArray{T, N}
    dims::NTuple{N, Int}
    DynamicHashNode{T, N}(data, dims) where {T, N} = new(data, dims)
    DynamicHashNode{T, N}(::UndefInitializer) where {T, N} = new()
end

function DynamicHashNode{T, N}(dims::Int...) where {T, N}
    data = MaskedHashArray{T}(undef, dims)
    DynamicHashNode{T, N}(data, dims)
end

Base.size(x::DynamicHashNode) = x.dims


text:
using Random
using PyCall
using ArviZ: attributes

try
    ArviZ.initialize_bokeh()
catch
    @info "bokeh backend not available. bokeh tests will be skipped."
end

py"""
class PyNullObject(object):
   def __init__(self):
       pass
"""

function create_model(seed=10)
    rng = MersenneTwister(seed)
    J = 8
    nchains = 4
    ndraws = 500
    data = Dict(
        "J" => J,
        "y" => [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0],
        "sigma" => [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0],
    )
    posterior = Dict(
        "mu" => randn(rng, nchains, ndraws),
        "tau" => abs.(randn(rng, nchains, ndraws)),
        "eta" => randn(rng, nchains, ndraws, J),
        "theta" => randn(rng, nchains, ndraws, J),
    )
    posterior_predictive = Dict("y" => randn(rng, nchains, ndraws, J))
    sample_stats = Dict(
        "energy" => randn(rng, nchains, ndraws),
        "diverging" => (randn(rng, nchains, ndraws) .> 0.90),
        "max_depth" => (randn(rng, nchains, ndraws) .> 0.90),
        "log_likelihood" => randn(rng, nchains, ndraws, J),
    )
    prior = Dict(
        "mu" => randn(rng, nchains, ndraws) / 2,
        "tau" => abs.(randn(rng, nchains, ndraws)) / 2,
        "eta" => randn(rng, nchains, ndraws, data["J"]) / 2,
        "theta" => randn(rng, nchains, ndraws, data["J"]) / 2,
    )
    prior_predictive = Dict("y" => randn(rng, nchains, ndraws, J) / 2)
    sample_stats_prior = Dict(
        "energy" => randn(rng, nchains, ndraws),
        "diverging" => Int.(randn(rng, nchains, ndraws) .> 0.95),
    )
    model = from_dict(;
        posterior=posterior,
        posterior_predictive=posterior_predictive,
        sample_stats=sample_stats,
        prior=prior,
        prior_predictive=prior_predictive,
        sample_stats_prior=sample_stats_prior,
        observed_data=Dict("y" => data["y"]),
        dims=Dict("y" => ["obs_dim"], "log_likelihood" => ["obs_dim"]),
        coords=Dict("obs_dim" => 1:J),
    )
    return model
end

models() = (model_1=create_model(10), model_2=create_model(11))

function noncentered_schools_data()
    return Dict(
        "J" => 8,
        "y" => [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0],
        "sigma" => [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0],
    )
end

function check_multiple_attrs(test_dict, parent)
    failed_attrs = []
    for (dataset_name, attributes) in test_dict
        if Symbol(dataset_name) ∈ propertynames(parent)
            dataset = getproperty(parent, Symbol(dataset_name))
            for attribute in attributes
                if Symbol(attribute) ∉ propertynames(dataset)
                    push!(failed_attrs, (dataset_name, attribute))
                end
            end
        else
            push!(failed_attrs, dataset_name)
        end
    end
    return failed_attrs
end

dimsizes(ds) = ds._dims
convertindex(x::AbstractArray) = x
convertindex(o::PyObject) = o.array.values
vardict(ds) = Dict(k => convertindex(v._data) for (k, v) in ds._variables)
dimdict(ds) = Dict(k => v._dims for (k, v) in ds._variables)

function test_namedtuple_data(
    idata, group, names, nchains, ndraws; library="MyLib", coords=Dict(), dims=Dict()
)
    @test idata isa InferenceData
    @test group in ArviZ.groupnames(idata)
    ds = getproperty(idata, group)
    sizes = dimsizes(ds)
    @test length(sizes) == 2 + length(coords)
    vars = vardict(ds)
    for name in string.(names)
        @test name in keys(vars)
        dim = get(dims, name, get(dims, Symbol(name), []))
        s = (x -> length(get(coords, x, get(coords, Symbol(x), [])))).(dim)
        @test size(vars[name]) == (nchains, ndraws, s...)
    end
    @test "inference_library" in keys(attributes(ds))
    @test attributes(ds)["inference_library"] == library
    return nothing
end